MITgcm_contrib/gmaze_pv/B_compute_relative_vorticity.m

%
% [] = B_COMPUTE_RELATIVE_VORTICITY(SNAPSHOT)
%
% For a time snapshot, this program computes the 
% 3D relative vorticity field from 3D 
% horizontal speed fields U,V (x,y,z) as:
% OMEGA = ( -dVdz ; dUdz ; dVdx - dUdy )
%       = (   Ox  ;  Oy  ;     ZETA    )
% 3 output files are created.
%
% 06/07/2006
% gmaze@mit.edu
%
  
function [] = B_compute_relative_vorticity(snapshot)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global sla netcdf_UVEL netcdf_VVEL netcdf_domain netcdf_suff
pv_checkpath


%% U,V files name:
filU = strcat(netcdf_UVEL,'.',netcdf_domain);
filV = strcat(netcdf_VVEL,'.',netcdf_domain);


%% Path and extension to find them:
pathname = strcat('netcdf-files',sla,snapshot,sla);
ext      = strcat('.',netcdf_suff);


%% Load files:
ferfile          = strcat(pathname,sla,filU,ext);
ncU              = netcdf(ferfile,'nowrite');
[Ulon Ulat Udpt] = coordfromnc(ncU);

ferfile          = strcat(pathname,sla,filV,ext);
ncV              = netcdf(ferfile,'nowrite');
[Vlon Vlat Vdpt] = coordfromnc(ncV);

clear ext ferfile

%% Optionnal flags
computeZETA = 1; % Compute ZETA or not ?
global toshow % Turn to 1 to follow the computing process


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VERTICAL COMPONENT: ZETA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% U field is on the zonal side of the c-grid and
% V field on the meridional one.
% So computing meridional gradient for U and 
% zonal gradient for V makes the relative vorticity
% zeta defined on the corner of the c-grid.

%%%%%%%%%%%%%%
%% Dimensions of ZETA field:
if toshow,disp('Dim'),end
  ny = length(Ulat)-1; 
  nx = length(Vlon)-1; 
  nz = length(Udpt); % Note that Udpt=Vdpt
  ZETA_lon = Ulon(1:nx);
  ZETA_lat = Vlat(1:ny);
  
%%%%%%%%%%%%%%
%% Pre-allocation:
if toshow,disp('Pre-allocate'),end
ZETA = zeros(nz,ny-1,nx-1).*NaN;
dx   = zeros(ny-1,nx-1);
dy   = zeros(ny-1,nx-1);

%%%%%%%%%%%%%%
%% Compute relative vorticity for each z-level:
if computeZETA
for iz=1:nz
  if toshow
    disp(strcat('Computing \zeta at depth : ',num2str(Udpt(iz)),...
                'm (',num2str(iz),'/',num2str(nz),')'   ));
  end
  
  % Get velocities:
  U = ncU{4}(iz,:,:);
  V = ncV{4}(iz,:,:);
  
  % And now compute the vertical component of relative vorticity:
  % (TO DO: m_lldist accepts tables as input, so this part may be
  % done without x,y loop ...)
  for iy = 1 : ny
    for ix = 1 : nx
      if iz==1 % It's more efficient to make this test each time than
              % recomputing distance each time. m_lldist is a slow routine.
         % ZETA axis and grid distance:
         dx(iy,ix) = m_lldist([Vlon(ix+1) Vlon(ix)],[1 1]*Vlat(iy));
         dy(iy,ix) = m_lldist([1 1]*Vlon(ix),[Ulat(iy+1) Ulat(iy)]);
      end %if 
      % Horizontal gradients and ZETA:
      dVdx        = ( V(iy,ix+1)-V(iy,ix) ) / dx(iy,ix) ;
      dUdy        = ( U(iy+1,ix)-U(iy,ix) ) / dy(iy,ix) ;
      ZETA(iz,iy,ix) = dVdx - dUdy;      
    end %for ix
  end %for iy

end %for iz

%%%%%%%%%%%%%%
%% Netcdf record:

% General informations: 
netfil     = strcat('ZETA','.',netcdf_domain,'.',netcdf_suff);
units      = '1/s';
ncid       = 'ZETA';
longname   = 'Vertical Component of the Relative Vorticity';
uniquename = 'vertical_relative_vorticity';

