%
% THIS IS NOT A FUNCTION !
%
% Plot time series of all variables
% 1 var / figure
% Outputs recording
%

clear
global sla netcdf_domain
pv_checkpath

netcdf_domain = 'climode';

% Date series:
TIME = ['20030201';'20030205';'20030212';'20030219';'20030226';...
        '20030301';'20030305';'20030312'];
TIME = ['20030201';'20030205';'20030212';'20030219';'20030226';...
        '20030301';'20030305'];
nt   = size(TIME,1);


iso    = 25.8; % Which sigma-theta surface ?
getiso = 1;    % We compute the isoST by default
outimg = 'img_tmp'; % Output directory
prtimg = 1; % Do we record figures as jpg files ?

% Plots type:
% 1: individual simple maps
% 2: meridional or zonal sections
% 3: meridional or zonal sections with Jz vectors at the surface
% 4: custom it
pl = 3;

if pl==2|pl==3, getiso = 0; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Video loop:
for it = 1 : nt
  
% Date:
snapshot = TIME(it,:);

%%%%%%%%%%%%%%%%
% NETCDF files name:
filPV   = 'splPV';
filST   = 'SIGMATHETA';
filT    = 'THETA';
filTx   = 'TAUX';
filTy   = 'TAUY';
filJFz  = 'JFz';
filJBz  = 'JBz';
filQnet = 'Qnet';
filQEk  = 'QEk';
filMLD  = 'KPPmld';

% Path and extension to find them:
pathname = strcat('netcdf-files',sla);
ext      = 'nc';

% Load fields:
disp('load fields...')
% (I keep proper axis for each variables in case of one day they would be different)
         ferfile = strcat(pathname,sla,snapshot,sla,filPV,'.',netcdf_domain,'.',ext);
             ncQ = netcdf(ferfile,'nowrite');
[Qlon Qlat Qdpt] = coordfromnc(ncQ);
               Q = ncQ{4}(:,:,:); clear ncQ ferfile
      [nz ny nx] = size(Q);

            ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext);
               ncST = netcdf(ferfile,'nowrite');
[STlon STlat STdpt] = coordfromnc(ncST);
                 ST = ncST{4}(:,:,:); clear ncST ferfile

         ferfile = strcat(pathname,sla,snapshot,sla,filT,'.',netcdf_domain,'.',ext);
             ncT = netcdf(ferfile,'nowrite');
[Tlon Tlat Tdpt] = coordfromnc(ncT);
               T = ncT{4}(:,:,:); clear ncT ferfile
	      
            ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext);
               ncTx = netcdf(ferfile,'nowrite');
[Txlon Txlat Txdpt] = coordfromnc(ncTx);
                 Tx = ncTx{4}(1,:,:); clear ncTx ferfile
            ferfile = strcat(pathname,sla,snapshot,sla,filTy,'.',netcdf_domain,'.',ext);
               ncTy = netcdf(ferfile,'nowrite');
[Tylon Tylat Tydpt] = coordfromnc(ncTy);
                 Ty = ncTy{4}(1,:,:); clear ncTy ferfile

               ferfile = strcat(pathname,sla,snapshot,sla,filJFz,'.',netcdf_domain,'.',ext);
                 ncJFz = netcdf(ferfile,'nowrite');
[JFzlon JFzlat JFzdpt] = coordfromnc(ncJFz);
                   JFz = ncJFz{4}(1,:,:);

               ferfile = strcat(pathname,sla,snapshot,sla,filJBz,'.',netcdf_domain,'.',ext);
                 ncJBz = netcdf(ferfile,'nowrite');
[JBzlon JBzlat JBzdpt] = coordfromnc(ncJBz);
                   JBz = ncJBz{4}(1,:,:);

               ferfile = strcat(pathname,sla,snapshot,sla,filQnet,'.',netcdf_domain,'.',ext);
                 ncQnet = netcdf(ferfile,'nowrite');
[Qnetlon Qnetlat Qnetdpt] = coordfromnc(ncQnet);
                   Qnet = ncQnet{4}(1,:,:);

               ferfile = strcat(pathname,sla,snapshot,sla,filQEk,'.',netcdf_domain,'.',ext);
                 ncQEk = netcdf(ferfile,'nowrite');
