gmaze_pv/visu/eg_view_Timeserie.m

%
% THIS IS NOT A FUNCTION !
%
% Plot time series of all variables in different ways
% Outputs recording possible
%

clear
global sla netcdf_domain
pv_checkpath

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

% Date series:
ID    = datenum(2000,12,31,12,0,0); % Start date
ID    = datenum(2000,12,31,0,0,0); % Start date
ID    = datenum(2001,1,1,12,0,0); % Start date
ID    = datenum(2001,4,1,0,0,0); % Start date
%IDend = datenum(2001,2,26,12,0,0); % End date
IDend = datenum(2001,7,4,0,0,0); % End date

dt = datenum(0,0,1,0,0,0); % Time step between input: 1 day
%dt = datenum(0,0,2,0,0,0); % Time step between input: 2 days
%dt = datenum(0,0,7,0,0,0); % Time step between input: 1 week
%dt = datenum(0,0,0,12,0,0); % Time step between input: 12 hours
IDend = ID + 1*dt; % 
nt = (IDend-ID)/dt;

% Create TIME table:
for it = 1 : nt
  ID = ID + dt;
  snapshot = datestr(ID,'yyyymmddHHMM'); % For twice-daily data
%  snapshot = datestr(ID,'yyyymmdd'); % For daily data
  TIME(it,:) = snapshot;
end %for it


% Some settings
iso    = 25.25; % Which sigma-theta surface ?
getiso = 0;    % We do not compute the isoST by default
outimg = 'img_tmp'; % Output directory
%outimg = 'img_tmp2'; % Output directory
%outimg = 'img_tmp3'; % Output directory
prtimg = 0; % Do we record figures as jpg files ?

% Plot modules available:
sub = get_plotlist('eg_view_Timeserie','.');
disp('Available plots:')
sub = get_plotlistdef('eg_view_Timeserie','.');
disp('Set the variable <pl> in view_Timeserie.m with wanted plots')

% Selected plots list:
pl = [7]; %getiso=1;

% Verif plots:
disp(char(2));disp('You have choosed to plot:')
for i = 1 : length(pl)
  disp(strcat(num2str(pl(i)),' -> ', sub(pl(i)).description ) )
end
s = input(' Are you sure ([y]/n) ?','s');
if ~isempty(s) & s == 'n'
    return
end

% To find a specific date
%find(str2num(TIME)==200103300000),break

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Video loop:
for it = 1 : nt
  snapshot = TIME(it,:);
  %titf='.section_32N';if ~exist(strcat(outimg,sla,'PV.',snapshot,titf,'.jpg'),'file')
  
%%%%%%%%%%%%%%%%
% NETCDF files name:
filPV   = 'PV';
filST   = 'SIGMATHETA';
filT    = 'THETA';
filTx   = 'TAUX';
filTy   = 'TAUY';
filJFz  = 'JFz';
filJBz  = 'JBz';
filQnet = 'TFLUX';
filQEk  = 'QEk';
%filMLD  = 'KPPmld';
filMLD  = 'MLD';
filOx   = 'OMEGAX';
filOy   = 'OMEGAY';
filZET  = 'ZETA';
filEKL  = 'EKL';


% 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);
      %Qdpt = -Qdpt;

                  ferfile = strcat(pathname,sla,snapshot,sla,filZET,'.',netcdf_domain,'.',ext);
                    ncZET = netcdf(ferfile,'nowrite');
[ZETAlon ZETAlat ZETAdpt] = coordfromnc(ncZET);
                     ZETA = ncZET{4}(:,:,:); clear ncZET ferfile
  % Move ZETA on the same grid as Q:
  ZETA = ( ZETA(:,:,2:nx-1) + ZETA(:,:,1:nx-2) )./2;
  ZETA = ( ZETA(:,2:ny-1,:) + ZETA(:,1:ny-2,:) )./2;
  ZETAlon = ( ZETAlon(2:nx-1) + ZETAlon(1:nx-2) )./2;   
  ZETAlat = ( ZETAlat(2:ny-1) + ZETAlat(1:ny-2) )./2;

            ferfile = strcat(pathname,sla,snapshot,sla,filOx,'.',netcdf_domain,'.',ext);
               ncOX = netcdf(ferfile,'nowrite');
