NS_PERIODIC boundary conditions

[isabelmonteiro]

I'm using ROMS in the western Iberia coast in an upwelling configuration. Because I have a rather small domain I've imposed NS_PERIODIC boundary conditions (with bathymetry and pn,pm adjusted).
Somehow when I consider NS_PERIODIC the history file changes the bathymetry and puts at the western boundary h=10m and at the eastern boundary the bathymetry of the western boundary!.
To better show my problem with NS_PERIODIC, I've 3 figures with h. "Grid_bathy.jpg" is the bathymetry from the grid file, "NS_periodic_his_bathy.jpg" is the bathymetry read from history file with NS_PERIODIC boundary conditions where the wrong h appears at Western and Eastern boundaries and finally "without_NS_periodic_his_bathy.jpg" is the bathymetry read from history file without NS_PERIODIC boundary conditions with no problem with the h at the boundaries
I know that NS_PERIODIC condition is not used by the majoroty of ROMS community, but is a possibility offered by ROMS...

without_NS_periodic_his_bathy.jpg (46.58 KiB) Viewed 22293 times

[jcwarner]

when roms reads in the grid, it will exchange h along the boundaries (see ROMS/Utility/get_grid.F near lines 267). Just as a test, try
#undef NS_PERIODIC
#define EW_PERIODIC

and see if it exchanges in the other direction. if so, load your grid into matlab and see that lon_rho(1,1) and lat_rho(1,1) are indeed the lower left corners.

[isabelmonteiro]

I did:
#undef NS_PERIODIC
#define EW_PERIODIC
and I get no changes in the bathymetry from history file.

I've checked lon_rho(1,1) and lat_rho(1,1) and corresponds to the domain southwestern corner (lower left)

[jcwarner]

i just did an svn update for roms and ran the shoreface test case. that has a NS_PERIODIC call in it. it worked fine. no weirdness like you are getting. can u try that test case? please check your cppdefs carefully to make sure you are not re-defining something twice. -j

[arango]

I also checked the NS_PERIODIC boundary conditions are they are working for me. I don't see anything anomalous. I think that the problem here is that you have a flaw configuration. Setting a north-south periodic boundary conditions in a realistic application of the Iberia peninsula doesn't make much sense to me. You need to approach modeling this application in a different way.

I periodic boundary applications we need to make sure that the bathymetry, land/sea mask, and metric factors (pm, pn, dmde, dndx, angler) are all periodic . That's the reason that you get different values between grid and history files. We need to apply periodic boundary conditions. You don't have any choice in the matter. Otherwise, your dynamical fluxes are no periodic and you will be incorporating bogus circulation in your simulations. If you think more about this, you will see why your configuration is not appropriate. This probably will work for flat bathymetry and symmetric north and south land/sea masking. Then, you need to worry about the western boundary. But what is the point? This will not give you any realistic circulation for the Iberia peninsula region.

The only realistic application with periodic boundary E-W conditions possible is modeling the ACC around Antarctica. However, we still need to fine tune all the above grid field to avoid bogus circulations.

You need to reconsider your configuration and explore the possibility of taking open boundary conditions from a larger model. I don't see any other choice for you to get realistic circulations... This is what mostly all users of ROMS do in realistic regional applications. You can have highly idealized applications, but you need to have an idealized expression for all the grid fields mentioned above.

[isabelmonteiro]

I will consider your suggestions.
I was considering doing something analogous to the work of Kuebel-Cervantes and Allen for northern California (Deep-Sea Research II 53 (2006) 2956-2984) they used ROMS with periodic NS, or some POM users like Oke et al. for Oregon coast (J. Phys. Oceanogr 32 1383-1403 (2002)).
But you are right I'm getting away from the Iberian circulation

[kate]

The three plots in the first post look the same to me. What is the difference you are seeing? Can you plot a difference field?

[isabelmonteiro]

Is just the first and last columns, see att difference with W and E zooms

[isabelmonteiro]

I did the shoreface test (ROMS version 3.2) I get the same problem. It probably is a problem in my grid parameters that make history file confused

