MITgcm_contrib/gmaze_pv/C_compute_potential_vorticity.m

%
% [Q] = C_compute_potential_vorticity(SNAPSHOT,[WANTSPLPV])
% [Q1,Q2,Q3] = C_compute_potential_vorticity(SNAPSHOT,[WANTSPLPV])
%
% This file computes the potential vorticity Q from
% netcdf files of relative vorticity (OMEGAX, OMEGAY, ZETA)
% and potential density (SIGMATHETA) as
% Q = OMEGAX . dSIGMATHETA/dx + OMEGAY . dSIGMATHETA/dy + (f+ZETA).dSIGMATHETA/dz
% 
% The optional flag WANTSPLPV is set to 0 by defaut. If turn to 1,
% then the program computes the simple PV defined by:
% splQ = f.dSIGMATHETA/dz
%
% Note that none of the fields are defined on the same grid points.
% So, I decided to compute Q on the same grid as SIGMATHETA, ie. the 
% center of the c-grid.
%
% Files names are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/OMEGAX.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/OMEGAY.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/ZETA.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/SIGMATHETA.<netcdf_domain>.<netcdf_suff>
% OUPUT:
% ./netcdf-files/<SNAPSHOT>/PV.<netcdf_domain>.<netcdf_suff>
% or 
% ./netcdf-files/<SNAPSHOT>/splPV.<netcdf_domain>.<netcdf_suff>
%
% 06/07/2006
% gmaze@mit.edu
%
  
function varargout = C_compute_potential_vorticity(snapshot,varargin)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global sla netcdf_domain netcdf_suff
pv_checkpath

%% Flags to choose which term to compute (by default, all):
FLpv1 = 1;
FLpv2 = 1;
FLpv3 = 1;
if nargin==2  % case of optional flag presents:
  if varargin{1}(1) == 1 % Case of the simple PV:
    FLpv1 = 0;
    FLpv2 = 0;
    FLpv3 = 2;
  end
end %if
%[FLpv1 FLpv2 FLpv3]


%% Optionnal flags:
global toshow % Turn to 1 to follow the computing process


%% NETCDF files:

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

% Names:
if FLpv3 ~= 2 % We don't need them for splPV
  filOx = strcat('OMEGAX'    ,'.',netcdf_domain);
  filOy = strcat('OMEGAY'    ,'.',netcdf_domain);
  filOz = strcat('ZETA'      ,'.',netcdf_domain);
end %if
  filST = strcat('SIGMATHETA','.',netcdf_domain);

% Load files and coordinates:
if FLpv3 ~= 2 % We don't need them for splPV
  ferfile             = strcat(pathname,sla,filOx,ext);
  ncOx                = netcdf(ferfile,'nowrite');
  [Oxlon Oxlat Oxdpt] = coordfromnc(ncOx);
  ferfile             = strcat(pathname,sla,filOy,ext);
  ncOy                = netcdf(ferfile,'nowrite');
  [Oylon Oylat Oydpt] = coordfromnc(ncOy);
  ferfile             = strcat(pathname,sla,filOz,ext);
  ncOz                = netcdf(ferfile,'nowrite');
  [Ozlon Ozlat Ozdpt] = coordfromnc(ncOz);
end %if
  ferfile             = strcat(pathname,sla,filST,ext);
  ncST                = netcdf(ferfile,'nowrite');
  [STlon STlat STdpt] = coordfromnc(ncST);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Then, compute the first term:  OMEGAX . dSIGMATHETA/dx  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if FLpv1
  
%%%%%  
%% 1: Compute zonal gradient of SIGMATHETA:

% Dim:
if toshow,disp('dim'),end
nx = length(STlon) - 1;
ny = length(STlat);
nz = length(STdpt);

% Pre-allocate:
if toshow,disp('pre-allocate'),end
dSIGMATHETAdx = zeros(nz,ny,nx-1)*NaN;
           dx = zeros(1,nx).*NaN;
         STup = zeros(nz,nx);
         STdw = zeros(nz,nx);

% Zonal gradient of SIGMATHETA:
if toshow,disp('grad'), end
for iy = 1 : ny
  if toshow
    disp(strcat('Computing dSIGMATHETA/dx at latitude : ',num2str(STlat(iy)),...
                '^o (',num2str(iy),'/',num2str(ny),')'   ));
  end  
  [dx b] = meshgrid( m_lldist(STlon(1:nx+1),[1 1]*STlat(iy)), STdpt ) ; clear b
  STup   = squeeze(ncST{4}(:,iy,2:nx+1));
  STdw   = squeeze(ncST{4}(:,iy,1:nx));
  dSTdx  = ( STup - STdw ) ./ dx;
  % Change horizontal grid point definition to fit with SIGMATHETA:
  dSTdx  = ( dSTdx(:,1:nx-1) + dSTdx(:,2:nx) )./2; 
  dSIGMATHETAdx(:,iy,:) = dSTdx;