[OXlon OXlat OXdpt] = coordfromnc(ncOX);
                 OX = ncOX{4}(:,:,:); clear ncOX ferfile
  % Move OMEGAx on the same grid as Q:
  OX = ( OX(:,2:ny-1,:) + OX(:,1:ny-2,:) )./2;
  OX = ( OX(2:nz-1,:,:) + OX(1:nz-2,:,:) )./2;
  OXlat = ( OXlat(2:ny-1) + OXlat(1:ny-2) )./2;
  OXdpt = ( OXdpt(2:nz-1) + OXdpt(1:nz-2) )./2;

            ferfile = strcat(pathname,sla,snapshot,sla,filOy,'.',netcdf_domain,'.',ext);
               ncOY = netcdf(ferfile,'nowrite');
[OYlon OYlat OYdpt] = coordfromnc(ncOY);
                 OY = ncOY{4}(:,:,:); clear ncOY ferfile
  % Move OMEGAy on the same grid as Q:
  OY = ( OY(2:nz-1,:,:) + OY(1:nz-2,:,:) )./2;
  OY = ( OY(:,:,2:nx-1) + OY(:,:,1:nx-2) )./2;
  OYdpt = ( OYdpt(2:nz-1) + OYdpt(1:nz-2) )./2;
  OYlon = ( OYlon(2:nx-1) + OYlon(1:nx-2) )./2;
  
  
            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,:,:);
  
                 ferfile = strcat(pathname,sla,snapshot,sla,filEKL,'.',netcdf_domain,'.',ext);
                   ncEKL = netcdf(ferfile,'nowrite');
  [EKLlon EKLlat EKLdpt] = coordfromnc(ncEKL);
                     EKL = ncEKL{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...
% In case of missing points, we add NaN.
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);
EKL    = squeeze(EKL(2:ny+1,2:nx+1));
EKLlon = EKLlon(2:nx+1);
EKLlat = EKLlat(2:ny+1);
ZETA = squeeze(ZETA(2:nz+1,:,:));
ZETA = cat(2,ZETA,ones(size(ZETA,1),1,size(ZETA,3)).*NaN);
ZETA = cat(2,ones(size(ZETA,1),1,size(ZETA,3)).*NaN,ZETA);
ZETA = cat(3,ZETA,ones(size(ZETA,1),size(ZETA,2),1).*NaN);
ZETA = cat(3,ones(size(ZETA,1),size(ZETA,2),1).*NaN,ZETA);
ZETAdpt = ZETAdpt(2:nz+1);
ZETAlon = STlon;
ZETAlat = STlat;
OX = squeeze(OX(:,:,2:nx+1));
OX = cat(1,OX,ones(1,size(OX,2),size(OX,3)).*NaN);
OX = cat(1,ones(1,size(OX,2),size(OX,3)).*NaN,OX);
OX = cat(2,OX,ones(size(OX,1),1,size(OX,3)).*NaN);
OX = cat(2,ones(size(OX,1),1,size(OX,3)).*NaN,OX);
OXlon = STlon;
OXlat = STlat;
OXdpt = STdpt;
OY = squeeze(OY(:,2:ny+1,:));
OY = cat(1,OY,ones(1,size(OY,2),size(OY,3)).*NaN);
OY = cat(1,ones(1,size(OY,2),size(OY,3)).*NaN,OY);
OY = cat(3,OY,ones(size(OY,1),size(OY,2),1).*NaN);
OY = cat(3,ones(size(OY,1),size(OY,2),1).*NaN,OY);
OYlon = STlon;
OYlat = STlat;
OYdpt = STdpt;


% Planetary vorticity:
  f = 2*(2*pi/86400)*sin(ZETAlat*pi/180);
  [a f c]=meshgrid(ZETAlon,f,ZETAdpt); clear a c
  f = permute(f,[3 1 2]);
  
% Apply mask:
MASK = ones(size(ST,1),size(ST,2),size(ST,3)); 
MASK(find(isnan(ST))) = NaN;
T = T.*MASK;
Qnet = Qnet.*squeeze(MASK(1,:,:));

  
% Grid:
global domain subdomain1 subdomain2 subdomain3
grid_setup
subdomain = subdomain1;