[QEklon QEklat QEkdpt] = coordfromnc(ncQEk);
                   QEk = ncQEk{4}(1,:,:);

               ferfile = strcat(pathname,sla,snapshot,sla,filMLD,'.',netcdf_domain,'.',ext);
                 ncMLD = netcdf(ferfile,'nowrite');
[MLDlon MLDlat MLDdpt] = coordfromnc(ncMLD);
                   MLD = -ncMLD{4}(1,:,:);

%%%%%%%%%%%%%%%%
% Q is defined on the same grid of ST but troncated by extrem 2 points, then here
% make all fields defined with same limits...
disp('Reshape them')
ST    = squeeze(ST(2:nz+1,2:ny+1,2:nx+1));
STdpt = STdpt(2:nz+1);
STlon = STlon(2:nx+1);
STlat = STlat(2:ny+1);
T    = squeeze(T(2:nz+1,2:ny+1,2:nx+1));
Tdpt = Tdpt(2:nz+1);
Tlon = Tlon(2:nx+1);
Tlat = Tlat(2:ny+1);
JBz    = squeeze(JBz(2:ny+1,2:nx+1));
JBzlon = JBzlon(2:nx+1);
JBzlat = JBzlat(2:ny+1);
Qnet    = squeeze(Qnet(2:ny+1,2:nx+1));
Qnetlon = Qnetlon(2:nx+1);
Qnetlat = Qnetlat(2:ny+1);
MLD    = squeeze(MLD(2:ny+1,2:nx+1));
MLDlon = MLDlon(2:nx+1);
MLDlat = MLDlat(2:ny+1);
  
% Grid:
global domain subdomain2 subdomain
grid_setup
%subdomain = subdomain2;


%%%%%%%%%%%%%%%%
% Here we determine the isosurface and its depth:
if getiso
  disp('Get iso-ST')
[Iiso mask] = subfct_getisoS(ST,iso);
Diso = ones(size(Iiso)).*NaN;
Qiso = ones(size(Iiso)).*NaN;
for ix = 1 : size(ST,3)
  for iy = 1 : size(ST,2) 
    if ~isnan(Iiso(iy,ix))
       Diso(iy,ix) = STdpt(Iiso(iy,ix));
       Qiso(iy,ix) =     Q(Iiso(iy,ix),iy,ix);
    end %if
end, end %for iy, ix
end %if

%%%%%%%%%%%%%%%%
% "Normalise" the PV:
fO  = 2*(2*pi/86400)*sin(32*pi/180);
dST = 27.6-25.4;
H   = -1000;
RHOo = 1000;
Qref = -fO/RHOo*dST/H;
if getiso, QisoN = Qiso./Qref; end


%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
% Plots:
disp('Plots ...')

