1 |
edhill |
1.1 |
t0=input('Enter temperature : ') |
2 |
|
|
|
3 |
|
|
% Interpolate Dw on T=T0 surface |
4 |
|
|
|
5 |
|
|
h=zeros(Nx,Ny); |
6 |
|
|
hh=zeros(Nx,Ny); |
7 |
|
|
hhh=zeros(Nx,Ny); |
8 |
|
|
Dh=zeros(Nx,Ny); |
9 |
|
|
|
10 |
|
|
Nz1=4; |
11 |
|
|
he=dz*(Nz1-1); |
12 |
|
|
|
13 |
|
|
meantheta(:,:,Nz)=0.; |
14 |
|
|
for i=1:Nx |
15 |
|
|
for j=1:Ny |
16 |
|
|
kk=find(meantheta(i,j,:)<t0); |
17 |
|
|
|
18 |
|
|
if kk(1)>1 |
19 |
|
|
h(i,j)=(kk(1)-1)*dz+dz*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1))); |
20 |
|
|
|
21 |
|
|
Dh(i,j)=Dw(i,j,kk(1)-1) ... |
22 |
|
|
+(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ... |
23 |
|
|
*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1))); |
24 |
|
|
else |
25 |
|
|
h(i,j)=0; |
26 |
|
|
Dh(i,j)=0; |
27 |
|
|
end |
28 |
|
|
|
29 |
|
|
end |
30 |
|
|
end |
31 |
|
|
|
32 |
|
|
hmax=max(max(h)) |
33 |
|
|
|
34 |
|
|
|
35 |
|
|
figure |
36 |
|
|
|
37 |
|
|
%pcolor(Dh.');shading flat; colorbar; axis square |
38 |
|
|
% title('mean Wstar'); |
39 |
|
|
|
40 |
|
|
%figure |
41 |
|
|
|
42 |
|
|
%pcolor(hhh.');shading flat; colorbar; axis square |
43 |
|
|
%pcolor(flipud(h.'));shading flat; colorbar; axis square |
44 |
|
|
%title(['mean H from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]); |
45 |
|
|
%figure |
46 |
|
|
|
47 |
|
|
|
48 |
|
|
we=0; |
49 |
|
|
ws=0; |
50 |
|
|
WE=zeros(Nx,Ny); |
51 |
|
|
WS=zeros(Nx,Ny); |
52 |
|
|
hh=zeros(Nx,Ny); |
53 |
|
|
|
54 |
|
|
%for i=Nx/2:Nx |
55 |
|
|
for i=1:Nx |
56 |
|
|
for j=1:Ny |
57 |
|
|
|
58 |
|
|
r=sqrt((i-3*Nx/4)*(i-3*Nx/4)+(j-Ny/2)*(j-Ny/2))*dx; |
59 |
|
|
WE(i,j)=meanw(i,j,Nz1); |
60 |
|
|
WS(i,j)=Dh(i,j); |
61 |
|
|
|
62 |
|
|
% if r<0.40 |
63 |
|
|
if h(i,j)>he |
64 |
|
|
we=meanw(i,j,Nz1)*dx*dx+we; |
65 |
|
|
ws=Dh(i,j)*dx*dx+ws; |
66 |
|
|
% end |
67 |
|
|
end |
68 |
|
|
end |
69 |
|
|
end |
70 |
|
|
|
71 |
|
|
hh=h; |
72 |
|
|
|
73 |
|
|
for k=1:3 |
74 |
|
|
WE1(:,:,k)=WE(:,:); |
75 |
|
|
WS1(:,:,k)=WS(:,:); |
76 |
|
|
hh1(:,:,k)=hh(:,:); |
77 |
|
|
end |
78 |
|
|
|
79 |
|
|
|
80 |
|
|
for k=1:5 |
81 |
|
|
WE1=smooth3(WE1); |
82 |
|
|
WS1=smooth3(WS1); |
83 |
|
|
hh1=smooth3(hh1); |
84 |
|
|
end |
85 |
|
|
|
86 |
|
|
hhh(2:Nx-1,2:Ny-1)=10000*(hh1(1:Nx-2,2:Ny-1,1)+hh1(3:Nx,2:Ny-1,1)+... |
87 |
|
|
hh1(2:Nx-1,1:Ny-2,1)+hh1(2:Nx-1,3:Ny,1)-4*hh1(2:Nx-1,2:Ny-1,1)); |
88 |
|
|
|
89 |
|
|
%hhh(2:Nx-1,2:Ny-1)=10000*(hh1(1:Nx-2,2:Ny-1,1)-hh1(3:Nx,2:Ny-1,1)).*... |
90 |
|
|
%(hh1(1:Nx-2,2:Ny-1,1)-hh1(3:Nx,2:Ny-1,1))+... |
91 |
|
|
%10000*(hh1(2:Nx-1,1:Ny-2,1)-hh1(2:Nx-1,3:Ny,1)).*... |
92 |
|
|
%(hh1(2:Nx-1,1:Ny-2,1)-hh1(2:Nx-1,3:Ny,1)); |
93 |
|
|
|
94 |
|
|
i=find(h<0.015); |
95 |
|
|
h(i)=0; |
96 |
|
|
% i=find(abs(hhh)>0.001); |
97 |
|
|
% hhh(i)=0; |
98 |
|
|
|
99 |
|
|
for i=1:Nx |
100 |
|
|
for j=1:Ny |
101 |
|
|
%r=sqrt((i-3*Nx/4)*(i-3*Nx/4)+(j-Ny/2)*(j-Ny/2))*dx; |
102 |
|
|
|
103 |
|
|
%if r>0.44 |
104 |
|
|
%hhh(i,j)=NaN; |
105 |
|
|
%hh(i,j)=NaN; |
106 |
|
|
%WS1(i,j,1)=NaN; |
107 |
|
|
%end |
108 |
|
|
|
109 |
|
|
if h(i,j)<he |
110 |
|
|
WE1(i,j,1)=0; |
111 |
|
|
WS1(i,j,1)=NaN; |
112 |
|
|
hh(i,j)=NaN; |
113 |
|
|
hhh(i,j)=NaN; |
114 |
|
|
end |
115 |
|
|
|
116 |
|
|
end |
117 |
|
|
end |
118 |
|
|
|
119 |
|
|
|
120 |
|
|
%pcolor(hh(1:Nx,:).');shading flat; colorbar; axis square |
121 |
|
|
% set(gca,'DataAspectRatio',[2,2,2]) |
122 |
|
|
% title('h','FontSize',16); |
123 |
|
|
%figure |
124 |
|
|
|
125 |
|
|
x=1:3:Nx-2; |
126 |
|
|
y=1:3:Ny; |
127 |
|
|
hh(:,:)=max(he,hh(:,:)); |
128 |
|
|
subplot('position',[0.30 0.7 0.4 0.3]) |
129 |
|
|
surf(x,y,200*(-hh(1:3:Nx-2,1:3:Ny)'+he)); |
130 |
|
|
axis([1 Nx-4 0 Ny-2 -10.0 0]) |
131 |
|
|
a=0:1:63; |
132 |
|
|
a=a/63; |
133 |
|
|
cmap=[a' a' a']; |
134 |
|
|
colormap(cmap); |
135 |
|
|
set(gca,'DataAspectRatio',[20,20,7]) |
136 |
|
|
set(gca,'XtickLabel','|') |
137 |
|
|
set(gca,'YtickLabel','|') |
138 |
|
|
set(gca,'ZtickLabel','|') |
139 |
|
|
xlabel('x','FontSize',20) |
140 |
|
|
ylabel('y','FontSize',20) |
141 |
|
|
zlabel('z','FontSize',20,'Rotation',[0]) |
142 |
|
|
%title('h','FontSize',19); |
143 |
|
|
view(-20,10); |
144 |
|
|
|
145 |
|
|
hhh(1:Nx,:)=-10*hhh(1:Nx,:); |
146 |
|
|
% hhh(Nx/2:Nx,:)=100*hhh(Nx/2:Nx,:); |
147 |
|
|
a=max(max(hhh(1:Nx,:))); |
148 |
|
|
|
149 |
|
|
%pcolor(hhh(1:Nx,:).');caxis([0 a]);shading flat; colorbar; axis square |
150 |
|
|
%set(gca,'DataAspectRatio',[2,2,2]) |
151 |
|
|
%title('-\nabla^2 h','FontSize',12); |
152 |
|
|
%xlabel('x','Fontsize',13); |
153 |
|
|
%ylabel('y','Fontsize',13); |
154 |
|
|
%text(50,-10,'x'); |
155 |
|
|
%text(-10,60,'y'); |
156 |
|
|
|
157 |
|
|
%pcolor(squeeze(WE1(Nx/2:Nx,:,1)).');shading flat; colorbar; axis square |
158 |
|
|
% set(gca,'DataAspectRatio',[2,2,2]) |
159 |
|
|
%title('Ek, smooth'); |
160 |
|
|
%figure |
161 |
|
|
|
162 |
|
|
%pcolor(meanw(Nx/2:Nx,:,3).');shading flat; colorbar; axis square |
163 |
|
|
%set(gca,'DataAspectRatio',[2,2,2]) |
164 |
|
|
%title('Ek'); |
165 |
|
|
figure |
166 |
|
|
|
167 |
|
|
subplot('position',[0.15 0.05 0.7 0.5]) |
168 |
|
|
WS1(1:Nx,:,1)=-WS1(1:Nx,:,1)*16667; |
169 |
|
|
a=max(max(WS1(1:Nx,:,1))); |
170 |
|
|
pcolor(squeeze(WS1(1:Nx,:,1)).');caxis([0 a]);shading flat;colormap('default'); |
171 |
|
|
set(gca,'DataAspectRatio',[2,2,2]) |
172 |
|
|
set(gca,'Xtick',[1 100],'XtickLabel','|') |
173 |
|
|
set(gca,'Ytick',[1 120],'YtickLabel','|') |
174 |
|
|
title('w^*','FontSize',20); |
175 |
|
|
xlabel('x','Fontsize',20); |
176 |
|
|
ylabel('y','Fontsize',20,'Rotation',[0]); |
177 |
|
|
|
178 |
|
|
|
179 |
|
|
we |
180 |
|
|
ws |
181 |
|
|
|