MITgcm_contrib/gmaze_pv/D_compute_potential_vorticity.m

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

  
function 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);

%% Error check:
if PV_dpt~RHO_dpt | PV_lon~RHO_lon | PV_lat~RHO_lat
  error('D_compute_potential_vorticity.m: RHO and PV fields must be defined on the same grid');
  return;
end
  
%% Flags:
global toshow % Turn to 1 to follow the computing process


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

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

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


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

% General informations: 
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);

1	%
2	% [] = D_COMPUTE_POTENTIAL_VORTICITY(SNAPSHOT,[WANTSPLPV])
3	%
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	% if the PV computed is the simple one (f.dSIGMATHETA/dz).
9	%
10	% 06/21/2006
11	% gmaze@mit.edu
12	%
13
14
15	function D_compute_potential_vorticity(snapshot,varargin)
16
17
18	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19	%% Setup
20	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21	global sla netcdf_RHO netcdf_PV netcdf_domain netcdf_suff
22	pv_checkpath
23
24	%% Flags to choose which term to compute (by default, all):
25	FLpv3 = 1;
26	if nargin==2 % case of optional flag presents:
27	if varargin{1}(1) == 1 % Case of the simple PV:
28	FLpv3 = 0;
29	end
30	end %if
31
32	%% PV and RHO netcdf-files:
33	filPV = strcat(netcdf_PV ,'.',netcdf_domain);
34	filRHO = strcat(netcdf_RHO,'.',netcdf_domain);
35
36	%% Path and extension to find them:
37	pathname = strcat('netcdf-files',sla,snapshot);
38	ext = strcat('.',netcdf_suff);
39
40	%% Load netcdf files:
41	ferfile = strcat(pathname,sla,filPV,ext);
42	ncPV = netcdf(ferfile,'nowrite');
43	[PV_lon PV_lat PV_dpt] = coordfromnc(ncPV);
44
45	ferfile = strcat(pathname,sla,filRHO,ext);
46	ncRHO = netcdf(ferfile,'nowrite');
47	[RHO_lon RHO_lat RHO_dpt] = coordfromnc(ncRHO);
48
49	%% Error check:
50	if PV_dpt~RHO_dpt \| PV_lon~RHO_lon \| PV_lat~RHO_lat
51	error('D_compute_potential_vorticity.m: RHO and PV fields must be defined on the same grid');
52	return;
53	end
54
55	%% Flags:
56	global toshow % Turn to 1 to follow the computing process
57
58
59
60	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61	%% Apply the coefficient
62	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63
64	%% Pre-allocate:
65	if toshow,disp('Pre-allocate');end
66	PV = zeros(length(PV_dpt),length(PV_lat),length(PV_dpt)).*NaN;
67
68	%% Apply:
69	if toshow,disp('Multiplying PV field by -1/RHO'),end
70	PV = - ncPV{4}(:,:,:) ./ ncRHO{4}(:,:,:) ;
71
72
73
74	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
75	% Record:
76	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77	if toshow,disp('Now reccording PV file ...'),end
78
79	% General informations:
80	if FLpv3 == 1
81	netfil = strcat('PV','.',netcdf_domain,'.',netcdf_suff);
82	units = '1/s/m';
83	ncid = 'PV';
84	longname = 'Potential vorticity';
85	uniquename = 'potential_vorticity';
86	else
87	netfil = strcat('splPV','.',netcdf_domain,'.',netcdf_suff);
88	units = '1/s/m';
89	ncid = 'splPV';
90	longname = 'Simple Potential vorticity';
91	uniquename = 'simple_potential_vorticity';
92	end %if
93
94	% Open output file:
95	nc = netcdf(strcat(pathname,sla,netfil),'clobber');
96
97	% Define axis:
98	nc('X') = length(PV_lon);
99	nc('Y') = length(PV_lat);
100	nc('Z') = length(PV_dpt);
101
102	nc{'X'} = 'X';
103	nc{'Y'} = 'Y';
104	nc{'Z'} = 'Z';
105
106	nc{'X'} = ncfloat('X');
107	nc{'X'}.uniquename = ncchar('X');
108	nc{'X'}.long_name = ncchar('longitude');
109	nc{'X'}.gridtype = nclong(0);
110	nc{'X'}.units = ncchar('degrees_east');
111	nc{'X'}(:) = PV_lon;
112
113	nc{'Y'} = ncfloat('Y');
114	nc{'Y'}.uniquename = ncchar('Y');
115	nc{'Y'}.long_name = ncchar('latitude');
116	nc{'Y'}.gridtype = nclong(0);
117	nc{'Y'}.units = ncchar('degrees_north');
118	nc{'Y'}(:) = PV_lat;
119
120	nc{'Z'} = ncfloat('Z');
121	nc{'Z'}.uniquename = ncchar('Z');
122	nc{'Z'}.long_name = ncchar('depth');
123	nc{'Z'}.gridtype = nclong(0);
124	nc{'Z'}.units = ncchar('m');
125	nc{'Z'}(:) = PV_dpt;
126
127	% And main field:
128	nc{ncid} = ncfloat('Z', 'Y', 'X');
129	nc{ncid}.units = ncchar(units);
130	nc{ncid}.missing_value = ncfloat(NaN);
131	nc{ncid}.FillValue_ = ncfloat(NaN);
132	nc{ncid}.longname = ncchar(longname);
133	nc{ncid}.uniquename = ncchar(uniquename);
134	nc{ncid}(:,:,:) = PV;
135
136	nc=close(nc);
137