if find(pl==1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Map projection:
m_proj('mercator','long',subdomain.limlon,'lat',subdomain.limlat);
%m_proj('mercator','long',subdomain.limlon,'lat',[25 40]);
m_proj('mercator','long',subdomain.limlon,'lat',[25 50]);

% Other settings:
CBAR = 'h'; % Colorbar orientation
Tcontour = [17 19]; Tcolor = [0 0 0]; % Theta contours
Hcontour = [0:200:600]; Hcolor = [0 0 0]; % MLD contours
unit = ''; % Default unit

% Which variables to plot:
%wvar = [1 2 3 4 6];
wvar = [1 2];

% Variables loop:
for ip=1:length(wvar)
  % Default:
  showT = 1; % Iso-Theta contours
  showW = 0; % Windstress arrows
  showH = 0; % Mixed Layer Depth
  CONT  = 0; % Contours instead of pcolor
  
  figure(ip);clf;hold on;iw=1;jw=1;
  if it==1, mini; end 
switch wvar(ip)
 case 1
  C     = Diso;
  Clon  = Qlon; Clat = Qlat;
  tit   = strcat(snapshot,' / Depth of \sigma_\theta=',num2str(iso),'kg.m^{-3}');
  showW = 0; % Windstress
  cx    = [-600 0];
  unit  = 'm';
  titf  = 'Diso';
 case 2
  C     = QisoN; % C = Qiso;
  Clon  = Qlon; Clat = Qlat;
  tit   = strcat(snapshot,'/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz) / q_{ref}');
  %tit   = strcat(snapshot,'/ Potential vorticity field: q = - f . d\sigma_\theta/dz / \rho');
  showW = 0; % Windstress
  cx    = [0 1]*10;
  unit  = char(1);
  titf  = 'QisoN';
 case 3
  C     = JFz;
  Clon  = JFzlon; Clat = JFzlat;
  tit   = strcat(snapshot,'/ Mechanical PV flux J^F_z and windstress');
  showW = 1; % Windstress
  cx    = [-1 1]*10^(-11);
  unit  = 'kg/m^3/s^2';
  titf  = 'JFz';
 case 4
  C     = JBz;
  Clon  = JBzlon; Clat = JBzlat;
  tit   = strcat(snapshot,'/ Diabatic PV flux J^B_z and windstress');
  showW = 1; % Windstress
  cx    = [-1 1]*10^(-11);
  unit  = 'kg/m^3/s^2';
  titf  = 'JBz';
 case 5
  C     = Qnet;
  Clon  = Qnetlon; Clat = Qnetlat;
  tit   = strcat(snapshot,'/ Net surface heat flux Q_{net} and MLD');
  showH = 1;
  showW = 1;
  cx    = [-1 1]*600;
  unit  = 'W/m^2';
  titf  = 'Qnet';
 case 6
  C     = JFz./JBz;
  Clon  = JFzlon; Clat = JFzlat;
  tit   = strcat(snapshot,'/ Ratio: J^F_z/J^B_z');
  cx    = [-1 1]*5;
  unit  = char(1);
  titf  = 'JFz_vs_JBz';
end %switch what to plot

% Draw variable:
if CONT ~= 1
  m_pcolor(Clon,Clat,C);shading flat
else
  m_contourf(Clon,Clat,C,CONTv);
end
caxis(cx);
ccol(ip) = colorbar(CBAR,'fontsize',10);
ctitle(ccol(ip),unit);
title(tit);

if showT
  [cs,h] = m_contour(Tlon,Tlat,squeeze(T(1,:,:)),Tcontour);
  clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200);
  for ih=1:length(h)
    set(h(ih),'edgecolor',Tcolor,'linewidth',1);
  end
end %if show THETA contours

