angle in grid file

[shrikantmp]

Hi ROMS users

In the ROMS grid file, the description of the variable 'angle' is the angle between XI-axis and EAST. My question is how do you determine this angle? If I want to run a simulation for Indian Ocean, how will I determine this angle? Also, should the units of angle be degrees or radians? Please help.

Thanks in advance

[kate]

How are you making your grid? Here's the code I use (radians):

Code: Select all

!  calculate the angle between the xi-eta grid and the lon-lat grid
!  at rho points.

        do j = 1, Mm
          do i = 1, Lm
            a1 = lat_u(i+1, j) - lat_u(i, j)
            a2 = lon_u(i+1, j) - lon_u(i, j)
            if (abs(a2) .gt. 180.) then
              if (a2 .lt. -180. ) then
                a2 = a2 + 360.
              else
                a2 = a2 - 360.
              endif
            endif
            a2 = a2 * cos(0.5*DTOR*(lat_u(i, j) + lat_u(i+1, j)))
            angle(i, j) = atan2(a1, a2) 
          enddo
        enddo

        do j = 1, Mm
          do i = 1, Lm
            a2 = lat_v(i, j) - lat_v(i, j+1)
            a1 = (lon_v(i, j) - lon_v(i, j+1))
            if (abs(a1) .gt. 180.) then
              if (a1 .lt. -180. ) then
                a1 = a1 + 360.
              else
                a1 = a1 - 360.
              endif
            endif
            a1 = a1 * cos(0.5*DTOR*(lat_v(i, j) + lat_v(i, j+1)))
            angle(i, j) = 0.5*(angle(i, j) + atan2(a1, -a2))
          enddo
        enddo

[shchepet]

See the bottom half of the following routine.

Angle can be either in degrees or radians -- it is just a matter of personal preference.
ROMS code itself does not use angle, but various pre- and post- processing routines use
it to rotate vectors -- velocities if one needs to prepare files for initial and
boundary conditions for ROMS, wind stress or velocity; also plotting packages obviously
need it for vector plots. So as long as there routines consistent with the grid file,
or can read attribute "units" and handle both degrees or radians, everything is just
fine. Degrees are easier to comprehend when looking at grid file by ncview.


In the code below: arrays lonr,latr are longitude and latitude at RHO-points defined
at full ranges of indices (interior + one boundary row at each side); lone,late are
longitude and latitude defined at VORTICITY points over EXTENDED range of indices, that
is having one point more in each direction that lonr,latr. The rest is more or less
self explanatory. A bit of care is needed around dateline (this routine fully capable
to handle it) north pole (this routine us known handle it as well, except the pole point itself, of course)

Code: Select all

function [lone,late, lonr,latr, pm,pn, ang,orterr] = compute_grid(  ...
                                         nx,ny, size_x,size_y, cent_lat, ...
                                         tra_lon,tra_lat,alpha, flip_xy, ...
                                                         compute_angle )

% Produces a logically rectangular curvilinear orthogonal grid with minimal
% distortion. At first it creates a patch of Mercator grid centered at a user
% specified latitude "cent_lat" on Greenwich meridian, sized as "size_y" [km]
% in north-south and "size_x" (also km] in east-west direction measured along
% the longitudinal line passing through the center of the patch.  To center
% the patch at Equator, one shouls set cent_lat=0; positive(negative) values
% move it north(sourth), which also makes it more tapered toward one of the
% poles. Thus, an infinitely small patch is almost rectangular if placed on
% Equator, while behaving more and more like a sector of polar coodinates if
% placed closer to one of the poles.  However in any case it is a perfectly
% Mercator grid with both dx,dy ~ cos(latitude) at any location.

% Then this patch is first moved toward Equator as a solid object, where it
% is asimuthally rotated around it center by angle "alpha", after which it is
% moved to the north or south to latitude "tra_lat", and finally moved east or
% west to longitude "tra_lon".

%   inputs: nx,ny: numbers of grid points in X- and Y-directions, counting
%                  only inner RHO-points (boundary rows are excluded) -- same
%                  as Lm,Mm in ROMS code, with the exception of the possibility
%                  to switch their roles if "flip_xy" is set to transpose the
%                  grid;
%   size_x,size_y: domain size in X- and Y-direction measured in [km];
%        cent_lat: latutude of of the initial location of the grid center;
%         tra_lat: desired latitude of grid center;
%         tra_lon: desired longitude of grid center;
%           alpha: asimuthal angle of grid direction (alpha=0 means that
%                  positive X-direction at the center of the grid points
%                  exactly to the east);
%         flip_xy: flag to transpose the grid: 0 do nothing; 1 flip X-direction
%                  and transpose; 2 flip Y-direction and transpose.  Thus 1 and
%                  2 corresponds to 90-degree azimutal rotation in counter- or
%                  or clockwise directions (respectively) however this rotation
%                  is "logical" rather than physical, i.e., the two are the
%                  same for would be perfectly rectangular grid, but are not
%                  equavalent in the case of spherical.
%  compute_angle:  switch to proceed with computation of east angle, =1 yes,
%                  =0 no. This is merely to speed up the responsiveness of GUI
%                  update button because for the plotting purposes all what is
%                  needed is lon,lat, but not pm,pn, ang,orterr.


 disp(['compute_grid :: nx=',num2str(nx),     ' ny=',num2str(ny), ...
               ' size_x=',num2str(size_x), ' size_y=',num2str(size_y),...
            ' cent_lat=',num2str(cent_lat) ])
 disp(['        tra_lon=',num2str(tra_lon),' tra_lat=',num2str(tra_lat),...
               ' alpha=',num2str(alpha),  ' flip_xy=',num2str(flip_xy) ])

