MITgcm_contrib/gmaze_pv/B_compute_relative_vorticity.m

%
% [OMEGA] = 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 outputs files are created.
%
% (U,V) must have same dimensions and by default are defined on
% a C-grid. 
% If (U,V) are defined on an A-grid (coming from a cube-sphere
% to lat/lon grid interpolation for example), ie at the same points
% as THETA, SALTanom, ... the global variable 'griddef' must
% be set to 'A-grid'. Then (U,V) are moved to a C-grid for the computation.
%
% ZETA is computed at the upper-right corner of the C-grid.
% OMEGAX and OMEGAY are computed at V and U locations but shifted downward
% by 1/2 grid. In case of a A-grid for (U,V), OMEGAX and OMEGAY are moved 
% to a C-grid according to the ZETA computation.
% 
%
% Files names are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/<netcdf_UVEL>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_VVEL>.<netcdf_domain>.<netcdf_suff>
% OUPUT:
% ./netcdf-files/<SNAPSHOT>/OMEGAX.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/OMEGAY.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/ZETA.<netcdf_domain>.<netcdf_suff>
%
% 2006/06/07
% gmaze@mit.edu
%
% Last update: 
% 2007/02/01 (gmaze) : Fix bug in ZETA grid and add compatibility with A-grid
%
  
% On the C-grid, U and V are supposed to have the same dimensions and are
% defined like this:
%
%  y
%  ^      -------------------------
%  |      |     |     |     |     |
%  | ny   U  *  U  *  U  *  U  *  |
%  |      |     |     |     |     |
%  |   ny -- V --- V --- V --- V --
%  |      |     |     |     |     |
%  |      U  *  U  *  U  *  U  *  |
%  |      |     |     |     |     |
%  |      -- V --- V --- V --- V --
%  |      |     |     |     |     |
%  |      U  *  U  *  U  *  U  *  |
%  |      |     |     |     |     |
%  |      -- V --- V --- V --- V --
%  |      |     |     |     |     |
%  |  1   U  *  U  *  U  *  U  *  |
%  |      |     |     |     |     |
%  |    1 -- V --- V --- V --- V --
%  |       
%  |      1                 nx
%  |         1                 nx
%--|-------------------------------------> x
%  | 
%
% On the A-grid, U and V are defined on *, so we simply shift U westward by 1/2 grid
% and V southward by 1/2 grid. New (U,V) have the same dimensions as original fields
% but with first col for U, and first row for V set to NaN. Values are computed by
% averaging two contiguous values.
%

function varargout = B_compute_relative_vorticity(snapshot)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global sla netcdf_UVEL netcdf_VVEL netcdf_domain netcdf_suff griddef
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 and axis:
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

%% Load grid definition:
global griddef
if length(griddef) == 0
  griddef = 'C-grid'; % By default
end
switch lower(griddef)
 case {'c-grid','cgrid','c'}
    % Nothing to do here
 case {'a-grid','agrid','a'}
    disp('Found (U,V) defined on A-grid')
    % Move Ulon westward by 1/2 grid point:
     Ulon = [Ulon(1)-abs(diff(Ulon(1:2))/2) ; (Ulon(1:end-1)+Ulon(2:end))/2];
    % Move V southward by 1/2 grid point:
     Vlat = [Vlat(1)-abs(diff(Vlat(1:2))/2); (Vlat(1:end-1)+Vlat(2:end))/2];
    % Now, (U,V) axis are defined as if they came from a C-grid
    % (U,V) fields are moved to a C-grid during computation...
 otherwise
    error('The grid must be: C-grid or A-grid');
    return
end %switch griddef
  
  
%% 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
  
%%%%%%%%%%%%%%
%% 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);

ZETA_lon = Ulon(2:nx+1);
ZETA_lat = Vlat(2:ny+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,:,:);
  switch lower(griddef)
   case {'a-grid','agrid','a'}
    % Move U westward by 1/2 grid point:
    % (1st col is set to nan, but axis defined)
    U = [ones(ny+1,1).*NaN  (U(:,1:end-1) + U(:,2:end))/2];
    % Move V southward by 1/2 grid point:
    % (1st row is set to nan but axis defined)
    V = [ones(1,nx+1).*NaN;  (V(1:end-1,:) + V(2:end,:))/2];
    % Now, U and V are defined as if they came from a C-grid
  end  
  
  % 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;

