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 |
% |
19 |
% Input files are supposed to be in a subdirectory called: |
20 |
% ./netcdf-files/<snapshot>/ |
21 |
% |
22 |
% File names id are stored in global variables: |
23 |
% netcdf_UVEL, netcdf_VVEL, netcdf_THETA, netcdf_SALTanom |
24 |
% with the format: |
25 |
% netcdf_<ID>.<netcdf_domain>.<netcdf_suff> |
26 |
% where netcdf_domain and netcdf_suff are also in global |
27 |
% THE DOT IS ADDED IN SUB-PROG, SO AVOID IT IN DEFINITIONS |
28 |
% |
29 |
% Note that Q is not initialy defined with the ratio by RHO (see below). |
30 |
% |
31 |
% A simple potential vorticity (splQ) computing is also available. |
32 |
% It is defined as: splQ = f. dSIGMATHETA/dz |
33 |
% |
34 |
% It's also possible to add a real last step 4 to compute PV as: |
35 |
% Q = -1/RHO * [Ox.dSTdx + Oy.dSTdy + (f+Oz).dSTdz] |
36 |
% Note that in this case, program loads the PV output from the |
37 |
% routine C_compute_potential_vorticity and simply multiply it |
38 |
% by: -1/RHO. |
39 |
% RHO may be computed with the routine compute_density.m |
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 |
|
55 |
% Files are looked for in subdirectory defined by: ./netcdf-files/<snapshot>/ |
56 |
% So let's define the snapshot: |
57 |
snapshot = ''; |
58 |
|
59 |
|
60 |
% File's name: |
61 |
global netcdf_UVEL netcdf_VVEL netcdf_THETA netcdf_SALTanom |
62 |
global netcdf_suff netcdf_domain |
63 |
netcdf_UVEL = 'UVEL'; |
64 |
netcdf_VVEL = 'VVEL'; |
65 |
netcdf_THETA = 'THETA'; |
66 |
netcdf_SALTanom = 'SALTanom'; |
67 |
netcdf_suff = 'nc'; |
68 |
netcdf_domain = 'north_atlantic'; % Must not be empty ! |
69 |
|
70 |
|
71 |
% FLAGS: |
72 |
% Turn 0/1 the following flag to determine which PV to compute: |
73 |
wantsplPV = 0; % (turn 1 for simple PV computing) |
74 |
% Turn 0/1 this flag to get online computing informations: |
75 |
global toshow |
76 |
toshow = 0; |
77 |
|
78 |
|
79 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMPUTING PV: |
80 |
if isempty('snapshot') |
81 |
disp(' THIS FILE IS AN EXAMPLE OF USE, COPY IT INTO YOUR') |
82 |
disp(' WORKING DIRECTORY WITH A DIFFERENT NAME AND') |
83 |
disp(' CUSTOMIZE IT !'); |
84 |
end |
85 |
|
86 |
% STEP 1: |
87 |
% Output netcdf file is: |
88 |
% ./netcdf-files/<snapshot>/SIGMATHETA.<netcdf_domain>.<netcdf_suff> |
89 |
A_compute_potential_density(snapshot) |
90 |
|
91 |
|
92 |
% STEP 2: |
93 |
% Output netcdf files are: |
94 |
% ./netcdf-files/<snapshot>/OMEGAX.<netcdf_domain>.<netcdf_suff> |
95 |
% ./netcdf-files/<snapshot>/OMEGAY.<netcdf_domain>.<netcdf_suff> |
96 |
% ./netcdf-files/<snapshot>/ZETA.<netcdf_domain>.<netcdf_suff> |
97 |
% No interest for the a splPV computing |
98 |
if ~wantsplPV |
99 |
B_compute_relative_vorticity(snapshot) |
100 |
end %if |
101 |
|
102 |
% STEP 3: |
103 |
% Output netcdf file is: |
104 |
% ./netcdf-files/<snapshot>/PV.<netcdf_domain>.<netcdf_suff> |
105 |
C_compute_potential_vorticity(snapshot,wantsplPV) |
106 |
|
107 |
% STEP 4: |
108 |
% Output netcdf file is (replace last one): |
109 |
% ./netcdf-files/<snapshot>/PV.<netcdf_domain>.<netcdf_suff> |
110 |
global netcdf_RHO netcdf_PV |
111 |
netcdf_RHO = 'RHO'; |
112 |
netcdf_PV = 'PV'; |
113 |
D_compute_potential_vorticity(snapshot) |
114 |
|
115 |
|
116 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THAT'S IT ! |