| 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))) |
Parent Directory
| 
Revision Log
| 
Revision Graph