NS_PERIODIC boundary conditions

Bug reports, work arounds and fixes

Moderators: arango, robertson

Post Reply
Message
Author
isabelmonteiro
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

NS_PERIODIC boundary conditions

#1 Unread post by 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
without_NS_periodic_his_bathy.jpg (46.58 KiB) Viewed 18997 times
Attachments
NS_periodic_his_bathy.jpg
NS_periodic_his_bathy.jpg (46.98 KiB) Viewed 18997 times
Grid_bathy.jpg
Grid_bathy.jpg (48.39 KiB) Viewed 18997 times

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

Re: NS_PERIODIC boundary conditions

#2 Unread post by 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
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

Re: NS_PERIODIC boundary conditions

#3 Unread post by 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
Posts: 1172
Joined: Wed Dec 31, 2003 6:16 pm
Location: USGS, USA

Re: NS_PERIODIC boundary conditions

#4 Unread post by 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

User avatar
arango
Site Admin
Posts: 1347
Joined: Wed Feb 26, 2003 4:41 pm
Location: DMCS, Rutgers University
Contact:

Re: NS_PERIODIC boundary conditions

#5 Unread post by 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.

:idea: 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
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

Re: NS_PERIODIC boundary conditions

#6 Unread post by 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

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

Re: NS_PERIODIC boundary conditions

#7 Unread post by 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
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

Re: NS_PERIODIC boundary conditions

#8 Unread post by isabelmonteiro »

Is just the first and last columns, see att difference with W and E zooms
Attachments
diff_NSper_vs_without_his_bathy_zoomW.jpg
diff_NSper_vs_without_his_bathy_zoomW.jpg (16.57 KiB) Viewed 18922 times
diff_NSper_vs_without_his_bathy_zoomE.jpg
diff_NSper_vs_without_his_bathy_zoomE.jpg (15.85 KiB) Viewed 18922 times
diff_NSper_vs_without_his_bathy.jpg
diff_NSper_vs_without_his_bathy.jpg (14.62 KiB) Viewed 18922 times

isabelmonteiro
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

Re: NS_PERIODIC boundary conditions

#9 Unread post by 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
Posts: 1172
Joined: Wed Dec 31, 2003 6:16 pm
Location: USGS, USA

Re: NS_PERIODIC boundary conditions

#10 Unread post by 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
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

Re: NS_PERIODIC boundary conditions

#11 Unread post by isabelmonteiro »

yes.

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

Re: NS_PERIODIC boundary conditions

#12 Unread post by 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
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

Re: NS_PERIODIC boundary conditions

#13 Unread post by isabelmonteiro »

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

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

Re: NS_PERIODIC boundary conditions

#14 Unread post by 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
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

Re: NS_PERIODIC boundary conditions

#15 Unread post by 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
Posts: 1172
Joined: Wed Dec 31, 2003 6:16 pm
Location: USGS, USA

Re: NS_PERIODIC boundary conditions

#16 Unread post by 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?

User avatar
arango
Site Admin
Posts: 1347
Joined: Wed Feb 26, 2003 4:41 pm
Location: DMCS, Rutgers University
Contact:

Re: NS_PERIODIC boundary conditions

#17 Unread post by 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.

:idea: 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
Posts: 87
Joined: Wed Mar 31, 2010 3:29 pm
Location: SKLEC,ECNU,Shanghai,China

Re: NS_PERIODIC boundary conditions

#18 Unread post by 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! :D

User avatar
wilkin
Posts: 875
Joined: Mon Apr 28, 2003 5:44 pm
Location: Rutgers University
Contact:

Re: NS_PERIODIC boundary conditions

#19 Unread post by 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.
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu

User avatar
arango
Site Admin
Posts: 1347
Joined: Wed Feb 26, 2003 4:41 pm
Location: DMCS, Rutgers University
Contact:

Re: NS_PERIODIC boundary conditions

#20 Unread post by 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
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

Re: NS_PERIODIC boundary conditions

#21 Unread post by 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
Posts: 39
Joined: Wed Feb 18, 2004 3:50 pm
Location: Instituto de Meteorologia

Re: NS_PERIODIC boundary conditions

#22 Unread post by 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

:mrgreen: Do I have to give up NS_PERIODIC with a realistic grid :?:

User avatar
arango
Site Admin
Posts: 1347
Joined: Wed Feb 26, 2003 4:41 pm
Location: DMCS, Rutgers University
Contact:

Re: NS_PERIODIC boundary conditions

#23 Unread post by 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 :arrow: 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 :arrow: 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.

Post Reply