MPI code bug in interpolating station locations

[longmtm]

Comrades,

There seems to be a bug in 2.2 version for the station location interpolation.

Here, in my case, I have 126 stations for the Chesapeake Bay (Red) in the following link. In the code, the station input file can take (lat, lon) as input data. And ROMS is supposed to find the (Ipos, Jpos) in the code. I have divided the domain into 8 partitions in the run. The triangles are the input locations. The yellow *'s are the interpolated locations found by ROMS. The interesting thing is that 3 of them (the red triangles) are missing. And all 3 of them are close to or on the dividing lines. I checked the interpolation output results, they missing ones are interpolated to infinity. (Ios, Jpos)= (1e35,1e35) for the pink ones.

Here is the matlab figure:

http://131.118.211.30/chesroms/retrospe ... nCheck.fig




Any quick fix can be done for the code?

Thanks

Wen

[jivica]

Yip,
there is a problem when using more than 1 tile (*MPI*) at the connecting nodes with Lon,Lat , so you have to find locations in I and J before (*NOT LON,LAT) and put them in the station file. You can test this with serial run and then grab i,j model coordinates and feed them into station file and run MPI version...
Ivica

[longmtm]

I tried to modify the grid_coords.F in directory Utility, and it seemed to be working.

The old and new results are located at:

http://131.118.211.30/chesroms/retrospe ... nCheck.fig

and

http://131.118.211.30/chesroms/retrospe ... ck_new.fig

(you may have to type the URL)

and the modifled code is at:

http://131.118.211.30/chesroms/LatestRO ... tationfix/

Ofcause, I think this needs further testing for I may not understand the code. I hope the station outputs are okay.

Wen

[longmtm]

I have to take it back. The fix doesn't realy work. The out put of (lat lon) are correct but the (Ipos and Jpos) are not.

Wen

[longmtm]

Finally here is the quickest fix the problem,,

We only need to change two places in file grid_coords.F in Utility

1) use a larger search box, LBi, UBi, LBj,UBj instead of Istr, Iend, Jstr, Jend
2) replace spv by -1.0_r8 in the call of hindices() and mp_collect() i as the following (only this call , not for the other hindices callin the same file)


CALL hindices (ng, LBi, UBi, LBj, UBj, &
!--commented by Wen Long---better use a larger search box
! & Istr, Iend, Jstr, Jend, &
!------
& LBi, UBi, LBj, UBj, &
& GRID(ng)%angler, &
& GRID(ng)%lonr, &
& GRID(ng)%latr, &
& -1.0_r8, mc, Slon, Slat, Ista, Jsta)

!--- (-1.0_r8) is used here--collect the value not equal to -1.0_r8 instead
! of summing up in mp_collect

# ifdef DISTRIBUTE
CALL mp_collect (ng, model, mc, -1.0_r8, Ista)
CALL mp_collect (ng, model, mc, -1.0_r8, Jsta)

!--- (-1.0_r8) is used here--collect the value not equal to -1.0_r8 instead
! of summing up in mp_collect

# endif

[jivica]

Well,
For my application (Adriatic Sea in spherical) it is not working, your patch I ment. Still hindices does not find corresponding I,J in tiles.
The simplest solution for me is to run serial code for 1 iteration with stations defined in station.in like 1 1 lon,lat then get from out_sta.nc Ipos and Jpos (FRACTAL!!! values) and feed them in station.in file (1 0 I J) for parallel run...
That's it in 15 sec.
Ugly but it's working

[longmtm]

Are you using ROMS2.2?

The fix I mentioned works for me. Would you please try replacing two files with the following link and post back?

http://131.118.211.30/chesroms/LatestRO ... tationfix/

Thanks,

Wen

[arango]

Yes, there is a bug in the way that stations grid locations are computed from (lon,lat) in grid_coords.F. The problem only occurs at the tile boundaries. The fix proposed above is incorrect We cannot have spv=-1 in a call to mp_collect because it will introduce a parallel bug. It is not a good idea to modify mp_collect. This routine is very special and tricky.

The fix is somewhat different:

(1) Remove lines 61 and 62 in routine grid_coords.F

Code: Select all

!!    IF (Itile.eq.(NtileI(ng)-1)) Iend=Lm(ng)+1
!!    IF (Jtile.eq.(NtileJ(ng)-1)) Jend=Mm(ng)+1

(2) Change Iend and Jend arguments to hindices to Iend+1 and Jend+1. Notice that you need to change it in two places. The first one is for the floats:

Code: Select all

          CALL hindices (ng, LBi, UBi, LBj, UBj,                        &
     &                   Istr, Iend+1, Jstr, Jend+1,                    &
     &                   GRID(ng)%angler,                               &
     &                   GRID(ng)%lonr,                                 &
     &                   GRID(ng)%latr,                                 &
     &                   spv, mc, FLT(ng)%Flon, FLT(ng)%Flat,           &
     &                   Iflt, Jflt)

and the second one is for the stations:

Code: Select all

          CALL hindices (ng, LBi, UBi, LBj, UBj,                        &
     &                   Istr, Iend+1, Jstr, Jend+1,                    &
     &                   GRID(ng)%angler,                               &
     &                   GRID(ng)%lonr,                                 &
     &                   GRID(ng)%latr,                                 &
     &                   spv, mc, Slon, Slat, Ista, Jsta)

Thank you for finding and reporting this bug. I corrected the tar file corrections to ROMS 2.2. This fix will be corrected in the next version of ROMS. We are in the process of releasing a new version.

[isabelmonteiro]

I'm using ROMS 3.2 version,
I have a problem with stations location:
For a given Ipos, Jpos, my station location is lat=38.5391230508468 and lon=-11.0589999999828
But in history file lat_rho(Jpos,Ipos)=38.5186762826554 and lon_rho(Jpos,Ipos)=-11.072999999983

What am I doing wrong?

[jfiechter]

I am experiencing a similar (and probably related) problem with ROMS 3.4:

When specifying the station position in grid units i,j (FLAG=0), the output values for that station (e.g., h, lon_rho, lat_rho) correspond to the cell i+1, j+1 in the grid file. For example, if I specify X-POS=10, Y-POS=20 in the station file, then the latitude, longitude, and depth information in sta.nc for that location match the grid file values for cell i=11, j=21.

Has anybody else encountered this issue?

[arango]

I have some time to check the issue about the values of array(Ipos,Jpos), reported recently in this thread. The code is correct as it is The problem is a matter of interpretation and how NetCDF files stores array data. The positions (Ipos, Jpos) specified in the input script file stations.in are all assumed to be located at RHO-points. That is, if Ipos and Jpos are whole numbers and not fractions, the station is located at the center of grid cell (RHO-point) defined by Ipos and Jpos in ROMS native index coordinates and not NetCDF indices in the file. Recall, that ROMS uses a staggered Arakawa C-grid discretization. This is a numerical artifact and nothing to do with geography. If you put an instrument in the ocean, you specify its location by the (lon,lat) coordinates and you never take into consideration if this is a B-grid or C-grid when deploying the instrument. Now, ROMS is a finite volume model and we assume that any data (scalar or vector) that it is measured in the real world is located at RHO-points (center of the finite volume cell). As a matter of fact, we use this assumption routinely in bathymetry, preparing initial conditions and climatology, data assimilation, model-data comparisons, etc.

If you check wikiROMS, you will find a diagram of the horizontal index discretization in ROMS. You will notice the following (I,J) ranges for ROMS variables:

Code: Select all

RHO-points  0:L  0:M
PSI-points  1:L  1:M
  U-points  1:L  0:M
  V-points  0:L  1:M

Now if you carefully examine the NetCDF file, you will discover that there is no zero neither negative indices in the NetCDF variable dimension. This doesn't means that we cannot store such arrays in a NetCDF file. There is always a linear re-map of indices in each dimension. The data is always offset by the amount required so the fist value of the array is always written at the (1,1,...) location in the NetCDF file. We don't have control about this.

Therefore, ROMS arrays at RHO-points are always shifted by 1 in the I- and J-indices when processing NetCDF data. Notice that at U-points the one-shift is in the J-index whereas at V-points the one-shift is in the I-index. This is the reason why we always recommend the users to use the software that it is distributed from our repositories. This is always a problem in Matlab since it cannot handle zero or negative indices. We are not responsible neither can check the software developed for ROMS from third parties.

So it makes a lot of sense that if we have Ipos=10 and Jpos=20, the information stored in a RHO-point variable in the NetCDF file is located at array(11,21) and not at array(10,20).

The C-preprocessing option STATIONS_CGRID can be used to write into the station NetCDF files the data at its C-grid native locations instead of the RHO-points default locations. In this case, the interpretation is different for U- and V-point variables and vectors are staggered.

I hope that this clarifies well this issue. It is important for each user to understand clearly this aspect about ROMS and NetCDF files. It is very important when pre- or post-processing ROMS data. If you are developing software to process ROMS data you need to take this into consideration.