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