MITgcm_contrib/gmaze_pv/D_compute_potential_vorticity.m

%
% [Q] = D_compute_potential_vorticity(SNAPSHOT,[WANTSPLPV])
%
% For a time snapshot, this program multiplies the potential
% vorticity computed with C_COMPUTE_POTENTIAL_VORTICITY by the
% coefficient: -1/RHO
% Optional flag WANTSPLPV is turn to 0 by default. Turn it to 1
% if the PV computed is the simple one (f.dSIGMATHETA/dz). It's 
% needed for the output netcdf file informations.
% 
% CAUTION:
%% If all the PV computing procedure has been performed with routines
%% from the package, the PV field has less points than the RHO one, exactly
%% first and last in all directions have to be removed from RHO.
%
% Files names are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/<netcdf_RHO>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_PV>.<netcdf_domain>.<netcdf_suff>
% OUPUT:
% ./netcdf-files/<SNAPSHOT>/PV.<netcdf_domain>.<netcdf_suff>
% or 
% ./netcdf-files/<SNAPSHOT>/splPV.<netcdf_domain>.<netcdf_suff>
%
% 06/21/2006
% gmaze@mit.edu
%

  
function varargout = D_compute_potential_vorticity(snapshot,varargin)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global sla netcdf_RHO netcdf_PV netcdf_domain netcdf_suff
pv_checkpath

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

%% PV and RHO netcdf-files:
filPV  = strcat(netcdf_PV ,'.',netcdf_domain);
filRHO = strcat(netcdf_RHO,'.',netcdf_domain);

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

%% Load netcdf files:
ferfile = strcat(pathname,sla,filPV,ext);
ncPV    = netcdf(ferfile,'nowrite');
[PV_lon PV_lat PV_dpt] = coordfromnc(ncPV);

ferfile = strcat(pathname,sla,filRHO,ext);
ncRHO   = netcdf(ferfile,'nowrite');
[RHO_lon RHO_lat RHO_dpt] = coordfromnc(ncRHO);

%% Flags:
global toshow % Turn to 1 to follow the computing process


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Apply the coefficient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Pre-allocate:
if toshow,disp('Pre-allocate');end
nx = length(PV_lon);
ny = length(PV_lat);
nz = length(PV_dpt);
PV = zeros(nz,ny,nx).*NaN;

%% Apply:
if toshow,disp('Multiplying PV field by -1/RHO'),end
PV =  - ncPV{4}(:,:,:) ./ ncRHO{4}(2:nz+1,2:ny+1,2:nx+1) ;


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

% General informations: 
%ncclose(ncPV);

if FLpv3 == 1
  netfil     = strcat('PV','.',netcdf_domain,'.',netcdf_suff);
  units      = '1/s/m';
  ncid       = 'PV';
  longname   = 'Potential vorticity';
  uniquename = 'potential_vorticity';
else
  netfil     = strcat('splPV','.',netcdf_domain,'.',netcdf_suff);
  units      = '1/s/m';
  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);
close(ncPV);
close(ncRHO);

% Outputs:
OUT = struct('PV',PV,'dpt',PV_dpt,'lat',PV_lat,'lon',PV_lon);
switch nargout
 case 1
  varargout(1) = {OUT};