%% Initialize the grid, then move and rotate it to proper location and
%% orientation.

 pi=3.14159265358979323846; deg2rad=pi/180; r_earth = 6371315.;

 [lone,late, lonr,latr, pm,pn] = init_grid(nx,ny, size_x,size_y, cent_lat);
 if (flip_xy > 0)
   lone = flipdim(lone,flip_xy)'; lonr = flipdim(lonr,flip_xy)';
   late = flipdim(late,flip_xy)'; latr = flipdim(latr,flip_xy)';
   tmp=pm;   pm=flipdim(pn,flip_xy)';  pn=flipdim(tmp,flip_xy)'; clear tmp;
 end

 [lonr,latr] = move_and_turn(tra_lon,tra_lat,alpha,cent_lat, lonr,latr);
 [lone,late] = move_and_turn(tra_lon,tra_lat,alpha,cent_lat, lone,late);

%% Compute angles of local grid positive x-axis relative to east

 [ii,jj]=size(latr) ; ang=zeros(ii,jj) ; orterr=zeros(ii,jj) ;

 if (compute_angle == 1)
   disp('Calculating angle ...')


%% ang=0.5*cos(latr);  %% <-- temporally
%%
%% a11 = ang .* ( lone(2:ii+1,2:jj+1) -lone(1:ii,2:jj+1) ... % cos(Lat)*dLon/dX
%%               +lone(2:ii+1,1:jj  ) -lone(1:ii,1:jj  ) ) ;
%%
%% a12 = ang .* ( lone(1:ii,2:jj+1) +lone(2:ii+1,2:jj+1) ... % cos(Lat)*dLon/dY
%%               -lone(1:ii,1:jj  ) -lone(2:ii+1,1:jj  ) ) ;
%%
%% a21 =  0.5 * ( late(2:ii+1,2:jj+1) -late(1:ii,2:jj+1) ... %     dLat/dX
%%               +late(2:ii+1,1:jj  ) -late(1:ii,1:jj  ) ) ;
%%
%% a22 =  0.5 * ( late(1:ii,2:jj+1) +late(2:ii+1,2:jj+1) ... %     dLat/dY
%%               -late(1:ii,1:jj  ) -late(2:ii+1,1:jj  ) ) ;

   for j=1:jj
     for i=1:ii
       dLnX1=lone(i+1,j+1)-lone(i,j+1);        % Most of the time "lone"
       if (dLnX1 > pi)                         % is expected to be a smooth
         dLnX1=dLnX1-2*pi;                     % function of its indices,
       elseif (dLnX1 < -pi)                    % however it may experience
         dLnX1=dLnX1+2*pi;                     % 2*pi jumps due to periodicity
       end                                     % resulting in large valies of
       dLnX=lone(i+1,j)-lone(i,j);             % differences dLnX,dLnY. If this
       if (dLnX > pi)                          % happes, 2*pi is added or
         dLnX=dLnX-2*pi;                       % subtracted as appropriate.
       elseif (dLnX < -pi)
         dLnX=dLnX+2*pi;
       end

       dLnY1=lone(i+1,j+1)-lone(i+1,j);        % Because of this possibility,
       if (dLnY1 > pi)                         % the elegant Matlab-style array
         dLnY1=dLnY1-2*pi;                     % syntax code above is commented
       elseif (dLnY1 < -pi)                    % out permanently, but still
         dLnY1=dLnY1+2*pi;                     % left here for Matlab lovers
       end                                     % to admire, while quasi-Fortran
       dLnY=lone(i,j+1)-lone(i,j);             % code of the left does the job.
       if (dLnY > pi)
         dLnY=dLnY-2*pi;
       elseif (dLnY < -pi)
         dLnY=dLnY+2*pi;
       end

       cff=0.5*cos(latr(i,j));
       a11=cff*(dLnX+dLnX1);                    % cos(Lat)*dLon/dX

       a12=cff*(dLnY+dLnY1);                    % cos(Lat)*dLon/dY

       a21=0.5*( late(i+1,j+1) -late(i,j+1) ... %     dLat/dX
                +late(i+1,j  ) -late(i,j  ) ) ;

       a22=0.5*( late(i,j+1) +late(i+1,j+1) ... %     dLat/dY
                -late(i,j  ) -late(i+1,j  ) ) ;
                                               % Note that this computation