% Open output file:
nc = netcdf(strcat(pathname,sla,netfil),'clobber');

% Define axis:
nc('X') = nx;
nc('Y') = ny;
nc('Z') = nz;

nc{'X'} = 'X';
nc{'Y'} = 'Y';
nc{'Z'} = 'Z';

nc{'X'}            = ncfloat('X');
nc{'X'}.uniquename = ncchar('X');
nc{'X'}.long_name  = ncchar('longitude');
nc{'X'}.gridtype   = nclong(0);
nc{'X'}.units      = ncchar('degrees_east');
nc{'X'}(:)         = ZETA_lon;

nc{'Y'}            = ncfloat('Y'); 
nc{'Y'}.uniquename = ncchar('Y');
nc{'Y'}.long_name  = ncchar('latitude');
nc{'Y'}.gridtype   = nclong(0);
nc{'Y'}.units      = ncchar('degrees_north');
nc{'Y'}(:)         = ZETA_lat;
 
nc{'Z'}            = ncfloat('Z');
nc{'Z'}.uniquename = ncchar('Z');
nc{'Z'}.long_name  = ncchar('depth');
nc{'Z'}.gridtype   = nclong(0);
nc{'Z'}.units      = ncchar('m');
nc{'Z'}(:)         = Udpt;

% And main field:
nc{ncid}               = ncfloat('Z', 'Y', 'X'); 
nc{ncid}.units         = ncchar(units);
nc{ncid}.missing_value = ncfloat(NaN);
nc{ncid}.FillValue_    = ncfloat(NaN);
nc{ncid}.longname      = ncchar(longname);
nc{ncid}.uniquename    = ncchar(uniquename);
nc{ncid}(:,:,:)        = ZETA;

nc=close(nc);

clear x y z U V dx dy nx ny nz DVdx dUdy

end %if compute ZETA


%%%%%%%%%%%%%%%%%%%%%%%%%
% HORIZONTAL COMPONENTS %
%%%%%%%%%%%%%%%%%%%%%%%%%
if toshow, disp('')
           disp('Now compute horizontal components of relative vorticity ...'); end

% U and V are defined on the same Z grid.

%%%%%%%%%%%%%%
%% Dimensions of OMEGA x and y fields:
if toshow,disp('Dim'),end
  O_nx = [length(Vlon) length(Ulon)];
  O_ny = [length(Vlat) length(Ulat)];
  O_nz = length(Udpt) - 1; % Idem Vdpt
  
%%%%%%%%%%%%%%
%% Pre-allocations:
if toshow,disp('Pre-allocate'),end
Ox = zeros(O_nz,O_ny(1),O_nx(1)).*NaN;
Oy = zeros(O_nz,O_ny(2),O_nx(2)).*NaN;

%%%%%%%%%%%%%%
%% Horizontal components:

%% Vertical grid differences:
dZ   = diff(Udpt); 
Odpt = Udpt(1:O_nz) + dZ/2;

%% Zonal component of OMEGA:
if toshow,disp('Zonal direction ...'); end
[a dZ_3D c] = meshgrid(Vlat,dZ,Vlon); clear a c
          V = ncV{4}(:,:,:);
         Ox = - ( V(2:O_nz+1,:,:) - V(1:O_nz,:,:) ) ./ dZ_3D;
clear V dZ_3D % For memory use

%% Meridional component of OMEGA:
if toshow,disp('Meridional direction ...'); end
[a dZ_3D c] = meshgrid(Ulat,dZ,Ulon); clear a c
          U = ncU{4}(:,:,:);
         Oy = ( U(2:O_nz+1,:,:) - U(1:O_nz,:,:) ) ./ dZ_3D;
clear U dZ_3D % For memory use

clear dZ

%%%%%%%%%%%%%%
%% Record Zonal component:
if toshow,disp('Records ...'); end

% General informations: 
netfil     = strcat('OMEGAX','.',netcdf_domain,'.',netcdf_suff);
units      = '1/s';
ncid       = 'OMEGAX';
longname   = 'Zonal Component of the Relative Vorticity';
uniquename = 'zonal_relative_vorticity';

% Open output file:
nc = netcdf(strcat(pathname,sla,netfil),'clobber');

% Define axis:
nc('X') = O_nx(1);
nc('Y') = O_ny(1);
nc('Z') = O_nz;

