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

Annotation of /MITgcm_contrib/timour_matlab/mscripts/meant1.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 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