Reading forcing data with DT < 1 second

[jcwarner]

I am using a NetCDF forcing data file with a baroclinic DT of 0.2 sec to drive a lab test case. But ROMS does not interpolate the data correctly because in set_ngfld (same for set_2dlfd and set_3dlfd) we have:

     fac1=ANINT(Tintrp(it2,ifield,ng)-time(ng),r8)
     fac2=ANINT(time(ng)-Tintrp(it1,ifield,ng),r8)

which truncates the time interpolation weights to a whole number for the nearest second towards zero.

We got it to work for smaller than a second baroclinic timestep by using:

      fac1=ANINT((Tintrp(it2,ifield,ng)-time(ng))*SecScale,r8)
      fac2=ANINT((time(ng)-Tintrp(it1,ifield,ng))*SecScale,r8)

where SecScale=1000. That is, the time interpolation weights are rounded to the nearest millisecond instead.

The following statements at the full precision did not work:

     fac1=Tintrp(it2,ifield,ng)-time(ng)
     fac2=time(ng)-Tintrp(it1,ifield,ng)

because there can be a small value of fac1 that is negative because of roundoff, and then the interpolation is stopped by

      ELSE IF (((fac1*fac2).ge.0.0_r8).and.(fac1+fac2).gt.0.0_r8) THEN
...

indicating unbounded interpolants.

---

WARNING:

Notice that we no longer will get identical solutions with previous versions due to very small differences in the time interpolated fields. It does not matter much because the differences are in the order of roundoff. However, users need to be aware of such fact.

[colucix]

I have had the same problem using a DT of 0.1 sec. The fix I adopted was

fac1=ANINT(Tintrp(it2,ifield,ng)*10.0_r8-time(ng)*10.0_r8,r8)
fac2=ANINT(time(ng)*10.0_r8-Tintrp(it1,ifield,ng)*10.0_r8,r8)

The division of the "ANINT number" by 10.0_r8 (or 100.0_r8 in your case) is not necessary since the factor are scaled correctly by

fac=1.0_r8/(fac1+fac2)
fac1=fac*fac1
fac2=fac*fac2

where all the terms of the fractions have been previously multiplied by a factor of 10 (or 100). Perhaps it is precisely the "unnecessary division" that fails in some cases.

Regards, alex

[arango]

Agree about division precision. I added the SecScale variable and set it to 1000, so the rounding is in the order of milliseconds. Notice that we need to divide by this scale when reporting interpolation fails because of unbounded time snapshots.