/[MITgcm]/MITgcm_contrib/gmaze_pv/B_compute_relative_vorticity.m

Diff of /MITgcm_contrib/gmaze_pv/B_compute_relative_vorticity.m

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

-revision 1.3 by gmaze,
Tue Jan 30 22:10:10 2007 UTC
+revision 1.4 by gmaze,
Thu Feb  1 17:02:02 2007 UTC
 Line 6
  % horizontal speed fields U,V (x,y,z) as:
  % OMEGA = ( -dVdz ; dUdz ; dVdx - dUdy )
  %       = (   Ox  ;  Oy  ;     ZETA    )
- % 3 output files are created.
+ % 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:
-Line 17
+Line 30
  % ./netcdf-files/<SNAPSHOT>/OMEGAY.<netcdf_domain>.<netcdf_suff>
  % ./netcdf-files/<SNAPSHOT>/ZETA.<netcdf_domain>.<netcdf_suff>
  %
- % 06/07/2006
+ % 2006/06/07
  % gmaze@mit.edu
  %
+ % Last update:
+ % 2007/02/01 (gmaze) : Fix bug in ZETA grid and add compatibility with A-grid
+ %
- function [] = B_compute_relative_vorticity(snapshot)
+ % 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
+ global sla netcdf_UVEL netcdf_VVEL netcdf_domain netcdf_suff griddef
  pv_checkpath
-Line 41 
 pathname = strcat('netcdf-files',sla,sna
+Line 90 
 pathname = strcat('netcdf-files',sla,sna
  ext      = strcat('.',netcdf_suff);
- %% Load files:
+ %% Load files and axis:
  ferfile          = strcat(pathname,sla,filU,ext);
  ncU              = netcdf(ferfile,'nowrite');
  [Ulon Ulat Udpt] = coordfromnc(ncU);
-Line 52 
 ncV              = netcdf(ferfile,'nowri
+Line 101 
 ncV              = netcdf(ferfile,'nowri
  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
-Line 73 
 if toshow,disp('Dim'),end
+Line 144 
 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:
-Line 83 
 ZETA = zeros(nz,ny-1,nx-1).*NaN;
+Line 152 
 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
+ for iz = 1 : nz
    if toshow
      disp(strcat('Computing \zeta at depth : ',num2str(Udpt(iz)),...
                  'm (',num2str(iz),'/',num2str(nz),')'   ));
-Line 95 
 for iz=1:nz
+Line 167 
 for iz=1:nz
    % 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
-Line 113 
 for iz=1:nz
+Line 195 
 for iz=1:nz
        ZETA(iz,iy,ix) = dVdx - dUdy;
      end %for ix
    end %for iy
  end %for iz
  %%%%%%%%%%%%%%
-Line 197 
 Ox = zeros(O_nz,O_ny(1),O_nx(1)).*NaN;
+Line 278 
 Ox = zeros(O_nz,O_ny(1),O_nx(1)).*NaN;
  Oy = zeros(O_nz,O_ny(2),O_nx(2)).*NaN;
  %%%%%%%%%%%%%%
- %% Horizontal components:
+ %% Computation:
  %% Vertical grid differences:
  dZ   = diff(Udpt);
-Line 206 
 Odpt = Udpt(1:O_nz) + dZ/2;
+Line 287 
 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}(:,:,:);
+ V = ncV{4}(:,:,:);
-          Ox = - ( V(2:O_nz+1,:,:) - V(1:O_nz,:,:) ) ./ dZ_3D;
+ 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}(:,:,:);
+ U = ncU{4}(:,:,:);
-          Oy = ( U(2:O_nz+1,:,:) - U(1:O_nz,:,:) ) ./ dZ_3D;
+ 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
-Line 326 
 nc{ncid}.uniquename    = ncchar(uniquena
+Line 422 
 nc{ncid}.uniquename    = ncchar(uniquena
  nc{ncid}(:,:,:)        = Oy;
  nc=close(nc);
+ close(ncU);
+ close(ncV);
  % Outputs:
  OMEGA = struct(...
-Line 336 
 OMEGA = struct(...
+Line 433 
 OMEGA = struct(...
      );
  switch nargout
   case 1
-   varargout(1) = OMEGA;
+   varargout(1) = {OMEGA};
  end

 Legend:



Removed from v.1.3
 


changed lines


 
Added in v.1.4
 Legend:



Removed from v.1.3
 


changed lines


 
Added in v.1.4
-Removed from v.1.3
+Added in v.1.4

	ViewVC Help
Powered by ViewVC 1.1.22