Ocean Modeling Discussion

ROMS/TOMS

Search for:
It is currently Tue Jul 23, 2019 3:44 pm




Post new topic Reply to topic  [ 12 posts ] 

All times are UTC

Author Message
 Post subject: Curvilinear Coordinates
PostPosted: Tue Feb 28, 2006 6:54 pm 
Offline
Site Admin
User avatar

Joined: Wed Feb 26, 2003 4:41 pm
Posts: 1079
Location: IMCS, Rutgers University
This is to remind you that ROMS/TOMS is formulated in general horizontal curvilinear coordinates. It is possible to formulate any application in Cartesian, spherical, or polar coordinates. You can activate a curvilnear application by turning on CPP option CURVILINEAR. Curvilinear applications include extra terms in the advection and any vector or tensor needs to be rotated to model generic coordinates (XI,ETA).

In curvilinear applications, the coordinate rotation angle is found in the ROMS grid and history NetCDF files. See NetCDF variable angle. It is defined as the counterclockwise angle (radians) between the XI-axis and true EAST at RHO-points.

Therefore, to transform between (XI,ETA) coordinates to (LON,LAT) coordinates, vectors and tensors need to be rotated according to

Code:
u(LON,LAT)=u(XI,ETA)*cos(angle(i,j))-v(XI,ETA)*sin(angle(i,j))
v(LON,LAT)=u(XI,ETA)*sin(angle(i,j))+v(XI,ETA)*cos(angle(i,j))


conversely,

Code:
u(XI,ETA)=u(LON,LAT)*cos(angle(i,j))+v(LON,LAT)*sin(angle(i,j))
v(XI,ETA)=u(LON,LAT)*sin(angle(i,j))-v(LON,LAT)*cos(angle(i,j))


Then, you need to consider that ROMS uses a staggered C-grid. That is, the vector components are not located in the same place:

Code:
                ------- V(i,j+1)-------
                |                     |
             U(i,j)     RHO(i,j)   U(i+1,j)
                |                     |
                -------- V(i,j)--------


If plotting, you have two choices:

