Bug with grid coordinates for station file

Bug reports, work arounds and fixes

Moderators: arango, robertson

Post Reply
Message
Author
User avatar
CharlesJames
Posts: 43
Joined: Thu May 24, 2007 12:12 pm
Location: South Australian Research and Development Institute

Bug with grid coordinates for station file

#1 Unread post by CharlesJames »

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

User avatar
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

#2 Unread post by CharlesJames »

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

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

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

Re: Bug with grid coordinates for station file

#3 Unread post by arango »

Thank you. I have to think about this. I need to reproduce your problem so I can look in the debugger.

Post Reply