angle in grid file

 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 XIaxis 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
In the ROMS grid file, the description of the variable 'angle' is the angle between XIaxis 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
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 xieta grid and the lonlat 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
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 RHOpoints 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)
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 RHOpoints 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 northsouth and "size_x" (also km] in eastwest 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 Ydirections, counting
% only inner RHOpoints (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 Ydirection 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 Xdirection at the center of the grid points
% exactly to the east);
% flip_xy: flag to transpose the grid: 0 do nothing; 1 flip Xdirection
% and transpose; 2 flip Ydirection and transpose. Thus 1 and
% 2 corresponds to 90degree 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 xaxis 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=dLnX12*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=dLnX2*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 Matlabstyle array
dLnY1=dLnY12*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 quasiFortran
dLnY=lone(i,j+1)lone(i,j); % code of the left does the job.
if (dLnY > pi)
dLnY=dLnY2*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*piatan(a11/a21); % For a perfectly othogonal
elseif (a21 > abs(a11)) % curvilinear grid "ang1" and
ang1 = 0.5*piatan(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=ang1pi; % 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) = ang1ang2 ; %% 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

 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.
Thank you very much.
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 xieta 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 xieta
is orthogonal curvilinear system, therefore it can be treated as "locally Cartesian".]
Alternatively, consider moving along eta coordinate by increment Deta: in nonrotated case
xi points to the east, eta to the north, alpha=0. Eastwest 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 northwest), and
in displacement in Latitide is (dLat/deta)*Dlat, so angle (this time between etadirection
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 finitedifference expressions to approximate these
derivatives (+ technicalities to avoid division by zero and handling 360degree 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.
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 xieta 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 xieta
is orthogonal curvilinear system, therefore it can be treated as "locally Cartesian".]
Alternatively, consider moving along eta coordinate by increment Deta: in nonrotated case
xi points to the east, eta to the north, alpha=0. Eastwest 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 northwest), and
in displacement in Latitide is (dLat/deta)*Dlat, so angle (this time between etadirection
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 finitedifference expressions to approximate these
derivatives (+ technicalities to avoid division by zero and handling 360degree 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.
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 wavecurrent 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 XIaxis and EAST at RHOpoints.
!
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
[ROMS/Utility/get_grid.F]
!
! Read in angle (radians) between XIaxis and EAST at RHOpoints.
!
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
Re: angle in grid file
True indeed, but one could teach ROMS to read it in degrees and convert if necessary.
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
j

 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.