A couple of times we have entered station coordinates at known wet point on a grid and found the station file contained NaN's for the location and all variables for that station.
This would occur for a particular tiling but the problem would go away if we changed the tiling slightly. But sometimes the problem came back if we made the tiling an integer multiple of the original problem tiling.
This suggested some sort of bug so I investigated further and found that sometimes the station is being found in multiple tiles by the hindicies subroutine when the coordinates lie very close to a tile boundary
The problem with this is that the parallel routine mp_collect in grids_coords.f90 collects and sums all Ista and Jsta variables found by all tiles and creates a fractional index 2 to 4 times the correct value (depending on how many tiles it is found in).
Often this point lies outside the grid which is when we get nan's in the station file but it is possible for the point to lie inside the grid for low Ista and Jsta indicies and the resulting station file will actually contain the incorrect location and data for the intended station. Then the problem is much harder to pick up.
I think this error can be replicated by setting the station location to any lon lat combination corresponding to IstrR,Iend+1,JstrR,Jend+1 as they occur in grid_coords
e.g.
station longitude=GRID(ng)%lonr(IstrR,JstrR)
station latitude=GRID(ng)%latr(IstrR,JstrR)
look at the value of Ista and Jsta before and after the mp_collect call, after the collect call the value will be the value of Ista and Jsta summed over every cell the station was found in.
If the values of Ista or Jsta exceed the size of the grid the station file will have NaN's in it.
In case it is a compiler issue we are using gfortran 4.4.7 on Linux 6.6 with mpich2 libraries
Bug with grid coordinates for station file
-
CharlesJames
- Posts: 43
- Joined: Thu May 24, 2007 12:12 pm
- Location: South Australian Research and Development Institute
-
CharlesJames
- Posts: 43
- Joined: Thu May 24, 2007 12:12 pm
- Location: South Australian Research and Development Institute
Re: Bug with grid coordinates for station file
This seems to work
I set a station count to 1 if the tile determines that it contains the station, 0 otherwise. Then I mp_collect that variable and the result is the number of times the station was found in different tiles (usually 1 but can be up to 4) - then divide the mp_collected fractional indicies by the station count and the proper fractional index should be recovered
search 'CJ fix' below to see what I tried to do in grid_coords.F
I set a station count to 1 if the tile determines that it contains the station, 0 otherwise. Then I mp_collect that variable and the result is the number of times the station was found in different tiles (usually 1 but can be up to 4) - then divide the mp_collected fractional indicies by the station count and the proper fractional index should be recovered
search 'CJ fix' below to see what I tried to do in grid_coords.F
Code: Select all
#include "cppdefs.h"
#if defined FLOATS || defined STATIONS
SUBROUTINE grid_coords (ng, model)
!
!svn $Id: grid_coords.F 751 2015-01-07 22:56:36Z arango $
!================================================== Hernan G. Arango ===
! Copyright (c) 2002-2015 The ROMS/TOMS Group !
! Licensed under a MIT/X style license !
! See License_ROMS.txt !
!=======================================================================
! !
! This routine converts initial locations to fractional grid (I,J) !
! coordinates, if appropriate. !
! !
!=======================================================================
!
USE mod_param
USE mod_parallel
# ifdef FLOATS
USE mod_floats
# endif
USE mod_grid
USE mod_scalars
!
# ifdef DISTRIBUTE
USE distribute_mod, ONLY : mp_collect
# endif
USE interpolate_mod
!
implicit none
!
! Imported variable declarations.
!
integer, intent(in) :: ng, model
!
! Local variable declarations.
!
integer :: IstrR, Iend, JstrR, Jend
integer :: LBi, UBi, LBj, UBj
integer :: i, j, k, l, mc
real(r8), parameter :: spv = 0.0_r8
# ifdef FLOATS
real(r8) :: Xstr, Xend, Ystr, Yend, zfloat
logical, dimension(Nfloats(ng)) :: my_thread
real(r8), dimension(Nfloats(ng)) :: Iflt, Jflt
# ifdef SOLVE3D
real(r8), dimension(Nfloats(ng)) :: Kflt
# endif
# endif
# ifdef STATIONS
real(r8), dimension(Nstation(ng)) :: Slon, Slat
real(r8), dimension(Nstation(ng)) :: Ista, Jsta
# ifdef DISTRIBUTE
! CJ fix create new variables to count number of occurances of station in tiles
real(r8), dimension(Nstation(ng)) :: STNcount
# endif
# endif
!
!-----------------------------------------------------------------------
! Determine searching model grid box and arrays bounds.
!-----------------------------------------------------------------------
!
# ifdef DISTRIBUTE
IstrR=BOUNDS(ng)%IstrR(MyRank)
Iend =BOUNDS(ng)%Iend (MyRank)
JstrR=BOUNDS(ng)%JstrR(MyRank)
Jend =BOUNDS(ng)%Jend (MyRank)
# else
IstrR=0
Iend =Lm(ng)
JstrR=0
Jend =Mm(ng)
# endif
!
LBi=LBOUND(GRID(ng)%h,DIM=1)
UBi=UBOUND(GRID(ng)%h,DIM=1)
LBj=LBOUND(GRID(ng)%h,DIM=2)
UBj=UBOUND(GRID(ng)%h,DIM=2)
# ifdef FLOATS
!
Xstr=REAL(BOUNDS(ng)%Istr(MyRank),r8)-0.5_r8
Xend=REAL(BOUNDS(ng)%Iend(MyRank),r8)+0.5_r8
Ystr=REAL(BOUNDS(ng)%Jstr(MyRank),r8)-0.5_r8
Yend=REAL(BOUNDS(ng)%Jend(MyRank),r8)+0.5_r8
!
!-----------------------------------------------------------------------
! If applicable, convert initial floats locations (Flon,Flat) to
! fractional grid coordinates.
!-----------------------------------------------------------------------
!
IF (spherical) THEN
IF (Lfloats(ng)) THEN
mc=DRIFTER(ng)%Findex(0)
IF (DRIFTER(ng)%Findex(0).gt.0) THEN
CALL hindices (ng, LBi, UBi, LBj, UBj, &
& IstrR, Iend+1, JstrR, Jend+1, &
& GRID(ng)%angler, &
& GRID(ng)%lonr, &
& GRID(ng)%latr, &
& 1, mc, 1, 1, &
& 1, mc, 1, 1, &
& DRIFTER(ng)%Flon, &
& DRIFTER(ng)%Flat, &
& Iflt, Jflt, spv, .FALSE.)
# ifdef DISTRIBUTE
CALL mp_collect (ng, model, mc, spv, Iflt)
CALL mp_collect (ng, model, mc, spv, Jflt)
# endif
DO i=1,mc
l=DRIFTER(ng)%Findex(i)
DRIFTER(ng)%Tinfo(ixgrd,l)=Iflt(i)
DRIFTER(ng)%Tinfo(iygrd,l)=Jflt(i)
END DO
END IF
END IF
END IF
# ifdef SOLVE3D
!
! Determine which node bounds the initial float location.
!
# ifdef DISTRIBUTE
IF (Lfloats(ng)) THEN
DO l=1,Nfloats(ng)
IF ((Xstr.le.DRIFTER(ng)%Tinfo(ixgrd,l)).and. &
& (DRIFTER(ng)%Tinfo(ixgrd,l).lt.Xend).and. &
& (Ystr.le.DRIFTER(ng)%Tinfo(iygrd,l)).and. &
& (DRIFTER(ng)%Tinfo(iygrd,l).lt.Yend)) THEN
my_thread(l)=.TRUE.
ELSE
my_thread(l)=.FALSE.
END IF
END DO
END IF
# else
DO l=1,Nfloats(ng)
my_thread(l)=.TRUE.
END DO
# endif
# endif
!
!-----------------------------------------------------------------------
! Set float initial vertical level position, if inside application
! grid. If the initial float depth (in meters) is not found, release
! float at the surface model level.
!-----------------------------------------------------------------------
!
DO l=1,Nfloats(ng)
IF (Lfloats(ng)) THEN
# ifdef SOLVE3D
DRIFTER(ng)%Fz0(l)=spv
IF (my_thread(l).and. &
& ((DRIFTER(ng)%Tinfo(ixgrd,l).ge.0.5_r8).and. &
& (DRIFTER(ng)%Tinfo(iygrd,l).ge.0.5_r8).and. &
& (DRIFTER(ng)%Tinfo(ixgrd,l).le. &
& REAL(Lm(ng),r8)+0.5_r8).and. &
& (DRIFTER(ng)%Tinfo(iygrd,l).le. &
& REAL(Mm(ng),r8)+0.5_r8))) THEN
zfloat=DRIFTER(ng)%Tinfo(izgrd,l)
DRIFTER(ng)%Fz0(l)=zfloat ! Save original value
Kflt(l)=zfloat
IF (zfloat.le.0.0_r8) THEN
i=INT(DRIFTER(ng)%Tinfo(ixgrd,l)) ! Fractional positions
j=INT(DRIFTER(ng)%Tinfo(iygrd,l)) ! are still in this cell
IF (zfloat.lt.GRID(ng)%z_w(i,j,0)) THEN
zfloat=GRID(ng)%z_w(i,j,0)+5.0_r8
DRIFTER(ng)%Fz0(l)=zfloat
END IF
DRIFTER(ng)%Tinfo(izgrd,l)=REAL(N(ng),r8)
DO k=N(ng),1,-1
IF ((GRID(ng)%z_w(i,j,k)-zfloat)* &
& (zfloat-GRID(ng)%z_w(i,j,k-1)).ge.0.0_r8) THEN
Kflt(l)=REAL(k-1,r8)+ &
& (zfloat-GRID(ng)%z_w(i,j,k-1))/ &
& GRID(ng)%Hz(i,j,k)
END IF
END DO
END IF
ELSE
Kflt(l)=spv
END IF
# else
DRIFTER(ng)%Tinfo(izgrd,l)=0.0_r8
# endif
END IF
END DO
# ifdef SOLVE3D
IF (Lfloats(ng)) THEN
# ifdef DISTRIBUTE
CALL mp_collect (ng, model, Nfloats(ng), spv, DRIFTER(ng)%Fz0)
CALL mp_collect (ng, model, Nfloats(ng), spv, Kflt)
# endif
DO l=1,Nfloats(ng)
DRIFTER(ng)%Tinfo(izgrd,l)=Kflt(l)
END DO
END IF
# endif
# endif
# ifdef STATIONS
!
!-----------------------------------------------------------------------
! If applicable, convert station locations (SposX,SposY) to fractional
! grid coordinates.
!-----------------------------------------------------------------------
!
IF (spherical) THEN
mc=0
DO l=1,Nstation(ng)
IF (SCALARS(ng)%Sflag(l).gt.0) THEN
mc=mc+1
Slon(mc)=SCALARS(ng)%SposX(l)
Slat(mc)=SCALARS(ng)%SposY(l)
END IF
END DO
IF (mc.gt.0) THEN
CALL hindices (ng, LBi, UBi, LBj, UBj, &
& IstrR, Iend+1, JstrR, Jend+1, &
& GRID(ng)%angler, &
& GRID(ng)%lonr, &
& GRID(ng)%latr, &
& 1, mc, 1, 1, &
& 1, mc, 1, 1, &
& Slon, Slat, &
& Ista, Jsta, &
& spv, .FALSE.)
# ifdef DISTRIBUTE
! CJ fix check to see if a station is in the tile (assume station index can not be 0.d0)
do l=1,Nstation(ng)
if ((Ista(l).gt.0.D0).OR.(Jsta(l).gt.0.D0)) then
STNcount(l)=1.D0
else
STNcount(l)=0.D0
end if
end do
! CJ fix collect the counts across all processors
call mp_collect(ng,model,mc,spv,STNcount)
CALL mp_collect (ng, model, mc, spv, Ista)
CALL mp_collect (ng, model, mc, spv, Jsta)
# endif
mc=0
DO l=1,Nstation(ng)
IF (SCALARS(ng)%Sflag(l).gt.0) THEN
mc=mc+1
# ifdef DISTRIBUTE
! CJ fix divide Ista by number of counts (1-4) to recover original value if necessary
Ista(mc)=Ista(mc)/STNcount(mc)
Jsta(mc)=Jsta(mc)/STNcount(mc)
# endif
SCALARS(ng)%SposX(l)=Ista(mc)
SCALARS(ng)%SposY(l)=Jsta(mc)
END IF
END DO
END IF
END IF
# endif
RETURN
END SUBROUTINE grid_coords
#else
SUBROUTINE grid_coords
RETURN
END SUBROUTINE grid_coords
#endif-
arango
- Site Admin
- Posts: 1368
- Joined: Wed Feb 26, 2003 4:41 pm
- Location: DMCS, Rutgers University
- Contact:
Re: Bug with grid coordinates for station file
Thank you. I have to think about this. I need to reproduce your problem so I can look in the debugger.