nc{'X'} = 'X';
nc{'Y'} = 'Y';
nc{'Z'} = 'Z';

nc{'X'}            = ncfloat('X');
nc{'X'}.uniquename = ncchar('X');
nc{'X'}.long_name  = ncchar('longitude');
nc{'X'}.gridtype   = nclong(0);
nc{'X'}.units      = ncchar('degrees_east');
nc{'X'}(:)         = Vlon;

nc{'Y'}            = ncfloat('Y'); 
nc{'Y'}.uniquename = ncchar('Y');
nc{'Y'}.long_name  = ncchar('latitude');
nc{'Y'}.gridtype   = nclong(0);
nc{'Y'}.units      = ncchar('degrees_north');
nc{'Y'}(:)         = Vlat;
 
nc{'Z'}            = ncfloat('Z');
nc{'Z'}.uniquename = ncchar('Z');
nc{'Z'}.long_name  = ncchar('depth');
nc{'Z'}.gridtype   = nclong(0);
nc{'Z'}.units      = ncchar('m');
nc{'Z'}(:)         = Odpt;

% And main field:
nc{ncid}               = ncfloat('Z', 'Y', 'X'); 
nc{ncid}.units         = ncchar(units);
nc{ncid}.missing_value = ncfloat(NaN);
nc{ncid}.FillValue_    = ncfloat(NaN);
nc{ncid}.longname      = ncchar(longname);
nc{ncid}.uniquename    = ncchar(uniquename);
nc{ncid}(:,:,:)        = Ox;

nc=close(nc);

%%%%%%%%%%%%%%
%% Record Meridional component:
% General informations: 
netfil     = strcat('OMEGAY','.',netcdf_domain,'.',netcdf_suff);
units      = '1/s';
ncid       = 'OMEGAY';
longname   = 'Meridional Component of the Relative Vorticity';
uniquename = 'meridional_relative_vorticity';

% Open output file:
nc = netcdf(strcat(pathname,sla,netfil),'clobber');

% Define axis:
nc('X') = O_nx(2);
nc('Y') = O_ny(2);
nc('Z') = O_nz;

nc{'X'} = 'X';
nc{'Y'} = 'Y';
nc{'Z'} = 'Z';

nc{'X'}            = ncfloat('X');
nc{'X'}.uniquename = ncchar('X');
nc{'X'}.long_name  = ncchar('longitude');
nc{'X'}.gridtype   = nclong(0);
nc{'X'}.units      = ncchar('degrees_east');
nc{'X'}(:)         = Ulon;

nc{'Y'}            = ncfloat('Y'); 
nc{'Y'}.uniquename = ncchar('Y');
nc{'Y'}.long_name  = ncchar('latitude');
nc{'Y'}.gridtype   = nclong(0);
nc{'Y'}.units      = ncchar('degrees_north');
nc{'Y'}(:)         = Ulat;
 
nc{'Z'}            = ncfloat('Z');
nc{'Z'}.uniquename = ncchar('Z');
nc{'Z'}.long_name  = ncchar('depth');
nc{'Z'}.gridtype   = nclong(0);
nc{'Z'}.units      = ncchar('m');
nc{'Z'}(:)         = Odpt;

% And main field:
nc{ncid}               = ncfloat('Z', 'Y', 'X'); 
nc{ncid}.units         = ncchar(units);
nc{ncid}.missing_value = ncfloat(NaN);
nc{ncid}.FillValue_    = ncfloat(NaN);
nc{ncid}.longname      = ncchar(longname);
nc{ncid}.uniquename    = ncchar(uniquename);
nc{ncid}(:,:,:)        = Oy;

nc=close(nc);

