MITgcm_contrib/gmaze_pv/eg_main_getPV.m

% 
% THIS IS NOT A FUNCTION !
%
% Here is the main program to compute the potential vorticity Q
% from the flow (UVEL,VVEL), potential temperature (THETA) and
% salinity (SALTanom), given snapshot fields.
% 3 steps to do it:
%   1- compute the potential density SIGMATHETA (also called ST)
%      from THETA and SALTanom: 
%      ST = SIGMA(S,THETA,p=0)
%   2- compute the 3D relative vorticity field OMEGA (called O)
%      without vertical velocity terms:
%      O = ( -dVdz ; dUdz ; dVdx - dUdy )
%   3- compute the potential vorticity Q:
%      Q = Ox.dSTdx + Oy.dSTdy + (f+Oz).dSTdz
%      (note that we only add the planetary vorticity at this last
%      step).
%      It's also possible to add a real last step 4 to compute PV as:
%      Q = -1/RHO * [Ox.dSTdx + Oy.dSTdy + (f+Oz).dSTdz]
%      Note that in this case, program loads the PV output from the
%      routine C_compute_potential_vorticity (step 3) and simply multiply 
%      it by: -1/RHO.
%      RHO may be computed with the routine compute_density.m
%
%
% Input files are supposed to be in a subdirectory called: 
% ./netcdf-files/<snapshot>/
%
% File names id are stored in global variables:
%    netcdf_UVEL, netcdf_VVEL, netcdf_THETA, netcdf_SALTanom
% with the format:
%    netcdf_<ID>.<netcdf_domain>.<netcdf_suff>
% where netcdf_domain and netcdf_suff are also in global
% THE DOT IS ADDED IN SUB-PROG, SO AVOID IT IN DEFINITIONS
%
% Note that Q is not initialy defined with the ratio by -RHO.
%
% A simple potential vorticity (splQ) computing is also available.
% It is defined as: splQ = f. dSIGMATHETA/dz
% 
% 30Jan/2007
% gmaze@mit.edu
%
clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   SETUP:
pv_checkpath


% File's name:
global netcdf_UVEL netcdf_VVEL netcdf_THETA 
global netcdf_SALTanom is_SALTanom
global netcdf_TAUX netcdf_TAUY netcdf_SIGMATHETA 
global netcdf_RHO netcdf_EKL netcdf_Qnet netcdf_MLD
global netcdf_JFz netcdf_JBz
global netcdf_suff netcdf_domain sla
netcdf_UVEL     = 'UVEL';
netcdf_VVEL     = 'VVEL';
netcdf_THETA    = 'THETA';
netcdf_SALTanom = 'SALTanom'; is_SALTanom = 1;
netcdf_TAUX     = 'TAUX';
netcdf_TAUY     = 'TAUY';
netcdf_SIGMATHETA = 'SIGMATHETA';
netcdf_RHO      = 'RHO';
netcdf_EKL      = 'EKL';
netcdf_MLD      = 'KPPmld'; %netcdf_MLD      = 'MLD';
netcdf_Qnet     = 'TFLUX';
netcdf_JFz      = 'JFz'; 
netcdf_JBz      = 'JBz'; 
netcdf_suff     = 'nc';
netcdf_domain   = 'north_atlantic'; % Must not be empty !


% FLAGS:
% Turn 0/1 the following flag to determine which PV to compute:
wantsplPV = 0; % (turn 1 for simple PV computing)
% Turn 0/1 this flag to get online computing informations:
global toshow
toshow = 0;

% Get date list:
ll = dir(strcat('netcdf-files',sla));
nt = 0;
for il = 1 : size(ll,1)
  if ll(il).isdir & findstr(ll(il).name,'00')
    nt = nt + 1;
    list(nt).name = ll(il).name;
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TIME LOOP
for it = 1 : nt
  % Files are looked for in subdirectory defined by: ./netcdf-files/<snapshot>/
  snapshot = list(it).name;
  disp('********************************************************')
  disp('********************************************************')
  disp(snapshot)
  disp('********************************************************')
  disp('********************************************************')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   COMPUTING PV:
% STEP 1:
% Output netcdf file is:
%       ./netcdf-files/<snapshot>/SIGMATHETA.<netcdf_domain>.<netcdf_suff>
A_compute_potential_density(snapshot)
compute_density(snapshot)


% STEP 2:
% Output netcdf files are:
%       ./netcdf-files/<snapshot>/OMEGAX.<netcdf_domain>.<netcdf_suff>
%       ./netcdf-files/<snapshot>/OMEGAY.<netcdf_domain>.<netcdf_suff>
%       ./netcdf-files/<snapshot>/ZETA.<netcdf_domain>.<netcdf_suff>
% No interest for the a splPV computing
if ~wantsplPV
   B_compute_relative_vorticity(snapshot)
end %if

% STEP 3:
% Output netcdf file is:
%       ./netcdf-files/<snapshot>/PV.<netcdf_domain>.<netcdf_suff>
C_compute_potential_vorticity(snapshot,wantsplPV)

% STEP 4:
% Output netcdf file is (replace last one):
%       ./netcdf-files/<snapshot>/PV.<netcdf_domain>.<netcdf_suff>
global netcdf_PV
if wantsplPV == 1
  netcdf_PV = 'splPV';
else
  netcdf_PV = 'PV';
end %if
D_compute_potential_vorticity(snapshot,wantsplPV)


% OTHER computations:
if 0
 compute_alpha(snapshot)
 compute_MLD(snapshot)
 compute_EKL(snapshot)
 compute_JFz(snapshot);
 compute_JBz(snapshot);
 compute_Qek(snapshot);
end %if 1/0


fclose('all');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THAT'S IT !
end %for it


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THAT'S IT !