%%%%%%%%%%%%%%%%
% 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)) & ~isnan( Q(Iiso(iy,ix),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 ...')


for i = 1 : length(pl)
  disp(strcat('Plotting module:',sub(pl(i)).name))
  eval(sub(pl(i)).name(1:end-2),'disp(''Oups scratch...'');return');
end


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

  %else,disp(strcat('Skip:',snapshot));end

fclose('all');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
end %for it  
1	gmaze	1.1	%
2			% THIS IS NOT A FUNCTION !
3			%
4	gmaze	1.2	% Plot time series of all variables in different ways
5			% Outputs recording possible
6	gmaze	1.1	%
7
8			clear
9			global sla netcdf_domain
10			pv_checkpath
11
12	gmaze	1.2	% 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	gmaze	1.1
19			% Date series:
20	gmaze	1.2	ID = datenum(2000,12,31,12,0,0); % Start date
21			ID = datenum(2000,12,31,0,0,0); % Start date
22			ID = datenum(2001,1,1,12,0,0); % Start date
23			ID = datenum(2001,4,1,0,0,0); % Start date
24			%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	gmaze	1.1
42
43	gmaze	1.2	% Some settings
44			iso = 25.25; % Which sigma-theta surface ?
45			getiso = 0; % We do not compute the isoST by default
46	gmaze	1.1	outimg = 'img_tmp'; % Output directory
47	gmaze	1.2	%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	gmaze	1.1
70	gmaze	1.2	% To find a specific date
71			%find(str2num(TIME)==200103300000),break
72	gmaze	1.1
73			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74			% Video loop:
75			for it = 1 : nt
76	gmaze	1.2	snapshot = TIME(it,:);
77			%titf='.section_32N';if ~exist(strcat(outimg,sla,'PV.',snapshot,titf,'.jpg'),'file')
78	gmaze	1.1
79			%%%%%%%%%%%%%%%%
80			% NETCDF files name:
81	gmaze	1.2	filPV = 'PV';
82	gmaze	1.1	filST = 'SIGMATHETA';
83			filT = 'THETA';
84			filTx = 'TAUX';
85			filTy = 'TAUY';
86			filJFz = 'JFz';
87			filJBz = 'JBz';
88	gmaze	1.2	filQnet = 'TFLUX';
89	gmaze	1.1	filQEk = 'QEk';
90	gmaze	1.2	%filMLD = 'KPPmld';
91			filMLD = 'MLD';
92			filOx = 'OMEGAX';
93			filOy = 'OMEGAY';
94			filZET = 'ZETA';
95			filEKL = 'EKL';
96	gmaze	1.1
97
98			% Load fields:
99			disp('load fields...')
100			% (I keep proper axis for each variables in case of one day they would be different)
101			ferfile = strcat(pathname,sla,snapshot,sla,filPV,'.',netcdf_domain,'.',ext);
102			ncQ = netcdf(ferfile,'nowrite');
103			[Qlon Qlat Qdpt] = coordfromnc(ncQ);
104			Q = ncQ{4}(:,:,:); clear ncQ ferfile
105			[nz ny nx] = size(Q);
106	gmaze	1.2	%Qdpt = -Qdpt;
107	gmaze	1.1
108	gmaze	1.2	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	gmaze	1.1	ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext);
140			ncST = netcdf(ferfile,'nowrite');
141			[STlon STlat STdpt] = coordfromnc(ncST);
142			ST = ncST{4}(:,:,:); clear ncST ferfile
143
144			ferfile = strcat(pathname,sla,snapshot,sla,filT,'.',netcdf_domain,'.',ext);
145			ncT = netcdf(ferfile,'nowrite');
146			[Tlon Tlat Tdpt] = coordfromnc(ncT);
147			T = ncT{4}(:,:,:); clear ncT ferfile
148
149	gmaze	1.2	ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext);
150			ncTx = netcdf(ferfile,'nowrite');
151			[Txlon Txlat Txdpt] = coordfromnc(ncTx);
152			Tx = ncTx{4}(1,:,:); clear ncTx ferfile
153			ferfile = strcat(pathname,sla,snapshot,sla,filTy,'.',netcdf_domain,'.',ext);
154			ncTy = netcdf(ferfile,'nowrite');
155			[Tylon Tylat Tydpt] = coordfromnc(ncTy);
156			Ty = ncTy{4}(1,:,:); clear ncTy ferfile
157
158			ferfile = strcat(pathname,sla,snapshot,sla,filJFz,'.',netcdf_domain,'.',ext);
159			ncJFz = netcdf(ferfile,'nowrite');
160			[JFzlon JFzlat JFzdpt] = coordfromnc(ncJFz);
161			JFz = ncJFz{4}(1,:,:);
162
163			ferfile = strcat(pathname,sla,snapshot,sla,filJBz,'.',netcdf_domain,'.',ext);
164			ncJBz = netcdf(ferfile,'nowrite');
165			[JBzlon JBzlat JBzdpt] = coordfromnc(ncJBz);
166			JBz = ncJBz{4}(1,:,:);
167
168			ferfile = strcat(pathname,sla,snapshot,sla,filQnet,'.',netcdf_domain,'.',ext);
169			ncQnet = netcdf(ferfile,'nowrite');
170	gmaze	1.1	[Qnetlon Qnetlat Qnetdpt] = coordfromnc(ncQnet);
171	gmaze	1.2	Qnet = ncQnet{4}(1,:,:);
172			% $$$
173			% $$$ ferfile = strcat(pathname,sla,snapshot,sla,filQEk,'.',netcdf_domain,'.',ext);
174			% $$$ ncQEk = netcdf(ferfile,'nowrite');
175			% $$$ [QEklon QEklat QEkdpt] = coordfromnc(ncQEk);
176			% $$$ QEk = ncQEk{4}(1,:,:);
177			% $$$
178			ferfile = strcat(pathname,sla,snapshot,sla,filMLD,'.',netcdf_domain,'.',ext);
179			ncMLD = netcdf(ferfile,'nowrite');
180			[MLDlon MLDlat MLDdpt] = coordfromnc(ncMLD);
181			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	gmaze	1.1
188	gmaze	1.2
189	gmaze	1.1	%%%%%%%%%%%%%%%%
190			% 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...
192	gmaze	1.2	% In case of missing points, we add NaN.
193	gmaze	1.1	disp('Reshape them')
194			ST = squeeze(ST(2:nz+1,2:ny+1,2:nx+1));
195			STdpt = STdpt(2:nz+1);
196			STlon = STlon(2:nx+1);
197			STlat = STlat(2:ny+1);
198			T = squeeze(T(2:nz+1,2:ny+1,2:nx+1));
199			Tdpt = Tdpt(2:nz+1);
200			Tlon = Tlon(2:nx+1);
201			Tlat = Tlat(2:ny+1);
202			JBz = squeeze(JBz(2:ny+1,2:nx+1));
203			JBzlon = JBzlon(2:nx+1);
204			JBzlat = JBzlat(2:ny+1);
205			Qnet = squeeze(Qnet(2:ny+1,2:nx+1));
206			Qnetlon = Qnetlon(2:nx+1);
207			Qnetlat = Qnetlat(2:ny+1);
208			MLD = squeeze(MLD(2:ny+1,2:nx+1));
209			MLDlon = MLDlon(2:nx+1);
210			MLDlat = MLDlat(2:ny+1);
211	gmaze	1.2	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(2pi/86400)sin(ZETAlatpi/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	gmaze	1.1
252			% Grid:
253	gmaze	1.2	global domain subdomain1 subdomain2 subdomain3
254	gmaze	1.1	grid_setup
255	gmaze	1.2	subdomain = subdomain1;
256	gmaze	1.1
257
258			%%%%%%%%%%%%%%%%
259			% Here we determine the isosurface and its depth:
260			if getiso
261			disp('Get iso-ST')
262			[Iiso mask] = subfct_getisoS(ST,iso);
263			Diso = ones(size(Iiso)).*NaN;
264			Qiso = ones(size(Iiso)).*NaN;
265			for ix = 1 : size(ST,3)
266			for iy = 1 : size(ST,2)
267	gmaze	1.2	if ~isnan(Iiso(iy,ix)) & ~isnan( Q(Iiso(iy,ix),iy,ix) )
268	gmaze	1.1	Diso(iy,ix) = STdpt(Iiso(iy,ix));
269			Qiso(iy,ix) = Q(Iiso(iy,ix),iy,ix);
270			end %if
271			end, end %for iy, ix
272			end %if
273
274	gmaze	1.2
275
276	gmaze	1.1	%%%%%%%%%%%%%%%%
277			% "Normalise" the PV:
278			fO = 2(2pi/86400)sin(32pi/180);
279			dST = 27.6-25.4;
280			H = -1000;
281			RHOo = 1000;
282			Qref = -fO/RHOo*dST/H;
283			if getiso, QisoN = Qiso./Qref; end
284
285
286			%%%%%%%%%%%%%%%%
287			%%%%%%%%%%%%%%%%
288			% Plots:
289			disp('Plots ...')
290
291	gmaze	1.2
292			for i = 1 : length(pl)
293			disp(strcat('Plotting module:',sub(pl(i)).name))
294			eval(sub(pl(i)).name(1:end-2),'disp(''Oups scratch...'');return');
295	gmaze	1.1	end
296
297
298			%%%%%%%%%%%%%%%%
299			%%%%%%%%%%%%%%%%
300
301	gmaze	1.2	%else,disp(strcat('Skip:',snapshot));end
302	gmaze	1.1
303	gmaze	1.2	fclose('all');
304	gmaze	1.1
305
306			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
307			end %for it