MITgcm_contrib/gmaze_pv/main_getPV_eg.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
% 
% 06/21/2006
% gmaze@mit.edu
%
clear

disp('')
disp('This program will compute the potential vorticity from horizontal flow,')
disp('potential temperature and salinity fields')
disp('(For a 50*800*700 grid this takes around 30 minutes to compute PV ...)')
disp('')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   SETUP:
pv-checkpath

% Files are looked for in subdirectory defined by: ./netcdf-files/<snapshot>/
% So let's define the snapshot:
snapshot = '';


% File's name:
global netcdf_UVEL netcdf_VVEL netcdf_THETA netcdf_SALTanom
global netcdf_suff netcdf_domain
netcdf_UVEL     = 'UVEL';
netcdf_VVEL     = 'VVEL';
netcdf_THETA    = 'THETA';
netcdf_SALTanom = 'SALTanom';
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;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   COMPUTING PV:
if isempty('snapshot')
  disp(' THIS FILE IS AN EXAMPLE OF USE, COPY IT INTO YOUR')
  disp(' WORKING DIRECTORY WITH A DIFFERENT NAME AND')
  disp(' CUSTOMIZE IT !');
end  

% STEP 1:
% Output netcdf file is:
%       ./netcdf-files/<snapshot>/SIGMATHETA.<netcdf_domain>.<netcdf_suff>
A_compute_potential_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_RHO netcdf_PV
netcdf_RHO = 'RHO';
netcdf_PV  = 'PV';
D_compute_potential_vorticity(snapshot,wantsplPV)


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

% Keep clean workspace:
clear wantsplPV toshow netcdf_*
clear global wantsplPV toshow netcdf_*
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	% 06/21/2006
42	% gmaze@mit.edu
43	%
44	clear
45
46	disp('')
47	disp('This program will compute the potential vorticity from horizontal flow,')
48	disp('potential temperature and salinity fields')
49	disp('(For a 50800700 grid this takes around 30 minutes to compute PV ...)')
50	disp('')
51
52
53	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SETUP:
54	pv-checkpath
55
56	% Files are looked for in subdirectory defined by: ./netcdf-files/<snapshot>/
57	% So let's define the snapshot:
58	snapshot = '';
59
60
61	% File's name:
62	global netcdf_UVEL netcdf_VVEL netcdf_THETA netcdf_SALTanom
63	global netcdf_suff netcdf_domain
64	netcdf_UVEL = 'UVEL';
65	netcdf_VVEL = 'VVEL';
66	netcdf_THETA = 'THETA';
67	netcdf_SALTanom = 'SALTanom';
68	netcdf_suff = 'nc';
69	netcdf_domain = 'north_atlantic'; % Must not be empty !
70
71
72	% FLAGS:
73	% Turn 0/1 the following flag to determine which PV to compute:
74	wantsplPV = 0; % (turn 1 for simple PV computing)
75	% Turn 0/1 this flag to get online computing informations:
76	global toshow
77	toshow = 0;
78
79
80	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMPUTING PV:
81	if isempty('snapshot')
82	disp(' THIS FILE IS AN EXAMPLE OF USE, COPY IT INTO YOUR')
83	disp(' WORKING DIRECTORY WITH A DIFFERENT NAME AND')
84	disp(' CUSTOMIZE IT !');
85	end
86
87	% STEP 1:
88	% Output netcdf file is:
89	% ./netcdf-files/<snapshot>/SIGMATHETA.<netcdf_domain>.<netcdf_suff>
90	A_compute_potential_density(snapshot)
91
92
93	% STEP 2:
94	% Output netcdf files are:
95	% ./netcdf-files/<snapshot>/OMEGAX.<netcdf_domain>.<netcdf_suff>
96	% ./netcdf-files/<snapshot>/OMEGAY.<netcdf_domain>.<netcdf_suff>
97	% ./netcdf-files/<snapshot>/ZETA.<netcdf_domain>.<netcdf_suff>
98	% No interest for the a splPV computing
99	if ~wantsplPV
100	B_compute_relative_vorticity(snapshot)
101	end %if
102
103	% STEP 3:
104	% Output netcdf file is:
105	% ./netcdf-files/<snapshot>/PV.<netcdf_domain>.<netcdf_suff>
106	C_compute_potential_vorticity(snapshot,wantsplPV)
107
108	% STEP 4:
109	% Output netcdf file is (replace last one):
110	% ./netcdf-files/<snapshot>/PV.<netcdf_domain>.<netcdf_suff>
111	global netcdf_RHO netcdf_PV
112	netcdf_RHO = 'RHO';
113	netcdf_PV = 'PV';
114	D_compute_potential_vorticity(snapshot,wantsplPV)
115
116
117	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THAT'S IT !
118
119	% Keep clean workspace:
120	clear wantsplPV toshow netcdf_*
121	clear global wantsplPV toshow netcdf_*