/[MITgcm]/MITgcm_contrib/timour_matlab/mscripts/DIAG4.m
ViewVC logotype

Contents of /MITgcm_contrib/timour_matlab/mscripts/DIAG4.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Wed Sep 3 21:22:22 2003 UTC (21 years, 10 months ago) by edhill
Branch: MAIN
CVS Tags: HEAD
Error occurred while calculating annotation data.
initial checkin of Timour's MatLAB scripts

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

  ViewVC Help
Powered by ViewVC 1.1.22