Changeset 3739

Timestamp:

2012-12-18T13:11:03+01:00 (11 years ago)

Author:

acc

Message:

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.

File:

• trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90 (modified) (7 diffs)

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90

@@ -102 +102 @@
          rcc(:,:) = zconvrad * glamt(:,:) - rpi
          ! time of midday
-         rtmd(:,:) = 0.5 - glamt(:,:) / 360.
-         rtmd(:,:) = MOD( (rtmd(:,:) + 1.), 1. )
+         rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp
+         rtmd(:,:) = MOD( (rtmd(:,:) + 1._wp) , 1._wp)
       ENDIF
 
@@ -118 +118 @@
          zdsws = REAL(11 + nday_year, wp)
          ! declination of the earths orbit
-         zdecrad = (-23.5 * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )
+         zdecrad = (-23.5_wp * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )
          ! Compute A and B needed to compute the time integral of the diurnal cycle
         
@@ -136 +136 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
+               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
          ! When is it night?
                   ztx = zinvtwopi * (ACOS(rab(ji,jj)) - rcc(ji,jj))
                   ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx )
          ! is it dawn or dusk?
-                  IF ( ztest > 0 ) THEN
+                  IF ( ztest > 0._wp ) THEN
                      rdawn(ji,jj) = ztx
                      rdusk(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) )
@@ -149 +149 @@
                   ENDIF
                ELSE
-                  rdawn(ji,jj) = rtmd(ji,jj) + 0.5
+                  rdawn(ji,jj) = rtmd(ji,jj) + 0.5_wp
                   rdusk(ji,jj) = rdawn(ji,jj)
                ENDIF
@@ -157 +157 @@
          rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp )
 
-         !     2.2 Compute the scalling function:
-         !         S* = the inverse of the time integral of the diurnal cycle from dawm to dusk
+         !     2.2 Compute the scaling function:
+         !         S* = the inverse of the time integral of the diurnal cycle from dawn to dusk
+         !         Avoid possible infinite scaling factor, associated with very short daylight
+         !         periods, by ignoring periods less than 1/1000th of a day (ticket #1040)
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
+               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
+                  rscal(ji,jj) = 0.0_wp
                   IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part
-                     rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
-                     rscal(ji,jj) = 1. / rscal(ji,jj)
+                     IF( (rdusk(ji,jj) - rdawn(ji,jj) ) .ge. 0.001_wp ) THEN
+                       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
                   ELSE                                         ! day time in two parts
-                     rscal(ji,jj) = fintegral(0., rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   &
-                        &         + fintegral(rdawn(ji,jj), 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
-                     rscal(ji,jj) = 1. / rscal(ji,jj)
+                     IF( (rdusk(ji,jj) + (1._wp - rdawn(ji,jj)) ) .ge. 0.001_wp ) THEN
+                       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
                   ENDIF
                ELSE
                   IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
-                     rscal(ji,jj) = fintegral(0., 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
-                     rscal(ji,jj) = 1. / rscal(ji,jj)
+                     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.e0
+                     rscal(ji,jj) = 0.0_wp
                   ENDIF
                ENDIF
@@ -191 +198 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
+            IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
                !
                IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN       ! day time in one part
@@ -218 +225 @@
                   !
                ELSE                                         ! No day
-                  zqsrout(ji,jj) = 0.e0
+                  zqsrout(ji,jj) = 0.0_wp
                ENDIF
             ENDIF