% Keep clean workspace:
clear wantsplPV toshow netcdf_*
clear global wantsplPV toshow netcdf_*
1	gmaze	1.1	%
2			% THIS IS NOT A FUNCTION !
3			%
4			% Here is the main program to compute the potential vorticity Q
5			% from the flow (UVEL,VVEL), potential temperature (THETA) and
6			% salinity (SALTanom), given snapshot fields.
7			% 3 steps to do it:
8			% 1- compute the potential density SIGMATHETA (also called ST)
9			% from THETA and SALTanom:
10			% ST = SIGMA(S,THETA,p=0)
11			% 2- compute the 3D relative vorticity field OMEGA (called O)
12			% without vertical velocity terms:
13			% O = ( -dVdz ; dUdz ; dVdx - dUdy )
14			% 3- compute the potential vorticity Q:
15			% Q = Ox.dSTdx + Oy.dSTdy + (f+Oz).dSTdz
16			% (note that we only add the planetary vorticity at this last
17			% step).
18			% It's also possible to add a real last step 4 to compute PV as:
19			% Q = -1/RHO * [Ox.dSTdx + Oy.dSTdy + (f+Oz).dSTdz]
20			% Note that in this case, program loads the PV output from the
21			% routine C_compute_potential_vorticity (step 3) and simply multiply
22			% it by: -1/RHO.
23			% RHO may be computed with the routine compute_density.m
24			%
25			%
26			% Input files are supposed to be in a subdirectory called:
27			% ./netcdf-files/<snapshot>/
28			%
29			% File names id are stored in global variables:
30			% netcdf_UVEL, netcdf_VVEL, netcdf_THETA, netcdf_SALTanom
31			% with the format:
32			% netcdf_<ID>.<netcdf_domain>.<netcdf_suff>
33			% where netcdf_domain and netcdf_suff are also in global
34			% THE DOT IS ADDED IN SUB-PROG, SO AVOID IT IN DEFINITIONS
35			%
36			% Note that Q is not initialy defined with the ratio by -RHO.
37			%
38			% A simple potential vorticity (splQ) computing is also available.
39			% It is defined as: splQ = f. dSIGMATHETA/dz
40			%
41			% 30Jan/2007
42			% gmaze@mit.edu
43			%
44			clear
45
46			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SETUP:
47			pv_checkpath
48
49
50			% File's name:
51			global netcdf_UVEL netcdf_VVEL netcdf_THETA
52			global netcdf_SALTanom is_SALTanom
53			global netcdf_TAUX netcdf_TAUY netcdf_SIGMATHETA
54			global netcdf_RHO netcdf_EKL netcdf_Qnet netcdf_MLD
55			global netcdf_JFz netcdf_JBz
56			global netcdf_suff netcdf_domain sla
57			netcdf_UVEL = 'UVEL';
58			netcdf_VVEL = 'VVEL';
59			netcdf_THETA = 'THETA';
60			netcdf_SALTanom = 'SALTanom'; is_SALTanom = 1;
61			netcdf_TAUX = 'TAUX';
62			netcdf_TAUY = 'TAUY';
63			netcdf_SIGMATHETA = 'SIGMATHETA';
64			netcdf_RHO = 'RHO';
65			netcdf_EKL = 'EKL';
66			netcdf_MLD = 'KPPmld'; %netcdf_MLD = 'MLD';
67			netcdf_Qnet = 'TFLUX';
68			netcdf_JFz = 'JFz';
69			netcdf_JBz = 'JBz';
70			netcdf_suff = 'nc';
71			netcdf_domain = 'north_atlantic'; % Must not be empty !
72
73
74
75			% FLAGS:
76			% Turn 0/1 the following flag to determine which PV to compute:
77			wantsplPV = 0; % (turn 1 for simple PV computing)
78			% Turn 0/1 this flag to get online computing informations:
79			global toshow
80			toshow = 0;
81
82			% Get date list:
83			ll = dir(strcat('netcdf-files',sla));
84			nt = 0;
85			for il = 1 : size(ll,1)
86			if ll(il).isdir & findstr(ll(il).name,'00')
87			nt = nt + 1;
88			list(nt).name = ll(il).name;
89			end
90			end
91
92
93			%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TIME LOOP
94			for it = 1 : nt
95			% Files are looked for in subdirectory defined by: ./netcdf-files/<snapshot>/
96			snapshot = list(it).name;
97			disp('********************************************************')
98			disp('********************************************************')
99			disp(snapshot)
100			disp('********************************************************')
101			disp('********************************************************')
102
103
104			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMPUTING PV:
105			% STEP 1:
106			% Output netcdf file is:
107			% ./netcdf-files/<snapshot>/SIGMATHETA.<netcdf_domain>.<netcdf_suff>
108			A_compute_potential_density(snapshot)
109			compute_density(snapshot)
110
111
112			% STEP 2:
113			% Output netcdf files are:
114			% ./netcdf-files/<snapshot>/OMEGAX.<netcdf_domain>.<netcdf_suff>
115			% ./netcdf-files/<snapshot>/OMEGAY.<netcdf_domain>.<netcdf_suff>
116			% ./netcdf-files/<snapshot>/ZETA.<netcdf_domain>.<netcdf_suff>
117			% No interest for the a splPV computing
118			if ~wantsplPV
119			B_compute_relative_vorticity(snapshot)
120			end %if
121
122			% STEP 3:
123			% Output netcdf file is:
124			% ./netcdf-files/<snapshot>/PV.<netcdf_domain>.<netcdf_suff>
125			C_compute_potential_vorticity(snapshot,wantsplPV)
126
127			% STEP 4:
128			% Output netcdf file is (replace last one):
129			% ./netcdf-files/<snapshot>/PV.<netcdf_domain>.<netcdf_suff>
130			global netcdf_PV
131			if wantsplPV == 1
132			netcdf_PV = 'splPV';
133			else
134			netcdf_PV = 'PV';
135			end %if
136			D_compute_potential_vorticity(snapshot,wantsplPV)
137
138
139			% OTHER computations:
140			if 0
141			compute_alpha(snapshot)
142			compute_MLD(snapshot)
143			compute_EKL(snapshot)
144			compute_JFz(snapshot);
145			compute_JBz(snapshot);
146			compute_Qek(snapshot);
147			end %if 1/0
148
149
150			fclose('all');
151			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THAT'S IT !
152			end %for it
153
154
155			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THAT'S IT !
156
157			% Keep clean workspace:
158			clear wantsplPV toshow netcdf_*
159			clear global wantsplPV toshow netcdf_*