| 1 |  | 
| 2 | colormap gray | 
| 3 | cmap=colormap; | 
| 4 | cmap=flipud(cmap); | 
| 5 | %cmap=sqrt(cmap); | 
| 6 | cmap=0.5*(1-cos(3.1415*cmap)); | 
| 7 | colormap(cmap) | 
| 8 |  | 
| 9 | sumM=zeros(Nx,Ny); | 
| 10 | sumD=zeros(Nx,Ny); | 
| 11 | sumW=zeros(Nx,Ny); | 
| 12 |  | 
| 13 | dz=0.005; | 
| 14 |  | 
| 15 | % DEFINE BOX OF INTEGRATION | 
| 16 | Nx1=4; | 
| 17 | Nx2=Nx-3; | 
| 18 | Ny1=4; | 
| 19 | Ny2=Ny-3; | 
| 20 | Nz1=5; | 
| 21 | Nz2=Nz-4; | 
| 22 |  | 
| 23 | pv1=pv; | 
| 24 | Mp1=Mp; | 
| 25 | Dp1=Dp; | 
| 26 |  | 
| 27 | %  pv1=smooth3(pv); | 
| 28 | %  Mp1=smooth3(Mp); | 
| 29 | %  Dp1=smooth3(Dp); | 
| 30 |  | 
| 31 | %  for kk=1:4 | 
| 32 | %  pv1=smooth3(pv1); | 
| 33 | %  Mp1=smooth3(Mp1); | 
| 34 | %  Dp1=smooth3(Dp1); | 
| 35 | %  end | 
| 36 |  | 
| 37 | for k=Nz1:Nz2 | 
| 38 | if k==Nz1 | 
| 39 | a=0.5; | 
| 40 | elseif k==Nz2 | 
| 41 | a=0.5; | 
| 42 | else | 
| 43 | a=1; | 
| 44 | end | 
| 45 |  | 
| 46 | sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+a*Mp1(3:Nx-2,3:Ny-2,k).*pv1(3:Nx-2,3:Ny-2,k)*dz; | 
| 47 | sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+a*Dp1(3:Nx-2,3:Ny-2,k).*pv1(3:Nx-2,3:Ny-2,k)*dz; | 
| 48 | end | 
| 49 | sumW(3:Nx-2,3:Ny-2)=0.5*meanw(3:Nx-2,3:Ny-2,Nz1).*pv1(3:Nx-2,3:Ny-2,Nz1).*pv1(3:Nx-2,3:Ny-2,Nz1); | 
| 50 |  | 
| 51 | for k=1:3 | 
| 52 | sumM1(:,:,k)=sumM(:,:); | 
| 53 | sumD1(:,:,k)=sumD(:,:); | 
| 54 | sumW1(:,:,k)=sumW(:,:); | 
| 55 | end | 
| 56 |  | 
| 57 | for k=1:25 | 
| 58 |  | 
| 59 | sumM1=smooth3(sumM1); | 
| 60 | sumD1=smooth3(sumD1); | 
| 61 | sumW1=smooth3(sumW1); | 
| 62 | end | 
| 63 |  | 
| 64 | V=[0.25*min(min(sumW(Nx1+2:Nx2,Ny1:Ny2-4))) 0]; | 
| 65 |  | 
| 66 | %       title='Vorticity input'; | 
| 67 | %        imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical'); | 
| 68 | %        set(gca,'ydir','norm') | 
| 69 | %        text(0,110,title); | 
| 70 | %figure | 
| 71 | v=zeros(10,1); | 
| 72 | v1=zeros(10,1); | 
| 73 | for i=1:10 | 
| 74 | %v(i)=-0.1*(i-0.5)*7; | 
| 75 | %v1(i)=0.1*(i-0.5)*7; | 
| 76 | v(i)=-0.1*(i)*10; | 
| 77 | v1(i)=0.1*(i)*10; | 
| 78 | end | 
| 79 |  | 
| 80 | %figure | 
| 81 |  | 
| 82 | % contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v) | 
| 83 | % hold on | 
| 84 | % contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--') | 
| 85 | % hold off | 
| 86 | %        text(20,0,'countour interval is 0.7 in the non-dimensional units'); | 
| 87 | %        text(20,95,'Vorticity input','Fontsize',16); | 
| 88 |  | 
| 89 | %figure | 
| 90 |  | 
| 91 | %V=[-0.5*max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4)))) 0.5*max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4))))]; | 
| 92 | %       title='Buoyancy diffusion integral'; | 
| 93 | %        imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical'); | 
| 94 | %        set(gca,'ydir','norm') | 
| 95 | %        text(0,110,title); | 
| 96 |  | 
| 97 |  | 
| 98 | figure | 
| 99 | x=Nx1:Nx2; | 
| 100 | y=Ny1:Ny2; | 
| 101 |  | 
| 102 | subplot(2,1,1) | 
| 103 | contourf(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',10) | 
| 104 | %       contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v) | 
| 105 | % hold on | 
| 106 | %       contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',10);colorbar; | 
| 107 | %       contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--') | 
| 108 | text(10,Ny2+5,'Buoyancy eddy-transfer integral','Fontsize',13); | 
| 109 | hold off | 
| 110 | xlabel('X ') | 
| 111 | ylabel('Y ') | 
| 112 | set(gca,'DataAspectRatio',[2,2,2]) | 
| 113 | set(gca,'XtickLabel','|||||||') | 
| 114 |  | 
| 115 | %figure | 
| 116 | %V=[-0.5*max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4)))) 0.5*max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4))))]; | 
| 117 | %       title='Momentum eddy-transfer integral'; | 
| 118 | %        imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical'); | 
| 119 | %        set(gca,'ydir','norm') | 
| 120 | %        text(0,110,title); | 
| 121 |  | 
| 122 | subplot(2,1,2) | 
| 123 | xx=Nx1:Nx2; | 
| 124 | yy=Ny1:Ny2; | 
| 125 | %       contour(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v) | 
| 126 | contourf(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',10) | 
| 127 |  | 
| 128 | text(10,Ny2+5,'Momentum eddy-transfer integral','Fontsize',13); | 
| 129 | %hold on | 
| 130 | %       contour(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',[0,0],' ') | 
| 131 | %       contour(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--') | 
| 132 | %       contour(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',10);colorbar; | 
| 133 | %hold off | 
| 134 | xlabel('X (gridpoints)') | 
| 135 | ylabel('Y (gridpoints)') | 
| 136 | set(gca,'DataAspectRatio',[2,2,2]) | 
| 137 | %        text(40,-25,'countour interval - 0.7 (non-dimensional units)','FontSize',7); | 
| 138 | %        text(40,-30,'negative values - solid line','FontSize',7); | 
| 139 | %        text(40,-35,'positive values - dashed line','FontSize',7); | 
| 140 |  | 
| 141 |  | 
| 142 |  | 
| 143 | %-----TEST------------------------------------------------- | 
| 144 |  | 
| 145 | sum=0; | 
| 146 | sumT=0; | 
| 147 | sumY=0; | 
| 148 | for i=Nx1:Nx2 | 
| 149 | for j=Ny1:Ny2 | 
| 150 |  | 
| 151 |  | 
| 152 | if j==Ny1 | 
| 153 | a=0.5; | 
| 154 | elseif j==Ny2 | 
| 155 | a=0.5; | 
| 156 | else | 
| 157 | a=1; | 
| 158 | end | 
| 159 |  | 
| 160 | if i==Nx1 | 
| 161 | b=0.5; | 
| 162 | elseif i==Nx2 | 
| 163 | b=0.5; | 
| 164 | else | 
| 165 | b=1; | 
| 166 | end | 
| 167 |  | 
| 168 |  | 
| 169 |  | 
| 170 | sum=sum+a*b*0.5*meanw(i,j,Nz1)*pv(i,j,Nz1)*pv(i,j,Nz1)*dx*dy; | 
| 171 | sumT=sum-0.5*meanw(i,j,Nz2)*pv(i,j,Nz2)*pv(i,j,Nz2)*dx*dy; | 
| 172 | end | 
| 173 | end | 
| 174 |  | 
| 175 | for j=Ny1:Ny2 | 
| 176 | for k=Nz1:Nz2 | 
| 177 | sumT=sumT+0.5*meanu(Nx2,j,k)*pv(Nx2,j,k)*pv(Nx2,j,k)*dz*dy; | 
| 178 | sumT=sumT-0.5*meanu(Nx1,j,k)*pv(Nx1,j,k)*pv(Nx1,j,k)*dz*dy; | 
| 179 | end | 
| 180 | end | 
| 181 |  | 
| 182 | for i=Nx1:Nx2 | 
| 183 | for k=Nz1:Nz2 | 
| 184 | sumT=sumT+0.5*meanv(i,Ny2,k)*pv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx; | 
| 185 | sumY=sumY+0.5*meanv(i,Ny2,k)*pv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx; | 
| 186 | sumT=sumT-0.5*meanv(i,Ny1,k)*pv(i,Ny1,k)*pv(i,Ny1,k)*dz*dx; | 
| 187 | end | 
| 188 | end | 
| 189 | 'VORTICITY GENERATION' | 
| 190 | sum | 
| 191 | sumY | 
| 192 | sumT | 
| 193 |  | 
| 194 |  | 
| 195 | sum=0; | 
| 196 | for i=Nx1:Nx2 | 
| 197 | for j=Ny1:Ny2 | 
| 198 |  | 
| 199 | if j==Ny1 | 
| 200 | a=0.5; | 
| 201 | elseif j==Ny2 | 
| 202 | a=0.5; | 
| 203 | else | 
| 204 | a=1; | 
| 205 | end | 
| 206 |  | 
| 207 | if i==Nx1 | 
| 208 | b=0.5; | 
| 209 | elseif i==Nx2 | 
| 210 | b=0.5; | 
| 211 | else | 
| 212 | b=1; | 
| 213 | end | 
| 214 |  | 
| 215 | sum=sum+a*b*sumD(i,j)*dx*dy; | 
| 216 | end | 
| 217 | end | 
| 218 | 'DISSIPATION BY EDDY-DIFFUSIVITY' | 
| 219 | sum | 
| 220 |  | 
| 221 | sum=0; | 
| 222 | for i=Nx1:Nx2 | 
| 223 | for j=Ny1:Ny2 | 
| 224 |  | 
| 225 | if j==Ny1 | 
| 226 | a=0.5; | 
| 227 | elseif j==Ny2 | 
| 228 | a=0.5; | 
| 229 | else | 
| 230 | a=1; | 
| 231 | end | 
| 232 |  | 
| 233 | if i==Nx1 | 
| 234 | b=0.5; | 
| 235 | elseif i==Nx2 | 
| 236 | b=0.5; | 
| 237 | else | 
| 238 | b=1; | 
| 239 | end | 
| 240 |  | 
| 241 | sum=sum+a*b*sumM(i,j)*dx*dy; | 
| 242 | end | 
| 243 | end | 
| 244 | 'DISSIPATION BY EDDY-VISCOSITY' | 
| 245 | sum |