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

Annotation of /MITgcm_contrib/timour_matlab/mscripts/DIAGdonut.m

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


Revision 1.1 - (hide annotations) (download)
Wed Sep 3 21:22:22 2003 UTC (21 years, 10 months ago) by edhill
Branch: MAIN
CVS Tags: HEAD
initial checkin of Timour's MatLAB scripts

1 edhill 1.1 clear path
2    
3     global Nx Ny Nz
4     global dx dy
5     global lat long dz dm mdep
6     global delt_su su_its t_su delt
7     global descriptor this_path
8     global f0 deltaf Q beta r_expt r_heat H
9     global time rots it
10     global g Cp rho_bar alpha
11     global u v t w
12     global iteration
13    
14    
15     param_file_name = ...
16     input(' Please enter the name of the m-file with the parameters for this run : ','s') ;
17     feval(param_file_name) ;
18    
19     dz=0.005;
20     beta=0;
21     f=f0;
22     iteration
23    
24     itstart = input(' Please enter start iteration : ','s')
25     itend = input(' Please enter end iteration : ','s')
26    
27    
28     sizeit=size(iteration);
29     for i=1:sizeit(1)
30     iter(i)=eval(iteration(i,1:10));
31     end
32     nitstart=find(iter==eval(itstart))
33     nitend=find(iter==eval(itend))
34    
35     path = this_path
36     cmdstr=['cd ' path ];
37     eval(cmdstr);
38     path=pwd
39    
40     sumtheta=zeros(Nx,Ny,Nz);
41     sumu=zeros(Nx,Ny,Nz);
42     sumv=zeros(Nx,Ny,Nz);
43     sumw=zeros(Nx,Ny,Nz);
44     counter=0;
45    
46     for i=nitstart:nitend
47    
48     path = this_path;
49     eval(['cd ' path '/it',int2str(iter(i))]);
50     path=pwd;
51    
52     i,0
53     tfn='theta.sun.b';
54     ufn='u.sun.b';
55     vfn='v.sun.b';
56    
57     % T (m/s)
58     t=Readmodel(tfn,[Nx Ny Nz],'float32');
59     t=Extract([3 Nx Ny Nz t(5:4+Nx*Ny*Nz)],[1 Nx 1 Ny 1 Nz]);
60     t=reshape(t,Nx,Ny,Nz);
61    
62     % U (m/s)
63     u=Readmodel(ufn,[Nx Ny Nz],'float32');
64     u=Extract([3 Nx Ny Nz u(5:4+Nx*Ny*Nz)],[1 Nx 1 Ny 1 Nz]);
65     u=reshape(u,Nx,Ny,Nz);
66    
67     % V (m/s)
68     v=Readmodel(vfn,[Nx Ny Nz],'float32');
69     v=Extract([3 Nx Ny Nz v(5:4+Nx*Ny*Nz)],[1 Nx 1 Ny 1 Nz]);
70     v=reshape(v,Nx,Ny,Nz);
71    
72     % W (m/s)
73    
74     w=zeros(Nx,Ny,Nz+1);
75    
76     for k=Nz:-1:1,
77     w(1:Nx-1,1:Ny-1,k)= w(1:Nx-1,1:Ny-1,k+1)...
78     -u(2:Nx,1:Ny-1,k)*dz/dx +u(1:Nx-1,1:Ny-1,k)*dz/dx ...
79     -v(1:Nx-1,2:Ny,k)*dz/dy +v(1:Nx-1,1:Ny-1,k)*dz/dy ;
80     end
81    
82    
83     %tfilename=(['T.' iterations((i),1:10) ]) ;
84     %t=rdmds(tfilename,'b');
85    
86     %ufilename=(['U.' iterations((i),1:10) ]) ;
87     %u=rdmds(ufilename,'b');
88     u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));
89    
90     %vfilename=(['V.' iterations((i),1:10) ]) ;
91     %v=rdmds(vfilename,'b');
92     v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));
93    
94     %wfilename=(['W.' iterations((i),1:10) ]) ;
95     %w=rdmds(wfilename,'b');
96    
97     sumtheta=sumtheta+t;
98     sumu=sumu+u;
99     sumv=sumv+v;
100     sumw=sumw+w(:,:,1:Nz);
101     counter=counter+1;
102     end
103    
104     meantheta=sumtheta/counter;
105     meanu=sumu/counter;
106     meanv=sumv/counter;
107     meanw=sumw/counter;
108    
109     sumut=zeros(Nx,Ny,Nz);
110     sumvt=zeros(Nx,Ny,Nz);
111     sumwt=zeros(Nx,Ny,Nz);
112    
113     sumuu=zeros(Nx,Ny,Nz);
114     sumvu=zeros(Nx,Ny,Nz);
115     sumwu=zeros(Nx,Ny,Nz);
116     sumvv=zeros(Nx,Ny,Nz);
117     sumwv=zeros(Nx,Ny,Nz);
118    
119     %meanu=smooth3(meanu);
120     %meanv=smooth3(meanv);
121     %meanw=smooth3(meanw);
122     %meantheta=smooth3(meantheta);
123    
124     for i=nitstart:nitend
125     i
126     path = this_path;
127     eval(['cd ' path '/it',int2str(iter(i))]);
128     path=pwd;
129    
130     % T (m/s)
131     t=Readmodel(tfn,[Nx Ny Nz],'float32');
132     t=Extract([3 Nx Ny Nz t(5:4+Nx*Ny*Nz)],[1 Nx 1 Ny 1 Nz]);
133     t=reshape(t,Nx,Ny,Nz);
134    
135     % U (m/s)
136     u=Readmodel(ufn,[Nx Ny Nz],'float32');
137     u=Extract([3 Nx Ny Nz u(5:4+Nx*Ny*Nz)],[1 Nx 1 Ny 1 Nz]);
138     u=reshape(u,Nx,Ny,Nz);
139    
140     % V (m/s)
141     v=Readmodel(vfn,[Nx Ny Nz],'float32');
142     v=Extract([3 Nx Ny Nz v(5:4+Nx*Ny*Nz)],[1 Nx 1 Ny 1 Nz]);
143     v=reshape(v,Nx,Ny,Nz);
144    
145     % W (m/s)
146    
147     w=zeros(Nx,Ny,Nz+1);
148    
149     for k=Nz:-1:1,
150     w(1:Nx-1,1:Ny-1,k)= w(1:Nx-1,1:Ny-1,k+1)...
151     -u(2:Nx,1:Ny-1,k)*dz/dx +u(1:Nx-1,1:Ny-1,k)*dz/dx ...
152     -v(1:Nx-1,2:Ny,k)*dz/dy +v(1:Nx-1,1:Ny-1,k)*dz/dy ;
153     end
154    
155    
156    
157     %tfilename=(['T.' iterations((i),1:10) ]) ;
158     %t=rdmds(tfilename,'b');
159    
160     %ufilename=(['U.' iterations((i),1:10) ]) ;
161     %u=rdmds(ufilename,'b');
162     u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));
163    
164     %vfilename=(['V.' iterations((i),1:10) ]) ;
165     %v=rdmds(vfilename,'b');
166     v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));
167    
168     %wfilename=(['W.' iterations((i),1:10) ]) ;
169     %w=rdmds(wfilename,'b');
170    
171     t=t-meantheta;
172     u=u-meanu;
173     v=v-meanv;
174     w(:,:,1:Nz)=w(:,:,1:Nz)-meanw;
175    
176     sumut=sumut+u.*t;
177     sumvt=sumvt+v.*t;
178     sumwt=sumwt+w(:,:,1:Nz).*t;
179    
180     sumuu=sumuu+u.*u;
181     sumvu=sumvu+v.*u;
182     sumwu=sumwu+w(:,:,1:Nz).*u;
183     sumvv=sumvv+v.*v;
184     sumwv=sumwv+w(:,:,1:Nz).*v;
185     end
186    
187     sumut=sumut/counter;
188     sumvt=sumvt/counter;
189     sumwt=sumwt/counter;
190    
191     sumuu=sumuu/counter;
192     sumvu=sumvu/counter;
193     sumwu=sumwu/counter;
194     sumvv=sumvv/counter;
195     sumwv=sumwv/counter;
196    
197     Mu=zeros(Nx,Ny,Nz);
198     Mv=zeros(Nx,Ny,Nz);
199     D=zeros(Nx,Ny,Nz);
200    
201     %dx=dm;
202     %dy=dm;
203     dz=-dz;
204    
205     Mu(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumuu(3:Nx,2:Ny-1,2:Nz-1)-sumuu(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
206     -(sumvu(2:Nx-1,3:Ny,2:Nz-1)-sumvu(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
207     -(sumwu(2:Nx-1,2:Ny-1,3:Nz)-sumwu(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
208    
209     Mv(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumvu(3:Nx,2:Ny-1,2:Nz-1)-sumvu(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
210     -(sumvv(2:Nx-1,3:Ny,2:Nz-1)-sumvv(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
211     -(sumwv(2:Nx-1,2:Ny-1,3:Nz)-sumwv(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
212    
213     D(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumut(3:Nx,2:Ny-1,2:Nz-1)-sumut(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
214     -(sumvt(2:Nx-1,3:Ny,2:Nz-1)-sumvt(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
215     -(sumwt(2:Nx-1,2:Ny-1,3:Nz)-sumwt(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
216    
217     pv=zeros(Nx,Ny,Nz);
218     ff=zeros(Nx,Ny,Nz);
219    
220    
221     for j=1:Ny
222     ff(:,j,:)=f+(j-1)*dy*beta;
223     end
224    
225     pv(2:Nx-1,2:Ny-1,2:Nz-1)=(meanv(3:Nx,2:Ny-1,2:Nz-1)/(2*dx)-meanv(1:Nx-2,2:Ny-1,2:Nz-1)/(2*dx) ...
226     -meanu(2:Nx-1,3:Ny,2:Nz-1)/(2*dy)+meanu(2:Nx-1,1:Ny-2,2:Nz-1)/(2*dy) ...
227     +ff(2:Nx-1,2:Ny-1,2:Nz-1)) ...
228     .*(meantheta(2:Nx-1,2:Ny-1,3:Nz)-meantheta(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz) ...
229     -(meanv(2:Nx-1,2:Ny-1,3:Nz)-meanv(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz) ...
230     .*(meantheta(3:Nx,2:Ny-1,2:Nz-1)-meantheta(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx)+ ...
231     (meanu(2:Nx-1,2:Ny-1,3:Nz)-meanu(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz) ...
232     .*(meantheta(2:Nx-1,3:Ny,2:Nz-1)-meantheta(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy);
233    
234     Mp=zeros(Nx,Ny,Nz);
235     Dp=zeros(Nx,Ny,Nz);
236    
237     Mp(3:Nx-2,3:Ny-2,3:Nz-2)=((Mv(4:Nx-1,3:Ny-2,3:Nz-2)-Mv(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx)- ...
238     (Mu(3:Nx-2,4:Ny-1,3:Nz-2)-Mu(3:Nx-2,2:Ny-3,3:Nz-2))/(2*dy)).* ...
239     (meantheta(3:Nx-2,3:Ny-2,4:Nz-1)-meantheta(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz) ...
240     -(Mv(3:Nx-2,3:Ny-2,4:Nz-1)-Mv(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz).* ...
241     (meantheta(4:Nx-1,3:Ny-2,3:Nz-2)-meantheta(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx) ...
242     +(Mu(3:Nx-2,3:Ny-2,4:Nz-1)-Mu(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz).* ...
243     (meantheta(3:Nx-2,4:Ny-1,3:Nz-2)-meantheta(3:Nx-2,2:Ny-3,3:Nz-2))/(2*dy);
244    
245     Dp(3:Nx-2,3:Ny-2,3:Nz-2)=(D(3:Nx-2,3:Ny-2,4:Nz-1)-D(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz).* ...
246     (meanv(4:Nx-1,3:Ny-2,3:Nz-2)/(2*dx)-meanv(2:Nx-3,3:Ny-2,3:Nz-2)/(2*dx) ...
247     -meanu(3:Nx-2,4:Ny-1,3:Nz-2)/(2*dy)+meanu(3:Nx-2,2:Ny-3,3:Nz-2)/(2*dy)+ff(3:Nx-2,3:Ny-2,3:Nz-2)) ...
248     -(D(4:Nx-1,3:Ny-2,3:Nz-2)-D(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx).* ...
249     (meanv(3:Nx-2,3:Ny-2,4:Nz-1)-meanv(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz)+ ...
250     (D(3:Nx-2,4:Ny-1,3:Nz-2)-D(3:Nx-2,2:Ny-3,3:Nz-2))/(2*dy).* ...
251     (meanu(3:Nx-2,3:Ny-2,4:Nz-1)-meanu(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz);
252    
253    
254     sumM=zeros(Nx,Ny);
255     sumD=zeros(Nx,Ny);
256     sumW=zeros(Nx,Ny);
257    
258     dz=-dz;
259    
260     % DEFINE BOX OF INTEGRATION
261     Nx1=3;
262     Nx2=Nx-2;
263     Ny1=3;
264     %Ny2=Ny-2;
265     Ny2=Ny/2-2;
266     Nz1=5;
267     Nz2=Nz-4;
268    
269    
270     for k=Nz1:Nz2
271     sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+Mp(3:Nx-2,3:Ny-2,k).*pv(3:Nx-2,3:Ny-2,k)*dz;
272     sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+Dp(3:Nx-2,3:Ny-2,k).*pv(3:Nx-2,3:Ny-2,k)*dz;
273     end
274     sumW(3:Nx-2,3:Ny-2)=0.5*meanw(3:Nx-2,3:Ny-2,Nz1).*pv(3:Nx-2,3:Ny-2,Nz1).*pv(3:Nx-2,3:Ny-2,Nz1);
275    
276     %contourf(sumM',20);colorbar;
277     %figure
278     %contourf(sumD',20);colorbar;
279    
280     V=[0.25*min(min(sumW(Nx1:Nx2,Ny1:Ny2-4))) 0];
281    
282     title='Vorticity input';
283     imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2-4)');shading flat;caxis(V);axis image;colorbar('vertical');
284     set(gca,'ydir','norm')
285     text(0,110,title);
286     figure
287     V=[-0.5*max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4)))) 0.5*max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4))))]
288     title='Buoyancy eddy-transfer integral';
289     imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2-4)');shading flat;caxis(V);axis image;colorbar('vertical');
290     set(gca,'ydir','norm')
291     text(0,110,title);
292     figure
293     V=[-0.5*max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4)))) 0.5*max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4))))]
294     title='Momentum eddy-transfer integral';
295     imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2-4)');shading flat;caxis(V);axis image;colorbar('vertical');
296     set(gca,'ydir','norm')
297     text(0,110,title);
298    
299    
300     %-----TEST-------------------------------------------------
301    
302     sum=0;
303     sumT=0;
304     for i=Nx1:Nx2
305     for j=Ny1:Ny2
306     sum=sum+0.5*meanw(i,j,Nz1)*pv(i,j,Nz1)*pv(i,j,Nz1)*dx*dy;
307     sumT=sum-0.5*meanw(i,j,Nz2)*pv(i,j,Nz2)*pv(i,j,Nz2)*dx*dy;
308     end
309     end
310    
311     for j=Ny1:Ny2
312     for k=Nz1:Nz2
313     sumT=sumT+0.5*meanu(Nx2,j,k)*pv(Nx2,j,k)*pv(Nx2,j,k)*dz*dy;
314     sumT=sumT-0.5*meanu(Nx1,j,k)*pv(Nx1,j,k)*pv(Nx1,j,k)*dz*dy;
315     end
316     end
317    
318     for i=Nx1:Nx2
319     for k=Nz1:Nz2
320     sumT=sumT+0.5*meanv(i,Ny2,k)*pv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
321     sumT=sumT-0.5*meanv(i,Ny1,k)*pv(i,Ny1,k)*pv(i,Ny1,k)*dz*dx;
322     end
323     end
324     'VORTICITY GENERATION'
325     sum
326     sumT
327    
328    
329     sum=0;
330     for i=Nx1:Nx2
331     for j=Ny1:Ny2
332     sum=sum+sumD(i,j)*dx*dy;
333     end
334     end
335     'DISSIPATION BY EDDY-DIFFUSIVITY'
336     sum
337    
338     sum=0;
339     for i=Nx1:Nx2
340     for j=Ny1:Ny2
341     sum=sum+sumM(i,j)*dx*dy;
342     end
343     end
344     'DISSIPATION BY EDDY-VISCOSITY'
345     sum

  ViewVC Help
Powered by ViewVC 1.1.22