# Changeset 3739

Ignore:
Timestamp:
2012-12-18T13:11:03+01:00 (9 years ago)
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

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

 r3294 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 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 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) ) 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 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 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 ! ELSE                                         ! No day zqsrout(ji,jj) = 0.e0 zqsrout(ji,jj) = 0.0_wp ENDIF ENDIF
Note: See TracChangeset for help on using the changeset viewer.