%%%%%%%%%%%%%%
%% Computation:

%% 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}(:,:,:);
switch lower(griddef)
   case {'a-grid','agrid','a'}
    % Move V southward by 1/2 grid point:
    % (1st row is set to nan but axis defined)
    V = cat(2,ones(O_nz+1,1,O_nx(1)).*NaN,(V(:,1:end-1,:) + V(:,2:end,:))/2);
    % Now, V is defined as if it came from a C-grid
end
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}(:,:,:);
switch lower(griddef)
   case {'a-grid','agrid','a'}
    % Move U westward by 1/2 grid point:
    % (1st col is set to nan, but axis defined)
    U = cat(3,ones(O_nz+1,O_ny(2),1).*NaN,(U(:,:,1:end-1) + U(:,:,2:end))/2);
    % Now, V is defined as if it came from a C-grid
end  
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);
close(ncU);
close(ncV);

% Outputs:
OMEGA = struct(...
    'Ox',struct('value',Ox,'dpt',Odpt,'lat',Vlat,'lon',Vlon),...
    'Oy',struct('value',Oy,'dpt',Odpt,'lat',Ulat,'lon',Vlon),...
    'Oz',struct('value',ZETA,'dpt',Udpt,'lat',ZETA_lat,'lon',ZETA_lon)...
    );
switch nargout
 case 1
  varargout(1) = {OMEGA};