end %for iy


%%%%%
%% 2: Move OMEGAX on the same grid:
if toshow,disp('Move OMEGAX on the same grid as dSIGMATHETA/dx'), end

% Change vertical gridding of OMEGAX:
Ox = ncOx{4}(:,:,:);
Ox = ( Ox(2:nz-1,:,:) + Ox(1:nz-2,:,:) )./2;
% And horizontal gridding:
Ox = ( Ox(:,2:ny-1,:) + Ox(:,1:ny-2,:) )./2;

%%%%%
%% 3: Make both fields having same limits:
%%    (Keep points where both fields are defined)
           Ox = squeeze(Ox(:,:,2:nx));
dSIGMATHETAdx = squeeze( dSIGMATHETAdx (2:nz-1,2:ny-1,:) );

%%%%%
%% 4: Last, compute first term of PV:
PV1 = Ox.*dSIGMATHETAdx ; 

% and define axis fron the ST grid:
PV1_lon = STlon(2:length(STlon)-1);
PV1_lat = STlat(2:length(STlat)-1);
PV1_dpt = STdpt(2:length(STdpt)-1);

clear nx ny nz dx STup STdw iy dSTdx Ox dSIGMATHETAdx
end %if FLpv1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the second term:  OMEGAY . dSIGMATHETA/dy  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if FLpv2
  
%%%%%  
%% 1: Compute meridional gradient of SIGMATHETA:

% Dim:
if toshow,disp('dim'), end
nx = length(STlon) ;
ny = length(STlat) - 1 ;
nz = length(STdpt) ;

% Pre-allocate:
if toshow,disp('pre-allocate'), end
dSIGMATHETAdy = zeros(nz,ny-1,nx).*NaN;
           dy = zeros(1,ny).*NaN;
         STup = zeros(nz,ny);
         STdw = zeros(nz,ny);

% Meridional gradient of SIGMATHETA:
% (Assuming the grid is regular, dy is independent of x)
[dy b] = meshgrid( m_lldist([1 1]*STlon(1),STlat(1:ny+1) ), STdpt ) ; clear b
for ix = 1 : nx
  if toshow
    disp(strcat('Computing dSIGMATHETA/dy at longitude : ',num2str(STlon(ix)),...
                '^o (',num2str(ix),'/',num2str(nx),')'   ));
  end
  STup  = squeeze(ncST{4}(:,2:ny+1,ix));
  STdw  = squeeze(ncST{4}(:,1:ny,ix));
  dSTdy = ( STup - STdw ) ./ dy;
  % Change horizontal grid point definition to fit with SIGMATHETA:
  dSTdy = ( dSTdy(:,1:ny-1) + dSTdy(:,2:ny) )./2; 
  dSIGMATHETAdy(:,:,ix) = dSTdy;
end %for iy

%%%%%
%% 2: Move OMEGAY on the same grid:
if toshow,disp('Move OMEGAY on the same grid as dSIGMATHETA/dy'), end

% Change vertical gridding of OMEGAY:
Oy = ncOy{4}(:,:,:);
Oy = ( Oy(2:nz-1,:,:) + Oy(1:nz-2,:,:) )./2;
% And horizontal gridding:
Oy = ( Oy(:,:,2:nx-1) + Oy(:,:,1:nx-2) )./2;

%%%%%
%% 3: Make them having same limits:
%%    (Keep points where both fields are defined)
           Oy = squeeze(Oy(:,2:ny,:));
dSIGMATHETAdy = squeeze( dSIGMATHETAdy (2:nz-1,:,2:nx-1) );

%%%%%
%% 4: Last, compute second term of PV:
PV2 = Oy.*dSIGMATHETAdy ; 

% and defined axis fron the ST grid:
PV2_lon = STlon(2:length(STlon)-1);
PV2_lat = STlat(2:length(STlat)-1);
PV2_dpt = STdpt(2:length(STdpt)-1);


clear nx ny nz dy STup STdw dy dSTdy Oy dSIGMATHETAdy
end %if FLpv2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the third term: ( f + ZETA ) . dSIGMATHETA/dz  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if FLpv3

