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

Contents of /MITgcm_contrib/timour_matlab/mscripts/DIAGW.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
initial checkin of Timour's MatLAB scripts

1 clear path
2
3 global Nx Ny Nz
4 global lat long dz dm mdep
5 global delt_su su_its t_su delt
6 global descriptor this_path
7 global f deltaf Q beta r_expt r_heat H
8 global time rots it
9 global g Cp rho_bar alpha
10 global u v t w
11 global iterations
12
13
14 param_file_name = ...
15 input(' Please enter the name of the m-file with the parameters for this run : ','s') ;
16 feval(param_file_name) ;
17
18 % iterations
19
20 itstart = input(' Please enter start iteration : ','s')
21 itend = input(' Please enter end iteration : ','s')
22
23
24 sizeit=size(iterations);
25 for i=1:sizeit(1)
26 iter(i)=eval(iterations(i,1:10));
27 end
28 nitstart=find(iter==eval(itstart))
29 nitend=find(iter==eval(itend))
30
31 path = this_path
32 cmdstr=['cd ' path ];
33 eval(cmdstr);
34 path=pwd
35
36 sumtheta=zeros(Nx,Ny,Nz);
37 sumu=zeros(Nx,Ny,Nz);
38 sumv=zeros(Nx,Ny,Nz);
39 sumw=zeros(Nx,Ny,Nz);
40 counter=0;
41
42 for i=nitstart:nitend
43 i,0
44 tfilename=(['T.' iterations((i),1:10) ]) ;
45 t=rdmds(tfilename,'b');
46
47 ufilename=(['U.' iterations((i),1:10) ]) ;
48 u=rdmds(ufilename,'b');
49 u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));
50
51 vfilename=(['V.' iterations((i),1:10) ]) ;
52 v=rdmds(vfilename,'b');
53 v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));
54
55 wfilename=(['W.' iterations((i),1:10) ]) ;
56 w=rdmds(wfilename,'b');
57
58 sumtheta=sumtheta+t;
59 sumu=sumu+u;
60 sumv=sumv+v;
61 sumw=sumw+w;
62 counter=counter+1;
63 end
64
65 meantheta=sumtheta/counter;
66 meanu=sumu/counter;
67 meanv=sumv/counter;
68 meanw=sumw/counter;
69
70 sumut=zeros(Nx,Ny,Nz);
71 sumvt=zeros(Nx,Ny,Nz);
72 sumwt=zeros(Nx,Ny,Nz);
73
74 meanu=smooth3(meanu);
75 meanv=smooth3(meanv);
76 meanw=smooth3(meanw);
77 meantheta=smooth3(meantheta);
78
79 for i=nitstart:nitend
80 i
81 tfilename=(['T.' iterations((i),1:10) ]) ;
82 t=rdmds(tfilename,'b');
83
84 ufilename=(['U.' iterations((i),1:10) ]) ;
85 u=rdmds(ufilename,'b');
86 u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));
87
88 vfilename=(['V.' iterations((i),1:10) ]) ;
89 v=rdmds(vfilename,'b');
90 v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));
91
92 wfilename=(['W.' iterations((i),1:10) ]) ;
93 w=rdmds(wfilename,'b');
94
95 t=t-meantheta;
96 u=u-meanu;
97 v=v-meanv;
98 w=w-meanw;
99
100 sumut=sumut+u.*t;
101 sumvt=sumvt+v.*t;
102 sumwt=sumwt+w.*t;
103
104 end
105
106 sumut=sumut/counter;
107 sumvt=sumvt/counter;
108 sumwt=sumwt/counter;
109
110
111 D=zeros(Nx,Ny,Nz);
112
113 dx=dm;
114 dy=dm;
115 dz=-dz;
116
117 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) ...
118 -(sumvt(2:Nx-1,3:Ny,2:Nz-1)-sumvt(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
119 -(sumwt(2:Nx-1,2:Ny-1,3:Nz)-sumwt(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
120
121 Dw=zeros(Nx,Ny,Nz);
122 ff=zeros(Nx,Ny,Nz);
123
124
125 for j=1:Ny
126 ff(:,j,:)=f+(j-1)*dy*beta;
127 end
128
129 Dw(:,:,2:Nz-1)=D(:,:,2:Nz-1)*2*dz./(meantheta(:,:,3:Nz)-meantheta(:,:,1:Nz-2));
130
131 i=find(meantheta(:,:,3:Nz)==meantheta(:,:,1:Nz-2));
132 Dw(i)=0;
133
134
135 t0=input('Enter temperature : ')
136
137 % Interpolate Dw on T=T0 surface
138
139 h=zeros(Nx,Ny);
140 Dh=zeros(Nx,Ny);
141 dz=-dz
142
143 meantheta(:,:,Nz)=0.;
144 for i=1:Nx
145 for j=1:Ny
146 kk=find(meantheta(i,j,:)<t0);
147
148 if kk(1)>1
149 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)));
150
151 Dh(i,j)=Dw(i,j,kk(1)-1) ...
152 +(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ...
153 *(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
154 else
155 h(i,j)=0;
156 Dh(i,j)=0;
157 end
158
159 end
160 end
161
162 hmax=max(max(h))
163
164 figure
165
166 pcolor(Dh.');shading flat; colorbar; axis square
167 title('mean Wstar');
168 % title(['mean Wstar from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]);
169
170 figure
171
172 pcolor(h.');shading flat; colorbar; axis square
173 %pcolor(flipud(h.'));shading flat; colorbar; axis square
174 %title(['mean H from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]);
175
176
177 we=0;
178 ws=0;
179 for i=1:Nx
180 for j=1:Ny
181
182 r=sqrt((i-3*Nx/4)*(i-3*Nx/4)+(j-Ny/2)*(j-Ny/2))*dx;
183
184 %if r<0.44
185
186 we=meanw(i,j,3)*dx*dx+we;
187 ws=Dh(i,j)*dx*dx+ws;
188 %end
189 end
190 end
191
192 we
193 ws
194
195
196
197
198
199

  ViewVC Help
Powered by ViewVC 1.1.22