/[MITgcm]/MITgcm_contrib/gmaze_pv/visu/eg_view_Timeserie.m
ViewVC logotype

Diff of /MITgcm_contrib/gmaze_pv/visu/eg_view_Timeserie.m

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

revision 1.1 by gmaze, Mon Jul 10 15:12:11 2006 UTC revision 1.2 by gmaze, Fri Oct 6 19:02:09 2006 UTC
# Line 1  Line 1 
1  %  %
2  % THIS IS NOT A FUNCTION !  % THIS IS NOT A FUNCTION !
3  %  %
4  % Plot time series of all variables  % Plot time series of all variables in different ways
5  % 1 var / figure  % Outputs recording possible
 % Outputs recording  
6  %  %
7    
8  clear  clear
9  global sla netcdf_domain  global sla netcdf_domain
10  pv_checkpath  pv_checkpath
11    
12  netcdf_domain = 'climode';  % Path and extension to find files:
13    pathname = strcat('netcdf-files',sla);
14    %pathname = strcat('netcdf-files-twice-daily',sla);
15    %pathname = strcat('netcdf-files-daily',sla);
16    ext      = 'nc';
17    netcdf_domain = 'western_north_atlantic';
18    
19  % Date series:  % Date series:
20  TIME = ['20030201';'20030205';'20030212';'20030219';'20030226';...  ID    = datenum(2000,12,31,12,0,0); % Start date
21          '20030301';'20030305';'20030312'];  ID    = datenum(2000,12,31,0,0,0); % Start date
22  TIME = ['20030201';'20030205';'20030212';'20030219';'20030226';...  ID    = datenum(2001,1,1,12,0,0); % Start date
23          '20030301';'20030305'];  ID    = datenum(2001,4,1,0,0,0); % Start date
24  nt   = size(TIME,1);  %IDend = datenum(2001,2,26,12,0,0); % End date
25    IDend = datenum(2001,7,4,0,0,0); % End date
26    
27    dt = datenum(0,0,1,0,0,0); % Time step between input: 1 day
28    %dt = datenum(0,0,2,0,0,0); % Time step between input: 2 days
29    %dt = datenum(0,0,7,0,0,0); % Time step between input: 1 week
30    %dt = datenum(0,0,0,12,0,0); % Time step between input: 12 hours
31    IDend = ID + 1*dt; %
32    nt = (IDend-ID)/dt;
33    
34    % Create TIME table:
35    for it = 1 : nt
36      ID = ID + dt;
37      snapshot = datestr(ID,'yyyymmddHHMM'); % For twice-daily data
38    %  snapshot = datestr(ID,'yyyymmdd'); % For daily data
39      TIME(it,:) = snapshot;
40    end %for it
41    
42    
43  iso    = 25.8; % Which sigma-theta surface ?  % Some settings
44  getiso = 1;    % We compute the isoST by default  iso    = 25.25; % Which sigma-theta surface ?
45    getiso = 0;    % We do not compute the isoST by default
46  outimg = 'img_tmp'; % Output directory  outimg = 'img_tmp'; % Output directory
47  prtimg = 1; % Do we record figures as jpg files ?  %outimg = 'img_tmp2'; % Output directory
48    %outimg = 'img_tmp3'; % Output directory
49    prtimg = 0; % Do we record figures as jpg files ?
50    
51    % Plot modules available:
52    sub = get_plotlist('eg_view_Timeserie','.');
53    disp('Available plots:')
54    sub = get_plotlistdef('eg_view_Timeserie','.');
55    disp('Set the variable <pl> in view_Timeserie.m with wanted plots')
56    
57    % Selected plots list:
58    pl = [7]; %getiso=1;
59    
60    % Verif plots:
61    disp(char(2));disp('You have choosed to plot:')
62    for i = 1 : length(pl)
63      disp(strcat(num2str(pl(i)),' -> ', sub(pl(i)).description ) )
64    end
65    s = input(' Are you sure ([y]/n) ?','s');
66    if ~isempty(s) & s == 'n'
67        return
68    end
69    
70  % Plots type:  % To find a specific date
71  % 1: individual simple maps  %find(str2num(TIME)==200103300000),break
 % 2: meridional or zonal sections  
 % 3: meridional or zonal sections with Jz vectors at the surface  
 % 4: custom it  
 pl = 3;  