%%%%%
%% 1: Compute vertical gradient of SIGMATHETA:

% Dim:
if toshow,disp('dim'), end
nx = length(STlon) ;
ny = length(STlat) ;
nz = length(STdpt) - 1 ;

% Pre-allocate:
if toshow,disp('pre-allocate'), end
dSIGMATHETAdz = zeros(nz-1,ny,nx).*NaN;
           ST = zeros(nz+1,ny,nx);
           dz = zeros(1,nz).*NaN;

% Vertical grid differences:
% STdpt contains negative values with STdpt(1) at the surface
% and STdpt(end) at the bottom of the ocean.
% So dz is positive with respect to z axis upward:
         dz = -diff(STdpt); 
[a dz_3D c] = meshgrid(STlat,dz,STlon); clear a c

% Vertical gradient:
if toshow,disp('Vertical gradient of SIGMATHETA'), end
           ST = ncST{4}(:,:,:);
           % Z axis upward, so vertical derivative is upper-part
           % minus lower-part:
dSIGMATHETAdz = ( ST(1:nz,:,:) - ST(2:nz+1,:,:) ) ./ dz_3D;
clear dz_3D ST

% Change vertical gridding:
dSIGMATHETAdz = ( dSIGMATHETAdz(1:nz-1,:,:) + dSIGMATHETAdz(2:nz,:,:) )./2;

if FLpv3 == 1 % Just for full PV
  
  %%%%%
  %% 2: Move ZETA on the same grid:
  if toshow,disp('Move ZETA on the same grid as dSIGMATHETA/dz'), end
  Oz = ncOz{4}(:,:,:);
  % Change horizontal gridding:
  Oz = ( Oz(:,:,2:nx-1) + Oz(:,:,1:nx-2) )./2;
  Oz = ( Oz(:,2:ny-1,:) + Oz(:,1:ny-2,:) )./2;

end %if FLpv3=1

%%%%%
%% 3: Make them having same limits:
%%    (Keep points where both fields are defined)
if FLpv3 == 1
           Oz = squeeze(Oz(2:nz,:,:));     
end %if    
dSIGMATHETAdz = squeeze( dSIGMATHETAdz (:,2:ny-1,2:nx-1) );


%%%%%
%% 4: Last, compute third term of PV:
% and defined axis fron the ST grid:
PV3_lon = STlon(2:length(STlon)-1);
PV3_lat = STlat(2:length(STlat)-1);
PV3_dpt = STdpt(2:length(STdpt)-1);

% Planetary vorticity:
f = 2*(2*pi/86400)*sin(PV3_lat*pi/180);
[a f c]=meshgrid(PV3_lon,f,PV3_dpt); clear a c
f = permute(f,[3 1 2]);

% Third term of PV:
if FLpv3 == 2
  % Compute simple PV, just with planetary vorticity:
  PV3 = f.*dSIGMATHETAdz ;
else
  % To compute full PV:
  PV3 = (f+Oz).*dSIGMATHETAdz ; 
end
 

clear nx ny nz dz ST Oz dSIGMATHETAdz f
end %if FLpv3


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Then, compute potential vorticity:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if toshow,disp('Summing terms to get PV:'),end
% If we had computed the first term:
if FLpv1
  if toshow,disp('First term alone'),end
  PV = PV1;
  PV_lon=PV1_lon;PV_lat=PV1_lat;PV_dpt=PV1_dpt;
end
% If we had computed the second term:
if FLpv2
  if exist('PV') % and the first one:
    if toshow,disp('Second term added to first one'),end
    PV = PV + PV2; 
  else           % or not:
    if toshow,disp('Second term alone'),end
    PV = PV2; 
    PV_lon=PV2_lon;PV_lat=PV2_lat;PV_dpt=PV2_dpt;  
  end
end
% If we had computed the third term:
if FLpv3
  if exist('PV') % and one of the first or second one:
    if toshow,disp('Third term added to first and/or second one(s)'),end
    PV = PV + PV3; 
  else           % or not:
    if toshow,disp('Third term alone'),end
    PV = PV3;
    PV_lon=PV3_lon;PV_lat=PV3_lat;PV_dpt=PV3_dpt;  
  end
end  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Record:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if toshow,disp('Now reccording PV file ...'),end

% General informations: 
if FLpv3 == 1
  netfil     = strcat('PV','.',netcdf_domain,'.',netcdf_suff);
  units      = 'kg/s/m^4';
  ncid       = 'PV';
  longname   = 'Potential vorticity';
  uniquename = 'potential_vorticity';
