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