Possible bug with stations in a nested grid configuration

[kakearney]

When I attempt to turn on stations in a nested grid configuration, I am encountering the following error:

Code: Select all

forrtl: severe (408): fort: (2): Subscript #1 of the array APOS has value 2 which is greater than the upper bound of 1

Image              PC                Routine            Line        Source
oceanG             00000000027E7E2E  Unknown               Unknown  Unknown
oceanG             00000000027E68C6  Unknown               Unknown  Unknown
oceanG             00000000027A4F92  Unknown               Unknown  Unknown
oceanG             0000000002764D1B  Unknown               Unknown  Unknown
oceanG             0000000002765231  Unknown               Unknown  Unknown
oceanG             000000000236AAA9  extract_sta_mod_m         145  extract_sta.f90
oceanG             00000000022F922E  wrt_info_                 416  wrt_info.f90
oceanG             0000000001EA864B  def_station_              663  def_station.f90
oceanG             0000000000A7F6E4  output_                   300  output.f90
oceanG             00000000007DD3CE  main3d_                   273  main3d.f90
oceanG             000000000041FA03  ocean_control_mod         164  ocean_control.f90
oceanG             000000000041E489  MAIN__                    108  master.f90
oceanG             000000000041DDEC  Unknown               Unknown  Unknown
libc.so.6          0000003B1C01ECDD  Unknown               Unknown  Unknown
oceanG             000000000041DCE9  Unknown               Unknown  Unknown

The relevant part of code in extract_sta.90 is:

Code: Select all

DO np=1,Npos
  Xgrd=Xpos(np)
  Ygrd=Ypos(np)
  bounded(np)=0.0_r8
  IF (((Xmin.le.Xgrd).and.(Xgrd.lt.Xmax)).and.                  &
&        ((Ymin.le.Ygrd).and.(Ygrd.lt.Ymax))) THEN
    i1=INT(Xgrd)
    j1=INT(Ygrd)
    i2=i1+1
    j2=j1+1
    IF (i2.gt.Lm(ng)+1) THEN
      i2=i1                   ! station at the eastern boundary
    END IF
    IF (j2.gt.Mm(ng)+1) THEN
      j2=j1                   ! station at the northern boundary
    END IF
    bounded(np)=1.0_r8
    p2=REAL(i2-i1,r8)*(Xgrd-REAL(i1,r8))
    q2=REAL(j2-j1,r8)*(Ygrd-REAL(j1,r8))
    p1=1.0_r8-p2
    q1=1.0_r8-q2
    w111=p1*q1
    w211=p2*q1
    w121=p1*q2
    w221=p2*q2
    w111=w111*GRID(ng)%rmask(i1,j1)
    w211=w211*GRID(ng)%rmask(i2,j1)
    w121=w121*GRID(ng)%rmask(i1,j2)
    w221=w221*GRID(ng)%rmask(i2,j2)
    wsum=w111+w211+w121+w221
    IF (wsum.gt.0.0_r8) THEN
      wsum=1.0_r8/wsum
      w111=w111*wsum
      w211=w211*wsum
      w121=w121*wsum
      w221=w221*wsum
    ELSE
      bounded(np)=0.0_r8
    ENDIF
    Apos(np)=Ascl*(w111*A(i1,j1)+                               &
&                     w211*A(i2,j1)+                               &
&                     w121*A(i1,j2)+                               &
&                     w221*A(i2,j2))
    IF (ABS(Apos(ng)).eq.0.0_r8) Apos(ng)=0.0_r8      ! <-- line 145
  ELSE
    Apos(np)=Aspv
  END IF
END DO

If I'm understanding this code correctly (and I may not be... I'm pretty new to ROMS), Apos is allocated with as many elements as stations in the grid, but then is being indexed by grid number. In my case, I was attempting 5 stations in the parent grid (grid 1) and only 1 in the refinement grid (grid 2), so the latter caused a dimension clash. If I add a second station to grid 2, the error is avoided. I'm not quite sure what the purpose of line 145 is... should it include Apos(np) instead of Apos(ng)?

[arango]

Indeed, this should be:

Code: Select all

            IF (ABS(Apos(np)).eq.0.0_r8) Apos(np)=0.0_r8      ! positive
                                                              ! zero

This is done to impose positive zero. It is harmless except when you have only one station and ng > 1, so you get the out-of-bounds error. It occurs in several places in extract_sta.F.

Starting F95 zero values can be signed (-0 or +0) following the IEEE 754 floating point standard. This can be advantageous in some computations but not in extract_sta.F when Ascl is negative and Apos is zero. This will produce different output files in serial and distributed memory applications. Since comparing serial and parallel output is essential for tracking parallel partition bugs, positive zero is enforced.

Thank you for bringing this to my attention.