Tidal constituents bug

[m.hadfield]

I've detected a bug in the loading of tidal constituent data in get_data.F & get_2dfld.F. I have a work-around for my application, but not sure of the best general fix.

I think it relates to the way get_2dfld was originally used for arrays in which the outer dimension is time, but has now been adapted for arrays in which the outer dimension is tidal constiuent.

In get_data.F, the code to read in tidal constituent data looks like this

Code: Select all

#ifdef SSH_TIDES
!
!  Tidal elevation amplitude and phase. In order to read data as a
!  function of tidal period, we need to reset the model time variables
!  temporarily.
!
      IF (iic(ng).eq.0) THEN
        time_save=time(ng)
        time(ng)=8640000.0_r8
        tdays(ng)=time(ng)*sec2day
        CALL get_2dfld (ng, idTzam, ncfrcid(idTzam,ng),                 &
     &                  nFfiles(ng), frcname(1,ng),                     &
     &                  update(1), LBi, UBi, LBj, UBj, MTC, NTC(ng),    &
# ifdef MASKING
     &                  GRID(ng) % rmask(LBi,LBj),                      &
# endif
     &                  TIDES(ng) % SSH_Tamp(LBi,LBj,1))
        CALL get_2dfld (ng, idTzph, ncfrcid(idTzph,ng),                 &
     &                  nFfiles(ng), frcname(1,ng),                     &
     &                  update(1), LBi, UBi, LBj, UBj, MTC, NTC(ng),    &
# ifdef MASKING
     &                  GRID(ng) % rmask(LBi,LBj),                      &
# endif
     &                  TIDES(ng) % SSH_Tphase(LBi,LBj,1))
        time(ng)=time_save
        tdays(ng)=time(ng)*sec2day
      END IF
#endif

Note that 11th argument passed to get_2dfld is MTC, which is set to 7 in mod_param.F.

Inside get_2dfld, the 11th argument is Iout, described as "Size of the outer dimension, if any. Otherwise, Iout must be set to one by the calling program." In loading tidal constituents, during execution of get_2dfld, Tindex is set to Iout and at line 262 the following statement is executed

Code: Select all

Vtime(Tindex,ifield,ng)=Finfo(1,ifield,ng)

When bounds checking is enabled, this statement triggers an out-of-bounds warning at this line (number 241 after preprocessing):

Code: Select all

Subscript -4 is out of range for dimension 1 for array
'VTIME' at line 300 in file 'get_2dfld.f90' with bounds 1:7

Vtime is declared in mod_ncparam.F as Vtime(2,NV,Ngrids) and described as "Latest two-time values of processed input data".

For my present application, with only one tidal constituent, my work-around is to reduce MTC to 1, but what if I want more than 2 tidal constituents?

I'm not sure I understand the logic in get_data & get_2dfld, but I guess the inner dimension of Vtime needs to be increased to max(2,MTC).

[arango]

Yes, the first dimension of Vtime must be MAX(2,MTC) when the tides are activated. In fact, if you check Version 2.1, this bug was corrected by making Vtime an allocatable array. In mod_ncparam.F you can find the following statement:

Code: Select all

#if defined SSH_TIDES || defined UV_TIDES
        allocate ( Vtime(MAX(2,MTC),NV,Ngrids) )
#else
        allocate ( Vtime(2,NV,Ngrids) )
#endif

I think that you were using Version 2.0. This bug was corrected in Version 2.1. The array Vtime is ignored in the tidal forcing. It is only used for interpolation between snapshots. This is irrelevant in tidal forcing because periods (spectral time) is used instead. This in get_2dfld for other generic use. I don't know what this bug does to your tidal constituents. I assume that it will still load their values correctly. Of course, the behavior may change on different computers.

[m.hadfield]

Oops My bug report was somewhat garbled because, having found a bug in ROMS 2.1, I saw that it resembled one in an earlier version and confused the two.

You are correct, Hernan, that the dimension of Vtime is set to MAX(2,MTC) when tides are activated.

So with tides activated and with MTC set to the default value of 7, the first dimension of Vtime has bounds 1:7. The bounds-checking error message I get is

Code: Select all

  Subscript -4 is out of range for dimension 1 for array
  'VTIME' at line 300 in file 'get_2dfld.f90' with bounds 1:7.

where line 300 in the preprocessed file is

Code: Select all

          Vtime(Tindex,ifield,ng)=Tval*Tscale

The offending index, Tindex, was set equal to Iout earlier in get_2dfld. Iout is one of the arguments of this subroutine and was set equal to MTC by the caller. Then in the following code, Tindex is recalculated as -4, which is out of bounds

Code: Select all

!
!  Set rolling index for two-time record storage of input data.  If
!  "Iout" is unity, input data is stored in timeless array by the
!  calling program.  If Irec > 1, this routine is used to read a 2D
!  field varying in another non-time dimension.
!
          IF (Irec.eq.1) THEN
            IF (Iout.eq.1) THEN
              Tindex=1
            ELSE
              Tindex=3-Tindex
            END IF
            Iinfo(5,ifield,ng)=Tindex
          END IF

[wmartin]

I've been working on what I think is a related problem I found that with more than one tidal component only the first one gets scaled, so for example, the angle of the ellipse, UV_Tangle, for the second component gets used in ROMS as degrees instead of radians. This is because the variable Npts in nf_fread is only equal to Ilen*Jlen unless you set INLINE_2DIO and DISTRIBUTED in ccpdefs (which I don't). If those options are set then there is additional code for Npts = Npts * Klen and all the data gets processed. I tried dropping back to just one constituent and got hit by the problem you reported. Now I'm running with two constituents, but the second is a dummy. Is it safe to just use that extra line of code even thought it's not on a distributed system?
thanks
Wayne

[arango]

Well, there is a CPP construct bug in routine nf_fread.F. Please change lines 148 and 169 of nf_fread.F from

# if !defined INLINE_2DIO && defined DISTRIBUTE

to

# if !(defined INLINE_2DIO && defined DISTRIBUTE)

Notice that now the not defined CPP directive is in parenthesis. It is very easy to confuse logical and/or operations with CPP negation. I hope that this solves your problem.

I did noticed this one before but I did not realized that it affects reading of tidal data. I think that John Warner noticed this problem before too. This bug has been already corrected in the new version of ROMS/TOMS which will be released sometime in June.

Thank you for reporting this bug