end
1	gmaze	1.1	%
2	gmaze	1.3	% [Q] = D_compute_potential_vorticity(SNAPSHOT,[WANTSPLPV])
3	gmaze	1.1	%
4			% For a time snapshot, this program multiplies the potential
5			% vorticity computed with C_COMPUTE_POTENTIAL_VORTICITY by the
6			% coefficient: -1/RHO
7			% Optional flag WANTSPLPV is turn to 0 by default. Turn it to 1
8	gmaze	1.2	% if the PV computed is the simple one (f.dSIGMATHETA/dz). It's
9			% needed for the output netcdf file informations.
10	gmaze	1.1	%
11	gmaze	1.2	% CAUTION:
12	gmaze	1.3	%% If all the PV computing procedure has been performed with routines
13	gmaze	1.2	%% from the package, the PV field has less points than the RHO one, exactly
14			%% first and last in all directions have to be removed from RHO.
15			%
16			% Files names are:
17			% INPUT:
18			% ./netcdf-files/<SNAPSHOT>/<netcdf_RHO>.<netcdf_domain>.<netcdf_suff>
19			% ./netcdf-files/<SNAPSHOT>/<netcdf_PV>.<netcdf_domain>.<netcdf_suff>
20			% OUPUT:
21			% ./netcdf-files/<SNAPSHOT>/PV.<netcdf_domain>.<netcdf_suff>
22			% or
23			% ./netcdf-files/<SNAPSHOT>/splPV.<netcdf_domain>.<netcdf_suff>
24			%
25	gmaze	1.1	% 06/21/2006
26			% gmaze@mit.edu
27			%
28
29
30	gmaze	1.4	function varargout = D_compute_potential_vorticity(snapshot,varargin)
31	gmaze	1.1
32
33			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34			%% Setup
35			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36			global sla netcdf_RHO netcdf_PV netcdf_domain netcdf_suff
37			pv_checkpath
38
39			%% Flags to choose which term to compute (by default, all):
40			FLpv3 = 1;
41			if nargin==2 % case of optional flag presents:
42			if varargin{1}(1) == 1 % Case of the simple PV:
43			FLpv3 = 0;
44			end
45			end %if
46
47			%% PV and RHO netcdf-files:
48			filPV = strcat(netcdf_PV ,'.',netcdf_domain);
49			filRHO = strcat(netcdf_RHO,'.',netcdf_domain);
50
51			%% Path and extension to find them:
52			pathname = strcat('netcdf-files',sla,snapshot);
53			ext = strcat('.',netcdf_suff);
54
55			%% Load netcdf files:
56			ferfile = strcat(pathname,sla,filPV,ext);
57			ncPV = netcdf(ferfile,'nowrite');
58			[PV_lon PV_lat PV_dpt] = coordfromnc(ncPV);
59
60			ferfile = strcat(pathname,sla,filRHO,ext);
61			ncRHO = netcdf(ferfile,'nowrite');
62			[RHO_lon RHO_lat RHO_dpt] = coordfromnc(ncRHO);
63
64			%% Flags:
65			global toshow % Turn to 1 to follow the computing process
66
67
68
69			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70			%% Apply the coefficient
71			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72
73			%% Pre-allocate:
74			if toshow,disp('Pre-allocate');end
75	gmaze	1.2	nx = length(PV_lon);
76			ny = length(PV_lat);
77			nz = length(PV_dpt);
78			PV = zeros(nz,ny,nx).*NaN;
79	gmaze	1.1
80			%% Apply:
81			if toshow,disp('Multiplying PV field by -1/RHO'),end
82	gmaze	1.2	PV = - ncPV{4}(:,:,:) ./ ncRHO{4}(2:nz+1,2:ny+1,2:nx+1) ;
83	gmaze	1.1
84
85			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86			% Record:
87			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88			if toshow,disp('Now reccording PV file ...'),end
89
90			% General informations:
91	gmaze	1.2	%ncclose(ncPV);
92
93	gmaze	1.1	if FLpv3 == 1
94			netfil = strcat('PV','.',netcdf_domain,'.',netcdf_suff);
95			units = '1/s/m';
96			ncid = 'PV';
97			longname = 'Potential vorticity';
98			uniquename = 'potential_vorticity';
99			else
100			netfil = strcat('splPV','.',netcdf_domain,'.',netcdf_suff);
101			units = '1/s/m';
102			ncid = 'splPV';
103			longname = 'Simple Potential vorticity';
104			uniquename = 'simple_potential_vorticity';
105			end %if
106
107			% Open output file:
108			nc = netcdf(strcat(pathname,sla,netfil),'clobber');
109
110			% Define axis:
111			nc('X') = length(PV_lon);
112			nc('Y') = length(PV_lat);
113			nc('Z') = length(PV_dpt);
114
115			nc{'X'} = 'X';
116			nc{'Y'} = 'Y';
117			nc{'Z'} = 'Z';
118
119			nc{'X'} = ncfloat('X');
120			nc{'X'}.uniquename = ncchar('X');
121			nc{'X'}.long_name = ncchar('longitude');
122			nc{'X'}.gridtype = nclong(0);
123			nc{'X'}.units = ncchar('degrees_east');
124			nc{'X'}(:) = PV_lon;
125
126			nc{'Y'} = ncfloat('Y');
127			nc{'Y'}.uniquename = ncchar('Y');
128			nc{'Y'}.long_name = ncchar('latitude');
129			nc{'Y'}.gridtype = nclong(0);
130			nc{'Y'}.units = ncchar('degrees_north');
131			nc{'Y'}(:) = PV_lat;
132
133			nc{'Z'} = ncfloat('Z');
134			nc{'Z'}.uniquename = ncchar('Z');
135			nc{'Z'}.long_name = ncchar('depth');
136			nc{'Z'}.gridtype = nclong(0);
137			nc{'Z'}.units = ncchar('m');
138			nc{'Z'}(:) = PV_dpt;
139
140			% And main field:
141			nc{ncid} = ncfloat('Z', 'Y', 'X');
142			nc{ncid}.units = ncchar(units);
143			nc{ncid}.missing_value = ncfloat(NaN);
144			nc{ncid}.FillValue_ = ncfloat(NaN);
145			nc{ncid}.longname = ncchar(longname);
146			nc{ncid}.uniquename = ncchar(uniquename);
147			nc{ncid}(:,:,:) = PV;
148
149			nc=close(nc);
150	gmaze	1.4	close(ncPV);
151			close(ncRHO);
152	gmaze	1.3
153			% Outputs:
154	gmaze	1.4	OUT = struct('PV',PV,'dpt',PV_dpt,'lat',PV_lat,'lon',PV_lon);
155	gmaze	1.3	switch nargout
156			case 1
157	gmaze	1.4	varargout(1) = {OUT};
158	gmaze	1.3	end