end
1	gmaze	1.1	%
2	gmaze	1.3	% [OMEGA] = B_compute_relative_vorticity(SNAPSHOT)
3	gmaze	1.1	%
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	gmaze	1.4	% 3 outputs files are created.
10			%
11			% (U,V) must have same dimensions and by default are defined on
12			% a C-grid.
13			% If (U,V) are defined on an A-grid (coming from a cube-sphere
14			% to lat/lon grid interpolation for example), ie at the same points
15			% as THETA, SALTanom, ... the global variable 'griddef' must
16			% be set to 'A-grid'. Then (U,V) are moved to a C-grid for the computation.
17			%
18			% ZETA is computed at the upper-right corner of the C-grid.
19			% OMEGAX and OMEGAY are computed at V and U locations but shifted downward
20			% by 1/2 grid. In case of a A-grid for (U,V), OMEGAX and OMEGAY are moved
21			% to a C-grid according to the ZETA computation.
22			%
23	gmaze	1.1	%
24	gmaze	1.2	% Files names are:
25			% INPUT:
26			% ./netcdf-files/<SNAPSHOT>/<netcdf_UVEL>.<netcdf_domain>.<netcdf_suff>
27			% ./netcdf-files/<SNAPSHOT>/<netcdf_VVEL>.<netcdf_domain>.<netcdf_suff>
28			% OUPUT:
29			% ./netcdf-files/<SNAPSHOT>/OMEGAX.<netcdf_domain>.<netcdf_suff>
30			% ./netcdf-files/<SNAPSHOT>/OMEGAY.<netcdf_domain>.<netcdf_suff>
31			% ./netcdf-files/<SNAPSHOT>/ZETA.<netcdf_domain>.<netcdf_suff>
32			%
33	gmaze	1.4	% 2006/06/07
34	gmaze	1.1	% gmaze@mit.edu
35			%
36	gmaze	1.4	% Last update:
37			% 2007/02/01 (gmaze) : Fix bug in ZETA grid and add compatibility with A-grid
38			%
39	gmaze	1.1
40	gmaze	1.4	% On the C-grid, U and V are supposed to have the same dimensions and are
41			% defined like this:
42			%
43			% y
44			% ^ -------------------------
45			% \| \| \| \| \| \|
46			% \| ny U * U * U * U * \|
47			% \| \| \| \| \| \|
48			% \| ny -- V --- V --- V --- V --
49			% \| \| \| \| \| \|
50			% \| U * U * U * U * \|
51			% \| \| \| \| \| \|
52			% \| -- V --- V --- V --- V --
53			% \| \| \| \| \| \|
54			% \| U * U * U * U * \|
55			% \| \| \| \| \| \|
56			% \| -- V --- V --- V --- V --
57			% \| \| \| \| \| \|
58			% \| 1 U * U * U * U * \|
59			% \| \| \| \| \| \|
60			% \| 1 -- V --- V --- V --- V --
61			% \|
62			% \| 1 nx
63			% \| 1 nx
64			%--\|-------------------------------------> x
65			% \|
66			%
67			% On the A-grid, U and V are defined on *, so we simply shift U westward by 1/2 grid
68			% and V southward by 1/2 grid. New (U,V) have the same dimensions as original fields
69			% but with first col for U, and first row for V set to NaN. Values are computed by
70			% averaging two contiguous values.
71			%
72
73			function varargout = B_compute_relative_vorticity(snapshot)
74	gmaze	1.1
75
76			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77			% Setup
78			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
79	gmaze	1.4	global sla netcdf_UVEL netcdf_VVEL netcdf_domain netcdf_suff griddef
80	gmaze	1.1	pv_checkpath
81
82
83			%% U,V files name:
84			filU = strcat(netcdf_UVEL,'.',netcdf_domain);
85			filV = strcat(netcdf_VVEL,'.',netcdf_domain);
86
87
88			%% Path and extension to find them:
89			pathname = strcat('netcdf-files',sla,snapshot,sla);
90			ext = strcat('.',netcdf_suff);
91
92
93	gmaze	1.4	%% Load files and axis:
94	gmaze	1.1	ferfile = strcat(pathname,sla,filU,ext);
95			ncU = netcdf(ferfile,'nowrite');
96			[Ulon Ulat Udpt] = coordfromnc(ncU);
97
98			ferfile = strcat(pathname,sla,filV,ext);
99			ncV = netcdf(ferfile,'nowrite');
100			[Vlon Vlat Vdpt] = coordfromnc(ncV);
101
102			clear ext ferfile
103
104	gmaze	1.4	%% Load grid definition:
105			global griddef
106			if length(griddef) == 0
107			griddef = 'C-grid'; % By default
108			end
109			switch lower(griddef)
110			case {'c-grid','cgrid','c'}
111			% Nothing to do here
112			case {'a-grid','agrid','a'}
113			disp('Found (U,V) defined on A-grid')
114			% Move Ulon westward by 1/2 grid point:
115			Ulon = [Ulon(1)-abs(diff(Ulon(1:2))/2) ; (Ulon(1:end-1)+Ulon(2:end))/2];
116			% Move V southward by 1/2 grid point:
117			Vlat = [Vlat(1)-abs(diff(Vlat(1:2))/2); (Vlat(1:end-1)+Vlat(2:end))/2];
118			% Now, (U,V) axis are defined as if they came from a C-grid
119			% (U,V) fields are moved to a C-grid during computation...
120			otherwise
121			error('The grid must be: C-grid or A-grid');
122			return
123			end %switch griddef
124
125
126	gmaze	1.1	%% Optionnal flags
127			computeZETA = 1; % Compute ZETA or not ?
128			global toshow % Turn to 1 to follow the computing process
129
130
131			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
132			% VERTICAL COMPONENT: ZETA %
133			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
134
135			% U field is on the zonal side of the c-grid and
136			% V field on the meridional one.
137			% So computing meridional gradient for U and
138			% zonal gradient for V makes the relative vorticity
139			% zeta defined on the corner of the c-grid.
140
141			%%%%%%%%%%%%%%
142			%% Dimensions of ZETA field:
143			if toshow,disp('Dim'),end
144			ny = length(Ulat)-1;
145			nx = length(Vlon)-1;
146			nz = length(Udpt); % Note that Udpt=Vdpt
147
148			%%%%%%%%%%%%%%
149			%% Pre-allocation:
150			if toshow,disp('Pre-allocate'),end
151			ZETA = zeros(nz,ny-1,nx-1).*NaN;
152			dx = zeros(ny-1,nx-1);
153			dy = zeros(ny-1,nx-1);
154
155	gmaze	1.4	ZETA_lon = Ulon(2:nx+1);
156			ZETA_lat = Vlat(2:ny+1);
157
158	gmaze	1.1	%%%%%%%%%%%%%%
159			%% Compute relative vorticity for each z-level:
160			if computeZETA
161	gmaze	1.4	for iz = 1 : nz
162	gmaze	1.1	if toshow
163			disp(strcat('Computing \zeta at depth : ',num2str(Udpt(iz)),...
164			'm (',num2str(iz),'/',num2str(nz),')' ));
165			end
166
167			% Get velocities:
168			U = ncU{4}(iz,:,:);
169			V = ncV{4}(iz,:,:);
170	gmaze	1.4	switch lower(griddef)
171			case {'a-grid','agrid','a'}
172			% Move U westward by 1/2 grid point:
173			% (1st col is set to nan, but axis defined)
174			U = [ones(ny+1,1).*NaN (U(:,1:end-1) + U(:,2:end))/2];
175			% Move V southward by 1/2 grid point:
176			% (1st row is set to nan but axis defined)
177			V = [ones(1,nx+1).*NaN; (V(1:end-1,:) + V(2:end,:))/2];
178			% Now, U and V are defined as if they came from a C-grid
179			end
180	gmaze	1.1
181			% And now compute the vertical component of relative vorticity:
182			% (TO DO: m_lldist accepts tables as input, so this part may be
183			% done without x,y loop ...)
184			for iy = 1 : ny
185			for ix = 1 : nx
186			if iz==1 % It's more efficient to make this test each time than
187			% recomputing distance each time. m_lldist is a slow routine.
188			% ZETA axis and grid distance:
189			dx(iy,ix) = m_lldist([Vlon(ix+1) Vlon(ix)],[1 1]*Vlat(iy));
190			dy(iy,ix) = m_lldist([1 1]*Vlon(ix),[Ulat(iy+1) Ulat(iy)]);
191			end %if
192			% Horizontal gradients and ZETA:
193			dVdx = ( V(iy,ix+1)-V(iy,ix) ) / dx(iy,ix) ;
194			dUdy = ( U(iy+1,ix)-U(iy,ix) ) / dy(iy,ix) ;
195			ZETA(iz,iy,ix) = dVdx - dUdy;
196			end %for ix
197			end %for iy
198			end %for iz
199
200			%%%%%%%%%%%%%%
201			%% Netcdf record:
202
203			% General informations:
204			netfil = strcat('ZETA','.',netcdf_domain,'.',netcdf_suff);
205			units = '1/s';
206			ncid = 'ZETA';
207			longname = 'Vertical Component of the Relative Vorticity';
208			uniquename = 'vertical_relative_vorticity';
209
210			% Open output file:
211			nc = netcdf(strcat(pathname,sla,netfil),'clobber');
212
213			% Define axis:
214			nc('X') = nx;
215			nc('Y') = ny;
216			nc('Z') = nz;
217
218			nc{'X'} = 'X';
219			nc{'Y'} = 'Y';
220			nc{'Z'} = 'Z';
221
222			nc{'X'} = ncfloat('X');
223			nc{'X'}.uniquename = ncchar('X');
224			nc{'X'}.long_name = ncchar('longitude');
225			nc{'X'}.gridtype = nclong(0);
226			nc{'X'}.units = ncchar('degrees_east');
227			nc{'X'}(:) = ZETA_lon;
228
229			nc{'Y'} = ncfloat('Y');
230			nc{'Y'}.uniquename = ncchar('Y');
231			nc{'Y'}.long_name = ncchar('latitude');
232			nc{'Y'}.gridtype = nclong(0);
233			nc{'Y'}.units = ncchar('degrees_north');
234			nc{'Y'}(:) = ZETA_lat;
235
236			nc{'Z'} = ncfloat('Z');
237			nc{'Z'}.uniquename = ncchar('Z');
238			nc{'Z'}.long_name = ncchar('depth');
239			nc{'Z'}.gridtype = nclong(0);
240			nc{'Z'}.units = ncchar('m');
241			nc{'Z'}(:) = Udpt;
242
243			% And main field:
244			nc{ncid} = ncfloat('Z', 'Y', 'X');
245			nc{ncid}.units = ncchar(units);
246			nc{ncid}.missing_value = ncfloat(NaN);
247			nc{ncid}.FillValue_ = ncfloat(NaN);
248			nc{ncid}.longname = ncchar(longname);
249			nc{ncid}.uniquename = ncchar(uniquename);
250			nc{ncid}(:,:,:) = ZETA;
251
252			nc=close(nc);
253
254			clear x y z U V dx dy nx ny nz DVdx dUdy
255
256			end %if compute ZETA
257
258
259			%%%%%%%%%%%%%%%%%%%%%%%%%
260			% HORIZONTAL COMPONENTS %
261			%%%%%%%%%%%%%%%%%%%%%%%%%
262			if toshow, disp('')
263			disp('Now compute horizontal components of relative vorticity ...'); end
264
265			% U and V are defined on the same Z grid.
266
267			%%%%%%%%%%%%%%
268			%% Dimensions of OMEGA x and y fields:
269			if toshow,disp('Dim'),end
270			O_nx = [length(Vlon) length(Ulon)];
271			O_ny = [length(Vlat) length(Ulat)];
272			O_nz = length(Udpt) - 1; % Idem Vdpt
273
274			%%%%%%%%%%%%%%
275			%% Pre-allocations:
276			if toshow,disp('Pre-allocate'),end
277			Ox = zeros(O_nz,O_ny(1),O_nx(1)).*NaN;
278			Oy = zeros(O_nz,O_ny(2),O_nx(2)).*NaN;
279
280			%%%%%%%%%%%%%%
281	gmaze	1.4	%% Computation:
282	gmaze	1.1
283			%% Vertical grid differences:
284			dZ = diff(Udpt);
285			Odpt = Udpt(1:O_nz) + dZ/2;
286
287			%% Zonal component of OMEGA:
288			if toshow,disp('Zonal direction ...'); end
289			[a dZ_3D c] = meshgrid(Vlat,dZ,Vlon); clear a c
290	gmaze	1.4	V = ncV{4}(:,:,:);
291			switch lower(griddef)
292			case {'a-grid','agrid','a'}
293			% Move V southward by 1/2 grid point:
294			% (1st row is set to nan but axis defined)
295			V = cat(2,ones(O_nz+1,1,O_nx(1)).*NaN,(V(:,1:end-1,:) + V(:,2:end,:))/2);
296			% Now, V is defined as if it came from a C-grid
297			end
298			Ox = - ( V(2:O_nz+1,:,:) - V(1:O_nz,:,:) ) ./ dZ_3D;
299	gmaze	1.1	clear V dZ_3D % For memory use
300
301			%% Meridional component of OMEGA:
302			if toshow,disp('Meridional direction ...'); end
303			[a dZ_3D c] = meshgrid(Ulat,dZ,Ulon); clear a c
304	gmaze	1.4	U = ncU{4}(:,:,:);
305			switch lower(griddef)
306			case {'a-grid','agrid','a'}
307			% Move U westward by 1/2 grid point:
308			% (1st col is set to nan, but axis defined)
309			U = cat(3,ones(O_nz+1,O_ny(2),1).*NaN,(U(:,:,1:end-1) + U(:,:,2:end))/2);
310			% Now, V is defined as if it came from a C-grid
311			end
312			Oy = ( U(2:O_nz+1,:,:) - U(1:O_nz,:,:) ) ./ dZ_3D;
313	gmaze	1.1	clear U dZ_3D % For memory use
314
315			clear dZ
316
317	gmaze	1.4
318	gmaze	1.1	%%%%%%%%%%%%%%
319			%% Record Zonal component:
320			if toshow,disp('Records ...'); end
321
322			% General informations:
323			netfil = strcat('OMEGAX','.',netcdf_domain,'.',netcdf_suff);
324			units = '1/s';
325			ncid = 'OMEGAX';
326			longname = 'Zonal Component of the Relative Vorticity';
327			uniquename = 'zonal_relative_vorticity';
328
329			% Open output file:
330			nc = netcdf(strcat(pathname,sla,netfil),'clobber');
331
332			% Define axis:
333			nc('X') = O_nx(1);
334			nc('Y') = O_ny(1);
335			nc('Z') = O_nz;
336
337			nc{'X'} = 'X';
338			nc{'Y'} = 'Y';
339			nc{'Z'} = 'Z';
340
341			nc{'X'} = ncfloat('X');
342			nc{'X'}.uniquename = ncchar('X');
343			nc{'X'}.long_name = ncchar('longitude');
344			nc{'X'}.gridtype = nclong(0);
345			nc{'X'}.units = ncchar('degrees_east');
346			nc{'X'}(:) = Vlon;
347
348			nc{'Y'} = ncfloat('Y');
349			nc{'Y'}.uniquename = ncchar('Y');
350			nc{'Y'}.long_name = ncchar('latitude');
351			nc{'Y'}.gridtype = nclong(0);
352			nc{'Y'}.units = ncchar('degrees_north');
353			nc{'Y'}(:) = Vlat;
354
355			nc{'Z'} = ncfloat('Z');
356			nc{'Z'}.uniquename = ncchar('Z');
357			nc{'Z'}.long_name = ncchar('depth');
358			nc{'Z'}.gridtype = nclong(0);
359			nc{'Z'}.units = ncchar('m');
360			nc{'Z'}(:) = Odpt;
361
362			% And main field:
363			nc{ncid} = ncfloat('Z', 'Y', 'X');
364			nc{ncid}.units = ncchar(units);
365			nc{ncid}.missing_value = ncfloat(NaN);
366			nc{ncid}.FillValue_ = ncfloat(NaN);
367			nc{ncid}.longname = ncchar(longname);
368			nc{ncid}.uniquename = ncchar(uniquename);
369			nc{ncid}(:,:,:) = Ox;
370
371			nc=close(nc);
372
373			%%%%%%%%%%%%%%
374			%% Record Meridional component:
375			% General informations:
376			netfil = strcat('OMEGAY','.',netcdf_domain,'.',netcdf_suff);
377			units = '1/s';
378			ncid = 'OMEGAY';
379			longname = 'Meridional Component of the Relative Vorticity';
380			uniquename = 'meridional_relative_vorticity';
381
382			% Open output file:
383			nc = netcdf(strcat(pathname,sla,netfil),'clobber');
384
385			% Define axis:
386			nc('X') = O_nx(2);
387			nc('Y') = O_ny(2);
388			nc('Z') = O_nz;
389
390			nc{'X'} = 'X';
391			nc{'Y'} = 'Y';
392			nc{'Z'} = 'Z';
393
394			nc{'X'} = ncfloat('X');
395			nc{'X'}.uniquename = ncchar('X');
396			nc{'X'}.long_name = ncchar('longitude');
397			nc{'X'}.gridtype = nclong(0);
398			nc{'X'}.units = ncchar('degrees_east');
399			nc{'X'}(:) = Ulon;
400
401			nc{'Y'} = ncfloat('Y');
402			nc{'Y'}.uniquename = ncchar('Y');
403			nc{'Y'}.long_name = ncchar('latitude');
404			nc{'Y'}.gridtype = nclong(0);
405			nc{'Y'}.units = ncchar('degrees_north');
406			nc{'Y'}(:) = Ulat;
407
408			nc{'Z'} = ncfloat('Z');
409			nc{'Z'}.uniquename = ncchar('Z');
410			nc{'Z'}.long_name = ncchar('depth');
411			nc{'Z'}.gridtype = nclong(0);
412			nc{'Z'}.units = ncchar('m');
413			nc{'Z'}(:) = Odpt;
414
415			% And main field:
416			nc{ncid} = ncfloat('Z', 'Y', 'X');
417			nc{ncid}.units = ncchar(units);
418			nc{ncid}.missing_value = ncfloat(NaN);
419			nc{ncid}.FillValue_ = ncfloat(NaN);
420			nc{ncid}.longname = ncchar(longname);
421			nc{ncid}.uniquename = ncchar(uniquename);
422			nc{ncid}(:,:,:) = Oy;
423
424			nc=close(nc);
425	gmaze	1.4	close(ncU);
426			close(ncV);
427	gmaze	1.3
428			% Outputs:
429			OMEGA = struct(...
430			'Ox',struct('value',Ox,'dpt',Odpt,'lat',Vlat,'lon',Vlon),...
431			'Oy',struct('value',Oy,'dpt',Odpt,'lat',Ulat,'lon',Vlon),...
432			'Oz',struct('value',ZETA,'dpt',Udpt,'lat',ZETA_lat,'lon',ZETA_lon)...
433			);
434			switch nargout
435			case 1
436	gmaze	1.4	varargout(1) = {OMEGA};
437	gmaze	1.3	end