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

Contents of /MITgcm_contrib/timour_matlab/mscripts/meant1.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 (20 years, 8 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='data11';
15 feval(param_file_name) ;
16
17 itstart='0000200000';
18 itend='0000240000';
19 iii=5;
20 istart=142;
21 iend=230;
22 t0=25;
23
24
25 iterations
26
27 sizeit=size(iterations);
28 for i=1:sizeit(1)
29 iter(i)=eval(iterations(i,1:10));
30 end
31
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 counter=0;
42
43 for i=nitstart:nitend
44 tfilename=(['T.' iterations((i),1:10) ]) ;
45 t=rdmds(tfilename,'b');
46 sumtheta=sumtheta+t;
47 counter=counter+1;
48 end
49 meantheta=sumtheta/counter;
50
51 W1=smooth3(meantheta);
52
53 for i=1:iii,
54 W1=smooth3(W1);
55 end
56
57 hh=zeros(Nx,1);
58 for i=1:Nx,
59 kk=find(W1(i,Ny/2,:)<t0);
60 if kk(1)>1
61 hh(i)=(kk(1)-1)+(W1(i,Ny/2,kk(1)-1)-t0)/(W1(i,Ny/2,kk(1)-1)-W1(i,Ny/2,kk(1)));
62 else
63 hh(i)=0;
64 end
65 end
66 hmax=max(hh);
67 I=find(hh==hmax);
68 x0=(I-istart)/(iend-istart)
69
70
71 contour(flipud(squeeze(W1(:,Ny/2,3:30))'),[25,25],'b')
72 hold on
73
74 itstart='0000240000';
75 itend='0000300000';
76 nitstart=find(iter==eval(itstart));
77 nitend=find(iter==eval(itend));
78
79 sumtheta=zeros(Nx,Ny,Nz);
80 counter=0;
81 for i=nitstart:nitend
82 tfilename=(['T.' iterations((i),1:10) ]) ;
83 t=rdmds(tfilename,'b');
84 sumtheta=sumtheta+t;
85 counter=counter+1;
86 end
87 meantheta=sumtheta/counter;
88 W2=smooth3(meantheta);
89 for i=1:iii,
90 W2=smooth3(W2);
91 end
92 hh=zeros(Nx,1);
93 for i=1:Nx,
94 kk=find(W2(i,Ny/2,:)<t0);
95 if kk(1)>1
96 hh(i)=(kk(1)-1)+(W2(i,Ny/2,kk(1)-1)-t0)/(W2(i,Ny/2,kk(1)-1)-W2(i,Ny/2,kk(1)));
97 else
98 hh(i)=0;
99 end
100 end
101 hmax=max(hh);
102 I=find(hh==hmax);
103 x0=(I-istart)/(iend-istart)
104 contour(flipud(squeeze(W2(:,Ny/2,3:30))'),[25,25],'b')
105 hold on
106
107 itstart='0000300000';
108 itend='0000350000';
109 nitstart=find(iter==eval(itstart));
110 nitend=find(iter==eval(itend));
111
112 sumtheta=zeros(Nx,Ny,Nz);
113 counter=0;
114 for i=nitstart:nitend
115 tfilename=(['T.' iterations((i),1:10) ]) ;
116 t=rdmds(tfilename,'b');
117 sumtheta=sumtheta+t;
118 counter=counter+1;
119 end
120 meantheta=sumtheta/counter;
121 W3=smooth3(meantheta);
122 for i=1:iii,
123 W3=smooth3(W3);
124 end
125 hh=zeros(Nx,1);
126 for i=1:Nx,
127 kk=find(W3(i,Ny/2,:)<t0);
128 if kk(1)>1
129 hh(i)=(kk(1)-1)+(W3(i,Ny/2,kk(1)-1)-t0)/(W3(i,Ny/2,kk(1)-1)-W3(i,Ny/2,kk(1)));
130 else
131 hh(i)=0;
132 end
133 end
134 hmax=max(hh);
135 I=find(hh==hmax);
136 x0=(I-istart)/(iend-istart)
137 contour(flipud(squeeze(W3(:,Ny/2,3:30))'),[25,25],'b')
138 hold on
139
140 itstart='0000350000';
141 itend='0000400000';
142 nitstart=find(iter==eval(itstart));
143 nitend=find(iter==eval(itend));
144
145 sumtheta=zeros(Nx,Ny,Nz);
146 counter=0;
147 for i=nitstart:nitend
148 tfilename=(['T.' iterations((i),1:10) ]) ;
149 t=rdmds(tfilename,'b');
150 sumtheta=sumtheta+t;
151 counter=counter+1;
152 end
153 meantheta=sumtheta/counter;
154 W4=smooth3(meantheta);
155 for i=1:iii,
156 W4=smooth3(W4);
157 end
158 hh=zeros(Nx,1);
159 for i=1:Nx,
160 kk=find(W4(i,Ny/2,:)<t0);
161 if kk(1)>1
162 hh(i)=(kk(1)-1)+(W4(i,Ny/2,kk(1)-1)-t0)/(W4(i,Ny/2,kk(1)-1)-W4(i,Ny/2,kk(1)));
163 else
164 hh(i)=0;
165 end
166 end
167 hmax=max(hh);
168 I=find(hh==hmax);
169 x0=(I-istart)/(iend-istart)
170 contour(flipud(squeeze(W4(:,Ny/2,3:30))'),[25,25],'b')
171 hold on
172
173 itstart='0000400000';
174 itend='0000500000';
175 nitstart=find(iter==eval(itstart));
176 nitend=find(iter==eval(itend));
177
178 sumtheta=zeros(Nx,Ny,Nz);
179 counter=0;
180 for i=nitstart:nitend
181 tfilename=(['T.' iterations((i),1:10) ]) ;
182 t=rdmds(tfilename,'b');
183 sumtheta=sumtheta+t;
184 counter=counter+1;
185 end
186 meantheta=sumtheta/counter;
187 W5=smooth3(meantheta);
188 for i=1:iii,
189 W5=smooth3(W5);
190 end
191 hh=zeros(Nx,1);
192 for i=1:Nx,
193 kk=find(W5(i,Ny/2,:)<t0);
194 if kk(1)>1
195 hh(i)=(kk(1)-1)+(W5(i,Ny/2,kk(1)-1)-t0)/(W5(i,Ny/2,kk(1)-1)-W5(i,Ny/2,kk(1)));
196 else
197 hh(i)=0;
198 end
199 end
200 hmax=max(hh);
201 I=find(hh==hmax);
202 x0=(I-istart)/(iend-istart)
203 contour(flipud(squeeze(W5(:,Ny/2,3:30))'),[25,25],'b')
204 hold off

  ViewVC Help
Powered by ViewVC 1.1.22