[jcwarner]

i am not clear. Are you saying that if you do an svn update, run shoreface, that the history file will have the grid edges exchanged E-W?

[isabelmonteiro]

yes.

[jcwarner]

this sounds very strange. not sure how to check on this, if i can not reproduce the problem. can u tar the *.f90 files and post them to me at:

If an local user, a member of the Woods Hole Science Center, wishes to receive files from an external user, the Public, the external user does the following:
C:\ ftp ftp.whoi.edu
login name: anonymous
password: email address
ftp> cd pub/usgs/filedrop
ftp> hash (this command simply writes # signs across your screen while ftp is working)
ftp> bin (this command allows binary image transfer)
ftp> put filename.ext (use mput if multiple files to deposit)
ftp> bye

[isabelmonteiro]

Did you find anything with the material that I put on the ftp site?

[jcwarner]

sorry. did not know u put anything there. I got your code. Compiled SHOREFACE and it ran fine.
Then i loaded your grid:
roms_grd_Roca_old_data_20100720_smth02_periodic.nc
It is flipped.
h(1,1) =10.
h(1,end)=10.
so the bottom row is all 10's. However, your plots show that the 10's (ie land) should be on the east side.
Do h=fliplr(h.');
and replace the h in your grid and try again. Let me know how that goes.
-j

[isabelmonteiro]

I did what you suggested (also for all grid parameters, lat_rho,lon_rho,pn,pm,f,etc, because they were all flipped ) and I still have the same problem land at west.
Meanwhile I also tested NS_PERIODIC with a flat bathymetry with land at eastern boundary and in history file I also have the western bathymetry changed with eastern bathymetry. It seems independent of the grid (it is me that constructed both grids, maybe I'm doing the same mistake).

I’ve putted my "upwelling.h" at the ftp site

isabel

[jcwarner]

i am not sure what i am supposed to do with the upwelling.h.

Take the code you sent.
edit makefile
change ROMS_APPLICATION?=SHOREFACE
make clean
make
./oceanS < ROMS/External/ocean_shoreface.in

when done,
ncload (or however) ocean_his.nc
pcolor(h)
and post that.

(i did these steps and all was fine.)
what ver matlab / netcdf are you using?

[arango]

Do you have a squared grid, Lm=Mm? I have mentioned this so many times in this forum and workshops to avoid at all cost having a squared grid application in ROMS. Some of the Matlab NetCDF interfaces and GUIs out there made the unfortunate decision of working with arrays dimensioned in row-major order, like array(Y,X), array(J,I), or array(lat,lon).

In computers, multidimensional data is stored linearly in memory. Fortran and Matlab stores data in column-major order and NOT in row-major order like the C-language. Numerical efficiency and page faulting is tied to this concept. Recall that ROMS is written in Fortran. If you have a square grid, you run into the danger that you may write tranposed arrays (order shown above) instead of array(X,Y), array(I,J), or array(lon,lat). This can give the wrong physics if you write the transpose of the Coriolis parameter.

The matlab scripts provided in the ROMS repository are all consistent with column-major order. These are the only ones that we support and processes the data as needed in the Fortran environment. You are free to use other non-official scripts provided that you know what are you doing and that sufficient guidelines are coded to check the tranpose problem.

Now, when creating a grid you need to be wise where to put the origin and need to be consistent with rotation. Recall that ROMS transform the grid to a rectangle in xi- and eta-space:

Code: Select all

               Mm  ________________
                  |                |
                  |                |
              eta |                |
                  |                |
                  |________________|
                O(1,1)             Lm
                          xi

Notice tha the origin is always in the left-bottom corner. This is going to govern the transformation from (xi,eta) to (lon,lat), curvilinear rotation angle, geographical true north and east, boundary conditons, normal and tagential flow, and vector and tensors rotations. I have seen so many users ignoring this basic geometric rule. It is very easy in SeaGrid or other grid generation packages to position the corners and the order in the wrong way.

[tony1230]

hi buddies
its a nice topic.i have questions with regard to arango‘s comments.
1.you said that if we have a squared grid the Fortran and Matlab would
transposed arrays (X,Y) into (Y,X).the point i
want to know is that the event you mentioned going to happen in any rate or by chance？
2.

This can give the wrong physics if you write the transpose of the Coriolis parameter.

i have no idea or concept about the mechanism between them.can you explains it to we all?
thank you in advance and any comments are appreciated from all of you！

[wilkin]

Hernan's comment is not that ROMS or Matlab automatically flip anything. It's just that if you have a square grid you won't immediately get a Matlab error if you unwittingly switch things by mistake. With a square grid you have to be more careful if, like me, you tend to rely on Matlab throwing the error...

Code: Select all

??? Error using ==> mtimes
Inner matrix dimensions must agree.  

... when you are careless.

[arango]

Then if you by mistake do a north-south flip the Coriolis parameter, f, in an squared grid, you will get the incorrect Coriolis force. Just think in terms of simple Newtonian physics to see what happen to a rotating frame of reference.

The Coriolis parameter parameter goes to zero at the Equator. It is positive in the northern hemisphere. Any change on the direction of its gradient will affect flow rotations, wave guides, and so on. You will get absurd circulations.

Users of ROMS have done this mistake in the past and waisted a lot of time finding the cause. Pretty scary stuff

[isabelmonteiro]

Sorry I was out of the office.
I did the shoreface example with my grid file and I still get the shift bathymetry at history file. Did you use my grid (from the ftp site) to do the shoreface example ?

[isabelmonteiro]

The NS_PERIODIC boundary conditions impose a exchange of h along east and west boundaries at history file.

I´ve got several suggestions:

1st suggestion: To impose EW_PERIODIC instead of NS_PERIODIC and see if h gets exchange.
result: No wrong h at history file=> no problem with EW_PERIODIC.
2nd suggestion: Try the upwelling test case with EW_PERIODIC changed to NS_PERIODIC and the domain walls at
the E and W boundaries
result: no problems but the E and W walls are the same anyway
3rd suggestion: Test shoreface test case (that has a NS_PERIODIC condition)
result: the same shift h at history file
4th suggestion: The square grid and grid manipulation that I do with my grid to be able to impose
NS periodic are a source of error column and row exchange and probably the reason why I get this kind of behaviour
result: I've tried all previous steps with a rectangular grid with 400 columns per 300 rows
directly from seagrid (that has a column major order consistency) and I did NO grid mannipulation at all and
I still get h exchange between E and W boundaries

Do I have to give up NS_PERIODIC with a realistic grid

[arango]

What function are you using to read/write NetCDF data in matlab I think that this is the source of confusion here. There are various NetCDF interface tools for matlab. You need to be fully aware that some of then transpose the arrays h(y,x), like in C-language.

The distributed scripts in the matlab repository to read (utility/nc_read.m) and write (utility/nc_write.m) NetCDF data in matlab uses the column-major order of Fortran, h(x,y). These routines are consistent with Matlab and Fortran order of storing multidimensional data linearly in memory. All the distributed scripts, except those associated with SeaGrid, uses the mexnc interface to process NetCDF data in Matlab.

I am suspecting that the problem is in your Matlab processing scripts and the source your miseries is in the NetCDF interface to Matlab that you are using. There are several of them: mexcdf, mexnc, snctools, NetCDF-java, and native Matlab (R2008b and newer). If you are using snctools, it will process the arrays in row-major order like the C-language This was an unwise choice, but it is too late to change it and you need to be fully aware of it. On the order hand, mexnc and native Matlab will give you the array in the same order that were defined and written to the NetCDF file. Please check the following PowerPoint presentation that discusses the issues with the NetCDF interface to Matlab.

Now if you are applying periodic boundary conditions in your processing scripts, you need to be cosistent with the ones applied in ROMS. Check either ROMS/Nonlinear/exchange_2d.F or ROMS/Nonlinear/exchange_3d.F.