1) Draw vectors at PSI-points (cell's corners):

Code:
 
          Upsi(i,j) = 0.5 * ( U(i,j-1) + U(i,j) )
          Vpsi(i,j) = 0.5 * ( V(i-1,j) + V(i,j) )


2) Draw vectors at interior RHO-points (cell's center):

Code:
           Urho(i,j) = 0.5 * ( U(i,j) + U(i+1,j) )
           Vrho(i,j) = 0.5 * ( V(i,j) + V(i,j+1) )


It is more physical to draw vectors at PSI-points. It will be consistent with the physical boundaries of the domain. It will also show the behavior around Land/Sea masking, if any.

Now, the actual world is assumed to be at RHO-points. This will allow a well posed interpolation of vectors/tensors to C-grid locations. For instance, if you are providing wind components (Uwind,Vwind) to ROMS, they must be at RHO-points. So we can use the rotation angle to convert from (LON,LAT) to (XI,ETA) or viceversa.

ROMS assumes that all the input vectors are already in (XI,ETA)coordinates. At output ROMS, writes vectors in (XI,ETA) coordinates at the appropriate C-grid location. To visualize vectors the User needs to interpolate to either PSI- or RHO-locations and rotate back to (LON,LAT) coordinates, using above formulas.


Last edited by arango on Fri Apr 07, 2006 6:53 pm, edited 1 time in total.

Top
 Profile  
Reply with quote  
 Post subject: Curvilinear grid errors
PostPosted: Mon Mar 06, 2006 4:26 pm 
Offline

Joined: Fri Oct 28, 2005 9:09 pm
Posts: 1
Location: Woods Hole Oceanographic Institution
I have been using the "SeaGrid" matlab code to generate a curvilinear grid, and it always returns a field of departure from orthogonality. I am wondering what sort of maximum departure is "good enough". After some web surfing, the only thing I have found is an exchange on the POM FAQ web site suggesting that this is not too big a problem.

Can someone point me to an appropriate reference or handy criterion on this?


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Mar 06, 2006 5:11 pm 
Offline
Site Admin
User avatar

Joined: Wed Feb 26, 2003 4:41 pm
Posts: 1079
Location: IMCS, Rutgers University
Well, this is a good question. Unfortunately, it depens on the application resolution and coastal topology. Personally, I always try to avoid big grid distorsions and I keep the orthogonality about 95 percent or better. Is it also important to keep the grid aspect ratio (dx/dy) as close to one as possible. The important issue here is that the solution should be independent of grid topology.

Nowadays, I try to avoid curvilinear grids. I prefer to rotate rectangles and use land/sea masking to define a coastline. This is still a curvilinear grid but the grid angle is near constant. Notice my emphasis on rectangles. I avoid having square applications at all cost. This is because many of us use matlab or similar software and is very easy to transpose the coordinates of a matrix. This is hard to detect in ROMS during reading because both Lm and Mm have the same value. I know of users transposing the Coriolis parameter to disastrous consequences.

Also notice that in boundary fitted coastlines there are usually big distortions on bays and very fine resolution on headlands. This kind of behavior is undesired in small coastal applications. I prefer to have land/sea masking instead.


Last edited by arango on Mon Mar 06, 2006 8:54 pm, edited 1 time in total.

Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Mar 06, 2006 7:27 pm 
Offline
User avatar

Joined: Wed Jul 02, 2003 5:29 pm
Posts: 3633
Location: IMS/UAF, USA
The Seagrid program only calls "rect" once while the Fortran version calls it more times. The errors should drop quite rapidly with each successive call, but you get to a point where the thing blows up if the starting errors are too small, so it's nice to be able to tune the number of calls to one less than "too many".

The circle test problem pitted a curvilinear circle against a staircase circle. The rectangularity errors didn't seem to be the limiting factor. Instead, it was the large dx in the "bays", as mentioned by Hernan. The staircases provide errors of their own, so there's no perfect solution as long as we're dealing with structured grids.

I'm not afraid of square grids - just debug your analysis code beforehand.


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Mar 13, 2006 2:11 pm 
Offline
User avatar

Joined: Mon Apr 28, 2003 5:44 pm
Posts: 491
Location: Rutgers University
When generating orthogonal grids you should be careful to ensure the geometry of the boundaries you input is such that you input square corners. If not, there is nothing the boundary-fitting algorithm can do to make the cell in the corner orthogonal, and it will dominate the errors and propagate a lack of orthogonality to many neighbouring cells.

_________________
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Tue Mar 14, 2006 1:50 am 
Offline
User avatar

Joined: Wed Jul 02, 2003 5:29 pm
Posts: 3633
Location: IMS/UAF, USA
Here's an example grid:
Image
In this case, the corners in the NE are not orthogonal, leading to some orthogonality error there. However, I believe the larger errors are due to the resulting jump in dx along that NE wall. This drastic change in dx is not a good thing and probably leads to just as much trouble as steep bathymetry.

Here's one case where Seagrid might be better than the Fortran code - the boundary line you draw in Seagrid will be through the outisde rho image points. The grid is computed at fine resolution, leaving the largest errors outside the grid. The fortran code is expecting a boundary through the wall psi points, putting the largest errors inside the grid and copying them to the outside image points.


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Mar 20, 2006 9:07 pm 
Offline
User avatar

Joined: Tue Jul 01, 2003 4:12 am
Posts: 515
Location: NIWA
Well someone has to ask the obvious, naive question: What do these orthogonality errors look like in terms of flow distortion or whatever?


Top
 Profile  
Reply with quote  
PostPosted: Thu Jul 24, 2014 7:56 am 
Offline

Joined: Mon Jul 22, 2013 4:54 am
Posts: 19
Location: Ocean University of China
I think may be the code for transforming real velocity to model velocity should be as follows, but I am not sure...

u(XI,ETA) = u(LON,LAT)*cos(angle(i,j))+v(LON,LAT)*sin(angle(i,j))
v(XI,ETA) = (u(LON,LAT)*sin(angle(i,j))-v(LON,LAT)*cos(angle(i,j)))*-1


Top
 Profile  
Reply with quote  
PostPosted: Mon Jul 28, 2014 3:20 pm 
Offline
User avatar

Joined: Mon Apr 28, 2003 5:44 pm
Posts: 491
Location: Rutgers University
I do the angle rotation this way (in Matlab but you get the idea) using complex numbers. The reverse transformation (geographic to roms) just needs the opposite sign on the angle.

% rotate coordinates from roms xi,eta grid directions to geographic east,north
% (but only after averaging/interpolating u,v to the rho points)
uveitheta = (u_roms+1i*v_roms).*exp(1i*angle); % 1i = sqrt(-1) in Matlab
u_east = real(uveitheta);
v_north = imag(uveitheta);

_________________
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu


Top
 Profile  
Reply with quote  
PostPosted: Tue Jul 29, 2014 2:37 am 
Offline

Joined: Mon Jul 22, 2013 4:54 am
Posts: 19
Location: Ocean University of China
wilkin wrote:
I do the angle rotation this way (in Matlab but you get the idea) using complex numbers. The reverse transformation (geographic to roms) just needs the opposite sign on the angle.

% rotate coordinates from roms xi,eta grid directions to geographic east,north
% (but only after averaging/interpolating u,v to the rho points)
uveitheta = (u_roms+1i*v_roms).*exp(1i*angle); % 1i = sqrt(-1) in Matlab
u_east = real(uveitheta);
v_north = imag(uveitheta);


Appreciation for wilkin's reply! I use Matlab to prepare NC files for ROMS as you do, your code is very simple and clear and I'll recode my program following your idea.

Thanks again!


Top
 Profile  
Reply with quote  
PostPosted: Sat Jul 02, 2016 8:53 pm 
Offline

Joined: Thu Oct 01, 2015 4:41 am
Posts: 5
Location: KAUST
I'm confused what should be the positive direction for the angle. Does "counterclockwise angle (radians) between the XI-axis and true EAST at RHO-points" mean that angle is calculated from the XI-axis to the the EAST direction with counterclockwise being positive direction?

In other words, if I were to use easygrid and specify rotangle=alpha, then I would save into the ocean grid netcdf file angle variable with opposite value, which is in matlab would be:
angle = zeros(size(lat_rho));
angle(:) = -1*alpha

Could someone clarify this please?


Top
 Profile  
Reply with quote  
PostPosted: Sun Jul 03, 2016 11:56 am 
Offline
User avatar

Joined: Mon Apr 28, 2003 5:44 pm
Posts: 491
Location: Rutgers University
These are the attributes I add to my ROMS grid files to remind me how angle is defined.

Quote:
double angle(eta_rho, xi_rho) ;
angle:long_name = "angle ROMS xi axis is rotated anticlockwise from due east" ;
angle:equation = "u(roms)=Real((Ueast+i*Vnorth)*exp(-i*angle)) and v(roms)=Imag((Ueast+i*Vnorth)*exp(-i*angle))" ;
angle:units = "radian" ;

_________________
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu


Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 12 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 2 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
Powered by phpBB® Forum Software © phpBB Group