## angle in grid file

General scientific issues regarding ROMS

Moderators: arango, robertson

Message
Author
shrikantmp
Posts: 26
Joined: Mon Jan 27, 2014 9:50 pm
Location: Indian Institute of Science

### angle in grid file

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
Posts: 3780
Joined: Wed Jul 02, 2003 5:29 pm
Location: CFOS/UAF, USA

### Re: angle in grid file

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
Posts: 185
Joined: Fri Nov 14, 2003 4:57 pm

### Re: angle in grid file

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.

[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
Posts: 26
Joined: Mon Jan 27, 2014 9:50 pm
Location: Indian Institute of Science

### Re: angle in grid file

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
Posts: 185
Joined: Fri Nov 14, 2003 4:57 pm

### Re: angle in grid file

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.
Last edited by shchepet on Fri Mar 11, 2016 6:00 pm, edited 2 times in total.

jcwarner
Posts: 855
Joined: Wed Dec 31, 2003 6:16 pm
Location: USGS, USA

### Re: angle in grid file

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]
!
!
CASE ('angle')
gtype=r2dvar
& var_name(i), var_id(i), &
& 0, gtype, Vsize, &
& LBi, UBi, LBj, UBj, &
& Fscl, Fmin, Fmax, &
# endif
& GRID(ng) % angler)

-john

kate
Posts: 3780
Joined: Wed Jul 02, 2003 5:29 pm
Location: CFOS/UAF, USA

### Re: angle in grid file

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

jcwarner
Posts: 855
Joined: Wed Dec 31, 2003 6:16 pm
Location: USGS, USA

### Re: angle in grid file

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
Posts: 26
Joined: Mon Jan 27, 2014 9:50 pm
Location: Indian Institute of Science

### Re: angle in grid file

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