| 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 |  |  |  |