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