| 1 |
% |
| 2 |
% THIS IS NOT A FUNCTION ! |
| 3 |
% |
| 4 |
% Plot time series of all variables |
| 5 |
% 1 var / figure |
| 6 |
% Outputs recording |
| 7 |
% |
| 8 |
|
| 9 |
clear |
| 10 |
global sla netcdf_domain |
| 11 |
pv_checkpath |
| 12 |
|
| 13 |
netcdf_domain = 'climode'; |
| 14 |
|
| 15 |
% Date series: |
| 16 |
TIME = ['20030201';'20030205';'20030212';'20030219';'20030226';... |
| 17 |
'20030301';'20030305';'20030312']; |
| 18 |
TIME = ['20030201';'20030205';'20030212';'20030219';'20030226';... |
| 19 |
'20030301';'20030305']; |
| 20 |
nt = size(TIME,1); |
| 21 |
|
| 22 |
|
| 23 |
iso = 25.8; % Which sigma-theta surface ? |
| 24 |
getiso = 1; % We compute the isoST by default |
| 25 |
outimg = 'img_tmp'; % Output directory |
| 26 |
prtimg = 1; % Do we record figures as jpg files ? |
| 27 |
|
| 28 |
% Plots type: |
| 29 |
% 1: individual simple maps |
| 30 |
% 2: meridional or zonal sections |
| 31 |
% 3: meridional or zonal sections with Jz vectors at the surface |
| 32 |
% 4: custom it |
| 33 |
pl = 3; |
| 34 |
|
| 35 |
if pl==2|pl==3, getiso = 0; end |
| 36 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 37 |
% Video loop: |
| 38 |
for it = 1 : nt |
| 39 |
|
| 40 |
% Date: |
| 41 |
snapshot = TIME(it,:); |
| 42 |
|
| 43 |
%%%%%%%%%%%%%%%% |
| 44 |
% NETCDF files name: |
| 45 |
filPV = 'splPV'; |
| 46 |
filST = 'SIGMATHETA'; |
| 47 |
filT = 'THETA'; |
| 48 |
filTx = 'TAUX'; |
| 49 |
filTy = 'TAUY'; |
| 50 |
filJFz = 'JFz'; |
| 51 |
filJBz = 'JBz'; |
| 52 |
filQnet = 'Qnet'; |
| 53 |
filQEk = 'QEk'; |
| 54 |
filMLD = 'KPPmld'; |
| 55 |
|
| 56 |
% Path and extension to find them: |
| 57 |
pathname = strcat('netcdf-files',sla); |
| 58 |
ext = 'nc'; |
| 59 |
|
| 60 |
% Load fields: |
| 61 |
disp('load fields...') |
| 62 |
% (I keep proper axis for each variables in case of one day they would be different) |
| 63 |
ferfile = strcat(pathname,sla,snapshot,sla,filPV,'.',netcdf_domain,'.',ext); |
| 64 |
ncQ = netcdf(ferfile,'nowrite'); |
| 65 |
[Qlon Qlat Qdpt] = coordfromnc(ncQ); |
| 66 |
Q = ncQ{4}(:,:,:); clear ncQ ferfile |
| 67 |
[nz ny nx] = size(Q); |
| 68 |
|
| 69 |
ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext); |
| 70 |
ncST = netcdf(ferfile,'nowrite'); |
| 71 |
[STlon STlat STdpt] = coordfromnc(ncST); |
| 72 |
ST = ncST{4}(:,:,:); clear ncST ferfile |
| 73 |
|
| 74 |
ferfile = strcat(pathname,sla,snapshot,sla,filT,'.',netcdf_domain,'.',ext); |
| 75 |
ncT = netcdf(ferfile,'nowrite'); |
| 76 |
[Tlon Tlat Tdpt] = coordfromnc(ncT); |
| 77 |
T = ncT{4}(:,:,:); clear ncT ferfile |
| 78 |
|
| 79 |
ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext); |
| 80 |
ncTx = netcdf(ferfile,'nowrite'); |
| 81 |
[Txlon Txlat Txdpt] = coordfromnc(ncTx); |
| 82 |
Tx = ncTx{4}(1,:,:); clear ncTx ferfile |
| 83 |
ferfile = strcat(pathname,sla,snapshot,sla,filTy,'.',netcdf_domain,'.',ext); |
| 84 |
ncTy = netcdf(ferfile,'nowrite'); |
| 85 |
[Tylon Tylat Tydpt] = coordfromnc(ncTy); |
| 86 |
Ty = ncTy{4}(1,:,:); clear ncTy ferfile |
| 87 |
|
| 88 |
ferfile = strcat(pathname,sla,snapshot,sla,filJFz,'.',netcdf_domain,'.',ext); |
| 89 |
ncJFz = netcdf(ferfile,'nowrite'); |
| 90 |
[JFzlon JFzlat JFzdpt] = coordfromnc(ncJFz); |
| 91 |
JFz = ncJFz{4}(1,:,:); |
| 92 |
|
| 93 |
ferfile = strcat(pathname,sla,snapshot,sla,filJBz,'.',netcdf_domain,'.',ext); |
| 94 |
ncJBz = netcdf(ferfile,'nowrite'); |
| 95 |
[JBzlon JBzlat JBzdpt] = coordfromnc(ncJBz); |
| 96 |
JBz = ncJBz{4}(1,:,:); |
| 97 |
|
| 98 |
ferfile = strcat(pathname,sla,snapshot,sla,filQnet,'.',netcdf_domain,'.',ext); |
| 99 |
ncQnet = netcdf(ferfile,'nowrite'); |
| 100 |
[Qnetlon Qnetlat Qnetdpt] = coordfromnc(ncQnet); |
| 101 |
Qnet = ncQnet{4}(1,:,:); |
| 102 |
|
| 103 |
ferfile = strcat(pathname,sla,snapshot,sla,filQEk,'.',netcdf_domain,'.',ext); |
| 104 |
ncQEk = netcdf(ferfile,'nowrite'); |
| 105 |
[QEklon QEklat QEkdpt] = coordfromnc(ncQEk); |
| 106 |
QEk = ncQEk{4}(1,:,:); |
| 107 |
|
| 108 |
ferfile = strcat(pathname,sla,snapshot,sla,filMLD,'.',netcdf_domain,'.',ext); |
| 109 |
ncMLD = netcdf(ferfile,'nowrite'); |
| 110 |
[MLDlon MLDlat MLDdpt] = coordfromnc(ncMLD); |
| 111 |
MLD = -ncMLD{4}(1,:,:); |
| 112 |
|
| 113 |
%%%%%%%%%%%%%%%% |
| 114 |
% Q is defined on the same grid of ST but troncated by extrem 2 points, then here |
| 115 |
% make all fields defined with same limits... |
| 116 |
disp('Reshape them') |
| 117 |
ST = squeeze(ST(2:nz+1,2:ny+1,2:nx+1)); |
| 118 |
STdpt = STdpt(2:nz+1); |
| 119 |
STlon = STlon(2:nx+1); |
| 120 |
STlat = STlat(2:ny+1); |
| 121 |
T = squeeze(T(2:nz+1,2:ny+1,2:nx+1)); |
| 122 |
Tdpt = Tdpt(2:nz+1); |
| 123 |
Tlon = Tlon(2:nx+1); |
| 124 |
Tlat = Tlat(2:ny+1); |
| 125 |
JBz = squeeze(JBz(2:ny+1,2:nx+1)); |
| 126 |
JBzlon = JBzlon(2:nx+1); |
| 127 |
JBzlat = JBzlat(2:ny+1); |
| 128 |
Qnet = squeeze(Qnet(2:ny+1,2:nx+1)); |
| 129 |
Qnetlon = Qnetlon(2:nx+1); |
| 130 |
Qnetlat = Qnetlat(2:ny+1); |
| 131 |
MLD = squeeze(MLD(2:ny+1,2:nx+1)); |
| 132 |
MLDlon = MLDlon(2:nx+1); |
| 133 |
MLDlat = MLDlat(2:ny+1); |
| 134 |
|
| 135 |
% Grid: |
| 136 |
global domain subdomain2 subdomain |
| 137 |
grid_setup |
| 138 |
%subdomain = subdomain2; |
| 139 |
|
| 140 |
|
| 141 |
%%%%%%%%%%%%%%%% |
| 142 |
% Here we determine the isosurface and its depth: |
| 143 |
if getiso |
| 144 |
disp('Get iso-ST') |
| 145 |
[Iiso mask] = subfct_getisoS(ST,iso); |
| 146 |
Diso = ones(size(Iiso)).*NaN; |
| 147 |
Qiso = ones(size(Iiso)).*NaN; |
| 148 |
for ix = 1 : size(ST,3) |
| 149 |
for iy = 1 : size(ST,2) |
| 150 |
if ~isnan(Iiso(iy,ix)) |
| 151 |
Diso(iy,ix) = STdpt(Iiso(iy,ix)); |
| 152 |
Qiso(iy,ix) = Q(Iiso(iy,ix),iy,ix); |
| 153 |
end %if |
| 154 |
end, end %for iy, ix |
| 155 |
end %if |
| 156 |
|
| 157 |
%%%%%%%%%%%%%%%% |
| 158 |
% "Normalise" the PV: |
| 159 |
fO = 2*(2*pi/86400)*sin(32*pi/180); |
| 160 |
dST = 27.6-25.4; |
| 161 |
H = -1000; |
| 162 |
RHOo = 1000; |
| 163 |
Qref = -fO/RHOo*dST/H; |
| 164 |
if getiso, QisoN = Qiso./Qref; end |
| 165 |
|
| 166 |
|
| 167 |
%%%%%%%%%%%%%%%% |
| 168 |
%%%%%%%%%%%%%%%% |
| 169 |
% Plots: |
| 170 |
disp('Plots ...') |
| 171 |
|
| 172 |
if find(pl==1) |
| 173 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 174 |
% Map projection: |
| 175 |
m_proj('mercator','long',subdomain.limlon,'lat',subdomain.limlat); |
| 176 |
%m_proj('mercator','long',subdomain.limlon,'lat',[25 40]); |
| 177 |
m_proj('mercator','long',subdomain.limlon,'lat',[25 50]); |
| 178 |
|
| 179 |
% Other settings: |
| 180 |
CBAR = 'h'; % Colorbar orientation |
| 181 |
Tcontour = [17 19]; Tcolor = [0 0 0]; % Theta contours |
| 182 |
Hcontour = [0:200:600]; Hcolor = [0 0 0]; % MLD contours |
| 183 |
unit = ''; % Default unit |
| 184 |
|
| 185 |
% Which variables to plot: |
| 186 |
%wvar = [1 2 3 4 6]; |
| 187 |
wvar = [1 2]; |
| 188 |
|
| 189 |
% Variables loop: |
| 190 |
for ip=1:length(wvar) |
| 191 |
% Default: |
| 192 |
showT = 1; % Iso-Theta contours |
| 193 |
showW = 0; % Windstress arrows |
| 194 |
showH = 0; % Mixed Layer Depth |
| 195 |
CONT = 0; % Contours instead of pcolor |
| 196 |
|
| 197 |
figure(ip);clf;hold on;iw=1;jw=1; |
| 198 |
if it==1, mini; end |
| 199 |
switch wvar(ip) |
| 200 |
case 1 |
| 201 |
C = Diso; |
| 202 |
Clon = Qlon; Clat = Qlat; |
| 203 |
tit = strcat(snapshot,' / Depth of \sigma_\theta=',num2str(iso),'kg.m^{-3}'); |
| 204 |
showW = 0; % Windstress |
| 205 |
cx = [-600 0]; |
| 206 |
unit = 'm'; |
| 207 |
titf = 'Diso'; |
| 208 |
case 2 |
| 209 |
C = QisoN; % C = Qiso; |
| 210 |
Clon = Qlon; Clat = Qlat; |
| 211 |
tit = strcat(snapshot,'/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz) / q_{ref}'); |
| 212 |
%tit = strcat(snapshot,'/ Potential vorticity field: q = - f . d\sigma_\theta/dz / \rho'); |
| 213 |
showW = 0; % Windstress |
| 214 |
cx = [0 1]*10; |
| 215 |
unit = char(1); |
| 216 |
titf = 'QisoN'; |
| 217 |
case 3 |
| 218 |
C = JFz; |
| 219 |
Clon = JFzlon; Clat = JFzlat; |
| 220 |
tit = strcat(snapshot,'/ Mechanical PV flux J^F_z and windstress'); |
| 221 |
showW = 1; % Windstress |
| 222 |
cx = [-1 1]*10^(-11); |
| 223 |
unit = 'kg/m^3/s^2'; |
| 224 |
titf = 'JFz'; |
| 225 |
case 4 |
| 226 |
C = JBz; |
| 227 |
Clon = JBzlon; Clat = JBzlat; |
| 228 |
tit = strcat(snapshot,'/ Diabatic PV flux J^B_z and windstress'); |
| 229 |
showW = 1; % Windstress |
| 230 |
cx = [-1 1]*10^(-11); |
| 231 |
unit = 'kg/m^3/s^2'; |
| 232 |
titf = 'JBz'; |
| 233 |
case 5 |
| 234 |
C = Qnet; |
| 235 |
Clon = Qnetlon; Clat = Qnetlat; |
| 236 |
tit = strcat(snapshot,'/ Net surface heat flux Q_{net} and MLD'); |
| 237 |
showH = 1; |
| 238 |
showW = 1; |
| 239 |
cx = [-1 1]*600; |
| 240 |
unit = 'W/m^2'; |
| 241 |
titf = 'Qnet'; |
| 242 |
case 6 |
| 243 |
C = JFz./JBz; |
| 244 |
Clon = JFzlon; Clat = JFzlat; |
| 245 |
tit = strcat(snapshot,'/ Ratio: J^F_z/J^B_z'); |
| 246 |
cx = [-1 1]*5; |
| 247 |
unit = char(1); |
| 248 |
titf = 'JFz_vs_JBz'; |
| 249 |
end %switch what to plot |
| 250 |
|
| 251 |
% Draw variable: |
| 252 |
if CONT ~= 1 |
| 253 |
m_pcolor(Clon,Clat,C);shading flat |
| 254 |
else |
| 255 |
m_contourf(Clon,Clat,C,CONTv); |
| 256 |
end |
| 257 |
caxis(cx); |
| 258 |
ccol(ip) = colorbar(CBAR,'fontsize',10); |
| 259 |
ctitle(ccol(ip),unit); |
| 260 |
title(tit); |
| 261 |
|
| 262 |
if showT |
| 263 |
[cs,h] = m_contour(Tlon,Tlat,squeeze(T(1,:,:)),Tcontour); |
| 264 |
clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200); |
| 265 |
for ih=1:length(h) |
| 266 |
set(h(ih),'edgecolor',Tcolor,'linewidth',1); |
| 267 |
end |
| 268 |
end %if show THETA contours |
| 269 |
|
| 270 |
if showW |
| 271 |
dx = 10*diff(Txlon(1:2)); dy = 8*diff(Txlat(1:2)); |
| 272 |
dx = 20*diff(Txlon(1:2)); dy = 10*diff(Txlat(1:2)); |
| 273 |
lo = [Txlon(1):dx:Txlon(length(Txlon))]; |
| 274 |
la = [Txlat(1):dy:Txlat(length(Txlat))]; |
| 275 |
[lo la] = meshgrid(lo,la); |
| 276 |
Txn = interp2(Txlat,Txlon,Tx',la,lo); |
| 277 |
Tyn = interp2(Txlat,Txlon,Ty',la,lo); |
| 278 |
s = 2; |
| 279 |
m_quiver(lo,la,Txn,Tyn,s,'w'); |
| 280 |
% m_quiver(lo,la,-(1+sin(la*pi/180)).*Txn,(1+sin(la*pi/180)).*Tyn,s,'w'); |
| 281 |
end %if show windstress |
| 282 |
|
| 283 |
if showH |
| 284 |
[cs,h] = m_contour(MLDlon,MLDlat,MLD,Hcontour); |
| 285 |
clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200); |
| 286 |
for ih=1:length(h) |
| 287 |
set(h(ih),'edgecolor',Hcolor,'linewidth',1); |
| 288 |
end |
| 289 |
end %if show Mixed Layer depth |
| 290 |
|
| 291 |
m_coast('patch',[0 0 0]);m_grid |
| 292 |
set(gcf,'name',titf); |
| 293 |
|
| 294 |
%%%%%%%%%%%%%%%% |
| 295 |
drawnow |
| 296 |
if prtimg |
| 297 |
set(gcf,'color','white') |
| 298 |
set(gcf,'paperposition',[0.6 6.5 25 14]); |
| 299 |
exportj(gcf,1,strcat(outimg,sla,titf,'.',snapshot)); |
| 300 |
end %if |
| 301 |
|
| 302 |
|
| 303 |
end %for ip |
| 304 |
|
| 305 |
|
| 306 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 307 |
end %if pl1 |
| 308 |
|
| 309 |
|
| 310 |
|
| 311 |
if find(pl==2) |
| 312 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 313 |
% Define sections: |
| 314 |
sectTYPE = 'z'; % m:meridional or z:zonal sections |
| 315 |
switch sectTYPE |
| 316 |
case 'm' % Fixed longitude |
| 317 |
lo = Qlon(max(find(Qlon<= 360-66 ))); sectNAME = '294E-66W'; |
| 318 |
case 'z' % Fixed latitude |
| 319 |
la = Qlat(max(find(Qlat<= 36 ))); sectNAME = '36N'; |
| 320 |
end %switch |
| 321 |
|
| 322 |
% Other settings: |
| 323 |
CBAR = 'h'; % Colorbar orientation |
| 324 |
Tcontour = [17 19]; Tcolor = [0 0 0]; % Theta contours |
| 325 |
STcontour = [20:.25:30]; STcolor = [1 1 1]; % Sigma-Theta contours |
| 326 |
unit = ''; % Default unit |
| 327 |
|
| 328 |
% Which variables to plot: |
| 329 |
wvar = [4]; |
| 330 |
|
| 331 |
% Variables loop: |
| 332 |
for ip=1:length(wvar) |
| 333 |
% Default: |
| 334 |
showT = 1; % Iso-Theta contours |
| 335 |
showST = 0; % Iso-Sigma-Theta contours |
| 336 |
showH = 0; % Mixed Layer Depth |
| 337 |
CONT = 0; % Contours instead of pcolor |
| 338 |
|
| 339 |
figure(ip);clf;hold on;iw=1;jw=1; |
| 340 |
%if it==1, mini; end |
| 341 |
switch wvar(ip) |
| 342 |
case 1 |
| 343 |
C = T; |
| 344 |
Clon = Tlon; Clat = Tlat; Cdpt = Tdpt; |
| 345 |
tit = strcat(snapshot,' / \theta'); |
| 346 |
cx = [0 30]; |
| 347 |
unit = 'K'; |
| 348 |
titf = 'THETA'; |
| 349 |
case 2 |
| 350 |
C = ST; |
| 351 |
Clon = STlon; Clat = STlat; Cdpt = STdpt; |
| 352 |
tit = strcat(snapshot,' / \sigma_\theta'); |
| 353 |
cx = [24 28]; |
| 354 |
unit = 'kg/m^3'; |
| 355 |
titf = 'SIGMATHETA'; |
| 356 |
case 3 |
| 357 |
C = Q/Qref; |
| 358 |
Clon = Qlon; Clat = Qlat; Cdpt = Qdpt; |
| 359 |
tit = strcat(snapshot,'/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz) / q_{ref}'); |
| 360 |
%tit = strcat(snapshot,'/ Potential vorticity field: q = - f . d\sigma_\theta/dz / \rho'); |
| 361 |
cx = [0 1]*2; |
| 362 |
unit = char(1); |
| 363 |
titf = 'splPVn'; |
| 364 |
CONTv = [0:1:5]; CONTv = [ [-1:.1:1] [1:1:10]]; |
| 365 |
CONT = 1; |
| 366 |
showT = 0; |
| 367 |
showST = 1; |
| 368 |
case 4 |
| 369 |
C = Q; |
| 370 |
Clon = Qlon; Clat = Qlat; Cdpt = Qdpt; |
| 371 |
tit = strcat(snapshot,'/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz)'); |
| 372 |
cx = [0 1]*10^(-9); |
| 373 |
unit = char(1); |
| 374 |
titf = 'splPV'; |
| 375 |
unit = '1/s/m'; |
| 376 |
CONTv = [ [0:.025:.5] [.5:.1:1] [1:1:10]]*10^(-9); |
| 377 |
CONT = 0; |
| 378 |
showT = 0; |
| 379 |
showST = 1; |
| 380 |
end %switch what to plot |
| 381 |
|
| 382 |
% Find the section: |
| 383 |
switch sectTYPE |
| 384 |
case 'm' % Fixed longitude |
| 385 |
s = max( find( Clon <= lo ) ); |
| 386 |
C = squeeze( C(:,:,s) ); |
| 387 |
Cx = Clat; Xlab = 'Latitude (^oN)'; |
| 388 |
Cy = Cdpt; Ylab = 'Depth(m)'; |
| 389 |
|
| 390 |
s = max( find( Tlon <= lo ) ); |
| 391 |
sectT = squeeze( T(:,:,s) ); sectTx = Tlat; sectTy = Tdpt; |
| 392 |
s = max( find(STlon <= lo ) ); |
| 393 |
sectST = squeeze(ST(:,:,s) ); sectSTx = STlat; sectSTy = STdpt; |
| 394 |
case 'z' % Fixed latitude |
| 395 |
s = max( find( Clat <= la ) ); |
| 396 |
C = squeeze( C(:,s,:) ); |
| 397 |
Cx = 360-Clon; Xlab = 'Longitude (^oW)'; |
| 398 |
Cy = Cdpt; Ylab = 'Depth(m)'; |
| 399 |
|
| 400 |
s = max( find( Tlat <= la ) ); |
| 401 |
sectT = squeeze( T(:,s,:) ); sectTx = 360-Tlon; sectTy = Tdpt; |
| 402 |
s = max( find(STlat <= la ) ); |
| 403 |
sectST = squeeze(ST(:,s,:) ); sectSTx = 360-STlon; sectSTy = STdpt; |
| 404 |
end %switch |
| 405 |
|
| 406 |
% Plots: |
| 407 |
clf; hold on |
| 408 |
|
| 409 |
if CONT ~= 1 |
| 410 |
pcolor(Cx,Cy,C); shading interp |
| 411 |
colormap(logcolormap(256,12)); |
| 412 |
else |
| 413 |
[cs,h]=contourf(Cx,Cy,C,CONTv); |
| 414 |
colormap(mycolormap(logcolormap(64,5),length(CONTv))); |
| 415 |
end |
| 416 |
caxis(cx); |
| 417 |
ccol(ip) = colorbar(CBAR,'fontsize',10); |
| 418 |
ctitle(ccol(ip),unit); |
| 419 |
title(strcat('Section:',sectNAME,'/',tit)); |
| 420 |
axis([min([Cx(1) Cx(length(Cx))]) max([Cx(1) Cx(length(Cx))]) ... |
| 421 |
min([Cy(1) Cy(length(Cy))]) max([Cy(1) Cy(length(Cy))]) ]); |
| 422 |
axis([55 75 -400 Cdpt(1)]); |
| 423 |
axis([55 75 min([Cy(1) Cy(length(Cy))]) max([Cy(1) Cy(length(Cy))]) ]); |
| 424 |
xlabel(Xlab); ylabel(Ylab); |
| 425 |
if sectTYPE=='z', set(gca,'xdir','reverse'); end |
| 426 |
box on |
| 427 |
|
| 428 |
if showT |
| 429 |
[cs,h]=contour(sectTx,sectTy,sectT,Tcontour); |
| 430 |
clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200); |
| 431 |
for ih=1:length(h) |
| 432 |
set(h(ih),'edgecolor',Tcolor,'linewidth',1); |
| 433 |
end |
| 434 |
end %if show THETA contours |
| 435 |
if showST |
| 436 |
STcontourE=[25.5 25.75]; |
| 437 |
[cs,h]=contour(sectSTx,sectSTy,sectST,STcontour); |
| 438 |
clabel(cs,h,'fontsize',8,'color',STcolor,'labelspacing',200); |
| 439 |
set(h,'edgecolor',STcolor); |
| 440 |
set(h,'linewidth',1); |
| 441 |
ch=get(h,'Children'); |
| 442 |
for ich=1:length(ch) |
| 443 |
%set(ch(ich),'edgecolor',STcolor); |
| 444 |
%set(ch(ich),'linewidth',1); |
| 445 |
for ii = 1 : length(STcontourE) |
| 446 |
if get(ch(ich),'UserData') == STcontourE(ii) |
| 447 |
set(ch(ich),'linewidth',2); |
| 448 |
end %if |
| 449 |
end %for ii |
| 450 |
end %for ich |
| 451 |
end %if show SIGMATHETA contours |
| 452 |
|
| 453 |
|
| 454 |
%%%%%%%%%%%%%%%% |
| 455 |
drawnow |
| 456 |
if prtimg |
| 457 |
set(gcf,'color','white') |
| 458 |
exportj(gcf,1,strcat(outimg,sla,titf,'.',snapshot,'.section_',sectNAME)); |
| 459 |
end %if |
| 460 |
|
| 461 |
end %for ip |
| 462 |
|
| 463 |
|
| 464 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 465 |
end %if pl2 |
| 466 |
|
| 467 |
|
| 468 |
|
| 469 |
|
| 470 |
|
| 471 |
|
| 472 |
|
| 473 |
|
| 474 |
|
| 475 |
if find(pl==3) |
| 476 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 477 |
% Define sections: |
| 478 |
sectTYPE = 'z'; % m:meridional or z:zonal sections |
| 479 |
|
| 480 |
switch sectTYPE |
| 481 |
case 'm' % Fixed longitude |
| 482 |
lo = 294; xlim=[20 42]; sectNAME = '294E'; |
| 483 |
|
| 484 |
case 'z' % Fixed latitude |
| 485 |
la = 36; xlim=[55 75]; sectNAME = '36N'; |
| 486 |
end %switch |
| 487 |
|
| 488 |
% Other settings: |
| 489 |
CBAR = 'v'; % Colorbar orientation |
| 490 |
Tcontour = [17 19]; Tcolor = [0 0 0]; % Theta contours |
| 491 |
STcontour = [20:.1:30]; STcolor = [1 1 1]*.5; % Sigma-Theta contours |
| 492 |
Hcolor = [0 0 0]; % Mixed layer depth color |
| 493 |
unit = ''; % Default unit |
| 494 |
|
| 495 |
% Which variables to plot: |
| 496 |
wvar = [4]; |
| 497 |
|
| 498 |
% Variables loop: |
| 499 |
for ip=1:length(wvar) |
| 500 |
% Default: |
| 501 |
showT = 1; % Iso-Theta contours |
| 502 |
showST = 0; % Iso-Sigma-Theta contours |
| 503 |
showH = 0; % Mixed Layer Depth |
| 504 |
CONT = 0; % Contours instead of pcolor |
| 505 |
|
| 506 |
figure(ip);clf;hold on;iw=1;jw=1; |
| 507 |
%if it==1, mini; end |
| 508 |
switch wvar(ip) |
| 509 |
case 1 |
| 510 |
C = T; |
| 511 |
Clon = Tlon; Clat = Tlat; Cdpt = Tdpt; |
| 512 |
tit = strcat(snapshot,' / \theta'); |
| 513 |
cx = [0 30]; |
| 514 |
unit = 'K'; |
| 515 |
titf = 'THETA'; |
| 516 |
case 2 |
| 517 |
C = ST; |
| 518 |
Clon = STlon; Clat = STlat; Cdpt = STdpt; |
| 519 |
tit = strcat(snapshot,' / \sigma_\theta'); |
| 520 |
cx = [24 28]; |
| 521 |
unit = 'kg/m^3'; |
| 522 |
titf = 'SIGMATHETA'; |
| 523 |
case 3 |
| 524 |
C = Q/Qref; |
| 525 |
Clon = Qlon; Clat = Qlat; Cdpt = Qdpt; |
| 526 |
tit = strcat('/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz) / q_{ref}'); |
| 527 |
%tit = strcat(snapshot,'/ Potential vorticity field: q = - f . d\sigma_\theta/dz / \rho'); |
| 528 |
cx = [0 1]*2; |
| 529 |
unit = char(1); |
| 530 |
titf = 'splPVn'; |
| 531 |
CONTv = [0:1:5]; CONTv = [ [-1:.1:1] [1:1:10]]; |
| 532 |
CONT = 1; |
| 533 |
showT = 0; |
| 534 |
showST = 1; |
| 535 |
case 4 |
| 536 |
C = Q; |
| 537 |
Clon = Qlon; Clat = Qlat; Cdpt = Qdpt; |
| 538 |
tit = strcat('/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz)'); |
| 539 |
cx = [0 1]*10^(-9); |
| 540 |
unit = char(1); |
| 541 |
titf = 'splPV'; |
| 542 |
unit = '1/s/m'; |
| 543 |
CONTv = [ [0:.025:.5] [.5:.1:1] [1:1:10]]*10^(-9); |
| 544 |
CONT = 0; |
| 545 |
showT = 0; |
| 546 |
showST = 1; |
| 547 |
showH = 1; |
| 548 |
end %switch what to plot |
| 549 |
|
| 550 |
% Find the section: |
| 551 |
switch sectTYPE |
| 552 |
case 'm' % Fixed longitude |
| 553 |
s = max( find( Clon <= lo ) ); |
| 554 |
C = squeeze( C(:,:,s) ); |
| 555 |
Cx = Clat; Xlab = 'Latitude (^oN)'; |
| 556 |
Cy = Cdpt; Ylab = 'Depth(m)'; |
| 557 |
|
| 558 |
s = max( find( Tlon <= lo ) ); |
| 559 |
sectT = squeeze( T(:,:,s) ); sectTx = Tlat; sectTy = Tdpt; |
| 560 |
s = max( find(STlon <= lo ) ); |
| 561 |
sectST = squeeze(ST(:,:,s) ); sectSTx = STlat; sectSTy = STdpt; |
| 562 |
s = max( find(MLDlon <= lo ) ); |
| 563 |
sectH = squeeze(MLD(:,s) ); sectHx = MLDlat; sectHy = MLDdpt; |
| 564 |
|
| 565 |
dx = 1; |
| 566 |
s = max( find( JFzlon <= lo ) ); |
| 567 |
sectJFz = squeeze( JFz(:,s) ); sectJFzx = JFzlat; sectJFzy = JFzdpt; |
| 568 |
s = max( find( JBzlon <= lo ) ); |
| 569 |
sectJBz = squeeze( JBz(:,s) ); sectJBzx = JBzlat; sectJBzy = JBzdpt; |
| 570 |
|
| 571 |
case 'z' % Fixed latitude |
| 572 |
s = max( find( Clat <= la ) ); |
| 573 |
C = squeeze( C(:,s,:) ); |
| 574 |
Cx = 360-Clon; Xlab = 'Longitude (^oW)'; |
| 575 |
Cy = Cdpt; Ylab = 'Depth(m)'; |
| 576 |
|
| 577 |
s = max( find( Tlat <= la ) ); |
| 578 |
sectT = squeeze( T(:,s,:) ); sectTx = 360-Tlon; sectTy = Tdpt; |
| 579 |
s = max( find(STlat <= la ) ); |
| 580 |
sectST = squeeze(ST(:,s,:) ); sectSTx = 360-STlon; sectSTy = STdpt; |
| 581 |
s = max( find(MLDlat <= la ) ); |
| 582 |
sectH = squeeze(MLD(s,:) ); sectHx = 360-MLDlon; sectHy = MLDdpt; |
| 583 |
|
| 584 |
dx = 1; |
| 585 |
s = max( find( JFzlat <= la ) ); |
| 586 |
sectJFz = squeeze( JFz(s,:) ); sectJFzx = 360-JFzlon; sectJFzy = JFzdpt; |
| 587 |
s = max( find( JBzlat <= la ) ); |
| 588 |
sectJBz = squeeze( JBz(s,:) ); sectJBzx = 360-JBzlon; sectJBzy = JBzdpt; |
| 589 |
end %switch |
| 590 |
|
| 591 |
% Plots: |
| 592 |
clf; hold on |
| 593 |
iw=2;jw=1; |
| 594 |
|
| 595 |
%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 596 |
sp(1)=subplot(iw,jw,1); hold on |
| 597 |
% Change grid: |
| 598 |
x = sectJBzx(1:dx:end); |
| 599 |
B = sectJBz(1:dx:end); |
| 600 |
F = sectJFz(1:dx:end); |
| 601 |
BpF = B+F; |
| 602 |
col=['b','g','r']; dl = 3; |
| 603 |
for ix = 1 : length(x) |
| 604 |
FLUX = [B(ix) F(ix) BpF(ix)]; |
| 605 |
FLUXposI = find(FLUX>=0); |
| 606 |
if ~isempty(FLUXposI) |
| 607 |
FLUXpos = FLUX(FLUXposI); |
| 608 |
co = col(FLUXposI); |
| 609 |
[L Li] = sort(FLUXpos); |
| 610 |
for ii = length(L):-1:1 |
| 611 |
line([1 1]*x(ix),[0 FLUXpos(Li(ii))],'color',co(Li(ii)),'linewidth',dl) |
| 612 |
end, end |
| 613 |
FLUXnegI = find(FLUX<0); |
| 614 |
if ~isempty(FLUXnegI) |
| 615 |
FLUXneg = FLUX(FLUXnegI); |
| 616 |
co = col(FLUXnegI); |
| 617 |
[L Li] = sort(FLUXneg); |
| 618 |
for ii = 1:length(L) |
| 619 |
line([1 1]*x(ix),[0 FLUXneg(Li(ii))],'color',co(Li(ii)),'linewidth',dl) |
| 620 |
end, end |
| 621 |
end |
| 622 |
line(get(gca,'xlim'),[0 0],'color','k') |
| 623 |
plot(x,BpF,'r'); |
| 624 |
grid on, box on |
| 625 |
ylabel('PV flux (kg/m^3/s^2)'); |
| 626 |
set(gca,'xlim',xlim); |
| 627 |
set(gca,'xtick',[xlim(1):2:xlim(2)]); |
| 628 |
set(gca,'xticklabel',[]); |
| 629 |
set(gca,'ylim',[-1 1]*1*10^(-11)); |
| 630 |
set(gca,'xdir','reverse'); |
| 631 |
title(strvcat(strcat('Section:',sectNAME,tit,... |
| 632 |
' and surface PV flux (blue: J^B_z, green: J^F_z, red: J^B_z+J^F_z)'),... |
| 633 |
'(gray: \sigma_\theta, black: Mixed layer depth)')); |
| 634 |
|
| 635 |
%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 636 |
sp(2)=subplot(iw,jw,2);hold on |
| 637 |
|
| 638 |
if CONT ~= 1 |
| 639 |
pcolor(Cx,Cy,C); shading interp |
| 640 |
colormap(logcolormap(256,12)); |
| 641 |
else |
| 642 |
[cs,h]=contourf(Cx,Cy,C,CONTv); |
| 643 |
colormap(mycolormap(logcolormap(64,5),length(CONTv))); |
| 644 |
end |
| 645 |
caxis(cx); ccol(ip) = colorbar(CBAR,'fontsize',10); ctitle(ccol(ip),unit); |
| 646 |
set(gca,'xlim',xlim); |
| 647 |
xlabel(Xlab); ylabel(Ylab); |
| 648 |
set(gca,'xdir','reverse'); |
| 649 |
box on |
| 650 |
set(gca,'xtick',[xlim(1):2:xlim(2)]); |
| 651 |
set(gca,'xticklabel',[xlim(1):2:xlim(2)]); |
| 652 |
set(gca,'ylim',[min([Cy(1) Cy(length(Cy))]) max([Cy(1) Cy(length(Cy))]) ]); |
| 653 |
set(gca,'ylim',[-600 max([Cy(1) Cy(length(Cy))])]); |
| 654 |
|
| 655 |
if showT |
| 656 |
[cs,h]=contour(sectTx,sectTy,sectT,Tcontour); |
| 657 |
clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200); |
| 658 |
for ih=1:length(h) |
| 659 |
set(h(ih),'edgecolor',Tcolor,'linewidth',1); |
| 660 |
end |
| 661 |
end %if show THETA contours |
| 662 |
|
| 663 |
if showST |
| 664 |
STcontourE=[25.5 25.75]; |
| 665 |
[cs,h]=contour(sectSTx,sectSTy,sectST,STcontour); |
| 666 |
clabel(cs,h,'fontsize',8,'color',STcolor,'labelspacing',200); |
| 667 |
set(h,'edgecolor',STcolor); |
| 668 |
set(h,'linewidth',1); |
| 669 |
ch=get(h,'Children'); |
| 670 |
for ich=1:length(ch) |
| 671 |
%set(ch(ich),'edgecolor',STcolor); |
| 672 |
%set(ch(ich),'linewidth',1); |
| 673 |
for ii = 1 : length(STcontourE) |
| 674 |
if get(ch(ich),'UserData') == STcontourE(ii) |
| 675 |
% set(ch(ich),'linewidth',1); |
| 676 |
end %if |
| 677 |
end %for ii |
| 678 |
end %for ich |
| 679 |
[cs,h]=contour(sectSTx,sectSTy,sectST,STcontourE); |
| 680 |
clabel(cs,h,'fontsize',8,'color',STcolor,'labelspacing',200); |
| 681 |
set(h,'edgecolor',[1 1 1]); |
| 682 |
set(h,'linewidth',2); |
| 683 |
|
| 684 |
end %if show SIGMATHETA contours |
| 685 |
|
| 686 |
if showH |
| 687 |
l=line(sectHx,sectH,'color',Hcolor,'linewidth',2); |
| 688 |
end %if show Mixed layer depth |
| 689 |
|
| 690 |
pos1=get(sp(1),'position'); |
| 691 |
pos2=get(sp(2),'position'); |
| 692 |
dy = .1; |
| 693 |
set(sp(1),'position',[pos1(1) pos1(2)+dy pos2(3) pos1(4)-dy]); |
| 694 |
%set(sp(2),'position',[pos2(1) pos2(2) pos2(3) pos2(4)+1.3*dy]); |
| 695 |
set(sp(2),'position',[pos2(1) pos2(2) pos2(3) pos1(2)+dy-1.01*pos2(2)]); |
| 696 |
|
| 697 |
|
| 698 |
|
| 699 |
%%%%%%%%%%%%%%%% |
| 700 |
drawnow |
| 701 |
set(gcf,'position',[4 48 888 430]); |
| 702 |
videotimeline(TIME,it,'b') |
| 703 |
if prtimg |
| 704 |
set(gcf,'color','white') |
| 705 |
set(gcf,'paperposition',[0.6 6.5 27.5 14]); |
| 706 |
exportj(gcf,1,strcat(outimg,sla,titf,'.',snapshot,'.withJz.section_',sectNAME)); |
| 707 |
end %if |
| 708 |
|
| 709 |
end %for ip |
| 710 |
|
| 711 |
|
| 712 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 713 |
end %if pl3 |
| 714 |
|
| 715 |
|
| 716 |
|
| 717 |
|
| 718 |
|
| 719 |
|
| 720 |
if find(pl==4) |
| 721 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 722 |
|
| 723 |
clf; |
| 724 |
m_proj('mercator','long',subdomain.limlon,'lat',subdomain.limlat); |
| 725 |
%m_proj('mercator','long',subdomain.limlon,'lat',[25 40]); |
| 726 |
%m_proj('mercator','long',subdomain.limlon,'lat',[25 50]); |
| 727 |
|
| 728 |
hold on |
| 729 |
|
| 730 |
m_pcolor(Qnetlon,Qnetlat,Qnet); |
| 731 |
caxis([-1 1]*500);canom |
| 732 |
|
| 733 |
dx = 10*diff(Txlon(1:2)); dy = 8*diff(Txlat(1:2)); |
| 734 |
dx = 20*diff(Txlon(1:2)); dy = 10*diff(Txlat(1:2)); |
| 735 |
lo = [Txlon(1):dx:Txlon(length(Txlon))]; |
| 736 |
la = [Txlat(1):dy:Txlat(length(Txlat))]; |
| 737 |
[lo la] = meshgrid(lo,la); |
| 738 |
Txn = interp2(Txlat,Txlon,Tx',la,lo); |
| 739 |
Tyn = interp2(Txlat,Txlon,Ty',la,lo); |
| 740 |
s = 3; |
| 741 |
m_quiver(lo,la,Txn,Tyn,s,'k'); |
| 742 |
|
| 743 |
m_coast('patch',[0 0 0]);m_grid |
| 744 |
ccol=colorbar('h');%ctitle(ccol,'m'); |
| 745 |
title(strcat(snapshot,'/ Net heat flux (W/m^2)')); |
| 746 |
[cs,h]=m_contour(MLDlon,MLDlat,MLD,[0:100:500]); |
| 747 |
clabel(cs,h,'labelspacing',600,'fontsize',8); |
| 748 |
for ih=1:length(h),set(h(ih),'edgecolor',[1 1 1],'linewidth',1);end; |
| 749 |
|
| 750 |
M(it) = getframe; |
| 751 |
|
| 752 |
drawnow |
| 753 |
if prtimg |
| 754 |
set(gcf,'color','white') |
| 755 |
set(gcf,'paperorientation','landscape'); |
| 756 |
%set(gcf,'paperposition',[0.6 6.5 25 14]); |
| 757 |
print(gcf,'-djpeg100',strcat(outimg,sla,snapshot,'.jpg')); |
| 758 |
%exportj(gcf,1,strcat(outimg,sla,snapshot)); |
| 759 |
end %if |
| 760 |
|
| 761 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 762 |
end %if pl4 |
| 763 |
|
| 764 |
|
| 765 |
%%%%%%%%%%%%%%%% |
| 766 |
%%%%%%%%%%%%%%%% |
| 767 |
|
| 768 |
|
| 769 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 770 |
end %for it |
Parent Directory
| 
Revision Log
| 
Revision Graph