72    
 if pl==2|pl==3, getiso = 0; end  
73  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74  % Video loop:  % Video loop:
75  for it = 1 : nt  for it = 1 : nt
76      snapshot = TIME(it,:);
77      %titf='.section_32N';if ~exist(strcat(outimg,sla,'PV.',snapshot,titf,'.jpg'),'file')
78        
 % Date:  
 snapshot = TIME(it,:);  
   
79  %%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%
80  % NETCDF files name:  % NETCDF files name:
81  filPV   = 'splPV';  filPV   = 'PV';
82  filST   = 'SIGMATHETA';  filST   = 'SIGMATHETA';
83  filT    = 'THETA';  filT    = 'THETA';
84  filTx   = 'TAUX';  filTx   = 'TAUX';
85  filTy   = 'TAUY';  filTy   = 'TAUY';
86  filJFz  = 'JFz';  filJFz  = 'JFz';
87  filJBz  = 'JBz';  filJBz  = 'JBz';
88  filQnet = 'Qnet';  filQnet = 'TFLUX';
89  filQEk  = 'QEk';  filQEk  = 'QEk';
90  filMLD  = 'KPPmld';  %filMLD  = 'KPPmld';
91    filMLD  = 'MLD';
92    filOx   = 'OMEGAX';
93    filOy   = 'OMEGAY';
94    filZET  = 'ZETA';
95    filEKL  = 'EKL';
96    
 % Path and extension to find them:  
 pathname = strcat('netcdf-files',sla);  
 ext      = 'nc';  