if showW
      dx = 10*diff(Txlon(1:2)); dy = 8*diff(Txlat(1:2));
      dx = 20*diff(Txlon(1:2)); dy = 10*diff(Txlat(1:2));
      lo = [Txlon(1):dx:Txlon(length(Txlon))];
      la = [Txlat(1):dy:Txlat(length(Txlat))];
      [lo la] = meshgrid(lo,la);
      Txn = interp2(Txlat,Txlon,Tx',la,lo);
      Tyn = interp2(Txlat,Txlon,Ty',la,lo);
      s = 2;
      m_quiver(lo,la,Txn,Tyn,s,'w');
%      m_quiver(lo,la,-(1+sin(la*pi/180)).*Txn,(1+sin(la*pi/180)).*Tyn,s,'w');
end %if show windstress

if showH
  [cs,h] = m_contour(MLDlon,MLDlat,MLD,Hcontour);
  clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200);
  for ih=1:length(h)
    set(h(ih),'edgecolor',Hcolor,'linewidth',1);
  end
end %if show Mixed Layer depth

m_coast('patch',[0 0 0]);m_grid
set(gcf,'name',titf);

%%%%%%%%%%%%%%%%
drawnow
if prtimg
set(gcf,'color','white') 
set(gcf,'paperposition',[0.6 6.5 25 14]);
exportj(gcf,1,strcat(outimg,sla,titf,'.',snapshot));
end %if


end %for ip


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end %if pl1



if find(pl==2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Define sections:
sectTYPE = 'z'; % m:meridional or z:zonal sections
switch sectTYPE
  case 'm' % Fixed longitude
    lo = Qlon(max(find(Qlon<= 360-66 ))); sectNAME = '294E-66W';
  case 'z' % Fixed latitude
    la = Qlat(max(find(Qlat<= 36 )));     sectNAME = '36N';
end %switch

% Other settings:
CBAR = 'h'; % Colorbar orientation
Tcontour  = [17 19];  Tcolor = [0 0 0]; % Theta contours
STcontour = [20:.25:30]; STcolor = [1 1 1]; % Sigma-Theta contours
unit = ''; % Default unit

% Which variables to plot:
wvar = [4];

% Variables loop:
for ip=1:length(wvar)
  % Default:
  showT  = 1; % Iso-Theta contours
  showST = 0; % Iso-Sigma-Theta contours
  showH  = 0; % Mixed Layer Depth
  CONT   = 0; % Contours instead of pcolor
  
  figure(ip);clf;hold on;iw=1;jw=1;
  %if it==1, mini; end 
switch wvar(ip)
 case 1
  C     = T;
  Clon  = Tlon; Clat = Tlat; Cdpt = Tdpt;
  tit   = strcat(snapshot,' / \theta');
  cx    = [0 30];
  unit  = 'K';
  titf  = 'THETA';
 case 2
  C     = ST;
  Clon  = STlon; Clat = STlat; Cdpt = STdpt;
  tit   = strcat(snapshot,' / \sigma_\theta');
  cx    = [24 28];
  unit  = 'kg/m^3';
  titf  = 'SIGMATHETA';
 case 3
  C      = Q/Qref;
  Clon   = Qlon; Clat = Qlat; Cdpt = Qdpt;
  tit    = strcat(snapshot,'/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz) / q_{ref}');
  %tit    = strcat(snapshot,'/ Potential vorticity field: q = - f . d\sigma_\theta/dz / \rho');
  cx     = [0 1]*2;
  unit   = char(1);
  titf   = 'splPVn';
  CONTv  = [0:1:5]; CONTv = [ [-1:.1:1] [1:1:10]];
  CONT   = 1;
  showT  = 0;
  showST = 1;
 case 4
  C      = Q;
  Clon   = Qlon; Clat = Qlat; Cdpt = Qdpt;
  tit    = strcat(snapshot,'/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz)');
  cx     = [0 1]*10^(-9);
  unit   = char(1);
  titf   = 'splPV';
  unit   = '1/s/m';
  CONTv = [ [0:.025:.5] [.5:.1:1] [1:1:10]]*10^(-9);
  CONT   = 0;
  showT  = 0;
  showST = 1;
end %switch what to plot

% Find the section:
switch sectTYPE
  case 'm' % Fixed longitude
    s = max( find( Clon <= lo ) );
    C = squeeze( C(:,:,s) );
    Cx = Clat; Xlab = 'Latitude (^oN)';
    Cy = Cdpt; Ylab = 'Depth(m)';
    
    s = max( find( Tlon <= lo ) );
    sectT  = squeeze( T(:,:,s) );  sectTx = Tlat;   sectTy = Tdpt;
    s = max( find(STlon <= lo ) );
    sectST = squeeze(ST(:,:,s) ); sectSTx = STlat; sectSTy = STdpt;
  case 'z' % Fixed latitude
    s = max( find( Clat <= la ) );
    C = squeeze( C(:,s,:) );
    Cx = 360-Clon; Xlab = 'Longitude (^oW)';
    Cy = Cdpt;     Ylab = 'Depth(m)';
    
    s = max( find( Tlat <= la ) );
    sectT  = squeeze( T(:,s,:) );  sectTx = 360-Tlon;   sectTy = Tdpt;
    s = max( find(STlat <= la ) );
    sectST = squeeze(ST(:,s,:) ); sectSTx = 360-STlon; sectSTy = STdpt;
end %switch

% Plots:
clf; hold on

if CONT ~= 1
  pcolor(Cx,Cy,C); shading interp
  colormap(logcolormap(256,12));
else
  [cs,h]=contourf(Cx,Cy,C,CONTv);
  colormap(mycolormap(logcolormap(64,5),length(CONTv)));
end
caxis(cx);
ccol(ip) = colorbar(CBAR,'fontsize',10);
ctitle(ccol(ip),unit);
title(strcat('Section:',sectNAME,'/',tit));
axis([min([Cx(1) Cx(length(Cx))]) max([Cx(1) Cx(length(Cx))]) ...
      min([Cy(1) Cy(length(Cy))]) max([Cy(1) Cy(length(Cy))]) ]);
axis([55 75 -400 Cdpt(1)]);
axis([55 75 min([Cy(1) Cy(length(Cy))]) max([Cy(1) Cy(length(Cy))]) ]);
xlabel(Xlab); ylabel(Ylab);
if sectTYPE=='z', set(gca,'xdir','reverse'); end
box on

if showT
  [cs,h]=contour(sectTx,sectTy,sectT,Tcontour);
  clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200);
  for ih=1:length(h)
    set(h(ih),'edgecolor',Tcolor,'linewidth',1);
  end
end %if show THETA contours
if showST
  STcontourE=[25.5 25.75];
  [cs,h]=contour(sectSTx,sectSTy,sectST,STcontour);
  clabel(cs,h,'fontsize',8,'color',STcolor,'labelspacing',200);
  set(h,'edgecolor',STcolor);
  set(h,'linewidth',1);
  ch=get(h,'Children');
  for ich=1:length(ch)
    %set(ch(ich),'edgecolor',STcolor);
    %set(ch(ich),'linewidth',1);
    for ii = 1 : length(STcontourE)
       if get(ch(ich),'UserData') == STcontourE(ii)
	 set(ch(ich),'linewidth',2);
       end %if
    end %for ii
  end %for ich
end %if show SIGMATHETA contours
  
 
%%%%%%%%%%%%%%%%
drawnow
if prtimg
set(gcf,'color','white') 
exportj(gcf,1,strcat(outimg,sla,titf,'.',snapshot,'.section_',sectNAME));
end %if

end %for ip


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end %if pl2









if find(pl==3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Define sections:
sectTYPE = 'z'; % m:meridional or z:zonal sections

switch sectTYPE
  case 'm' % Fixed longitude
    lo = 294; xlim=[20 42]; sectNAME = '294E';
 
  case 'z' % Fixed latitude
    la = 36;  xlim=[55 75]; sectNAME = '36N';
end %switch

% Other settings:
CBAR = 'v'; % Colorbar orientation
Tcontour  = [17 19];  Tcolor = [0 0 0]; % Theta contours
STcontour = [20:.1:30]; STcolor = [1 1 1]*.5; % Sigma-Theta contours
Hcolor = [0 0 0]; % Mixed layer depth color
unit = ''; % Default unit

% Which variables to plot:
wvar = [4];

% Variables loop:
for ip=1:length(wvar)
  % Default:
  showT  = 1; % Iso-Theta contours
  showST = 0; % Iso-Sigma-Theta contours
  showH  = 0; % Mixed Layer Depth
  CONT   = 0; % Contours instead of pcolor
  
  figure(ip);clf;hold on;iw=1;jw=1;
  %if it==1, mini; end 
switch wvar(ip)
 case 1
  C     = T;
  Clon  = Tlon; Clat = Tlat; Cdpt = Tdpt;
  tit   = strcat(snapshot,' / \theta');
  cx    = [0 30];
  unit  = 'K';
  titf  = 'THETA';
 case 2
  C     = ST;
  Clon  = STlon; Clat = STlat; Cdpt = STdpt;
  tit   = strcat(snapshot,' / \sigma_\theta');
  cx    = [24 28];
  unit  = 'kg/m^3';
  titf  = 'SIGMATHETA';
 case 3
  C      = Q/Qref;
  Clon   = Qlon; Clat = Qlat; Cdpt = Qdpt;
  tit    = strcat('/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz) / q_{ref}');
  %tit    = strcat(snapshot,'/ Potential vorticity field: q = - f . d\sigma_\theta/dz / \rho');
  cx     = [0 1]*2;
  unit   = char(1);
  titf   = 'splPVn';
  CONTv  = [0:1:5]; CONTv = [ [-1:.1:1] [1:1:10]];
  CONT   = 1;
  showT  = 0;
  showST = 1;
 case 4
  C      = Q;
  Clon   = Qlon; Clat = Qlat; Cdpt = Qdpt;
  tit    = strcat('/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz)');
  cx     = [0 1]*10^(-9);
  unit   = char(1);
  titf   = 'splPV';
  unit   = '1/s/m';
  CONTv = [ [0:.025:.5] [.5:.1:1] [1:1:10]]*10^(-9);
  CONT   = 0;
  showT  = 0;
  showST = 1;
  showH  = 1;
end %switch what to plot

% Find the section:
switch sectTYPE
  case 'm' % Fixed longitude
    s = max( find( Clon <= lo ) );
    C = squeeze( C(:,:,s) );
    Cx = Clat; Xlab = 'Latitude (^oN)';
    Cy = Cdpt; Ylab = 'Depth(m)';
    
    s = max( find( Tlon <= lo ) );
    sectT  = squeeze( T(:,:,s) );  sectTx = Tlat;   sectTy = Tdpt;
    s = max( find(STlon <= lo ) );
    sectST = squeeze(ST(:,:,s) ); sectSTx = STlat; sectSTy = STdpt;
    s = max( find(MLDlon <= lo ) );
    sectH = squeeze(MLD(:,s) ); sectHx = MLDlat; sectHy = MLDdpt;
    
    dx = 1;
    s = max( find( JFzlon <= lo ) );
    sectJFz  = squeeze( JFz(:,s) );  sectJFzx = JFzlat;   sectJFzy = JFzdpt;
    s = max( find( JBzlon <= lo ) );
    sectJBz  = squeeze( JBz(:,s) );  sectJBzx = JBzlat;   sectJBzy = JBzdpt;
    
  case 'z' % Fixed latitude
    s = max( find( Clat <= la ) );
    C = squeeze( C(:,s,:) );
    Cx = 360-Clon; Xlab = 'Longitude (^oW)';
    Cy = Cdpt;     Ylab = 'Depth(m)';
    
    s = max( find( Tlat <= la ) );
    sectT  = squeeze( T(:,s,:) );  sectTx = 360-Tlon;   sectTy = Tdpt;
    s = max( find(STlat <= la ) );
    sectST = squeeze(ST(:,s,:) ); sectSTx = 360-STlon; sectSTy = STdpt;
    s = max( find(MLDlat <= la ) );
    sectH = squeeze(MLD(s,:) ); sectHx = 360-MLDlon; sectHy = MLDdpt;
       
    dx = 1;
    s = max( find( JFzlat <= la ) ); 
    sectJFz  = squeeze( JFz(s,:) );  sectJFzx = 360-JFzlon;   sectJFzy = JFzdpt;
    s = max( find( JBzlat <= la ) );
    sectJBz  = squeeze( JBz(s,:) );  sectJBzx = 360-JBzlon;   sectJBzy = JBzdpt;
end %switch

% Plots:
clf; hold on
iw=2;jw=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
sp(1)=subplot(iw,jw,1);  hold on
% Change grid:
x = sectJBzx(1:dx:end);
B = sectJBz(1:dx:end);
F = sectJFz(1:dx:end);
BpF = B+F;
col=['b','g','r']; dl = 3;
for ix = 1 : length(x)
  FLUX = [B(ix) F(ix) BpF(ix)];  
  FLUXposI = find(FLUX>=0);
  if ~isempty(FLUXposI)
  FLUXpos  = FLUX(FLUXposI);
  co = col(FLUXposI);
  [L Li] = sort(FLUXpos);
  for ii = length(L):-1:1
    line([1 1]*x(ix),[0 FLUXpos(Li(ii))],'color',co(Li(ii)),'linewidth',dl)
  end, end  
  FLUXnegI = find(FLUX<0);
  if ~isempty(FLUXnegI)
  FLUXneg  = FLUX(FLUXnegI);
  co = col(FLUXnegI);
  [L Li] = sort(FLUXneg);
  for ii = 1:length(L)
    line([1 1]*x(ix),[0 FLUXneg(Li(ii))],'color',co(Li(ii)),'linewidth',dl)
  end,  end
end
line(get(gca,'xlim'),[0 0],'color','k')
plot(x,BpF,'r');
grid on, box on
ylabel('PV flux (kg/m^3/s^2)');
set(gca,'xlim',xlim);
set(gca,'xtick',[xlim(1):2:xlim(2)]);
set(gca,'xticklabel',[]);
set(gca,'ylim',[-1 1]*1*10^(-11));
set(gca,'xdir','reverse');
title(strvcat(strcat('Section:',sectNAME,tit,...
	             ' and surface PV flux (blue: J^B_z, green: J^F_z, red: J^B_z+J^F_z)'),...
              '(gray: \sigma_\theta, black: Mixed layer depth)'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%
sp(2)=subplot(iw,jw,2);hold on

if CONT ~= 1
  pcolor(Cx,Cy,C); shading interp
  colormap(logcolormap(256,12));
else
  [cs,h]=contourf(Cx,Cy,C,CONTv);
  colormap(mycolormap(logcolormap(64,5),length(CONTv)));
end
caxis(cx); ccol(ip) = colorbar(CBAR,'fontsize',10); ctitle(ccol(ip),unit);
set(gca,'xlim',xlim);
xlabel(Xlab); ylabel(Ylab);
set(gca,'xdir','reverse');
box on
set(gca,'xtick',[xlim(1):2:xlim(2)]);
set(gca,'xticklabel',[xlim(1):2:xlim(2)]);
set(gca,'ylim',[min([Cy(1) Cy(length(Cy))]) max([Cy(1) Cy(length(Cy))]) ]);
set(gca,'ylim',[-600 max([Cy(1) Cy(length(Cy))])]);

if showT
  [cs,h]=contour(sectTx,sectTy,sectT,Tcontour);
  clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200);
  for ih=1:length(h)
    set(h(ih),'edgecolor',Tcolor,'linewidth',1);
  end
end %if show THETA contours

if showST
  STcontourE=[25.5 25.75];
  [cs,h]=contour(sectSTx,sectSTy,sectST,STcontour);
  clabel(cs,h,'fontsize',8,'color',STcolor,'labelspacing',200);
  set(h,'edgecolor',STcolor);
  set(h,'linewidth',1);
  ch=get(h,'Children');
  for ich=1:length(ch)
    %set(ch(ich),'edgecolor',STcolor);
    %set(ch(ich),'linewidth',1);
    for ii = 1 : length(STcontourE)
       if get(ch(ich),'UserData') == STcontourE(ii)
	% set(ch(ich),'linewidth',1);
       end %if
    end %for ii
  end %for ich
  [cs,h]=contour(sectSTx,sectSTy,sectST,STcontourE);
  clabel(cs,h,'fontsize',8,'color',STcolor,'labelspacing',200);
  set(h,'edgecolor',[1 1 1]);
  set(h,'linewidth',2);
  
end %if show SIGMATHETA contours

if showH
  l=line(sectHx,sectH,'color',Hcolor,'linewidth',2);
end %if show Mixed layer depth
  
pos1=get(sp(1),'position');
pos2=get(sp(2),'position');
dy = .1;
set(sp(1),'position',[pos1(1) pos1(2)+dy pos2(3) pos1(4)-dy]);
%set(sp(2),'position',[pos2(1) pos2(2)    pos2(3) pos2(4)+1.3*dy]);
set(sp(2),'position',[pos2(1) pos2(2)    pos2(3) pos1(2)+dy-1.01*pos2(2)]);



%%%%%%%%%%%%%%%%
drawnow
set(gcf,'position',[4 48 888 430]);
videotimeline(TIME,it,'b')
if prtimg
set(gcf,'color','white') 
set(gcf,'paperposition',[0.6 6.5 27.5 14]);
exportj(gcf,1,strcat(outimg,sla,titf,'.',snapshot,'.withJz.section_',sectNAME));
end %if

end %for ip


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end %if pl3






if find(pl==4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clf;
m_proj('mercator','long',subdomain.limlon,'lat',subdomain.limlat);
%m_proj('mercator','long',subdomain.limlon,'lat',[25 40]);
%m_proj('mercator','long',subdomain.limlon,'lat',[25 50]);

hold on

m_pcolor(Qnetlon,Qnetlat,Qnet);
caxis([-1 1]*500);canom

      dx = 10*diff(Txlon(1:2)); dy = 8*diff(Txlat(1:2));
      dx = 20*diff(Txlon(1:2)); dy = 10*diff(Txlat(1:2));
      lo = [Txlon(1):dx:Txlon(length(Txlon))];
      la = [Txlat(1):dy:Txlat(length(Txlat))];
      [lo la] = meshgrid(lo,la);
      Txn = interp2(Txlat,Txlon,Tx',la,lo);
      Tyn = interp2(Txlat,Txlon,Ty',la,lo);
      s = 3;
      m_quiver(lo,la,Txn,Tyn,s,'k');
      
m_coast('patch',[0 0 0]);m_grid
ccol=colorbar('h');%ctitle(ccol,'m');
title(strcat(snapshot,'/ Net heat flux (W/m^2)'));
[cs,h]=m_contour(MLDlon,MLDlat,MLD,[0:100:500]);
clabel(cs,h,'labelspacing',600,'fontsize',8);
for ih=1:length(h),set(h(ih),'edgecolor',[1 1 1],'linewidth',1);end;

M(it) = getframe;

drawnow
if prtimg
set(gcf,'color','white') 
set(gcf,'paperorientation','landscape');
%set(gcf,'paperposition',[0.6 6.5 25 14]);
print(gcf,'-djpeg100',strcat(outimg,sla,snapshot,'.jpg'));
%exportj(gcf,1,strcat(outimg,sla,snapshot));
end %if

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end %if pl4


%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
end %for it