/[MITgcm]/MITgcm_contrib/stephd_matlab/plotall.m
ViewVC logotype

Annotation of /MITgcm_contrib/stephd_matlab/plotall.m

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


Revision 1.1 - (hide annotations) (download)
Thu Sep 4 02:49:26 2003 UTC (21 years, 10 months ago) by edhill
Branch: MAIN
CVS Tags: HEAD
initial checkin

1 edhill 1.1 % PLOTALL
2     %%%%%%%%%%%%%%%
3     % plots horizontal maps, zonal transects, streamfunction, overturning
4     % etc.
5     % cleaned, updated: 9 Aug 01 swd
6     % change from plotresults: 8 May 01 swd
7     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8     % Set up for several grids...(via tags: bathynum and ocmipornot)
9     % USE AT YOUR OWN RISK
10     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
11     %
12     clear all
13     path(path,'/data25/stephd/Global/Functions');
14     %
15     %%%%%%%%%% CHOICES %%%%%%%%%%%%%%%%
16     %%%%%%%%%%
17     titstr=['0204.03'];
18     yeard=365;
19     run1str=['/data25/stephd/c46c3/newerice/0204.03'];
20     timestr=['0003139000']; dt=43200;
21     bathynum=10; % bathymetric set up
22     ocmipornot=0; %1=ocmip setup
23     %%%%%%%%%
24     kn=2; % level to plot
25     knsp=3; % level to plot for subplt
26     choice=2; % 1=snapshot; 2=mean
27     %%%%%%%%%
28     if (ocmipornot==1), bathynum=20; end
29     atlanornot=0; %1=turn fields to atlantic centric
30     %
31     if (bathynum==3), direc=['/data25/stephd/Newdata/input_thirdtry']; end
32     if (bathynum==9),
33     direc=['/data25/stephd/Newdata/input_thirdtry_nomed']; end
34     if (bathynum==10),
35     direc=['/data25/stephd/Nesdata/input_tenthtry_ts']; end
36     if (bathynum==11),
37     direc=['/data25/stephd/Nesdata/input_tenthtry_noarctic']; end
38     if (bathynum==12),
39     direc=['/data25/stephd/Nesdata/input_tenthtry_ts_ideal_sill2']; end
40     if (bathynum==13),
41     direc=['/data25/stephd/Nesdata/input_tenthtry_ts_ideal_sill2_var2']; end
42     if (bathynum==14),
43     direc=['/data25/stephd/Nesdata/input_tenthtry_ts_newideal_new']; end
44     if (ocmipornot==1), direc=['/data25/stephd/c44/ocmip/data_ocmip']; end
45     if (bathynum==12 | bathynum==13 | bathynum==14),
46     dzfile=['dz_ideal'];
47     else
48     dzfile=['dz_vero'];
49     end % if
50     %%%%%%%%%%%%%%%%%%% END CHOICES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51     %
52     % some parameters
53     dtday=dt/(60*60*24); degm=112*1e3; degrad=pi/180; kgm=1;
54     %%%%%%%%%%%%%%%%%%%%%%%% GRID %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55     % horizontal grid
56     if ((bathynum>0 & bathynum<13) | bathynum==14),
57     dx=4; dy=4;
58     latg=[-90,-84:dy:84,90]; longg=[0:dx:360];
59     latc=[-87,-82:dy:82,87]; longc=2:dx:358;
60     end
61     if (bathynum==13),
62     dx=4; dy=4;
63     latg=[-90,-84:dy:84,90]; latc=[-87,-82:dy:82,87];
64     i1=24; i2=136; i3=300;
65     longg=[0:dx:i1 i1+3.6 i1+6.8 i1+9.6 ...
66     i1+12 i1+14 i1+16 ...
67     i1+18 i1+20 i1+22.4 ...
68     i1+25.2 i1+28.4 i1+32 ...
69     i1+36:dx:i2 i2+3.6 i2+6.8 i2+9.6 ...
70     i2+12 i2+14 i2+16 ...
71     i2+18 i2+20 i2+22.4 ...
72     i2+25.2 i2+28.4 i2+32 ...
73     i2+36:dx:i3 i3+3.6 i3+6.8 i3+9.6 ...
74     i3+12 i3+14 i3+16 ...
75     i3+18 i3+20 i3+22.4 ...
76     i3+25.2 i3+28.4 i3+32 ...
77     i3+36:4:360];
78     for i=1:length(longg)-1
79     longc(i)=longg(i)+( longg(i+1)-longg(i) ) /2;
80     end % for i
81     end % if
82     if (bathynum==20),
83     dx=2.8125; dy=2.8125;
84     latc=[-88.59375:2.8125:88.59375]; longc=1.40625:2.8125:358.59375;
85     latg=[-90:2.8125:90]; longg=0:2.8125:360;
86     end % if
87     latu=latc; longu=longg; latv=latg; longv=longc;
88     nx=length(longc); ny=length(latc); nx1=nx-1; ny1=ny-1;
89     %
90     % find length between grid points and center points
91     for i=1:nx, for j=1:ny
92     dyg(i,j)=(latg(j+1)-latg(j))*degm;
93     if (j<ny), dyc(i,j)=(latc(j+1)-latc(j))*degm; end
94     end, end % for i,j
95     dyc(:,ny)=dyg(:,ny)/2+dyg(:,1)/2;
96     for i=1:nx, for j=1:ny
97     dxg(i,j)=(longg(i+1)-longg(i))*cos(latg(j)*degrad)*degm;
98     if (i<nx) dxc(i,j)=(longc(i+1)-longc(i))*cos(latc(j)*degrad)*degm; end
99     end, end % for i,j
100     dxc(nx,:)=dxg(nx1,:)/2+dxg(1,:)/2;
101     % dxcg at center of grid
102     for i=1:nx, for j=1:ny
103     dxcg(i,j)=(longg(i+1)-longg(i))*cos(latc(j)*degrad)*degm;
104     end, end % for i j
105     % vertical resolution
106     load ([dzfile,'.dat']); eval(['dz=',dzfile,';']); dz=dz';
107     if (bathynum==14), clear tmp, tmp=dz(1:length(dz)-1); clear dz, dz=tmp; end
108     depth=cumsum(dz); nz=length(dz); depthm(1)=0+dz(1)/2;
109     for k=2:nz, depthm(k)=depthm(k-1)+dz(k-1)/2+dz(k)/2; end
110     %%%%%%
111     % face areas
112     for k=1:nz,
113     area_xyc(:,:,k)=dxcg(:,:).*dyg(:,:);
114     area_xzg(:,:,k)=dxg(:,:)*dz(k);
115     area_yzg(:,:,k)=dyg(:,:)*dz(k);
116     volc(:,:,k)=area_xyc(:,:,k)*dz(k);
117     end % for k
118     % now that areas are found, convert all nx+1,ny+1 arrays to nx,ny
119     clear tmp, tmp=longg(1:nx); clear longg, longg=tmp;
120     clear tmp, tmp=longu(1:nx); clear longu, longu=tmp;
121     clear tmp, tmp=dxg(1:nx,:); clear dxg, dxg=tmp;
122     clear tmp, tmp=latg(1:ny); clear latg, latg=tmp;
123     clear tmp, tmp=latv(1:ny); clear latv, latv=tmp;
124     clear tmp, tmp=dyg(:,1:ny); clear dyg, dyg=tmp;
125     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
126     %%%%
127     if (atlanornot==1)
128     olongg=longg; olongc=longc; olongu=longu; olongv=longv;
129     longg=turnlong(olongg,longg); longc=turnlong(olongc,longc);
130     longu=turnlong(olongu,longu); longv=turnlong(olongv,longv);
131     dyg=turn2(olongg,dyg); dyc=turn2(olongc,dyc);
132     dxg=turn2(olongg,dxg); dxc=turn2(olongc,dxc);
133     area_xyc=turn3(olongg,area_xyc); area_xzg=turn3(olongg,area_xzg);
134     area_yzg=turn3(olongg,area_yzg); volc=turn3(olongg,volc);
135     end % if atlanornot 1
136     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
137     % read bathymetry and masks
138     if (ocmipornot==1),
139     fid=fopen([direc,'/bathymetry.bin'],'r','b');
140     else
141     fid=fopen([direc,'/bathy.bin'],'r','b');
142     end % if
143     bathy=fread(fid,'float32'); bathy=reshape(bathy,nx,ny);
144     % masks
145     fid=fopen([direc,'/hFacC.bin'],'r','b');
146     if (bathynum==11), fid=fopen([direc,'/pmask.bin'],'r','b'); end
147     pmask=fread(fid,'float32'); pmask=reshape(pmask,nx,ny,nz);
148     % find u,v mask
149     for k=1:nz, for i=1:nx, for j=2:ny,
150     vmask(i,j,k)=min(pmask(i,j,k),pmask(i,j-1,k));
151     end, end, end % for k,i,j
152     for k=1:nz, for i=2:nx, for j=1:ny,
153     umask(i,j,k)=min(pmask(i,j,k),pmask(i-1,j,k));
154     if (i==2), umask(1,j,k)=min(pmask(1,j,k),pmask(nx,j,k)); end
155     end, end, end % for k,i,j
156     % masks for each ocean (down to antarctica)
157     fid=fopen([direc,'/tocean.bin'],'r','b');
158     tocean=fread(fid,'float32'); tocean=reshape(tocean,nx,ny);
159     % masks for each ocean, including southern ocean
160     fid=fopen([direc,'/cocean.bin'],'r','b');
161     cocean=fread(fid,'float32'); cocean=reshape(cocean,nx,ny);
162     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
163     % original T and S fields (ie Levitus)
164     if (ocmipornot==1),
165     fid=fopen([direc,'/lev_clim_temp.bin'],'r','b'); monnum=1;
166     theta=fread(fid,'float32');
167     theta=reshape(theta,nx,ny,nz); theta_clim=theta;
168     else
169     fid=fopen([direc,'/theta.bin'],'r','b'); monnum=12;
170     theta=fread(fid,'float32');
171     theta=reshape(theta,nx,ny,nz,monnum); theta_clim=mean(theta,4);
172     end
173     % clear fi, fi=find(isnan(theta_clim)==1); theta_clim(fi)=900;
174     % clear fi, fi=find(theta_clim==0); theta_clim(fi)=900;
175     %
176     if (ocmipornot==1),
177     fid=fopen([direc,'/lev_clim_salt.bin'],'r','b'); monnum=1;
178     salt=fread(fid,'float32');
179     salt=reshape(salt,nx,ny,nz); salt_clim=salt;
180     else
181     fid=fopen([direc,'/sal.bin'],'r','b'); monnum=12;
182     salt=fread(fid,'float32');
183     salt=reshape(salt,nx,ny,nz,monnum); salt_clim=mean(salt,4);
184     end
185     %%%%%
186     if (atlanornot==1),
187     bathy=turn2(olongc,bathy); pmask=turn3(olongc,pmask);
188     umask=turn3(olongu,umask); vmask=turn3(olongv,vmask);
189     tocean=turn2(olongc,tocean); cocean=turn2(olongc,cocean);
190     theta_clim=turn3(olongc,theta_clim);
191     salt_clim=turn3(olongc,salt_clim);
192     end % if atlanornot
193     %%%%%
194     % land
195     filand=find(pmask==0); filand1=find(pmask(:,:,1)==0);
196     if (kn>1); depthkn1=depth(kn-1); else; depthkn1=0; end
197     % salt_clim(filand)=NaN; theta_clim(filand)=NaN;
198     % time
199     time=str2num(timestr)*dt; timeyear=time/(yeard*24*60*60);
200     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
201     %%%%%%%%%%%%%%%%%%%
202     % SETUP PLOTS %
203     %%%%%%%%%%%%%%%%%%%
204     figure(1), hold off, subplot(1,1,1), pcolor(longc,latc,bathy');
205     colorbar, title([titstr,': bathymetry']);
206     figure(2), hold off, subplot(1,1,1), pcolor(longc,latc,tocean');
207     title([titstr,': Ocean Masks']);
208     figure(3), hold off, subplot(1,1,1), subplot(3,1,1)
209     pcolor(longc,latc,theta_clim(:,:,kn)'); shading flat, colorbar,
210     title([titstr,': Levitus Theta, ',num2str(depthm(kn)),'m'])
211     figure(4), hold off, subplot(1,1,1), subplot(3,1,1)
212     pcolor(longc,latc,salt_clim(:,:,kn)'); shading flat, colorbar,
213     title([titstr,': Levitus Salinity, ',num2str(depthm(kn)),'m'])
214     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
215     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
216     % initialize summary plots %
217     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
218     fgn=115;
219     figure(fgn), hold off, subplot(1,1,1);
220     figure(fgn+1), hold off, subplot(1,1,1);
221     figure(fgn+2), hold off, subplot(1,1,1);
222     figure(fgn+3), hold off, subplot(1,1,1);
223     % global
224     ax1=1; ax2=360; ax3=-90; ax4=90; inton= 20; vmult=150;
225     if (atlanornot==1), ax1=-180; ax2=180; end
226     % north atlantic for subplt (note could also be another level);
227     ax1_sp=280; ax2_sp=360; ax3_sp=-10; ax4_sp=90; inton_sp=5; vmult_sp=200;
228     if (atlanornot==1), ax1_sp=-100; ax2_sp=20; end
229     %%%%%%%%%%%%
230     %
231     ['finished set up, now starting reading fields'],
232     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
233     % MAIN FIELDS %
234     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
235     if (choice==1), maxfld=6; else, maxfld=7; end
236     if (choice>0),
237     figure(10), hold off, subplot(1,1,1),
238     for ifld=1:maxfld,
239     if (ifld==1), if (choice==1), str=['T']; else, str=['Ttave']; end
240     lat=latc; long=longc; cax1=-2; cax2=28;
241     cint1=-20; cint2=50; cint=5; end
242     if (ifld==2), if (choice==1), str=['S']; else, str=['Stave']; end
243     lat=latc; long=longc; cax1=28; cax2=40;
244     cint1=20; cint2=50; cint=.5; end
245     if (ifld==3), if (choice==1), str=['U']; else, str=['uVeltave']; end
246     lat=latu; long=longu; end
247     if (ifld==4), if (choice==1), str=['V']; else, str=['vVeltave']; end
248     lat=latv; long=longv; end
249     if (ifld==5), if (choice==1), str=['W']; else, str=['wVeltave']; end
250     lat=latc; long=longc; cax1=-9e-6; cax2=6e-6; cint=3e-6; end
251     if (ifld==6), if (choice==1), str=['Eta']; else, str=['ETAtave']; end
252     lat=latc; long=longc; cax1=-2; cax2=1.2;
253     cint1=-5; cint2=5; cint=.25; end
254     if (ifld==7 & choice>1 ), str=['Convtave']; lat=latc; long=longc;
255     cax1=0; cax2=1;
256     cint1=0; cint2=1; cint=.2; end
257     % read fields
258     fid=fopen([run1str,'/',str,'.',timestr,'.data'],'r','b');
259     clear field, field=fread(fid,'float32');
260     if (ifld==6),
261     field=reshape(field,nx,ny);
262     if (atlanornot==1), field=turn2(olongc,field); end
263     else
264     field=reshape(field,nx,ny,nz);
265     if (atlanornot==1),
266     if (ifld==1 |ifld==2| ifld==5 | ifld==6); olong=olongc; end
267     if (ifld==3), olong=olongu; end, if (ifld==4), olong=olongv; end
268     field=turn3(olong,field); end
269     end % ifld
270     %%%
271     if (ifld==1), temp=field; clear filand_temp,
272     filand_temp=find(temp==0);
273     filand1_temp=find(temp(:,:,1)==0);
274     temp_pmask=ones(size(temp)); temp_pmask(filand_temp)=0; end
275     if (ifld==2), salt=field; end
276     if (ifld==3), uvel=field; end
277     if (ifld==4), vvel=field; end
278     if (ifld==5), wvel=field; end
279     if (ifld==6), eta=field; end
280     if (ifld==7), conv=field; end
281     %%%%%%%%%%%%%%%%%%%%
282     % calculate max/min of each field
283     clear tmp, tmp=field;
284     if (ifld~=6), tmp(filand_temp)=NaN; else, tmp(filand1_temp)=NaN; end
285     totmax(ifld)=max(max(max(tmp)));
286     totmin(ifld)=min(min(min(tmp)));
287     % calculate mean
288     if (ifld~=6), volume=volc; volume(filand_temp)=0;
289     else, volume=volc(:,:,1); volume(filand1_temp)=0; end
290     tmp=field.*volume; totmean(ifld)=sum(sum(sum(tmp)))/(sum(sum(sum(tmp))));
291     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
292     % calculate zonal transports
293     if (ifld==3),
294     % for ACC
295     % j1=2; j2=10; i1=73; if (atlanornot==1), i1=i1-45; end
296     lat1=-83; lat2=-50; if (atlanornot==0), long1=290; else; long1=-70; end
297     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
298     clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
299     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
300     transp_accd=sum(sum(field(i1,j1:j2,:).*area_yzg(i1,j1:j2,:) ));
301     %
302     % j1=2; j2=15; i1=6; if (atlanornot==1), i1=i1+45; end
303     lat1=-83; lat2=-30; long1=22;
304     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
305     clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
306     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
307     transp_accg=sum(sum(field(i1,j1:j2,:).*area_yzg(i1,j1:j2,:) ));
308     %
309     % j1=2; j2=13; i1=37; if (atlanornot==1), i1=i1+45; end
310     lat1=-83; lat2=-38; long1=146;
311     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
312     clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
313     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
314     transp_accf=sum(sum(field(i1,j1:j2,:).*area_yzg(i1,j1:j2,:) ));
315     % for mediterranean
316     % j1=31; j2=33; i1=nx; if (atlanornot==1), i1=i1-45; end
317     lat1= 34; lat2= 42; long1=178;
318     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
319     clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
320     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
321     clear tmp2 fi, tmp2=field(i1,j1:j2,:).*area_xzg(i1,j1:j2,:);
322     fi=find(tmp2<0); transp_med=sum(tmp2(fi));
323     % indonesian throughflow
324     if (bathynum==12 | bathynum==13),
325     lat1=-6; lat2=6; long1=142;
326     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
327     clear fi, fi=find(latc>lat2-dy/2 & latc<lat2+dy/2); j2=fi(1);
328     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
329     clear tmp2 fi, tmp2=field(i1,j1:j2,:).*area_xzg(i1,j1:j2,:);
330     fi=find(tmp2<0); transp_indtf=sum(tmp2(fi));
331     end % if bathynum 12
332     end % if
333     % meridional transports
334     if (ifld==4),
335     % equatorial atlantic
336     % j1=22; i1=80; i2=85; if (atlanornot==1), i1=i1-45; i2=i2-45; end
337     lat1=-2; if (atlanornot==0), long1=318; long2=338; else
338     long1=-42; long2=-22; end
339     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
340     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
341     clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
342     transp_eqatl=sum(sum(field(i1:i2,j1,1:8).*area_xzg(i1:i2,j1,1:8) ));
343     % gulf stream
344     % j1=30; i1=70; i2=75; if (atlanornot==1),i1=i1-45; i2=i2-45; end
345     lat1=30; if (atlanornot==0), long1=278; long2=298; else
346     long1=-82; long2=-62; end
347     if (bathynum==12 | bathynum==13), if (atlanornot==0), long1=310; long2=318;
348     else, long1=-50; long2=-32; end, end
349     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
350     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
351     clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
352     transp_gstre=sum(sum(field(i1:i2,j1,1:8).*area_xzg(i1:i2,j1,1:8) ));
353     % kuroshio
354     % j1=29; i1=31; i2=36; if (atlanornot==1), i1=i1+45; i2=i2+45; end
355     lat1=26; long1=122; long2=142;
356     if (bathynum==12 | bathynum==13), long1=142; long2=162; end
357     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
358     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
359     clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
360     transp_kuros=sum(sum(field(i1:i2,j1,1:8).*area_xzg(i1:i2,j1,1:8) ));
361     % indonesian throughflow
362     % j1=22; i1=30; i2=35; if (atlanornot==1), i1=i1+45; i2=i2+45; end
363     if (bathynum~=12 & bathynum~=13)
364     lat1=-2; long1=118; long2=138;
365     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
366     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
367     clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
368     transp_indtf=sum(sum(field(i1:i2,j1,:).*area_xzg(i1:i2,j1,:) ));
369     end % if bathynum 12
370     % arctic exchange
371     % j1=39; i1=1; i2=7; if (atlanornot==1), i1=83-45; i2=7+45; end
372     lat1=66; long2=26; if (atlanornot==0), long1=2; else, long1=-30; end
373     if (bathynum==12 | bathynum==13),
374     if (atlanornot==0), long1=2; else, long1=-50; end, end
375     clear fi, fi=find(latc>lat1-dy/2 & latc<lat1+dy/2); j1=fi(1);
376     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
377     clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2); i2=fi(1);
378     clear tmp, tmp=field.*area_xzg;
379     clear tmp2 fi, tmp2=tmp(i1:i2,j1,:); fi=find(tmp2>0);
380     transp_arcex=sum(tmp2(fi));
381     if (atlanornot==0), long1=330; long2=360;
382     if (bathynum==12 | bathynum==13), long1=310; end
383     clear fi, fi=find(longc>long1-dx/2 & longc<long1+dx/2); i1=fi(1);
384     clear fi, fi=find(longc>long2-dx/2 & longc<long2+dx/2);
385     if (isempty(fi)==1), i2=nx; else, i2=fi(1); end
386     clear tmp2 fi, tmp2=tmp(i1:i2,j1,:); fi=find(tmp2>0);
387     transp_arcex= transp_arcex+sum(tmp2(fi));
388     end % if
389     end % if
390     %%%%%%%%%%%%%%%%%%%%
391     % XY FIGURES %
392     %%%%%%%%%%%%%%%%%%%%
393     clear tmp, if (ifld==6); tmp=field; else, tmp(:,:)=field(:,:,kn); end
394     clear fi, fi=find(pmask(:,:,kn)~=1); tmp(fi)=NaN;
395     if (ifld<7),
396     figure(10), subplot(3,2,ifld), pcolor(long,lat,tmp');
397     else
398     figure(11), hold off, subplot(1,1,1), pcolor(long,lat,tmp');
399     end
400     shading flat, colorbar
401     hold on, [c,han]=contour(long,lat,tmp',[-1000 0 1000],'k-');
402     set(han,'LineWidth',1.5);
403     if (ifld==6), title([direc,': ',str]);
404     else, title([str,' at ',num2str(depthm(kn)),'m']); end
405     if (ifld==5), title([str,' at ',num2str(depthkn1),'m']); end
406     %
407     % plot difference from Levitus (for salt and temp)
408     if (ifld==1 | ifld==2), clear tmp2
409     if (ifld==1), tmp2=tmp-theta_clim(:,:,kn); end
410     if (ifld==2), tmp2=tmp-salt_clim(:,:,kn); end
411     figure(30+ifld), hold off, subplot(1,1,1),
412     subplot(3,1,1), pcolor(long,lat,tmp2'); shading flat, colorbar
413     hold on, contour(long,lat,tmp2',[-1000 0 1000],'k-');
414     title([titstr,' minus ',str,'_{lev}, at ',num2str(depthm(kn)),'m']);
415     end % if
416     % **********************PLOT FGN*******************************************
417     if (ifld==1 | ifld==2), figure(fgn+2), subplot(3,2,ifld),
418     pcolor(long,lat,tmp'); caxis([cax1 cax2]); shading flat, colorbar
419     hold on, [c,han]=contour(long,lat,tmp',[-1000 0 1000],'k-');
420     set(han,'LineWidth',1.5);
421     contour(long,lat,tmp',[cint1:cint:cint2],'k-');
422     title([str,' at ',num2str(depthm(kn)),'m']);
423     axis([ax1 ax2 ax3 ax4]);
424     mintr=min(min(tmp')); maxtr=max(max(tmp'));
425     han=text(ax1,ax3-inton,[num2str(mintr),' to ',num2str(maxtr)]);
426     set(han,'FontSize',6);
427     % title for page fgn+2
428     if (ifld==2),
429     han=text(ax1-8*inton,ax4+3*inton,[titstr]); set(han,'FontSize',12);
430     text(ax1-4*inton,ax4+2*inton,['year ', num2str(timeyear)]);
431     end
432     end
433     if (ifld==6)
434     figure(fgn), subplot(3,2,6),pcolor(long,lat,tmp');
435     caxis([cax1 cax2]); shading flat, colorbar
436     hold on, [c,han]=contour(long,lat,tmp',[-1000 0 1000],'k-');
437     set(han,'LineWidth',1.5);
438     contour(long,lat,tmp',[cint1:cint:cint2],'k-');
439     title([str]); axis([ax1 ax2 ax3 ax4]);
440     mintr=min(min(tmp')); maxtr=max(max(tmp'));
441     han=text(ax1,ax3-inton,[num2str(mintr),' to ',num2str(maxtr)]);
442     set(han,'FontSize',6);
443     end % if ifld
444     % subplot
445     if (ifld==1 | ifld==2 | ifld==5);
446     inum=ifld+2; if (ifld==5), inum=2; end
447     clear tmp_sp, tmp_sp=field(:,:,knsp);
448     clear fi, fi=find(pmask(:,:,knsp)~=1); tmp_sp(fi)=NaN;
449     figure(fgn+3), subplot(2,2,inum), pcolor(long,lat,tmp_sp');
450     shading flat, colorbar
451     hold on, [c,han]=contour(long,lat,tmp',[-1000 0 1000],'k-');
452     set(han,'LineWidth',1.5);
453     % title for page
454     if (ifld==5),
455     han=text(ax1_sp-8*inton_sp,ax4_sp+3*inton_sp,[titstr]);
456     set(han,'FontSize',12);
457     text(ax1_sp-4*inton_sp,ax4_sp+2*inton_sp,['year ', num2str(timeyear)]);
458     end
459     if (ifld==5), title([str,' at ',num2str(depthm(knsp)),'m']);
460     else, title([str,' at ',num2str(depth(knsp)),'m']); end
461     axis([ax1_sp ax2_sp ax3_sp ax4_sp]);
462     mintr=min(min(tmp')); maxtr=max(max(tmp'));
463     han=text(ax1_sp,ax3_sp-inton,[num2str(mintr),' to ',num2str(maxtr)]);
464     set(han,'FontSize',6);
465     end % if ifld
466     % ***************************************************************************
467     % plot velocity vectors
468     if (ifld==4),
469     for ip=1:3,
470     if (ip==1), figure(12), hold off, subplot(1,1,1); knn=kn; vmultn=vmult; end
471     if (ip==2), figure(fgn), subplot(2,1,1); knn=kn; vmultn=vmult;end
472     if (ip==3), figure(fgn+3), subplot(2,2,1); knn=knsp; vmultn=vmult_sp; end
473     clear tmp, tmp=temp_pmask(:,:,knn);
474     contour(longc,latc,tmp',1,'k-'); hold on % draw land
475     if (nx>105), nv=2; else, nv=1; end
476     for i=1:nv:nx, for j=1:nv:ny1,
477     if (i<nx1), ix=i+1; else ix=1; end
478     u(i,j)=(uvel(i,j,knn)+uvel(ix,j,knn))/2;
479     v(i,j)=(vvel(i,j,knn)+vvel(i,j+1,knn))/2;
480     spd(i,j)=sqrt( u(i,j)^2 + v(i,j)^2 );
481     if (spd(i,j)~=0 )
482     if (u(i,j)~=0),
483     dir=(180*atan(v(i,j)/u(i,j))/pi) ;
484     if (u(i,j) < 0), dir = 180 + dir ; end,
485     else
486     if (v(i,j)>0) dir=90; else, dir=270; end
487     end
488     if (spd(i,j)>0),
489     velvct(longc(i),latc(j),spd(i,j),dir,vmultn,'k'),
490     end,
491     end % if
492     end, end %for i j
493     title([titstr,' Velocity at ',num2str(depthm(knn)),'m']);
494     velvct(50,-75,.2,0,vmultn,'k');
495     text(50,-70,'20cm/s');
496     %*********************** PLOT FGN *********************************************
497     if (ip==2),
498     axis([ax1 ax2 ax3 ax4]);
499     title(['Velocity at ',num2str(depthm(kn)),'m']);
500     han=text(ax1+(ax2-ax1)/8,ax3-2*inton,[titstr]); set(han,'FontSize',18);
501     han=text(ax1+(ax2-ax1)/8,ax3-3*inton,['year ', num2str(timeyear)]);
502     set(han,'FontSize',14);
503     text(ax1+(ax2-ax1)/8,ax3-3.5*inton,['dt=',num2str(dtday),' day']);
504     maxtr=max(max(spd(:,:))); han=text(ax1,ax3-inton,['max ',num2str(maxtr)]);
505     set(han,'FontSize',6);
506     c=axis; inton1=c(2)-c(1); inton2=c(4)-c(3);
507     % write transport numbers
508     text(ax1+2*(ax2-ax1)/4,ax3-1.2*inton,'TRANSPORTS (Sv)')
509     text(ax1+2*(ax2-ax1)/4,ax3-1.5*inton,'---------------')
510     text(ax1+2*(ax2-ax1)/4,ax3-2*inton, ...
511     ['ACC: ',num2str(transp_accd/1e6)]);
512     text(ax1+2*(ax2-ax1)/4,ax3-2.5*inton, ...
513     ['Arctic Exchange: ',num2str(transp_arcex/1e6)]);
514     text(ax1+2*(ax2-ax1)/4,ax3-3*inton, ...
515     ['Equatorial Crossing (Atlantic): ',num2str(transp_eqatl/1e6)]);
516     text(ax1+2*(ax2-ax1)/4,ax3-3.5*inton, ...
517     ['Gulf Stream (24N): ',num2str(transp_gstre/1e6)]);
518     text(ax1+2*(ax2-ax1)/4,ax3-4*inton, ...
519     ['Kuroshio (20N): ',num2str(transp_kuros/1e6)]);
520     text(ax1+2*(ax2-ax1)/4,ax3-4.5*inton, ...
521     ['Indonesian Throughflow: ',num2str(transp_indtf/1e6)]);
522     text(ax1+2*(ax2-ax1)/4,ax3-5*inton, ...
523     ['Mediterranean Exchange: ',num2str(transp_med/1e6)]);
524     end % if %ip 2
525     if (ip==3),
526     axis([ax1_sp ax2_sp ax3_sp ax4_sp]);
527     title(['Velocity at ',num2str(depthm(knsp)),'m']);
528     maxtr=max(max(spd(:,:)));
529     han=text(ax1_sp,ax3_sp-inton_sp,['max ',num2str(maxtr)]);
530     set(han,'FontSize',6);
531     end % if
532     %*******************************************************************************
533     end % for ip
534     end % if ifld 4
535     %
536     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
537     % meridional transects
538     if (ifld<6 | ifld==7),
539     num=0; if (ocmipornot==1), nocean=5; else, nocean=6; end
540     if (bathynum==11), nocean=5; end
541     clear tmp2 tmp3 tmp4
542     tmp2=field.*volc; tmp2(filand_temp)=0;
543     field(filand_temp)=NaN;
544     theta_clim(filand_temp)=NaN;
545     salt_clim(filand_temp)=NaN;
546     for iocean=1:nocean,
547     if (iocean==1), ioc=0; tstr=['Global']; end
548     if (iocean==2), ioc=3; tstr=['Atlantic']; end
549     if (iocean==3), ioc=2; tstr=['Pacific']; end
550     if (iocean==4), ioc=1; tstr=['Indian']; end
551     if (iocean==5), ioc=4; tstr=['Southern']; end
552     if (iocean==6), ioc=5; tstr=['Arctic']; end
553     if (iocean==7), ioc=6; tstr=['Med']; end
554     for j=1:ny, for k=1:nz,
555     clear fi,
556     if (iocean==1), fi=find(tocean(:,j)>ioc & temp_pmask(:,j,k)==1); end;
557     if (iocean>1 & iocean <5),
558     fi=find(tocean(:,j)==ioc & temp_pmask(:,j,k)==1); end
559     if (iocean>4) fi=find(cocean(:,j)==ioc & temp_pmask(:,j,k)==1); end
560     if (isempty(fi)==0),
561     trfield(j,k)=sum(field(fi,j,k).*volc(fi,j,k))/sum(volc(fi,j,k));
562     else,
563     trfield(j,k)=NaN;
564     end
565     end, end % for j k
566     figure(20+ifld),
567     if (iocean<5), subplot(3,2,iocean); else, subplot(3,3,iocean+2); end
568     pcolor(lat,-depth,trfield'), shading flat, colorbar,
569     hold on, c=contour(lat,-depth,trfield','k-'); clabel(c);
570     [c,han]=contour(lat,-depth,trfield',[-1000 0 1000],'k-');
571     set(han,'LineWidth',1.5);
572     if (iocean<5), axis([latg(1) latg(ny) -5300 0]); end
573     if (iocean==5), axis([latg(1) 0 -5300 0]); end
574     if (iocean==6), axis([60 latg(ny) -5300 0]); end
575     if (iocean==7), axis([30 50 -5300 0]); end
576     title([tstr,': ',str]);
577     % plot Levitus and differences from Levitus
578     if (ifld==1 | ifld==2), clear tmp2 lev,
579     for j=1:ny, for k=1:nz,
580     clear fi,
581     if (iocean==1), fi=find(tocean(:,j)>ioc & temp_pmask(:,j,k)==1); end;
582     if (iocean>1 & iocean <5),
583     fi=find(tocean(:,j)==ioc & temp_pmask(:,j,k)==1); end
584     if (iocean>4) fi=find(cocean(:,j)==ioc & temp_pmask(:,j,k)==1); end
585     if (isempty(fi)==0),
586     if (ifld==1),
587     lev(j,k)=sum(theta_clim(fi,j,k).*volc(fi,j,k))/sum(volc(fi,j,k));
588     else
589     lev(j,k)=sum(salt_clim(fi,j,k).*volc(fi,j,k))/sum(volc(fi,j,k));
590     end
591     else,
592     lev(j,k)=NaN;
593     end
594     end, end % for j k
595     figure(2+ifld), subplot(3,4,iocean+4);
596     pcolor(lat,-depth,lev'), shading flat, colorbar,
597     hold on, [c,han]=contour(lat,-depth,lev',[-1000 0 1000],'k-');
598     set(han,'LineWidth',1.5);
599     if (iocean<5), axis([latg(1) latg(ny) -5300 0]); end
600     if (iocean==5), axis([latg(1) 0 -5300 0]); end
601     if (iocean==6), axis([60 latg(ny) -5300 0]); end
602     if (iocean==7), axis([30 50 -5300 0]); end
603     title([tstr,': ',str,'_{lev}']);
604     %
605     figure(30+ifld), subplot(3,4,iocean+4);
606     pcolor(lat,-depth,(trfield-lev)'), shading flat, colorbar,
607     hold on, c=contour(lat,-depth,(trfield-lev)','k-'); clabel(c);
608     [c,han]=contour(lat,-depth,(trfield-lev)',[-1000 0 1000],'k-');
609     set(han,'LineWidth',1.5);
610     if (iocean<5), axis([latg(1) latg(ny) -5300 0]); end
611     if (iocean==5), axis([latg(1) 0 -5300 0]); end
612     if (iocean==6), axis([60 latg(ny) -5300 0]); end
613     if (iocean==7), axis([30 50 -5300 0]); end
614     title([tstr]);
615     end % if
616     % find single profile of differences
617     if (ifld==1 | ifld==2),
618     clear prof_field profclim
619     for k=1:nz, clear tmp3 tmp4 fi,
620     field(filand_temp)=NaN;
621     tmp3=volc(:,:,k); tmp4=field(:,:,k).*volc(:,:,k);
622     if (ifld==1), tmp5=theta_clim(:,:,k).*volc(:,:,k); end
623     if (ifld==2), tmp5=salt_clim(:,:,k).*volc(:,:,k); end
624     if (iocean==1), fi=find(tocean>ioc & isnan(tmp4)==0); else,
625     fi=find(cocean==ioc & isnan(tmp4)==0); end
626     prof_field(ifld,iocean,k)=sum(tmp4(fi))/sum(tmp3(fi));
627     prof_clim(ifld,iocean,k)=sum(tmp5(fi))/sum(tmp3(fi));
628     if (k==nz),
629     figure(32+ifld), if (iocean==1), hold off, subplot(1,1,1); end
630     subplot(2,7,iocean), clear tmp, tmp=squeeze(prof_clim(ifld,iocean,:));
631     plot(tmp,-depthm,'k--');
632     hold on, clear tmp, tmp=squeeze(prof_field(ifld,iocean,:));
633     plot(tmp,-depthm,'r-');
634     title([tstr]);
635     clear tmp,
636     tmp=squeeze(prof_field(ifld,iocean,:)-prof_clim(ifld,iocean,:));
637     subplot(2,7,iocean+7), plot(tmp,-depthm,'k-'); hold on,
638     plot([0 0],[-5300 0],':');
639     %****************************** PLOTFGN ******************************
640     if (iocean==1),
641     figure(fgn+2), subplot(3,2,4+ifld);
642     plot(tmp,-depthm,'k-'); hold on,
643     plot([0 0],[-5300 0],':');
644     title([tstr,': ',str,'-',str,'_{lev}']);
645     end % if iocean
646     %*********************************************************************
647     end % if k nz
648     end, % for k
649     end % if ifld
650    
651     %%******************PLOT FGN******************************
652     if (iocean==1 & (ifld==1 | ifld==2) ),
653     if (ifld==1), figure(fgn+2), pnum=3; end
654     if (ifld==2), figure(fgn+2), pnum=4; end
655     subplot(3,2,pnum), pcolor(lat,-depth',trfield');
656     if (ifld==2), caxis([32 36]); else
657     caxis([cax1 cax2]); end
658     shading flat, colorbar, hold on
659     [c,han]=contour(lat,-depth',trfield',[-1000 0 1000],'k-');
660     set(han,'LineWidth',1.5);
661     contour(lat,-depth',trfield',[cint1:cint:cint2],'k');
662     axis([latg(1) latg(ny) -5300 0]);
663     mintr=min(min(trfield')); maxtr=max(max(trfield'));
664     han=text(latg(1),-5500,[num2str(mintr),' to ',num2str(maxtr)]);
665     set(han,'FontSize',6);
666     title(['Global, zonal ',str]);
667     end % if
668     %***************************************************************
669     end % for iocean
670     end % if ifld
671     end % for ifld
672     end % if choice
673     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
674     % merdional overturning
675     dxelk=dxcg(1,:); dyelk=dyg(:,1);
676     if (choice==1),
677     str=['V'];
678     else,
679     str=['vVeltave'];
680     str2=['GM_Kwy-T'];
681     end
682     field1=vvel;
683     fid=fopen([run1str,'/',str2,'.',timestr,'.data'],'r','b');
684     clear field3, field3=fread(fid,'float32');
685     field3=reshape(field3,nx,ny,nz);
686     if (atlanornot==1), field3=turn3(olongv,field3); end
687     %
688     % calculate bolus velocity
689     psiby=field3*kgm/2;
690     for k=1:nz-1,
691     vbc(:,:,k)=-(psiby(:,:,k)-psiby(:,:,k+1))/dz(k);
692     end % for k
693     vbc(:,:,nz)=(psiby(:,:,nz)-0)/dz(nz);
694     for i=1:nx,
695     for j=1:ny-1,
696     vbg(i,j,:)=(vbc(i,j,:)+vbc(i,j+1,:))/2;
697     end, % for j
698     vbg(i,ny,:)=(vbc(i,ny,:)+vbc(i,1,:))/2;
699     end % for i
700     field3=vbg;
701     if (atlanornot==0), splocean=9; else, splocean=9; end
702     for iv=1:2,
703     for k=1:nz, vocean(:,:,k)=cocean.*vmask(:,:,k); end
704     for iocean=1:4,
705     clear dxtab field
706     if (iocean==1), ioc=0; tstr=['Global']; cax1=-30; cax2=30; end
707     if (iocean==2), ioc=3; tstr=['Atlantic']; end
708     if (iocean==3), ioc=2; tstr=['Pacific']; end
709     if (iocean==4), ioc=1; tstr=['Indian']; end
710     clear fi fj,
711     if (iocean==1), [fi]=find(vocean>ioc); else,
712     [fi]=find(vocean==ioc); end
713     mask=zeros(size(vocean)); mask(fi)=1;
714     if (iv==1), clear field2, field2=field1.*mask; end
715     if (iv==2), clear field2, field2=field3.*mask; end
716     clear fi fj fi1 fj1
717     clear minfi maxfi minfj maxfj minfi1 maxfi1 mfi mfj mfi1 mfj1
718     if (ioc~=splocean)
719     [fi,fj]=find(vocean(:,:,1)==ioc);
720     minfi=min(fi); maxfi=max(fi); mfi=maxfi-minfi+1;
721     minfj=min(fj); maxfj=max(fj); mfj=maxfj-minfj+1;
722     field=field2(minfi:maxfi,minfj:maxfj,:);
723     lat=latc(minfj:maxfj);
724     for k=1:nz,
725     dxtab(:,:,k) = dxcg(minfi:maxfi,minfj:maxfj);
726     end % for k
727     dzglob2 = ones(size(1:mfj))' * dz;
728     else
729     hnx=ceil(nx/2);
730     [fi,fj]=find(vocean(hnx:nx,:,1)==ioc);
731     [fi1,fj1]=find(vocean(1:hnx,:,1)==ioc);
732     if (length(fi)>0),
733     minfi=min(fi)+hnx-1; maxfi=max(fi)+hnx-1; mfi=maxfi-minfi+1;
734     minfj=min(fj); maxfj=max(fj); mfj=maxfj-minfj+1;
735     else
736     minfi=NaN; maxfi==NaN; mfi=0;
737     end % if
738     if (length(fi1)>0),
739     minfi1=min(fi1); maxfi1=max(fi1); mfi1=maxfi1-minfi1+1;
740     minfj=min(min(fj1),minfj); maxfj=max(max(fj1),maxfj);
741     mfj=maxfj-minfj+1;
742     end % if
743     field=field2(minfi:maxfi,minfj:maxfj,:);
744     for k=1:nz,
745     dxtab(:,:,k) = dxcg(minfi:maxfi,minfj:maxfj);
746     end % for k
747     if (length(fi1)>0),
748     field(mfi+1:mfi+mfi1,1:mfj,:)= ...
749     field2(minfi1:maxfi1,minfj:maxfj,:);
750     for k=1:nz,
751     dxtab(mfi+1:mfi+mfi1,1:mfj,k) = dxcg(minfi1:maxfi1,minfj:maxfj);
752     end % for k
753     end % if isempty
754     lat=latc(minfj:maxfj);
755     dzglob2 = ones(size(1:mfj))' * dz;
756     end % ioc 3
757     % calculate streamfunction
758     %elk% =============================
759     %elk% 1. Integrate along longitude
760     %elk% =============================
761     % tabx = dxtab.* squeeze( sum(field, 1) );
762     tabx = squeeze( sum(field.*dxtab,1) );
763     %elk% ===================================
764     %elk% 2. Integrate along depth
765     %elk% ===================================
766     psi = NaN*zeros(size(tabx));
767     for iz = nz:-1:1
768     psi(:, iz) = ...
769     -sum( dzglob2(:,nz:-1:iz) .* tabx(:,nz:-1:iz),2 );
770     end
771     psi(find(psi==0))=psi(find(psi==0))*NaN;
772     psi=psi/1e6;
773     %%%%%%%
774     figure(40+iv), if (iocean==1), hold off, subplot(1,1,1), end
775     subplot(2,2,iocean);
776     pcolor(lat,-depth',psi'); caxis([-30 30]), shading flat,
777     colorbar, hold on
778     c=contour(lat,-depth',psi',[-50:5:0],'k--'); clabel(c);
779     hold on
780     c=contour(lat,-depth',psi',[0:5:50],'k-'); clabel(c);
781     [c,han]=contour(lat,-depth',psi',[-1000 0 1000],'k-');
782     set(han,'LineWidth',1.5);
783     axis([latg(1) latg(ny) -5300 0]);
784     if (iv==1), title([tstr,' Eulerian Overturning (Sv)']); end
785     if (iv==2), title([tstr,' Bolus Overturning (Sv)']); end
786     if (iv==1), eval(['psi_eu',num2str(iocean),'=psi;']); end
787     if (iv==2), psi_bo=psi; end
788     %*******************PLOT FGN*****************************
789     if (iocean==1),
790     if (iv==1),ipl=1; end
791     if (iv==2),ipl=3; end
792     figure(fgn+1), subplot(3,2,ipl), pcolor(lat,-depth',psi');
793     caxis([cax1 cax2]), shading flat,
794     colorbar, hold on
795     contour(lat,-depth',psi',[-200:5:0],'k--');
796     hold on
797     contour(lat,-depth',psi',[0:5:200],'k-');
798     [c,han]=contour(lat,-depth',psi',[-1000 0 1000],'k-');
799     set(han,'LineWidth',1.5);
800     axis([latg(1) latg(ny) -5300 0]);
801     mintr=min(min(psi')); maxtr=max(max(psi'));
802     han=text(latg(1),-5500,[num2str(mintr),' to ',num2str(maxtr)]);
803     set(han,'FontSize',6);
804     if (iv==1), title([tstr,' Eulerian Overturning (Sv)']); end
805     if (iv==2), title([tstr,' Bolus Overturning (Sv)']); end
806     end % if icean 1
807     if (iv==2), eval(['psi=psi_bo+psi_eu',num2str(iocean),';']);
808     if (iocean==1), ipl=5; end
809     if (iocean==2), ipl=2; end
810     if (iocean==3), ipl=4; end
811     if (iocean==4), ipl=6; end
812     figure(fgn+1), subplot(3,2,ipl), pcolor(lat,-depth',psi');
813     caxis([cax1 cax2]), shading flat,
814     colorbar, hold on
815     contour(lat,-depth',psi',[-200:5:0],'k--');
816     hold on
817     contour(lat,-depth',psi',[0:5:200],'k-');
818     [c,han]=contour(lat,-depth',psi',[-1000 0 1000],'k-');
819     set(han,'LineWidth',1.5);
820     axis([latg(1) latg(ny) -5300 0]);
821     mintr=min(min(psi')); maxtr=max(max(psi'));
822     han=text(latg(1),-5500,[num2str(mintr),' to ',num2str(maxtr)]);
823     set(han,'FontSize',6);
824     title([tstr,' Total Overturning (Sv)']);
825     end %if iv 2
826    
827     %********************************************************
828     %
829     end % for iocean
830     end % for iv
831     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
832     % horizontal streamfunctions
833     dyelk=dyg(:,2);
834     field1=uvel;
835     % include Bolus transport
836     str2=['GM_Kwx-T'];
837     fid=fopen([run1str,'/',str2,'.',timestr,'.data'],'r','b');
838     clear field3, field3=fread(fid,'float32');
839     field3=reshape(field3,nx,ny,nz);
840     if (atlanornot==1), field3=turn3(olongu,field3); end
841     psibx=field3*kgm/2;
842     for k=1:nz-1,
843     ubc(:,:,k)=-(psibx(:,:,k)-psibx(:,:,k+1))/dz(k);
844     end % for k
845     ubc(:,:,nz)=(psibx(:,:,nz)-0)/dz(nz);
846     for j=1:ny,
847     for i=1:nx-1,
848     ubg(i,j,:)=(ubc(i,j,:)+ubc(i+1,j,:))/2;
849     end, % for i
850     ubg(nx,j,:)=(ubc(nx,j,:)+ubc(1,j,:))/2;
851     end % for j
852     field3=ubg;
853     field1=field1+field3;
854     % end % if
855     clear field2, field2=field1.*umask;
856     % calculate barotropic streamfunction
857     dztab = dyelk* ones(size(1:ny));
858     %elk% =============================
859     %elk% 1. Integrate along depth
860     %elk% =============================
861     tabz=zeros(nx,ny);
862     for k=1:nz, tabz=tabz+field2(:,:,k)*dz(k); end
863     %elk% ===================================
864     %elk% 2. Integrate along meridian
865     %elk% ===================================
866     psi = NaN*zeros(size(tabz));
867     for iy = ny:-1:1
868     psi(:, iy) = ...
869     sum( dyg(:,ny:-1:iy) .* tabz(:,ny:-1:iy),2 );
870     end
871     clear fi, fi=find(umask(:,:,1)==0); psi(fi)=NaN;
872     psi= psi/1e6; %psi(1,:)=psi(nx,:);
873     %%%%%%
874     figure(43); pcolor(longc,latc,psi'); shading flat,
875     colorbar, hold on
876     c=contour(longc,latc,psi',[-200:10:0],'k--'); clabel(c);
877     hold on
878     c=contour(longc,latc,psi',[0:10:200],'k-'); clabel(c);
879     [c,han]=contour(longc,latc,psi',[1000 0 1000],'k-');
880     set(han,'LineWidth',1.5);
881     title(['Barotropic Streamfunction (Sv)']);
882     % *****************PLOTFGN*****************************
883     figure(fgn), subplot(3,2,5), pcolor(longc,latc,psi');
884     caxis([cax1 cax2]); shading flat,
885     colorbar, hold on
886     contour(longc,latc,psi',[-200:10:0],'k--'); hold on
887     contour(longc,latc,psi',[0:10:200],'k-');
888     [c,han]=contour(longc,latc,psi',[1000 0 1000],'k-');
889     set(han,'LineWidth',1.5);
890     title(['Barotropic Streamfunction ']);
891     axis([ax1 ax2 ax3 ax4]);
892     mintr=min(min(psi')); maxtr=max(max(psi'));
893     han=text(ax1,ax3-inton,[num2str(mintr),' to ',num2str(maxtr)]);
894     set(han,'FontSize',6);
895     % *****************************************************
896     % heat flux
897     %*******************************************************************
898     figure(44), hold off, subplot(1,1,1);
899     for iocean=1:4,
900     if (iocean==1), ioc=0; tstr=['Global']; end
901     if (iocean==2), ioc=3; tstr=['Atlantic']; end
902     if (iocean==3), ioc=2; tstr=['Pacific']; end
903     if (iocean==4), ioc=1; tstr=['Indian']; end
904     % from adcroft
905     for j=1:ny;
906     jm1=max([j-1 1]);
907     Fadv(:,j,:)=vvel(:,j,:).*(273.15+(temp(:,jm1,:)+temp(:,j,:))/2);
908     end
909     for j=1:ny;
910     jm1=max([j-1 1]);
911     Fpar(:,j,:)=vbg(:,j,:).*(temp(:,j,:)-temp(:,jm1,:));
912     end
913     % end from adcroft
914     clear fi,
915     tmp=zeros(size(temp_pmask)); for k=1:nz, tmp(:,:,k)=cocean; end
916     if (iocean==1), fi=find(tmp>ioc); else, fi=find(tmp==ioc); end
917     clear tmp2, tmp2=zeros(size(tmp)); tmp2(fi)=1;
918     area=area_xzg.*temp_pmask.*tmp2;
919     for j=1:ny-1, for k=1:nz,
920     trfield(j,k)=sum((Fadv(:,j,k).*area(:,j,k)),1);
921     % trfield(j,k)=sum((Fpar(:,j,k)+Fadv(:,j,k)).*area(:,j,k),1);
922     end, end % for j k
923     %%%%
924     clear fi, fi=find(isnan(trfield)==1); trfield(fi)=0;
925     rho=1024; cp=4200.;
926     tottrfield=sum(trfield,2)*rho*cp/1e15;
927     clear fi, fi=find(tottrfield==0); tottrfield(fi)=NaN;
928     subplot(2,2,iocean); plot(latc,tottrfield,'k-');
929     title([tstr,': meridional heat transport']);
930     xlabel('latitude'); ylabel('PW');
931     end % for iocean
932    
933     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
934     % some statistics
935     clear tmp, tmp=temp-theta_clim;
936     max(max(max(tmp))), min(min(min(tmp))),
937     clear fi, fi=find(tmp<0); tmp(fi)=-tmp(fi);
938     clear fi, fi=find(tmp>0);
939     clear tmp1, tmp1=tmp(fi).*volc(fi); stmp1=sum(tmp1);
940     clear tmp2, tmp2=volc(fi); stmp2=sum(tmp2);
941     tdiff=stmp1/stmp2
942     clear tmp, tmp=salt-salt_clim;
943     max(max(max(tmp))), min(min(min(tmp))),
944     clear fi, fi=find(tmp<0); tmp(fi)=-tmp(fi);
945     clear fi, fi=find(tmp>0);
946     clear tmp1, tmp1=tmp(fi).*volc(fi); stmp1=sum(tmp1);
947     clear tmp2, tmp2=volc(fi); stmp2=sum(tmp2);
948     sdiff=stmp1/stmp2
949     %
950     clear tmp, tmp=temp.*volc;
951     clear fi, fi =find(temp==0); tmp(fi)=NaN;
952     clear fi, fi=find(isnan(tmp)==0);
953     clear tmp2 tmp3, tmp2=tmp(fi); tmp3=volc(fi);
954     theta_mean=sum(tmp2)/sum(tmp3)
955     %
956     clear tmp, tmp=salt.*volc;
957     clear fi, fi =find(salt==0); tmp(fi)=NaN;
958     clear fi, fi=find(isnan(tmp)==0);
959     clear tmp2 tmp3, tmp2=tmp(fi); tmp3=volc(fi);
960     salt_mean=sum(tmp2)/sum(tmp3)
961     %
962     temp_max=max(max(max(temp)))
963     temp_min=min(min(min(temp)))
964     salt_max=max(max(max(salt)))
965     clear tmp fi, tmp=salt; fi=find(salt==0); tmp(fi)=NaN;
966     salt_min=min(min(min(tmp)))

  ViewVC Help
Powered by ViewVC 1.1.22