else
  netfil     = strcat('splPV','.',netcdf_domain,'.',netcdf_suff);
  units      = 'kg/s/m^4';
  ncid       = 'splPV';
  longname   = 'Simple Potential vorticity';
  uniquename = 'simple_potential_vorticity';
end %if  

% Open output file:
nc = netcdf(strcat(pathname,sla,netfil),'clobber');

% Define axis:
nc('X') = length(PV_lon);
nc('Y') = length(PV_lat);
nc('Z') = length(PV_dpt);

nc{'X'} = 'X';
nc{'Y'} = 'Y';
nc{'Z'} = 'Z';

nc{'X'}            = ncfloat('X');
nc{'X'}.uniquename = ncchar('X');
nc{'X'}.long_name  = ncchar('longitude');
nc{'X'}.gridtype   = nclong(0);
nc{'X'}.units      = ncchar('degrees_east');
nc{'X'}(:)         = PV_lon;

nc{'Y'}            = ncfloat('Y'); 
nc{'Y'}.uniquename = ncchar('Y');
nc{'Y'}.long_name  = ncchar('latitude');
nc{'Y'}.gridtype   = nclong(0);
nc{'Y'}.units      = ncchar('degrees_north');
nc{'Y'}(:)         = PV_lat;
 
nc{'Z'}            = ncfloat('Z');
nc{'Z'}.uniquename = ncchar('Z');
nc{'Z'}.long_name  = ncchar('depth');
nc{'Z'}.gridtype   = nclong(0);
nc{'Z'}.units      = ncchar('m');
nc{'Z'}(:)         = PV_dpt;

% And main field:
nc{ncid}               = ncfloat('Z', 'Y', 'X'); 
nc{ncid}.units         = ncchar(units);
nc{ncid}.missing_value = ncfloat(NaN);
nc{ncid}.FillValue_    = ncfloat(NaN);
nc{ncid}.longname      = ncchar(longname);
nc{ncid}.uniquename    = ncchar(uniquename);
nc{ncid}(:,:,:)        = PV;

nc=close(nc);
if FLpv3 ~= 2
   close(ncOx);
   close(ncOy);
   close(ncOz);
end
close(ncST);

% Outputs:
OUT = struct('PV',PV,'dpt',PV_dpt,'lat',PV_lat,'lon',PV_lon);
switch nargout
 case 1
  varargout(1) = {OUT};
 case 2
  varargout(1) = {struct('PV1',PV1,'dpt',PV1_dpt,'lat',PV1_lat,'lon',PV1_lon)};
  varargout(2) = {struct('PV2',PV2,'dpt',PV2_dpt,'lat',PV2_lat,'lon',PV2_lon)};
 case 3
  varargout(1) = {struct('PV1',PV1,'dpt',PV1_dpt,'lat',PV1_lat,'lon',PV1_lon)};
  varargout(2) = {struct('PV2',PV2,'dpt',PV2_dpt,'lat',PV2_lat,'lon',PV2_lon)};
  varargout(3) = {struct('PV3',PV3,'dpt',PV3_dpt,'lat',PV3_lat,'lon',PV3_lon)};
