Ignore:
Timestamp:
2012-12-18T13:11:03+01:00 (8 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:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90

    r3294 r3739  
    102102         rcc(:,:) = zconvrad * glamt(:,:) - rpi 
    103103         ! time of midday 
    104          rtmd(:,:) = 0.5 - glamt(:,:) / 360. 
    105          rtmd(:,:) = MOD( (rtmd(:,:) + 1.), 1. ) 
     104         rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp 
     105         rtmd(:,:) = MOD( (rtmd(:,:) + 1._wp) , 1._wp) 
    106106      ENDIF 
    107107 
     
    118118         zdsws = REAL(11 + nday_year, wp) 
    119119         ! declination of the earths orbit 
    120          zdecrad = (-23.5 * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) ) 
     120         zdecrad = (-23.5_wp * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) ) 
    121121         ! Compute A and B needed to compute the time integral of the diurnal cycle 
    122122         
     
    136136         DO jj = 1, jpj 
    137137            DO ji = 1, jpi 
    138                IF ( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h 
     138               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    139139         ! When is it night? 
    140140                  ztx = zinvtwopi * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 
    141141                  ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx ) 
    142142         ! is it dawn or dusk? 
    143                   IF ( ztest > 0 ) THEN 
     143                  IF ( ztest > 0._wp ) THEN 
    144144                     rdawn(ji,jj) = ztx 
    145145                     rdusk(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) ) 
     
    149149                  ENDIF 
    150150               ELSE 
    151                   rdawn(ji,jj) = rtmd(ji,jj) + 0.5 
     151                  rdawn(ji,jj) = rtmd(ji,jj) + 0.5_wp 
    152152                  rdusk(ji,jj) = rdawn(ji,jj) 
    153153               ENDIF 
     
    157157         rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp ) 
    158158 
    159          !     2.2 Compute the scalling function: 
    160          !         S* = the inverse of the time integral of the diurnal cycle from dawm to dusk 
     159         !     2.2 Compute the scaling function: 
     160         !         S* = the inverse of the time integral of the diurnal cycle from dawn to dusk 
     161         !         Avoid possible infinite scaling factor, associated with very short daylight 
     162         !         periods, by ignoring periods less than 1/1000th of a day (ticket #1040) 
    161163         DO jj = 1, jpj 
    162164            DO ji = 1, jpi 
    163                IF ( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h 
     165               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
     166                  rscal(ji,jj) = 0.0_wp 
    164167                  IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part 
    165                      rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    166                      rscal(ji,jj) = 1. / rscal(ji,jj) 
     168                     IF( (rdusk(ji,jj) - rdawn(ji,jj) ) .ge. 0.001_wp ) THEN 
     169                       rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
     170                       rscal(ji,jj) = 1._wp / rscal(ji,jj) 
     171                     ENDIF 
    167172                  ELSE                                         ! day time in two parts 
    168                      rscal(ji,jj) = fintegral(0., rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   & 
    169                         &         + fintegral(rdawn(ji,jj), 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    170                      rscal(ji,jj) = 1. / rscal(ji,jj) 
     173                     IF( (rdusk(ji,jj) + (1._wp - rdawn(ji,jj)) ) .ge. 0.001_wp ) THEN 
     174                       rscal(ji,jj) = fintegral(0._wp, rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   & 
     175                          &         + fintegral(rdawn(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
     176                       rscal(ji,jj) = 1. / rscal(ji,jj) 
     177                     ENDIF 
    171178                  ENDIF 
    172179               ELSE 
    173180                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day 
    174                      rscal(ji,jj) = fintegral(0., 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    175                      rscal(ji,jj) = 1. / rscal(ji,jj) 
     181                     rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
     182                     rscal(ji,jj) = 1._wp / rscal(ji,jj) 
    176183                  ELSE                                          ! No day 
    177                      rscal(ji,jj) = 0.e0 
     184                     rscal(ji,jj) = 0.0_wp 
    178185                  ENDIF 
    179186               ENDIF 
     
    191198      DO jj = 1, jpj 
    192199         DO ji = 1, jpi 
    193             IF( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h 
     200            IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    194201               ! 
    195202               IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN       ! day time in one part 
     
    218225                  ! 
    219226               ELSE                                         ! No day 
    220                   zqsrout(ji,jj) = 0.e0 
     227                  zqsrout(ji,jj) = 0.0_wp 
    221228               ENDIF 
    222229            ENDIF 
Note: See TracChangeset for help on using the changeset viewer.