| 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_* |
Parent Directory
| 
Revision Log
| 
Revision Graph