end
1	gmaze	1.1	%
2	gmaze	1.4	% [Q] = C_compute_potential_vorticity(SNAPSHOT,[WANTSPLPV])
3	gmaze	1.6	% [Q1,Q2,Q3] = C_compute_potential_vorticity(SNAPSHOT,[WANTSPLPV])
4	gmaze	1.1	%
5			% This file computes the potential vorticity Q from
6			% netcdf files of relative vorticity (OMEGAX, OMEGAY, ZETA)
7			% and potential density (SIGMATHETA) as
8			% Q = OMEGAX . dSIGMATHETA/dx + OMEGAY . dSIGMATHETA/dy + (f+ZETA).dSIGMATHETA/dz
9			%
10			% The optional flag WANTSPLPV is set to 0 by defaut. If turn to 1,
11			% then the program computes the simple PV defined by:
12			% splQ = f.dSIGMATHETA/dz
13			%
14			% Note that none of the fields are defined on the same grid points.
15			% So, I decided to compute Q on the same grid as SIGMATHETA, ie. the
16			% center of the c-grid.
17			%
18	gmaze	1.3	% Files names are:
19			% INPUT:
20			% ./netcdf-files/<SNAPSHOT>/OMEGAX.<netcdf_domain>.<netcdf_suff>
21			% ./netcdf-files/<SNAPSHOT>/OMEGAY.<netcdf_domain>.<netcdf_suff>
22			% ./netcdf-files/<SNAPSHOT>/ZETA.<netcdf_domain>.<netcdf_suff>
23			% ./netcdf-files/<SNAPSHOT>/SIGMATHETA.<netcdf_domain>.<netcdf_suff>
24			% OUPUT:
25			% ./netcdf-files/<SNAPSHOT>/PV.<netcdf_domain>.<netcdf_suff>
26			% or
27			% ./netcdf-files/<SNAPSHOT>/splPV.<netcdf_domain>.<netcdf_suff>
28			%
29	gmaze	1.1	% 06/07/2006
30			% gmaze@mit.edu
31			%
32
33	gmaze	1.5	function varargout = C_compute_potential_vorticity(snapshot,varargin)
34	gmaze	1.1
35	gmaze	1.2
36	gmaze	1.1	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37			%% Setup
38			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39			global sla netcdf_domain netcdf_suff
40			pv_checkpath
41
42			%% Flags to choose which term to compute (by default, all):
43			FLpv1 = 1;
44			FLpv2 = 1;
45			FLpv3 = 1;
46			if nargin==2 % case of optional flag presents:
47			if varargin{1}(1) == 1 % Case of the simple PV:
48			FLpv1 = 0;
49			FLpv2 = 0;
50			FLpv3 = 2;
51			end
52			end %if
53	gmaze	1.6	%[FLpv1 FLpv2 FLpv3]
54
55	gmaze	1.1
56			%% Optionnal flags:
57			global toshow % Turn to 1 to follow the computing process
58
59
60			%% NETCDF files:
61
62			% Path and extension to find them:
63			pathname = strcat('netcdf-files',sla,snapshot,sla);
64	gmaze	1.6	%pathname = '.';
65	gmaze	1.1	ext = strcat('.',netcdf_suff);
66
67			% Names:
68			if FLpv3 ~= 2 % We don't need them for splPV
69			filOx = strcat('OMEGAX' ,'.',netcdf_domain);
70			filOy = strcat('OMEGAY' ,'.',netcdf_domain);
71			filOz = strcat('ZETA' ,'.',netcdf_domain);
72			end %if
73			filST = strcat('SIGMATHETA','.',netcdf_domain);
74
75			% Load files and coordinates:
76			if FLpv3 ~= 2 % We don't need them for splPV
77			ferfile = strcat(pathname,sla,filOx,ext);
78			ncOx = netcdf(ferfile,'nowrite');
79			[Oxlon Oxlat Oxdpt] = coordfromnc(ncOx);
80			ferfile = strcat(pathname,sla,filOy,ext);
81			ncOy = netcdf(ferfile,'nowrite');
82			[Oylon Oylat Oydpt] = coordfromnc(ncOy);
83			ferfile = strcat(pathname,sla,filOz,ext);
84			ncOz = netcdf(ferfile,'nowrite');
85			[Ozlon Ozlat Ozdpt] = coordfromnc(ncOz);
86			end %if
87			ferfile = strcat(pathname,sla,filST,ext);
88			ncST = netcdf(ferfile,'nowrite');
89			[STlon STlat STdpt] = coordfromnc(ncST);
90
91
92			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
93			% Then, compute the first term: OMEGAX . dSIGMATHETA/dx %
94			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
95			if FLpv1
96
97			%%%%%
98			%% 1: Compute zonal gradient of SIGMATHETA:
99
100			% Dim:
101			if toshow,disp('dim'),end
102			nx = length(STlon) - 1;
103			ny = length(STlat);
104			nz = length(STdpt);
105
106			% Pre-allocate:
107			if toshow,disp('pre-allocate'),end
108			dSIGMATHETAdx = zeros(nz,ny,nx-1)*NaN;
109			dx = zeros(1,nx).*NaN;
110			STup = zeros(nz,nx);
111			STdw = zeros(nz,nx);
112
113			% Zonal gradient of SIGMATHETA:
114			if toshow,disp('grad'), end
115			for iy = 1 : ny
116			if toshow
117			disp(strcat('Computing dSIGMATHETA/dx at latitude : ',num2str(STlat(iy)),...
118			'^o (',num2str(iy),'/',num2str(ny),')' ));
119			end
120			[dx b] = meshgrid( m_lldist(STlon(1:nx+1),[1 1]*STlat(iy)), STdpt ) ; clear b
121			STup = squeeze(ncST{4}(:,iy,2:nx+1));
122			STdw = squeeze(ncST{4}(:,iy,1:nx));
123			dSTdx = ( STup - STdw ) ./ dx;
124			% Change horizontal grid point definition to fit with SIGMATHETA:
125			dSTdx = ( dSTdx(:,1:nx-1) + dSTdx(:,2:nx) )./2;
126			dSIGMATHETAdx(:,iy,:) = dSTdx;
127			end %for iy
128
129
130			%%%%%
131			%% 2: Move OMEGAX on the same grid:
132			if toshow,disp('Move OMEGAX on the same grid as dSIGMATHETA/dx'), end
133
134			% Change vertical gridding of OMEGAX:
135			Ox = ncOx{4}(:,:,:);
136			Ox = ( Ox(2:nz-1,:,:) + Ox(1:nz-2,:,:) )./2;
137			% And horizontal gridding:
138			Ox = ( Ox(:,2:ny-1,:) + Ox(:,1:ny-2,:) )./2;
139
140			%%%%%
141			%% 3: Make both fields having same limits:
142			%% (Keep points where both fields are defined)
143			Ox = squeeze(Ox(:,:,2:nx));
144			dSIGMATHETAdx = squeeze( dSIGMATHETAdx (2:nz-1,2:ny-1,:) );
145
146			%%%%%
147			%% 4: Last, compute first term of PV:
148			PV1 = Ox.*dSIGMATHETAdx ;
149
150			% and define axis fron the ST grid:
151			PV1_lon = STlon(2:length(STlon)-1);
152			PV1_lat = STlat(2:length(STlat)-1);
153			PV1_dpt = STdpt(2:length(STdpt)-1);
154
155			clear nx ny nz dx STup STdw iy dSTdx Ox dSIGMATHETAdx
156			end %if FLpv1
157
158
159
160
161			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
162			% Compute the second term: OMEGAY . dSIGMATHETA/dy %
163			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
164			if FLpv2
165
166			%%%%%
167			%% 1: Compute meridional gradient of SIGMATHETA:
168
169			% Dim:
170			if toshow,disp('dim'), end
171			nx = length(STlon) ;
172			ny = length(STlat) - 1 ;
173			nz = length(STdpt) ;
174
175			% Pre-allocate:
176			if toshow,disp('pre-allocate'), end
177			dSIGMATHETAdy = zeros(nz,ny-1,nx).*NaN;
178			dy = zeros(1,ny).*NaN;
179			STup = zeros(nz,ny);
180			STdw = zeros(nz,ny);
181
182			% Meridional gradient of SIGMATHETA:
183			% (Assuming the grid is regular, dy is independent of x)
184			[dy b] = meshgrid( m_lldist([1 1]*STlon(1),STlat(1:ny+1) ), STdpt ) ; clear b
185			for ix = 1 : nx
186			if toshow
187			disp(strcat('Computing dSIGMATHETA/dy at longitude : ',num2str(STlon(ix)),...
188			'^o (',num2str(ix),'/',num2str(nx),')' ));
189			end
190			STup = squeeze(ncST{4}(:,2:ny+1,ix));
191			STdw = squeeze(ncST{4}(:,1:ny,ix));
192			dSTdy = ( STup - STdw ) ./ dy;
193			% Change horizontal grid point definition to fit with SIGMATHETA:
194			dSTdy = ( dSTdy(:,1:ny-1) + dSTdy(:,2:ny) )./2;
195			dSIGMATHETAdy(:,:,ix) = dSTdy;
196			end %for iy
197
198			%%%%%
199			%% 2: Move OMEGAY on the same grid:
200			if toshow,disp('Move OMEGAY on the same grid as dSIGMATHETA/dy'), end
201
202			% Change vertical gridding of OMEGAY:
203			Oy = ncOy{4}(:,:,:);
204			Oy = ( Oy(2:nz-1,:,:) + Oy(1:nz-2,:,:) )./2;
205			% And horizontal gridding:
206			Oy = ( Oy(:,:,2:nx-1) + Oy(:,:,1:nx-2) )./2;
207
208			%%%%%
209			%% 3: Make them having same limits:
210			%% (Keep points where both fields are defined)
211			Oy = squeeze(Oy(:,2:ny,:));
212			dSIGMATHETAdy = squeeze( dSIGMATHETAdy (2:nz-1,:,2:nx-1) );
213
214			%%%%%
215			%% 4: Last, compute second term of PV:
216			PV2 = Oy.*dSIGMATHETAdy ;
217
218			% and defined axis fron the ST grid:
219			PV2_lon = STlon(2:length(STlon)-1);
220			PV2_lat = STlat(2:length(STlat)-1);
221			PV2_dpt = STdpt(2:length(STdpt)-1);
222
223
224			clear nx ny nz dy STup STdw dy dSTdy Oy dSIGMATHETAdy
225			end %if FLpv2
226
227
228
229
230
231			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
232			% Compute the third term: ( f + ZETA ) . dSIGMATHETA/dz %
233			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
234			if FLpv3
235
236			%%%%%
237			%% 1: Compute vertical gradient of SIGMATHETA:
238
239			% Dim:
240			if toshow,disp('dim'), end
241			nx = length(STlon) ;
242			ny = length(STlat) ;
243			nz = length(STdpt) - 1 ;
244
245			% Pre-allocate:
246			if toshow,disp('pre-allocate'), end
247			dSIGMATHETAdz = zeros(nz-1,ny,nx).*NaN;
248			ST = zeros(nz+1,ny,nx);
249			dz = zeros(1,nz).*NaN;
250
251			% Vertical grid differences:
252	gmaze	1.6	% STdpt contains negative values with STdpt(1) at the surface
253			% and STdpt(end) at the bottom of the ocean.
254			% So dz is positive with respect to z axis upward:
255			dz = -diff(STdpt);
256	gmaze	1.1	[a dz_3D c] = meshgrid(STlat,dz,STlon); clear a c
257
258			% Vertical gradient:
259			if toshow,disp('Vertical gradient of SIGMATHETA'), end
260			ST = ncST{4}(:,:,:);
261	gmaze	1.6	% Z axis upward, so vertical derivative is upper-part
262			% minus lower-part:
263			dSIGMATHETAdz = ( ST(1:nz,:,:) - ST(2:nz+1,:,:) ) ./ dz_3D;
264	gmaze	1.1	clear dz_3D ST
265
266			% Change vertical gridding:
267			dSIGMATHETAdz = ( dSIGMATHETAdz(1:nz-1,:,:) + dSIGMATHETAdz(2:nz,:,:) )./2;
268
269			if FLpv3 == 1 % Just for full PV
270
271			%%%%%
272			%% 2: Move ZETA on the same grid:
273			if toshow,disp('Move ZETA on the same grid as dSIGMATHETA/dz'), end
274			Oz = ncOz{4}(:,:,:);
275			% Change horizontal gridding:
276			Oz = ( Oz(:,:,2:nx-1) + Oz(:,:,1:nx-2) )./2;
277			Oz = ( Oz(:,2:ny-1,:) + Oz(:,1:ny-2,:) )./2;
278
279			end %if FLpv3=1
280
281			%%%%%
282			%% 3: Make them having same limits:
283			%% (Keep points where both fields are defined)
284			if FLpv3 == 1
285			Oz = squeeze(Oz(2:nz,:,:));
286			end %if
287			dSIGMATHETAdz = squeeze( dSIGMATHETAdz (:,2:ny-1,2:nx-1) );
288
289
290			%%%%%
291			%% 4: Last, compute third term of PV:
292			% and defined axis fron the ST grid:
293			PV3_lon = STlon(2:length(STlon)-1);
294			PV3_lat = STlat(2:length(STlat)-1);
295			PV3_dpt = STdpt(2:length(STdpt)-1);
296
297			% Planetary vorticity:
298			f = 2(2pi/86400)sin(PV3_latpi/180);
299			[a f c]=meshgrid(PV3_lon,f,PV3_dpt); clear a c
300			f = permute(f,[3 1 2]);
301
302			% Third term of PV:
303			if FLpv3 == 2
304			% Compute simple PV, just with planetary vorticity:
305			PV3 = f.*dSIGMATHETAdz ;
306			else
307			% To compute full PV:
308			PV3 = (f+Oz).*dSIGMATHETAdz ;
309			end
310
311
312
313			clear nx ny nz dz ST Oz dSIGMATHETAdz f
314			end %if FLpv3
315
316
317
318			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
319			% Then, compute potential vorticity:
320			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
321			if toshow,disp('Summing terms to get PV:'),end
322			% If we had computed the first term:
323			if FLpv1
324			if toshow,disp('First term alone'),end
325			PV = PV1;
326			PV_lon=PV1_lon;PV_lat=PV1_lat;PV_dpt=PV1_dpt;
327			end
328			% If we had computed the second term:
329			if FLpv2
330			if exist('PV') % and the first one:
331			if toshow,disp('Second term added to first one'),end
332			PV = PV + PV2;
333			else % or not:
334			if toshow,disp('Second term alone'),end
335			PV = PV2;
336			PV_lon=PV2_lon;PV_lat=PV2_lat;PV_dpt=PV2_dpt;
337			end
338			end
339			% If we had computed the third term:
340			if FLpv3
341			if exist('PV') % and one of the first or second one:
342			if toshow,disp('Third term added to first and/or second one(s)'),end
343			PV = PV + PV3;
344			else % or not:
345			if toshow,disp('Third term alone'),end
346			PV = PV3;
347			PV_lon=PV3_lon;PV_lat=PV3_lat;PV_dpt=PV3_dpt;
348			end
349			end
350
351
352			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
353			% Record:
354			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
355			if toshow,disp('Now reccording PV file ...'),end
356
357			% General informations:
358			if FLpv3 == 1
359			netfil = strcat('PV','.',netcdf_domain,'.',netcdf_suff);
360			units = 'kg/s/m^4';
361			ncid = 'PV';
362			longname = 'Potential vorticity';
363			uniquename = 'potential_vorticity';
364			else
365			netfil = strcat('splPV','.',netcdf_domain,'.',netcdf_suff);
366			units = 'kg/s/m^4';
367			ncid = 'splPV';
368			longname = 'Simple Potential vorticity';
369			uniquename = 'simple_potential_vorticity';
370			end %if
371
372			% Open output file:
373			nc = netcdf(strcat(pathname,sla,netfil),'clobber');
374
375			% Define axis:
376			nc('X') = length(PV_lon);
377			nc('Y') = length(PV_lat);
378			nc('Z') = length(PV_dpt);
379
380			nc{'X'} = 'X';
381			nc{'Y'} = 'Y';
382			nc{'Z'} = 'Z';
383
384			nc{'X'} = ncfloat('X');
385			nc{'X'}.uniquename = ncchar('X');
386			nc{'X'}.long_name = ncchar('longitude');
387			nc{'X'}.gridtype = nclong(0);
388			nc{'X'}.units = ncchar('degrees_east');
389			nc{'X'}(:) = PV_lon;
390
391			nc{'Y'} = ncfloat('Y');
392			nc{'Y'}.uniquename = ncchar('Y');
393			nc{'Y'}.long_name = ncchar('latitude');
394			nc{'Y'}.gridtype = nclong(0);
395			nc{'Y'}.units = ncchar('degrees_north');
396			nc{'Y'}(:) = PV_lat;
397
398			nc{'Z'} = ncfloat('Z');
399			nc{'Z'}.uniquename = ncchar('Z');
400			nc{'Z'}.long_name = ncchar('depth');
401			nc{'Z'}.gridtype = nclong(0);
402			nc{'Z'}.units = ncchar('m');
403			nc{'Z'}(:) = PV_dpt;
404
405			% And main field:
406			nc{ncid} = ncfloat('Z', 'Y', 'X');
407			nc{ncid}.units = ncchar(units);
408			nc{ncid}.missing_value = ncfloat(NaN);
409			nc{ncid}.FillValue_ = ncfloat(NaN);
410			nc{ncid}.longname = ncchar(longname);
411			nc{ncid}.uniquename = ncchar(uniquename);
412			nc{ncid}(:,:,:) = PV;
413
414			nc=close(nc);
415	gmaze	1.6	if FLpv3 ~= 2
416			close(ncOx);
417			close(ncOy);
418			close(ncOz);
419			end
420	gmaze	1.5	close(ncST);
421	gmaze	1.4
422			% Outputs:
423	gmaze	1.5	OUT = struct('PV',PV,'dpt',PV_dpt,'lat',PV_lat,'lon',PV_lon);
424	gmaze	1.4	switch nargout
425			case 1
426	gmaze	1.5	varargout(1) = {OUT};
427	gmaze	1.6	case 2
428			varargout(1) = {struct('PV1',PV1,'dpt',PV1_dpt,'lat',PV1_lat,'lon',PV1_lon)};
429			varargout(2) = {struct('PV2',PV2,'dpt',PV2_dpt,'lat',PV2_lat,'lon',PV2_lon)};
430			case 3
431			varargout(1) = {struct('PV1',PV1,'dpt',PV1_dpt,'lat',PV1_lat,'lon',PV1_lon)};
432			varargout(2) = {struct('PV2',PV2,'dpt',PV2_dpt,'lat',PV2_lat,'lon',PV2_lon)};
433			varargout(3) = {struct('PV3',PV3,'dpt',PV3_dpt,'lat',PV3_lat,'lon',PV3_lon)};
434	gmaze	1.4	end