97    
98  % Load fields:  % Load fields:
99  disp('load fields...')  disp('load fields...')
# Line 65  disp('load fields...') Line 103  disp('load fields...')
103  [Qlon Qlat Qdpt] = coordfromnc(ncQ);  [Qlon Qlat Qdpt] = coordfromnc(ncQ);
104                 Q = ncQ{4}(:,:,:); clear ncQ ferfile                 Q = ncQ{4}(:,:,:); clear ncQ ferfile
105        [nz ny nx] = size(Q);        [nz ny nx] = size(Q);
106          %Qdpt = -Qdpt;
107    
108                      ferfile = strcat(pathname,sla,snapshot,sla,filZET,'.',netcdf_domain,'.',ext);
109                        ncZET = netcdf(ferfile,'nowrite');
110    [ZETAlon ZETAlat ZETAdpt] = coordfromnc(ncZET);
111                         ZETA = ncZET{4}(:,:,:); clear ncZET ferfile
112      % Move ZETA on the same grid as Q:
113      ZETA = ( ZETA(:,:,2:nx-1) + ZETA(:,:,1:nx-2) )./2;
114      ZETA = ( ZETA(:,2:ny-1,:) + ZETA(:,1:ny-2,:) )./2;
115      ZETAlon = ( ZETAlon(2:nx-1) + ZETAlon(1:nx-2) )./2;  
116      ZETAlat = ( ZETAlat(2:ny-1) + ZETAlat(1:ny-2) )./2;
117    
118                ferfile = strcat(pathname,sla,snapshot,sla,filOx,'.',netcdf_domain,'.',ext);
119                   ncOX = netcdf(ferfile,'nowrite');
120    [OXlon OXlat OXdpt] = coordfromnc(ncOX);
121                     OX = ncOX{4}(:,:,:); clear ncOX ferfile
122      % Move OMEGAx on the same grid as Q:
123      OX = ( OX(:,2:ny-1,:) + OX(:,1:ny-2,:) )./2;
124      OX = ( OX(2:nz-1,:,:) + OX(1:nz-2,:,:) )./2;
125      OXlat = ( OXlat(2:ny-1) + OXlat(1:ny-2) )./2;
126      OXdpt = ( OXdpt(2:nz-1) + OXdpt(1:nz-2) )./2;
127    
128                ferfile = strcat(pathname,sla,snapshot,sla,filOy,'.',netcdf_domain,'.',ext);
129                   ncOY = netcdf(ferfile,'nowrite');
130    [OYlon OYlat OYdpt] = coordfromnc(ncOY);
131                     OY = ncOY{4}(:,:,:); clear ncOY ferfile
132      % Move OMEGAy on the same grid as Q:
133      OY = ( OY(2:nz-1,:,:) + OY(1:nz-2,:,:) )./2;
134      OY = ( OY(:,:,2:nx-1) + OY(:,:,1:nx-2) )./2;
135      OYdpt = ( OYdpt(2:nz-1) + OYdpt(1:nz-2) )./2;
136      OYlon = ( OYlon(2:nx-1) + OYlon(1:nx-2) )./2;
137      
138      
139              ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext);              ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext);
140                 ncST = netcdf(ferfile,'nowrite');                 ncST = netcdf(ferfile,'nowrite');
141  [STlon STlat STdpt] = coordfromnc(ncST);  [STlon STlat STdpt] = coordfromnc(ncST);
# Line 76  disp('load fields...') Line 146  disp('load fields...')
146  [Tlon Tlat Tdpt] = coordfromnc(ncT);  [Tlon Tlat Tdpt] = coordfromnc(ncT);
147                 T = ncT{4}(:,:,:); clear ncT ferfile                 T = ncT{4}(:,:,:); clear ncT ferfile
148                                
149              ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext);                ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext);
150                 ncTx = netcdf(ferfile,'nowrite');                   ncTx = netcdf(ferfile,'nowrite');
151  [Txlon Txlat Txdpt] = coordfromnc(ncTx);    [Txlon Txlat Txdpt] = coordfromnc(ncTx);
152                   Tx = ncTx{4}(1,:,:); clear ncTx ferfile                     Tx = ncTx{4}(1,:,:); clear ncTx ferfile
153              ferfile = strcat(pathname,sla,snapshot,sla,filTy,'.',netcdf_domain,'.',ext);                ferfile = strcat(pathname,sla,snapshot,sla,filTy,'.',netcdf_domain,'.',ext);
154                 ncTy = netcdf(ferfile,'nowrite');                   ncTy = netcdf(ferfile,'nowrite');
155  [Tylon Tylat Tydpt] = coordfromnc(ncTy);    [Tylon Tylat Tydpt] = coordfromnc(ncTy);
156                   Ty = ncTy{4}(1,:,:); clear ncTy ferfile                     Ty = ncTy{4}(1,:,:); clear ncTy ferfile
157      
158                 ferfile = strcat(pathname,sla,snapshot,sla,filJFz,'.',netcdf_domain,'.',ext);                   ferfile = strcat(pathname,sla,snapshot,sla,filJFz,'.',netcdf_domain,'.',ext);
159                   ncJFz = netcdf(ferfile,'nowrite');                     ncJFz = netcdf(ferfile,'nowrite');
160  [JFzlon JFzlat JFzdpt] = coordfromnc(ncJFz);    [JFzlon JFzlat JFzdpt] = coordfromnc(ncJFz);
161                     JFz = ncJFz{4}(1,:,:);                       JFz = ncJFz{4}(1,:,:);
162      
163                 ferfile = strcat(pathname,sla,snapshot,sla,filJBz,'.',netcdf_domain,'.',ext);                   ferfile = strcat(pathname,sla,snapshot,sla,filJBz,'.',netcdf_domain,'.',ext);
164                   ncJBz = netcdf(ferfile,'nowrite');                     ncJBz = netcdf(ferfile,'nowrite');
165  [JBzlon JBzlat JBzdpt] = coordfromnc(ncJBz);    [JBzlon JBzlat JBzdpt] = coordfromnc(ncJBz);
166                     JBz = ncJBz{4}(1,:,:);                       JBz = ncJBz{4}(1,:,:);
167      
168                 ferfile = strcat(pathname,sla,snapshot,sla,filQnet,'.',netcdf_domain,'.',ext);                    ferfile = strcat(pathname,sla,snapshot,sla,filQnet,'.',netcdf_domain,'.',ext);
169                   ncQnet = netcdf(ferfile,'nowrite');                     ncQnet = netcdf(ferfile,'nowrite');
170  [Qnetlon Qnetlat Qnetdpt] = coordfromnc(ncQnet);  [Qnetlon Qnetlat Qnetdpt] = coordfromnc(ncQnet);
171                     Qnet = ncQnet{4}(1,:,:);                       Qnet = ncQnet{4}(1,:,:);
172    % $$$
173                 ferfile = strcat(pathname,sla,snapshot,sla,filQEk,'.',netcdf_domain,'.',ext);  % $$$                ferfile = strcat(pathname,sla,snapshot,sla,filQEk,'.',netcdf_domain,'.',ext);
174                   ncQEk = netcdf(ferfile,'nowrite');  % $$$                  ncQEk = netcdf(ferfile,'nowrite');
175  [QEklon QEklat QEkdpt] = coordfromnc(ncQEk);  % $$$ [QEklon QEklat QEkdpt] = coordfromnc(ncQEk);
176                     QEk = ncQEk{4}(1,:,:);  % $$$                    QEk = ncQEk{4}(1,:,:);
177    % $$$
178                 ferfile = strcat(pathname,sla,snapshot,sla,filMLD,'.',netcdf_domain,'.',ext);                   ferfile = strcat(pathname,sla,snapshot,sla,filMLD,'.',netcdf_domain,'.',ext);
179                   ncMLD = netcdf(ferfile,'nowrite');                     ncMLD = netcdf(ferfile,'nowrite');
180  [MLDlon MLDlat MLDdpt] = coordfromnc(ncMLD);    [MLDlon MLDlat MLDdpt] = coordfromnc(ncMLD);
181                     MLD = -ncMLD{4}(1,:,:);                       MLD = ncMLD{4}(1,:,:);
182      
183                     ferfile = strcat(pathname,sla,snapshot,sla,filEKL,'.',netcdf_domain,'.',ext);
184                       ncEKL = netcdf(ferfile,'nowrite');
185      [EKLlon EKLlat EKLdpt] = coordfromnc(ncEKL);
186                         EKL = ncEKL{4}(1,:,:);
187    
188                  
189  %%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%
190  % Q is defined on the same grid of ST but troncated by extrem 2 points, then here  % Q is defined on the same grid of ST but troncated by extrem 2 points, then here
191  % make all fields defined with same limits...  % make all fields defined with same limits...
192    % In case of missing points, we add NaN.
193  disp('Reshape them')  disp('Reshape them')
194  ST    = squeeze(ST(2:nz+1,2:ny+1,2:nx+1));  ST    = squeeze(ST(2:nz+1,2:ny+1,2:nx+1));
195  STdpt = STdpt(2:nz+1);  STdpt = STdpt(2:nz+1);
# Line 131  Qnetlat = Qnetlat(2:ny+1); Line 208  Qnetlat = Qnetlat(2:ny+1);
208  MLD    = squeeze(MLD(2:ny+1,2:nx+1));  MLD    = squeeze(MLD(2:ny+1,2:nx+1));
209  MLDlon = MLDlon(2:nx+1);  MLDlon = MLDlon(2:nx+1);
210  MLDlat = MLDlat(2:ny+1);  MLDlat = MLDlat(2:ny+1);
211    EKL    = squeeze(EKL(2:ny+1,2:nx+1));
212    EKLlon = EKLlon(2:nx+1);
213    EKLlat = EKLlat(2:ny+1);
214    ZETA = squeeze(ZETA(2:nz+1,:,:));
215    ZETA = cat(2,ZETA,ones(size(ZETA,1),1,size(ZETA,3)).*NaN);
216    ZETA = cat(2,ones(size(ZETA,1),1,size(ZETA,3)).*NaN,ZETA);
217    ZETA = cat(3,ZETA,ones(size(ZETA,1),size(ZETA,2),1).*NaN);
218    ZETA = cat(3,ones(size(ZETA,1),size(ZETA,2),1).*NaN,ZETA);
219    ZETAdpt = ZETAdpt(2:nz+1);
220    ZETAlon = STlon;
221    ZETAlat = STlat;
222    OX = squeeze(OX(:,:,2:nx+1));
223    OX = cat(1,OX,ones(1,size(OX,2),size(OX,3)).*NaN);
224    OX = cat(1,ones(1,size(OX,2),size(OX,3)).*NaN,OX);
225    OX = cat(2,OX,ones(size(OX,1),1,size(OX,3)).*NaN);
226    OX = cat(2,ones(size(OX,1),1,size(OX,3)).*NaN,OX);
227    OXlon = STlon;
228    OXlat = STlat;
229    OXdpt = STdpt;
230    OY = squeeze(OY(:,2:ny+1,:));
231    OY = cat(1,OY,ones(1,size(OY,2),size(OY,3)).*NaN);
232    OY = cat(1,ones(1,size(OY,2),size(OY,3)).*NaN,OY);
233    OY = cat(3,OY,ones(size(OY,1),size(OY,2),1).*NaN);
234    OY = cat(3,ones(size(OY,1),size(OY,2),1).*NaN,OY);
235    OYlon = STlon;
236    OYlat = STlat;
237    OYdpt = STdpt;
238    
239    
240    % Planetary vorticity:
241      f = 2*(2*pi/86400)*sin(ZETAlat*pi/180);
242      [a f c]=meshgrid(ZETAlon,f,ZETAdpt); clear a c
243      f = permute(f,[3 1 2]);
244      
245    % Apply mask:
246    MASK = ones(size(ST,1),size(ST,2),size(ST,3));
247    MASK(find(isnan(ST))) = NaN;
248    T = T.*MASK;
249    Qnet = Qnet.*squeeze(MASK(1,:,:));
250    
251        
252  % Grid:  % Grid:
253  global domain subdomain2 subdomain  global domain subdomain1 subdomain2 subdomain3
254  grid_setup  grid_setup
255  %subdomain = subdomain2;  subdomain = subdomain1;
256    
257    
258  %%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%
# Line 147  Diso = ones(size(Iiso)).*NaN; Line 264  Diso = ones(size(Iiso)).*NaN;
264  Qiso = ones(size(Iiso)).*NaN;  Qiso = ones(size(Iiso)).*NaN;
265  for ix = 1 : size(ST,3)  for ix = 1 : size(ST,3)
266    for iy = 1 : size(ST,2)    for iy = 1 : size(ST,2)
267      if ~isnan(Iiso(iy,ix))      if ~isnan(Iiso(iy,ix)) & ~isnan( Q(Iiso(iy,ix),iy,ix) )
268         Diso(iy,ix) = STdpt(Iiso(iy,ix));         Diso(iy,ix) = STdpt(Iiso(iy,ix));
269         Qiso(iy,ix) =     Q(Iiso(iy,ix),iy,ix);         Qiso(iy,ix) =     Q(Iiso(iy,ix),iy,ix);
270      end %if      end %if
271  end, end %for iy, ix  end, end %for iy, ix
272  end %if  end %if
273    
274    
275    
276  %%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%
277  % "Normalise" the PV:  % "Normalise" the PV:
278  fO  = 2*(2*pi/86400)*sin(32*pi/180);  fO  = 2*(2*pi/86400)*sin(32*pi/180);
# Line 169  if getiso, QisoN = Qiso./Qref; end Line 288  if getiso, QisoN = Qiso./Qref; end
288  % Plots:  % Plots:
289  disp('Plots ...')  disp('Plots ...')
290    
 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  
   
   
291    
292  if find(pl==2)  for i = 1 : length(pl)
293  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    disp(strcat('Plotting module:',sub(pl(i)).name))
294  % Define sections:    eval(sub(pl(i)).name(1:end-2),'disp(''Oups scratch...'');return');
 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)));  
295  end  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)]);  
   
296    
297    
298  %%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%
299  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  
300    
301  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    %else,disp(strcat('Skip:',snapshot));end
 end %if pl4  
302    
303    fclose('all');
304    
 %%%%%%%%%%%%%%%%  
 %%%%%%%%%%%%%%%%  
305    
     
306  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
307  end %for it    end %for it  
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22