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))) |