/[MITgcm]/MITgcm_contrib/gmaze_pv/C_compute_potential_vorticity.m
ViewVC logotype

Contents of /MITgcm_contrib/gmaze_pv/C_compute_potential_vorticity.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.8 - (show annotations) (download)
Wed Sep 19 15:37:38 2007 UTC (16 years, 7 months ago) by gmaze
Branch: MAIN
CVS Tags: HEAD
Changes since 1.7: +0 -0 lines
General Update

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

  ViewVC Help
Powered by ViewVC 1.1.22