1	gmaze	1.1	%
2			% [] = B_COMPUTE_RELATIVE_VORTICITY(SNAPSHOT)
3			%
4			% For a time snapshot, this program computes the
5			% 3D relative vorticity field from 3D
6			% horizontal speed fields U,V (x,y,z) as:
7			% OMEGA = ( -dVdz ; dUdz ; dVdx - dUdy )
8			% = ( Ox ; Oy ; ZETA )
9			% 3 output files are created.
10			%
11			% 06/07/2006
12			% gmaze@mit.edu
13			%
14
15			function [] = B_compute_relative_vorticity(snapshot)
16
17
18			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19			% Setup
20			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21			global sla netcdf_UVEL netcdf_VVEL netcdf_domain netcdf_suff
22			pv_checkpath
23
24
25			%% U,V files name:
26			filU = strcat(netcdf_UVEL,'.',netcdf_domain);
27			filV = strcat(netcdf_VVEL,'.',netcdf_domain);
28
29
30			%% Path and extension to find them:
31			pathname = strcat('netcdf-files',sla,snapshot,sla);
32			ext = strcat('.',netcdf_suff);
33
34
35			%% Load files:
36			ferfile = strcat(pathname,sla,filU,ext);
37			ncU = netcdf(ferfile,'nowrite');
38			[Ulon Ulat Udpt] = coordfromnc(ncU);
39
40			ferfile = strcat(pathname,sla,filV,ext);
41			ncV = netcdf(ferfile,'nowrite');
42			[Vlon Vlat Vdpt] = coordfromnc(ncV);
43
44			clear ext ferfile
45
46			%% Optionnal flags
47			computeZETA = 1; % Compute ZETA or not ?
48			global toshow % Turn to 1 to follow the computing process
49
50
51			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52			% VERTICAL COMPONENT: ZETA %
53			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
54
55			% U field is on the zonal side of the c-grid and
56			% V field on the meridional one.
57			% So computing meridional gradient for U and
58			% zonal gradient for V makes the relative vorticity
59			% zeta defined on the corner of the c-grid.
60
61			%%%%%%%%%%%%%%
62			%% Dimensions of ZETA field:
63			if toshow,disp('Dim'),end
64			ny = length(Ulat)-1;
65			nx = length(Vlon)-1;
66			nz = length(Udpt); % Note that Udpt=Vdpt
67			ZETA_lon = Ulon(1:nx);
68			ZETA_lat = Vlat(1:ny);
69
70			%%%%%%%%%%%%%%
71			%% Pre-allocation:
72			if toshow,disp('Pre-allocate'),end
73			ZETA = zeros(nz,ny-1,nx-1).*NaN;
74			dx = zeros(ny-1,nx-1);
75			dy = zeros(ny-1,nx-1);
76
77			%%%%%%%%%%%%%%
78			%% Compute relative vorticity for each z-level:
79			if computeZETA
80			for iz=1:nz
81			if toshow
82			disp(strcat('Computing \zeta at depth : ',num2str(Udpt(iz)),...
83			'm (',num2str(iz),'/',num2str(nz),')' ));
84			end
85
86			% Get velocities:
87			U = ncU{4}(iz,:,:);
88			V = ncV{4}(iz,:,:);
89
90			% And now compute the vertical component of relative vorticity:
91			% (TO DO: m_lldist accepts tables as input, so this part may be
92			% done without x,y loop ...)
93			for iy = 1 : ny
94			for ix = 1 : nx
95			if iz==1 % It's more efficient to make this test each time than
96			% recomputing distance each time. m_lldist is a slow routine.
97			% ZETA axis and grid distance:
98			dx(iy,ix) = m_lldist([Vlon(ix+1) Vlon(ix)],[1 1]*Vlat(iy));
99			dy(iy,ix) = m_lldist([1 1]*Vlon(ix),[Ulat(iy+1) Ulat(iy)]);
100			end %if
101			% Horizontal gradients and ZETA:
102			dVdx = ( V(iy,ix+1)-V(iy,ix) ) / dx(iy,ix) ;
103			dUdy = ( U(iy+1,ix)-U(iy,ix) ) / dy(iy,ix) ;
104			ZETA(iz,iy,ix) = dVdx - dUdy;
105			end %for ix
106			end %for iy
107
108			end %for iz
109
110			%%%%%%%%%%%%%%
111			%% Netcdf record:
112
113			% General informations:
114			netfil = strcat('ZETA','.',netcdf_domain,'.',netcdf_suff);
115			units = '1/s';
116			ncid = 'ZETA';
117			longname = 'Vertical Component of the Relative Vorticity';
118			uniquename = 'vertical_relative_vorticity';
119
120			% Open output file:
121			nc = netcdf(strcat(pathname,sla,netfil),'clobber');
122
123			% Define axis:
124			nc('X') = nx;
125			nc('Y') = ny;
126			nc('Z') = nz;
127
128			nc{'X'} = 'X';
129			nc{'Y'} = 'Y';
130			nc{'Z'} = 'Z';
131
132			nc{'X'} = ncfloat('X');
133			nc{'X'}.uniquename = ncchar('X');
134			nc{'X'}.long_name = ncchar('longitude');
135			nc{'X'}.gridtype = nclong(0);
136			nc{'X'}.units = ncchar('degrees_east');
137			nc{'X'}(:) = ZETA_lon;
138
139			nc{'Y'} = ncfloat('Y');
140			nc{'Y'}.uniquename = ncchar('Y');
141			nc{'Y'}.long_name = ncchar('latitude');
142			nc{'Y'}.gridtype = nclong(0);
143			nc{'Y'}.units = ncchar('degrees_north');
144			nc{'Y'}(:) = ZETA_lat;
145
146			nc{'Z'} = ncfloat('Z');
147			nc{'Z'}.uniquename = ncchar('Z');
148			nc{'Z'}.long_name = ncchar('depth');
149			nc{'Z'}.gridtype = nclong(0);
150			nc{'Z'}.units = ncchar('m');
151			nc{'Z'}(:) = Udpt;
152
153			% And main field:
154			nc{ncid} = ncfloat('Z', 'Y', 'X');
155			nc{ncid}.units = ncchar(units);
156			nc{ncid}.missing_value = ncfloat(NaN);
157			nc{ncid}.FillValue_ = ncfloat(NaN);
158			nc{ncid}.longname = ncchar(longname);
159			nc{ncid}.uniquename = ncchar(uniquename);
160			nc{ncid}(:,:,:) = ZETA;
161
162			nc=close(nc);
163
164			clear x y z U V dx dy nx ny nz DVdx dUdy
165
166			end %if compute ZETA
167
168
169			%%%%%%%%%%%%%%%%%%%%%%%%%
170			% HORIZONTAL COMPONENTS %
171			%%%%%%%%%%%%%%%%%%%%%%%%%
172			if toshow, disp('')
173			disp('Now compute horizontal components of relative vorticity ...'); end
174
175			% U and V are defined on the same Z grid.
176
177			%%%%%%%%%%%%%%
178			%% Dimensions of OMEGA x and y fields:
179			if toshow,disp('Dim'),end
180			O_nx = [length(Vlon) length(Ulon)];
181			O_ny = [length(Vlat) length(Ulat)];
182			O_nz = length(Udpt) - 1; % Idem Vdpt
183
184			%%%%%%%%%%%%%%
185			%% Pre-allocations:
186			if toshow,disp('Pre-allocate'),end
187			Ox = zeros(O_nz,O_ny(1),O_nx(1)).*NaN;
188			Oy = zeros(O_nz,O_ny(2),O_nx(2)).*NaN;
189
190			%%%%%%%%%%%%%%
191			%% Horizontal components:
192
193			%% Vertical grid differences:
194			dZ = diff(Udpt);
195			Odpt = Udpt(1:O_nz) + dZ/2;
196
197			%% Zonal component of OMEGA:
198			if toshow,disp('Zonal direction ...'); end
199			[a dZ_3D c] = meshgrid(Vlat,dZ,Vlon); clear a c
200			V = ncV{4}(:,:,:);
201			Ox = - ( V(2:O_nz+1,:,:) - V(1:O_nz,:,:) ) ./ dZ_3D;
202			clear V dZ_3D % For memory use
203
204			%% Meridional component of OMEGA:
205			if toshow,disp('Meridional direction ...'); end
206			[a dZ_3D c] = meshgrid(Ulat,dZ,Ulon); clear a c
207			U = ncU{4}(:,:,:);
208			Oy = ( U(2:O_nz+1,:,:) - U(1:O_nz,:,:) ) ./ dZ_3D;
209			clear U dZ_3D % For memory use
210
211			clear dZ
212
213			%%%%%%%%%%%%%%
214			%% Record Zonal component:
215			if toshow,disp('Records ...'); end
216
217			% General informations:
218			netfil = strcat('OMEGAX','.',netcdf_domain,'.',netcdf_suff);
219			units = '1/s';
220			ncid = 'OMEGAX';
221			longname = 'Zonal Component of the Relative Vorticity';
222			uniquename = 'zonal_relative_vorticity';
223
224			% Open output file:
225			nc = netcdf(strcat(pathname,sla,netfil),'clobber');
226
227			% Define axis:
228			nc('X') = O_nx(1);
229			nc('Y') = O_ny(1);
230			nc('Z') = O_nz;
231
232			nc{'X'} = 'X';
233			nc{'Y'} = 'Y';
234			nc{'Z'} = 'Z';
235
236			nc{'X'} = ncfloat('X');
237			nc{'X'}.uniquename = ncchar('X');
238			nc{'X'}.long_name = ncchar('longitude');
239			nc{'X'}.gridtype = nclong(0);
240			nc{'X'}.units = ncchar('degrees_east');
241			nc{'X'}(:) = Vlon;
242
243			nc{'Y'} = ncfloat('Y');
244			nc{'Y'}.uniquename = ncchar('Y');
245			nc{'Y'}.long_name = ncchar('latitude');
246			nc{'Y'}.gridtype = nclong(0);
247			nc{'Y'}.units = ncchar('degrees_north');
248			nc{'Y'}(:) = Vlat;
249
250			nc{'Z'} = ncfloat('Z');
251			nc{'Z'}.uniquename = ncchar('Z');
252			nc{'Z'}.long_name = ncchar('depth');
253			nc{'Z'}.gridtype = nclong(0);
254			nc{'Z'}.units = ncchar('m');
255			nc{'Z'}(:) = Odpt;
256
257			% And main field:
258			nc{ncid} = ncfloat('Z', 'Y', 'X');
259			nc{ncid}.units = ncchar(units);
260			nc{ncid}.missing_value = ncfloat(NaN);
261			nc{ncid}.FillValue_ = ncfloat(NaN);
262			nc{ncid}.longname = ncchar(longname);
263			nc{ncid}.uniquename = ncchar(uniquename);
264			nc{ncid}(:,:,:) = Ox;
265
266			nc=close(nc);
267
268			%%%%%%%%%%%%%%
269			%% Record Meridional component:
270			% General informations:
271			netfil = strcat('OMEGAY','.',netcdf_domain,'.',netcdf_suff);
272			units = '1/s';
273			ncid = 'OMEGAY';
274			longname = 'Meridional Component of the Relative Vorticity';
275			uniquename = 'meridional_relative_vorticity';
276
277			% Open output file:
278			nc = netcdf(strcat(pathname,sla,netfil),'clobber');
279
280			% Define axis:
281			nc('X') = O_nx(2);
282			nc('Y') = O_ny(2);
283			nc('Z') = O_nz;
284
285			nc{'X'} = 'X';
286			nc{'Y'} = 'Y';
287			nc{'Z'} = 'Z';
288
289			nc{'X'} = ncfloat('X');
290			nc{'X'}.uniquename = ncchar('X');
291			nc{'X'}.long_name = ncchar('longitude');
292			nc{'X'}.gridtype = nclong(0);
293			nc{'X'}.units = ncchar('degrees_east');
294			nc{'X'}(:) = Ulon;
295
296			nc{'Y'} = ncfloat('Y');
297			nc{'Y'}.uniquename = ncchar('Y');
298			nc{'Y'}.long_name = ncchar('latitude');
299			nc{'Y'}.gridtype = nclong(0);
300			nc{'Y'}.units = ncchar('degrees_north');
301			nc{'Y'}(:) = Ulat;
302
303			nc{'Z'} = ncfloat('Z');
304			nc{'Z'}.uniquename = ncchar('Z');
305			nc{'Z'}.long_name = ncchar('depth');
306			nc{'Z'}.gridtype = nclong(0);
307			nc{'Z'}.units = ncchar('m');
308			nc{'Z'}(:) = Odpt;
309
310			% And main field:
311			nc{ncid} = ncfloat('Z', 'Y', 'X');
312			nc{ncid}.units = ncchar(units);
313			nc{ncid}.missing_value = ncfloat(NaN);
314			nc{ncid}.FillValue_ = ncfloat(NaN);
315			nc{ncid}.longname = ncchar(longname);
316			nc{ncid}.uniquename = ncchar(uniquename);
317			nc{ncid}(:,:,:) = Oy;
318
319			nc=close(nc);
320