/[lmdze]/trunk/phylmd/zenang.f90

Diff of /trunk/phylmd/zenang.f90

Parent Directory | Revision Log | View Patch Patch

-revision 210 by guez,
Tue Dec 13 16:02:23 2016 UTC
+revision 213 by guez,
Mon Feb 27 15:44:55 2017 UTC
 Line 4 
 module zenang_m
  contains
-   SUBROUTINE zenang(longi, gmtime, pdtrad, mu0, frac)
+   SUBROUTINE zenang(longi, gmtime, pdtrad, mu0, fract)
      ! Author: O. Boucher (LMD/CNRS), d'après les routines "zenith" et
      ! "angle" de Z. X. Li
 Line 12 
 contains
      ! Date : première version le 13 décembre 1994, revu pour GCM le 30
      ! septembre 1996
-     ! Calcule les valeurs moyennes du cos de l'angle zénithal et
+     ! Calcule les valeurs moyennes du cosinus de l'angle zénithal et
      ! l'ensoleillement moyen entre "gmtime" et "gmtime + pdtrad"
      ! connaissant la déclinaison, la latitude et la longitude.
-     ! Différent de la routine "angle" en ce sens que "zenang" fournit
+     ! Différent de la routine "angle" parce que "zenang" fournit des
-     ! des moyennes de "mu0" et non des valeurs instantanées. Du coup
+     ! moyennes de "mu0" et non des valeurs instantanées. Du coup
-     ! "frac" prend toutes les valeurs entre 0 et 1. Cf. Capderou (2003
+     ! "fract" prend toutes les valeurs entre 0 et 1. Cf. Capderou (2003
      ! 784, equation 9.11).
      USE dimphy, ONLY: klon
 Line 33 
 contains
      REAL, INTENT(IN):: pdtrad ! pas de temps du rayonnement (s)
      REAL, INTENT(OUT):: mu0(:) ! (klon)
-     ! cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad"
+     ! cosine of mean zenith angle between "gmtime" and "gmtime + pdtrad"
-     REAL, INTENT(OUT), OPTIONAL:: frac(:) ! (klon)
+     REAL, INTENT(OUT), OPTIONAL:: fract(:) ! (klon)
-     ! ensoleillement moyen entre gmtime et gmtime+pdtrad
+     ! ensoleillement moyen entre gmtime et gmtime + pdtrad
      ! Local:
 Line 45 
 contains
      REAL omega1, omega2 ! temps 1 et 2 exprimés en radians avec 0 à midi
      REAL omega ! heure en rad du coucher de soleil
-     ! - omega est donc l'heure en rad de lever du soleil
+     ! "- omega" est donc l'heure en rad de lever du soleil.
      REAL omegadeb, omegafin
      REAL zfrac1, zfrac2, z1_mu, z2_mu
 Line 54 
 contains
      !----------------------------------------------------------------------
-     if (present(frac)) call assert((/size(mu0), size(frac)/) == klon, "zenang")
+     if (present(fract)) call assert((/size(mu0), size(fract)/) == klon, &
+          "zenang")
      lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
      ! Capderou (2003 784, equation 4.49)
      gmtime1 = gmtime * 86400.
-     gmtime2 = gmtime * 86400. + pdtrad
+     gmtime2 = gmtime1 + pdtrad
      DO i = 1, klon
         latr = rlat(i) * pi / 180.
-        omega = 0. ! nuit polaire
-        IF (latr>=(pi / 2.-lat_sun) .OR. latr<=(-pi / 2.-lat_sun)) THEN
+        IF (latr >= pi / 2. - lat_sun .OR. latr <= - pi / 2. - lat_sun) then
            omega = pi ! journée polaire
-        END IF
+        else IF (latr < pi / 2. + lat_sun .AND. latr > - pi / 2. + lat_sun) THEN
-        IF (latr<(pi / 2.+lat_sun) .AND. latr>(-pi / 2.+lat_sun) .AND. &
+           omega = acos(- tan(latr) * tan(lat_sun))
-             latr<(pi / 2.-lat_sun) .AND. latr>(-pi / 2.-lat_sun)) THEN
+        else
-           omega = -tan(latr) * tan(lat_sun)
+           omega = 0. ! nuit polaire
-           omega = acos(omega)
         END IF
-        omega1 = gmtime1 + rlon(i) * 86400. / 360.
+        omega1 = mod((gmtime1 + rlon(i) * 86400. / 360.) / 86400. * twopi &
-        omega1 = omega1 / 86400. * twopi
+             + twopi, twopi) - pi
-        omega1 = mod(omega1+twopi, twopi)
+        omega2 = mod((gmtime2 + rlon(i) * 86400. / 360.) / 86400. * twopi &
-        omega1 = omega1 - pi
+             + twopi, twopi) - pi
-        omega2 = gmtime2 + rlon(i) * 86400. / 360.
-        omega2 = omega2 / 86400. * twopi
-        omega2 = mod(omega2+twopi, twopi)
-        omega2 = omega2 - pi
-        IF (omega1<=omega2) THEN
+        IF (omega1 <= omega2) THEN
            ! on est dans la meme journee locale
-           IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN
+           IF (omega2 <= - omega .OR. omega1 >= omega .OR. omega < 1E-5) THEN
               ! nuit
-              IF (present(frac)) frac(i) = 0.
+              IF (present(fract)) fract(i) = 0.
               mu0(i) = 0.
            ELSE
               ! jour + nuit / jour
-              omegadeb = max(-omega, omega1)
+              omegadeb = max(- omega, omega1)
               omegafin = min(omega, omega2)
-              IF (present(frac)) frac(i) = (omegafin-omegadeb) / (omega2-omega1)
+              IF (present(fract)) fract(i) = (omegafin - omegadeb) &
+                   / (omega2 - omega1)
               mu0(i) = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
                    * (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
            END IF
         ELSE
            ! omega1 > omega2, à cheval sur deux journées
            ! entre omega1 et pi
-           IF (omega1>=omega) THEN ! nuit
+           IF (omega1 >= omega) THEN ! nuit
               zfrac1 = 0.
               z1_mu = 0.
-           ELSE ! jour+nuit
+           ELSE ! jour + nuit
-              omegadeb = max(-omega, omega1)
+              omegadeb = max(- omega, omega1)
               omegafin = omega
               zfrac1 = omegafin - omegadeb
               z1_mu = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
                    * (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
            END IF
-           ! entre -pi et omega2
+           ! entre - pi et omega2
-           IF (omega2<=-omega) THEN ! nuit
+           IF (omega2 <= - omega) THEN ! nuit
               zfrac2 = 0.
               z2_mu = 0.
-           ELSE ! jour+nuit
+           ELSE ! jour + nuit
-              omegadeb = -omega
+              omegadeb = - omega
               omegafin = min(omega, omega2)
               zfrac2 = omegafin - omegadeb
               z2_mu = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
                    * (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
            END IF
            ! moyenne
-           IF (present(frac)) frac(i) = (zfrac1+zfrac2) / (omega2+twopi-omega1)
+           IF (present(fract)) fract(i) = (zfrac1 + zfrac2) &
-           mu0(i) = (zfrac1 * z1_mu+zfrac2 * z2_mu) / max(zfrac1+zfrac2, 1e-10)
+                / (omega2 + twopi - omega1)
+           mu0(i) = (zfrac1 * z1_mu + zfrac2 * z2_mu) &
+                / max(zfrac1 + zfrac2, 1e-10)
         END IF
      END DO

 Legend:



Removed from v.210
 


changed lines


 
Added in v.213
 Legend:



Removed from v.210
 


changed lines


 
Added in v.213
-Removed from v.210
+Added in v.213

	ViewVC Help
Powered by ViewVC 1.1.21