%%     pm(i,j)=1./(r_earth*sqrt(a11^2+a21^2)); % of pm and pn overwrites their
%%     pn(i,j)=1./(r_earth*sqrt(a12^2+a22^2)); % previously computed analytical
                                               % counterparts from init_grid.m
       if (a21 < -abs(a11))
         ang1 = -0.5*pi-atan(a11/a21);         % For a perfectly othogonal
       elseif (a21 > abs(a11))                 % curvilinear grid "ang1" and
         ang1 =  0.5*pi-atan(a11/a21);         % "ang2" should be the same.
       else                                    % Here we chose to compute both
         ang1 = atan(a21/a11);                 % of them, use their average as
         if (a11 < 0.)                         % the final accepted value for
           if (a21 < 0.)                       % east angle and also compute
             ang1=ang1-pi;                     % and save their difference as
           else                                % the measure for orthogonality
             ang1=ang1+pi;                     % error.
           end
         end
       end

       if (a12 < -abs(a22))
         ang2 = 0.5*pi +atan(a22/a12);
       elseif (a12 > abs(a22))
         ang2 = -0.5*pi +atan(a22/a12);
       else
         ang2 = -atan(a12/a22);
         if (a22 < 0.)
           if (a12 < 0.)
             ang2 = ang2 +pi;
           else
             ang2 = ang2 -pi;
           end
         end
       end

       ang(i,j)=0.5*(ang1+ang2);  %% the two should be same or very close
       orterr(i,j) = ang1-ang2 ;  %% orthogonality error

     end  %% <-- for i=1:ii
   end  %% <-- for j=1:jj
   disp('fully updated compute_grid')
 else
   disp('updated lon,lat in compute_grid')
 end  %% <-- if compute_angle

[shrikantmp]

Thank you very much kate and shchepet. I used to use roms_tools but the algorithm it uses to calculate the angle was difficult to understand. Now I use matlab routines by Arango, but couldn't find an algorithm to compute angle.

Thank you very much.

[shchepet]

The rationale for computing angle is very simple:

Consider a small Dxi increment along the direction of ROMS grid curvilinear coordinate
xi: the associated displacement in longitudinal direction measured in km is
R * cos(Lat) * dLon/dxi * Dxi
where R is Earth radius; displacement in latitudinal direction is
R * dLat/dxi * Dxi.
Correspondingly the angle of rotation alpha is such that
tan(alpha) = (dLat/dxi)/(cos(Lat) * dLon/dxi)
is the ratio of the two.

[To avoid confusion about metric term -- cos(Lat) here and there -- think in terms of
local Cartesian coordinate in tangential plane touching the Earth at point Lon,Lat and
expressed in real distances measured in km. So local "x" direction points to the east,
"y" to the north, and xi-eta axes are rotated relative to them; then the increment in "x"
is R*cos(Lat)*Dlon; increment in "Y" is R*Dlat, and you are kind of rotating one
Cartesian system of coordinates relatively to another Cartesian system -- as xi-eta
is orthogonal curvilinear system, therefore it can be treated as "locally Cartesian".]

Alternatively, consider moving along eta coordinate by increment Deta: in non-rotated case
xi points to the east, eta to the north, alpha=0. East-west displacement associated with
increment Deta is -R * cos(Lat) * dLon/deta * Deta (minus sign is due to the fact that
eta direction rotated by a positive angle from the north points to the north-west), and
in displacement in Latitide is (dLat/deta)*Dlat, so angle (this time between eta-direction
and the north) is tan(alpha) = -[cos(Lat) * dLon/deta]/(dLat/deta).

Because grids are orthogonal, the these two angles should be the same,
-[cos(Lat) * dLon/deta]/(dLat/deta) = (dLat/dxi)/[cos(Lat) * dLon/dxi]
which is a way to check for orthogonality error.

All what you see above is different finite-difference expressions to approximate these
derivatives (+ technicalities to avoid division by zero and handling 360-degree jumps).
Note that we actually need not the derivatives, but their ratios, so dxi
and deta cancel out, leaving behind just ratios of finite differences.

[jcwarner]

I don't want to cause any confusion on this thread, but the Rutgers version of ROMS does internally use the variable 'angle' at several instances, so the units do matter. For example, if a user has a forcing file of Uwind_east and Vwind_north, and they impose that on a rotated ocean grid, then the code will use the grid angle to rotate the winds to the ocean grid orientation for the computation of surface stresses. We also use angle for the wave-current interactions, and sediment bottom stress computations. The angle needs to be in radians in the grid file.

[ROMS/Utility/get_grid.F]
!
! Read in angle (radians) between XI-axis and EAST at RHO-points.
!
CASE ('angle')
gtype=r2dvar
status=nf_fread2d(ng, model, ncname, GRD(ng)%ncid, &
& var_name(i), var_id(i), &
& 0, gtype, Vsize, &
& LBi, UBi, LBj, UBj, &
& Fscl, Fmin, Fmax, &
# ifdef MASKING
& GRID(ng) % rmask, &
# endif
& GRID(ng) % angler)



-john

[kate]

True indeed, but one could teach ROMS to read it in degrees and convert if necessary.

[jcwarner]

Kate- yes I agree. I think the more important part of my post was that ROMS uses the angle, not the units of it.
-j

[shrikantmp]

Thanks everyone. The rationale provided by shchepet is indeed easy to understand and implement. Thank you very much.