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_* |