1 |
% |
2 |
% [] = C_COMPUTE_POTENTIAL_VORTICITY(SNAPSHOT,[WANTSPLPV]) |
3 |
% |
4 |
% This file computes the potential vorticity Q from |
5 |
% netcdf files of relative vorticity (OMEGAX, OMEGAY, ZETA) |
6 |
% and potential density (SIGMATHETA) as |
7 |
% Q = OMEGAX . dSIGMATHETA/dx + OMEGAY . dSIGMATHETA/dy + (f+ZETA).dSIGMATHETA/dz |
8 |
% |
9 |
% The optional flag WANTSPLPV is set to 0 by defaut. If turn to 1, |
10 |
% then the program computes the simple PV defined by: |
11 |
% splQ = f.dSIGMATHETA/dz |
12 |
% |
13 |
% Note that none of the fields are defined on the same grid points. |
14 |
% So, I decided to compute Q on the same grid as SIGMATHETA, ie. the |
15 |
% center of the c-grid. |
16 |
% |
17 |
% 06/07/2006 |
18 |
% gmaze@mit.edu |
19 |
% |
20 |
|
21 |
function [] = C_compute_potential_vorticity(snapshot,varargin) |
22 |
|
23 |
|
24 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
25 |
%% Setup |
26 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
27 |
global sla netcdf_domain netcdf_suff |
28 |
pv_checkpath |
29 |
|
30 |
%% Flags to choose which term to compute (by default, all): |
31 |
FLpv1 = 1; |
32 |
FLpv2 = 1; |
33 |
FLpv3 = 1; |
34 |
if nargin==2 % case of optional flag presents: |
35 |
if varargin{1}(1) == 1 % Case of the simple PV: |
36 |
FLpv1 = 0; |
37 |
FLpv2 = 0; |
38 |
FLpv3 = 2; |
39 |
end |
40 |
end %if |
41 |
|
42 |
%% Optionnal flags: |
43 |
global toshow % Turn to 1 to follow the computing process |
44 |
|
45 |
|
46 |
%% NETCDF files: |
47 |
|
48 |
% Path and extension to find them: |
49 |
pathname = strcat('netcdf-files',sla,snapshot,sla); |
50 |
ext = strcat('.',netcdf_suff); |
51 |
|
52 |
% Names: |
53 |
if FLpv3 ~= 2 % We don't need them for splPV |
54 |
filOx = strcat('OMEGAX' ,'.',netcdf_domain); |
55 |
filOy = strcat('OMEGAY' ,'.',netcdf_domain); |
56 |
filOz = strcat('ZETA' ,'.',netcdf_domain); |
57 |
end %if |
58 |
filST = strcat('SIGMATHETA','.',netcdf_domain); |
59 |
|
60 |
% Load files and coordinates: |
61 |
if FLpv3 ~= 2 % We don't need them for splPV |
62 |
ferfile = strcat(pathname,sla,filOx,ext); |
63 |
ncOx = netcdf(ferfile,'nowrite'); |
64 |
[Oxlon Oxlat Oxdpt] = coordfromnc(ncOx); |
65 |
ferfile = strcat(pathname,sla,filOy,ext); |
66 |
ncOy = netcdf(ferfile,'nowrite'); |
67 |
[Oylon Oylat Oydpt] = coordfromnc(ncOy); |
68 |
ferfile = strcat(pathname,sla,filOz,ext); |
69 |
ncOz = netcdf(ferfile,'nowrite'); |
70 |
[Ozlon Ozlat Ozdpt] = coordfromnc(ncOz); |
71 |
end %if |
72 |
ferfile = strcat(pathname,sla,filST,ext); |
73 |
ncST = netcdf(ferfile,'nowrite'); |
74 |
[STlon STlat STdpt] = coordfromnc(ncST); |
75 |
|
76 |
|
77 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
78 |
% Then, compute the first term: OMEGAX . dSIGMATHETA/dx % |
79 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
80 |
if FLpv1 |
81 |
|
82 |
%%%%% |
83 |
%% 1: Compute zonal gradient of SIGMATHETA: |
84 |
|
85 |
% Dim: |
86 |
if toshow,disp('dim'),end |
87 |
nx = length(STlon) - 1; |
88 |
ny = length(STlat); |
89 |
nz = length(STdpt); |
90 |
|
91 |
% Pre-allocate: |
92 |
if toshow,disp('pre-allocate'),end |
93 |
dSIGMATHETAdx = zeros(nz,ny,nx-1)*NaN; |
94 |
dx = zeros(1,nx).*NaN; |
95 |
STup = zeros(nz,nx); |
96 |
STdw = zeros(nz,nx); |
97 |
|
98 |
% Zonal gradient of SIGMATHETA: |
99 |
if toshow,disp('grad'), end |
100 |
for iy = 1 : ny |
101 |
if toshow |
102 |
disp(strcat('Computing dSIGMATHETA/dx at latitude : ',num2str(STlat(iy)),... |
103 |
'^o (',num2str(iy),'/',num2str(ny),')' )); |
104 |
end |
105 |
[dx b] = meshgrid( m_lldist(STlon(1:nx+1),[1 1]*STlat(iy)), STdpt ) ; clear b |
106 |
STup = squeeze(ncST{4}(:,iy,2:nx+1)); |
107 |
STdw = squeeze(ncST{4}(:,iy,1:nx)); |
108 |
dSTdx = ( STup - STdw ) ./ dx; |
109 |
% Change horizontal grid point definition to fit with SIGMATHETA: |
110 |
dSTdx = ( dSTdx(:,1:nx-1) + dSTdx(:,2:nx) )./2; |
111 |
dSIGMATHETAdx(:,iy,:) = dSTdx; |
112 |
end %for iy |
113 |
|
114 |
|
115 |
%%%%% |
116 |
%% 2: Move OMEGAX on the same grid: |
117 |
if toshow,disp('Move OMEGAX on the same grid as dSIGMATHETA/dx'), end |
118 |
|
119 |
% Change vertical gridding of OMEGAX: |
120 |
Ox = ncOx{4}(:,:,:); |
121 |
Ox = ( Ox(2:nz-1,:,:) + Ox(1:nz-2,:,:) )./2; |
122 |
% And horizontal gridding: |
123 |
Ox = ( Ox(:,2:ny-1,:) + Ox(:,1:ny-2,:) )./2; |
124 |
|
125 |
%%%%% |
126 |
%% 3: Make both fields having same limits: |
127 |
%% (Keep points where both fields are defined) |
128 |
Ox = squeeze(Ox(:,:,2:nx)); |
129 |
dSIGMATHETAdx = squeeze( dSIGMATHETAdx (2:nz-1,2:ny-1,:) ); |
130 |
|
131 |
%%%%% |
132 |
%% 4: Last, compute first term of PV: |
133 |
PV1 = Ox.*dSIGMATHETAdx ; |
134 |
|
135 |
% and define axis fron the ST grid: |
136 |
PV1_lon = STlon(2:length(STlon)-1); |
137 |
PV1_lat = STlat(2:length(STlat)-1); |
138 |
PV1_dpt = STdpt(2:length(STdpt)-1); |
139 |
|
140 |
clear nx ny nz dx STup STdw iy dSTdx Ox dSIGMATHETAdx |
141 |
end %if FLpv1 |
142 |
|
143 |
|
144 |
|
145 |
|
146 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
147 |
% Compute the second term: OMEGAY . dSIGMATHETA/dy % |
148 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
149 |
if FLpv2 |
150 |
|
151 |
%%%%% |
152 |
%% 1: Compute meridional gradient of SIGMATHETA: |
153 |
|
154 |
% Dim: |
155 |
if toshow,disp('dim'), end |
156 |
nx = length(STlon) ; |
157 |
ny = length(STlat) - 1 ; |
158 |
nz = length(STdpt) ; |
159 |
|
160 |
% Pre-allocate: |
161 |
if toshow,disp('pre-allocate'), end |
162 |
dSIGMATHETAdy = zeros(nz,ny-1,nx).*NaN; |
163 |
dy = zeros(1,ny).*NaN; |
164 |
STup = zeros(nz,ny); |
165 |
STdw = zeros(nz,ny); |
166 |
|
167 |
% Meridional gradient of SIGMATHETA: |
168 |
% (Assuming the grid is regular, dy is independent of x) |
169 |
[dy b] = meshgrid( m_lldist([1 1]*STlon(1),STlat(1:ny+1) ), STdpt ) ; clear b |
170 |
for ix = 1 : nx |
171 |
if toshow |
172 |
disp(strcat('Computing dSIGMATHETA/dy at longitude : ',num2str(STlon(ix)),... |
173 |
'^o (',num2str(ix),'/',num2str(nx),')' )); |
174 |
end |
175 |
STup = squeeze(ncST{4}(:,2:ny+1,ix)); |
176 |
STdw = squeeze(ncST{4}(:,1:ny,ix)); |
177 |
dSTdy = ( STup - STdw ) ./ dy; |
178 |
% Change horizontal grid point definition to fit with SIGMATHETA: |
179 |
dSTdy = ( dSTdy(:,1:ny-1) + dSTdy(:,2:ny) )./2; |
180 |
dSIGMATHETAdy(:,:,ix) = dSTdy; |
181 |
end %for iy |
182 |
|
183 |
%%%%% |
184 |
%% 2: Move OMEGAY on the same grid: |
185 |
if toshow,disp('Move OMEGAY on the same grid as dSIGMATHETA/dy'), end |
186 |
|
187 |
% Change vertical gridding of OMEGAY: |
188 |
Oy = ncOy{4}(:,:,:); |
189 |
Oy = ( Oy(2:nz-1,:,:) + Oy(1:nz-2,:,:) )./2; |
190 |
% And horizontal gridding: |
191 |
Oy = ( Oy(:,:,2:nx-1) + Oy(:,:,1:nx-2) )./2; |
192 |
|
193 |
%%%%% |
194 |
%% 3: Make them having same limits: |
195 |
%% (Keep points where both fields are defined) |
196 |
Oy = squeeze(Oy(:,2:ny,:)); |
197 |
dSIGMATHETAdy = squeeze( dSIGMATHETAdy (2:nz-1,:,2:nx-1) ); |
198 |
|
199 |
%%%%% |
200 |
%% 4: Last, compute second term of PV: |
201 |
PV2 = Oy.*dSIGMATHETAdy ; |
202 |
|
203 |
% and defined axis fron the ST grid: |
204 |
PV2_lon = STlon(2:length(STlon)-1); |
205 |
PV2_lat = STlat(2:length(STlat)-1); |
206 |
PV2_dpt = STdpt(2:length(STdpt)-1); |
207 |
|
208 |
|
209 |
clear nx ny nz dy STup STdw dy dSTdy Oy dSIGMATHETAdy |
210 |
end %if FLpv2 |
211 |
|
212 |
|
213 |
|
214 |
|
215 |
|
216 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
217 |
% Compute the third term: ( f + ZETA ) . dSIGMATHETA/dz % |
218 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
219 |
if FLpv3 |
220 |
|
221 |
%%%%% |
222 |
%% 1: Compute vertical gradient of SIGMATHETA: |
223 |
|
224 |
% Dim: |
225 |
if toshow,disp('dim'), end |
226 |
nx = length(STlon) ; |
227 |
ny = length(STlat) ; |
228 |
nz = length(STdpt) - 1 ; |
229 |
|
230 |
% Pre-allocate: |
231 |
if toshow,disp('pre-allocate'), end |
232 |
dSIGMATHETAdz = zeros(nz-1,ny,nx).*NaN; |
233 |
ST = zeros(nz+1,ny,nx); |
234 |
dz = zeros(1,nz).*NaN; |
235 |
|
236 |
% Vertical grid differences: |
237 |
dz = diff(STdpt); |
238 |
[a dz_3D c] = meshgrid(STlat,dz,STlon); clear a c |
239 |
|
240 |
% Vertical gradient: |
241 |
if toshow,disp('Vertical gradient of SIGMATHETA'), end |
242 |
ST = ncST{4}(:,:,:); |
243 |
dSIGMATHETAdz = ( ST(2:nz+1,:,:) - ST(1:nz,:,:) ) ./ dz_3D; |
244 |
clear dz_3D ST |
245 |
|
246 |
% Change vertical gridding: |
247 |
dSIGMATHETAdz = ( dSIGMATHETAdz(1:nz-1,:,:) + dSIGMATHETAdz(2:nz,:,:) )./2; |
248 |
|
249 |
if FLpv3 == 1 % Just for full PV |
250 |
|
251 |
%%%%% |
252 |
%% 2: Move ZETA on the same grid: |
253 |
if toshow,disp('Move ZETA on the same grid as dSIGMATHETA/dz'), end |
254 |
Oz = ncOz{4}(:,:,:); |
255 |
% Change horizontal gridding: |
256 |
Oz = ( Oz(:,:,2:nx-1) + Oz(:,:,1:nx-2) )./2; |
257 |
Oz = ( Oz(:,2:ny-1,:) + Oz(:,1:ny-2,:) )./2; |
258 |
|
259 |
end %if FLpv3=1 |
260 |
|
261 |
%%%%% |
262 |
%% 3: Make them having same limits: |
263 |
%% (Keep points where both fields are defined) |
264 |
if FLpv3 == 1 |
265 |
Oz = squeeze(Oz(2:nz,:,:)); |
266 |
end %if |
267 |
dSIGMATHETAdz = squeeze( dSIGMATHETAdz (:,2:ny-1,2:nx-1) ); |
268 |
|
269 |
|
270 |
%%%%% |
271 |
%% 4: Last, compute third term of PV: |
272 |
% and defined axis fron the ST grid: |
273 |
PV3_lon = STlon(2:length(STlon)-1); |
274 |
PV3_lat = STlat(2:length(STlat)-1); |
275 |
PV3_dpt = STdpt(2:length(STdpt)-1); |
276 |
|
277 |
% Planetary vorticity: |
278 |
f = 2*(2*pi/86400)*sin(PV3_lat*pi/180); |
279 |
[a f c]=meshgrid(PV3_lon,f,PV3_dpt); clear a c |
280 |
f = permute(f,[3 1 2]); |
281 |
|
282 |
% Third term of PV: |
283 |
if FLpv3 == 2 |
284 |
% Compute simple PV, just with planetary vorticity: |
285 |
PV3 = f.*dSIGMATHETAdz ; |
286 |
else |
287 |
% To compute full PV: |
288 |
PV3 = (f+Oz).*dSIGMATHETAdz ; |
289 |
end |
290 |
|
291 |
|
292 |
|
293 |
clear nx ny nz dz ST Oz dSIGMATHETAdz f |
294 |
end %if FLpv3 |
295 |
|
296 |
|
297 |
|
298 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
299 |
% Then, compute potential vorticity: |
300 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
301 |
if toshow,disp('Summing terms to get PV:'),end |
302 |
% If we had computed the first term: |
303 |
if FLpv1 |
304 |
if toshow,disp('First term alone'),end |
305 |
PV = PV1; |
306 |
PV_lon=PV1_lon;PV_lat=PV1_lat;PV_dpt=PV1_dpt; |
307 |
end |
308 |
% If we had computed the second term: |
309 |
if FLpv2 |
310 |
if exist('PV') % and the first one: |
311 |
if toshow,disp('Second term added to first one'),end |
312 |
PV = PV + PV2; |
313 |
else % or not: |
314 |
if toshow,disp('Second term alone'),end |
315 |
PV = PV2; |
316 |
PV_lon=PV2_lon;PV_lat=PV2_lat;PV_dpt=PV2_dpt; |
317 |
end |
318 |
end |
319 |
% If we had computed the third term: |
320 |
if FLpv3 |
321 |
if exist('PV') % and one of the first or second one: |
322 |
if toshow,disp('Third term added to first and/or second one(s)'),end |
323 |
PV = PV + PV3; |
324 |
else % or not: |
325 |
if toshow,disp('Third term alone'),end |
326 |
PV = PV3; |
327 |
PV_lon=PV3_lon;PV_lat=PV3_lat;PV_dpt=PV3_dpt; |
328 |
end |
329 |
end |
330 |
|
331 |
|
332 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
333 |
% Record: |
334 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
335 |
if toshow,disp('Now reccording PV file ...'),end |
336 |
|
337 |
% General informations: |
338 |
if FLpv3 == 1 |
339 |
netfil = strcat('PV','.',netcdf_domain,'.',netcdf_suff); |
340 |
units = 'kg/s/m^4'; |
341 |
ncid = 'PV'; |
342 |
longname = 'Potential vorticity'; |
343 |
uniquename = 'potential_vorticity'; |
344 |
else |
345 |
netfil = strcat('splPV','.',netcdf_domain,'.',netcdf_suff); |
346 |
units = 'kg/s/m^4'; |
347 |
ncid = 'splPV'; |
348 |
longname = 'Simple Potential vorticity'; |
349 |
uniquename = 'simple_potential_vorticity'; |
350 |
end %if |
351 |
|
352 |
% Open output file: |
353 |
nc = netcdf(strcat(pathname,sla,netfil),'clobber'); |
354 |
|
355 |
% Define axis: |
356 |
nc('X') = length(PV_lon); |
357 |
nc('Y') = length(PV_lat); |
358 |
nc('Z') = length(PV_dpt); |
359 |
|
360 |
nc{'X'} = 'X'; |
361 |
nc{'Y'} = 'Y'; |
362 |
nc{'Z'} = 'Z'; |
363 |
|
364 |
nc{'X'} = ncfloat('X'); |
365 |
nc{'X'}.uniquename = ncchar('X'); |
366 |
nc{'X'}.long_name = ncchar('longitude'); |
367 |
nc{'X'}.gridtype = nclong(0); |
368 |
nc{'X'}.units = ncchar('degrees_east'); |
369 |
nc{'X'}(:) = PV_lon; |
370 |
|
371 |
nc{'Y'} = ncfloat('Y'); |
372 |
nc{'Y'}.uniquename = ncchar('Y'); |
373 |
nc{'Y'}.long_name = ncchar('latitude'); |
374 |
nc{'Y'}.gridtype = nclong(0); |
375 |
nc{'Y'}.units = ncchar('degrees_north'); |
376 |
nc{'Y'}(:) = PV_lat; |
377 |
|
378 |
nc{'Z'} = ncfloat('Z'); |
379 |
nc{'Z'}.uniquename = ncchar('Z'); |
380 |
nc{'Z'}.long_name = ncchar('depth'); |
381 |
nc{'Z'}.gridtype = nclong(0); |
382 |
nc{'Z'}.units = ncchar('m'); |
383 |
nc{'Z'}(:) = PV_dpt; |
384 |
|
385 |
% And main field: |
386 |
nc{ncid} = ncfloat('Z', 'Y', 'X'); |
387 |
nc{ncid}.units = ncchar(units); |
388 |
nc{ncid}.missing_value = ncfloat(NaN); |
389 |
nc{ncid}.FillValue_ = ncfloat(NaN); |
390 |
nc{ncid}.longname = ncchar(longname); |
391 |
nc{ncid}.uniquename = ncchar(uniquename); |
392 |
nc{ncid}(:,:,:) = PV; |
393 |
|
394 |
nc=close(nc); |
395 |
|