MITgcm_contrib/gmaze_pv/C_compute_potential_vorticity.m

%
% [] = 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.
%
% 06/07/2006
% gmaze@mit.edu
%
  
function [] = 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

%% 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);
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:
         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}(:,:,:);
dSIGMATHETAdz = ( ST(2:nz+1,:,:) - ST(1:nz,:,:) ) ./ 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);

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