[MITgcm-support] Bug in external_fields_load when using periodicExternalForcing

[Jean-Michel Campin]

Hi Christopher,

Thanks for your detailed bug report.
Going to fix this in a way similar to what you 
proposed (but without the "save" instruction, because
of TAF + not safe with multi-threading).

Cheers,
Jean-Michel

On Tue, Apr 12, 2011 at 12:47:03PM -0700, Christopher L. Wolfe wrote:
> Hi all,
>
> I appear to have found a bug in the part of external_fields_load which  
> calculates interpolation weights and which records to load when using  
> periodicExternalForcing. The bug arises due to round-off errors incurred  
> when the variables are converted to integers using INT and NINT. It can  
> cause the model to crash when it goes looking for "impossible" record  
> numbers and also result in the interpolation weights being calculated  
> incorrectly. As an example of what can go wrong, consider
>
> deltaTclock = 86400 (one day)
> externForcingCycle = 31536000 (one 365-day year)
> externForcingPeriod = 2628000 (1/12 of a 365-day year)
>
> For mod(myIter,365) = 10, 11, ..., 14 this gives intime0 = 13, which is  
> "impossible" given that there should be only 12 records. This causes the  
> model to crash since Iftm(n)-Ifprd*(intime0(n)-1) .eq. 0 when  
> mod(myIter,365) = 10, thus causing the model to go looking for record  
> number 13. Simply duplicating record 1 as record 13 will prevent the  
> crash, but now the interpolation weights will be calculated incorrectly.  
> In particular, for mod(myIter,365) = 10, 11, ..., 14 aWght < .5 when it  
> should be greater than 0.5.
>
> I've attached a diff which appears to fix these problems for the case  
> above, but I haven't extensively tested it with different values of  
> deltaTclock, externForcingCycle, etc. Also, I've introduced a SAVE  
> statement to make it easier to figure out whether a new record should be  
> loaded. This works in my case, but it will no doubt break TAF. It  
> shouldn't be hard to rewrite the routine to avoid the use of the SAVE  
> statement.
>
> Cheers,
> Christopher
>
> -----------------------------------------------------------
> Dr. Christopher L. Wolfe                  858-534-4560
> Climate, Atmospheric Science, and Physical Oceanography
> Scripps Institution of Oceanography, UCSD  clwolfe at ucsd.edu
> -----------------------------------------------------------
>
>
>
>

> --- export/home/cwolfe/docs/research/globe/global1x1_tot/code/external_fields_load.F	2011-04-12 12:14:03.000000000 -0700
> +++ home/cwolfe/data/output/MITgcm/model/src/external_fields_load.F	2010-03-29 13:05:49.000000000 -0700
> @@ -57,32 +57,40 @@
>  C     !LOCAL VARIABLES:
>  C     === Local arrays ===
>  C     aWght, bWght :: Interpolation weights
> -      INTEGER bi,bj,i,j,intime0,intime1,intimem1
> -      _RL aWght,bWght,rdt,mytm,ftm
> +      INTEGER bi,bj,i,j,intime0,intime1
> +      _RL aWght,bWght,rdt
>        _RL tmp1Wght, tmp2Wght
>        INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm
> -      save intimem1
> -      data intimem1 /0/
>  CEOP
>  
>        IF ( periodicExternalForcing ) THEN
> +
> +C First call requires that we initialize everything to zero for safety
> +cph    has been shifted to ini_forcing.F
> +
>  C Now calculate whether it is time to update the forcing arrays
> -       nForcingPeriods = int(externForcingCycle/externForcingPeriod)
> -       mytm = myTime
> -     &        + externForcingCycle*(1-nint(myTime/externForcingCycle))
> -       ftm = mod(mytm+externForcingCycle-externForcingPeriod/2,
> -     &        externForcingCycle)
> +      rdt = 1. _d 0 / deltaTclock
> +      nForcingPeriods = NINT(externForcingCycle/externForcingPeriod)
> +      Imytm = NINT(myTime*rdt)
> +      Ifprd = NINT(externForcingPeriod*rdt)
> +      Ifcyc = NINT(externForcingCycle*rdt)
> +      Imytm = Imytm + Ifcyc*( 1 - NINT(myTime/externForcingCycle) )
> +      Iftm  = MOD( Imytm+Ifcyc-Ifprd/2, Ifcyc)
>       
> -       intime0 = 1 + int(ftm/externForcingPeriod)
> +      intime0 = 1 + INT(Iftm/Ifprd)
>         intime1 = 1 + mod(intime0,nForcingPeriods)
> -       
> -       tmp1Wght = ftm-externForcingPeriod*(intime0 - 1)
> -       aWght = tmp1Wght/externForcingPeriod
> +C-jmc: with some option of g77, FLOAT results in real*4 evaluation
> +C      of aWght; using DFLOAT always force real*8 computation:
> +c     aWght = DFLOAT( Iftm-Ifprd*(intime0 - 1) ) / DFLOAT( Ifprd )
> +C-ph: however, TAF doesnt recognize DFLOAT,
> +C-jmc: so let me try this:
> +      tmp1Wght = FLOAT( Iftm-Ifprd*(intime0 - 1) )
> +      tmp2Wght = FLOAT( Ifprd )
> +      aWght =  tmp1Wght / tmp2Wght
>         bWght = 1. _d 0 - aWght
>  
>        IF (
> -     &  intime0 .EQ. intimem1
> -     &  .OR. myTime .Eq. 0
> +     &  Iftm-Ifprd*(intime0-1) .EQ. 0
>       &  .OR. myIter .EQ. nIter0
>       & ) THEN
>  
> @@ -208,8 +216,6 @@
>  
>        ENDIF
>  
> -      intimem1 = intime1
> -      
>  C--   Interpolate fu,fv,Qnet,EmPmR,SST,SSS,Qsw
>        DO bj = myByLo(myThid), myByHi(myThid)
>         DO bi = myBxLo(myThid), myBxHi(myThid)

> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support