Opened 8 years ago

Closed 8 years ago

#1040 closed Bug (fixed)

NaNs originating in sbcdcy (diurnal cycle) at some points in space and time

Reported by: acc Owned by: nemo
Priority: normal Milestone:
Component: OCE Version: release-3.4
Severity: Keywords:


Affects runs with ln_dm2dc = .true.

It seems the algorithm in sbcdcy can fail in some circumstances. In particular I have points at Arctic circle latitudes which return -Infinity values for qsr despite valid inputs at the winter solstice. At the point and time in question dusk and dawn are almost coincident which leads to infinite values of rscal. This program emulatess the problem using values written out from sbcdcy:

      program tests
      REAL(KIND=8) :: fintegral, pt1, pt2, paaa, pbbb, pccc, ztwopi, zinvtwopi, zconvrad
      REAL(KIND=8) :: raa, rbb, rcc, rdawn, rdusk
      REAL(KIND=8) :: gphit, zdecrad,zdsws
      INTEGER      :: nday
      fintegral( pt1, pt2, paaa, pbbb, pccc ) =                         &
         &   paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2)   &
         & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1)
       ztwopi    = 2.d0 * 4.d0*ATAN(1.d0)
       zinvtwopi = 1.d0 / ztwopi
       zconvrad  = ztwopi / 360.d0
! Point near the Arctic circle in December
! values written out from sbcdcy.F90 for this point
       gphit = 6.6500000000000000E+01
       zdecrad = -4.1015237421866746E-01
       zdsws = 3.6500000000000000E+02
       nday = 20
       rdawn = 4.8254282104903723E-01
       rdusk = 4.8254282579222396E-01
       raa = -3.6567685080958523E-01
       rbb = 3.6567685080958529E-01
       rcc = -3.0319059782014595E+00
       write(6,*) 'RSCAL is ' , 1.d0/fintegral(rdawn, rdusk, raa, rbb, rcc)
! Infinity!!
       end program tests

I can avoid the problem by checking for very short days (sbcdcy.F90) :

         !     2.2 Compute the scaling function:
         !         S* = the inverse of the time integral of the diurnal cycle from dawn to dusk
         DO jj = 1, jpj
            DO ji = 1, jpi
               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
                  rscal(ji,jj) = 0.0d0
                  IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part
                     IF( (rdusk(ji,jj) - rdawn(ji,jj) ) .ge. 0.001_wp ) THEN                                   ! bugfix
                       rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
                       rscal(ji,jj) = 1._wp / rscal(ji,jj)
                     ENDIF                                                                                     ! bugfix
                  ELSE                                         ! day time in two parts
                     IF( (rdusk(ji,jj) + (1._wp - rdawn(ji,jj)) ) .ge. 0.001_wp ) THEN                         ! bugfix
                       rscal(ji,jj) = fintegral(0._wp, rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   &
                          &         + fintegral(rdawn(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
                       rscal(ji,jj) = 1. / rscal(ji,jj)
                     ENDIF                                                                                     ! bugfix
                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
                     rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
                     rscal(ji,jj) = 1._wp / rscal(ji,jj)
                  ELSE                                          ! No day
                     rscal(ji,jj) = 0.0_wp
            END DO
         END DO

but this may not be the best solution?

Commit History (1)


ticket #1040. Fix to sbcdcy.F90 to avoid possible infinite scaling factor associated with very short daylight periods. Also upgraded all real constants to full precision.

Change History (1)

comment:1 Changed 8 years ago by acc

  • Resolution set to fixed
  • Status changed from new to closed

In trying to ascertain why this hasn't been a problem before I can only conclude that the behaviour is compiler dependent. In my case a switch from the PGI compiler to the CRAY compiler on a XT6 system is the cause. I can confirm different results with the different compilers:

crayftn -o tests tests.F90

RSCAL is Infinity

module swap Prg Env?-cray Prg Env?-pgi
pgf90 -o tests tests.F90

RSCAL is -3.1737137883082534E+017

so the PGI compiler didn't flag the result as an Inf. Although RSCAL is large I guess the duration and flux are both so small that the code copes. The Inf however eventually leads to a NaN which then propagates everywhere.

Tests with ifort mimic the CRAY compiler results. A solution has been submitted in Changeset: 3739

Note: See TracTickets for help on using tickets.