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 |
|
|
% 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 50*800*700 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 |
gmaze |
1.2 |
D_compute_potential_vorticity(snapshot,wantsplPV) |
115 |
gmaze |
1.1 |
|
116 |
|
|
|
117 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THAT'S IT ! |
118 |
gmaze |
1.2 |
|
119 |
|
|
% Keep clean workspace: |
120 |
|
|
clear wantsplPV toshow netcdf_* |
121 |
|
|
clear global wantsplPV toshow netcdf_* |