Sources/phylmd/zenang.f

module zenang_m

  IMPLICIT NONE

contains

  SUBROUTINE zenang(longi, gmtime, pdtrad, mu0, frac)

    ! Author: O. Boucher (LMD/CNRS), d'après les routines "zenith" et
    ! "angle" de Z.X. Li

    ! 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
    ! 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
    ! des moyennes de "mu0" et non des valeurs instantanées. Du coup
    ! "frac" prend toutes les valeurs entre 0 et 1. Cf. Capderou (2003
    ! 784, equation 9.11).

    USE dimphy, ONLY: klon
    USE yomcst, ONLY: r_incl
    USE phyetat0_m, ONLY: rlat, rlon
    use nr_util, only: assert, pi, twopi

    REAL, INTENT(IN):: longi
    ! longitude vraie de la terre dans son plan solaire à partir de
    ! l'équinoxe de printemps (in degrees)

    REAL, INTENT(IN):: gmtime ! temps universel en fraction de jour
    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"

    REAL, INTENT(OUT), OPTIONAL:: frac(:) ! (klon)
    ! ensoleillement moyen entre gmtime et gmtime+pdtrad

    ! Local:

    INTEGER i
    REAL gmtime1, gmtime2
    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

    REAL omegadeb, omegafin
    REAL zfrac1, zfrac2, z1_mu, z2_mu
    REAL lat_sun ! déclinaison en radians
    REAL latr ! latitude du point de grille en radians

    !----------------------------------------------------------------------

    if (present(frac)) call assert((/size(mu0), size(frac)/) == 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

    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
          omega = pi ! journée polaire
       END IF
       IF (latr<(pi / 2.+lat_sun) .AND. latr>(-pi / 2.+lat_sun) .AND. &
            latr<(pi / 2.-lat_sun) .AND. latr>(-pi / 2.-lat_sun)) THEN
          omega = -tan(latr) * tan(lat_sun)
          omega = acos(omega)
       END IF

       omega1 = gmtime1 + rlon(i) * 86400. / 360.
       omega1 = omega1 / 86400. * twopi
       omega1 = mod(omega1+twopi, twopi)
       omega1 = omega1 - pi

       omega2 = gmtime2 + rlon(i) * 86400. / 360.
       omega2 = omega2 / 86400. * twopi
       omega2 = mod(omega2+twopi, twopi)
       omega2 = omega2 - pi

       IF (omega1<=omega2) THEN
          ! on est dans la meme journee locale
          IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN
             ! nuit
             IF (present(frac)) frac(i) = 0.
             mu0(i) = 0.
          ELSE
             ! jour + nuit / jour
             omegadeb = max(-omega, omega1)
             omegafin = min(omega, omega2)
             IF (present(frac)) frac(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
             zfrac1 = 0.
             z1_mu = 0.
          ELSE ! jour+nuit
             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
          IF (omega2<=-omega) THEN ! nuit
             zfrac2 = 0.
             z2_mu = 0.
          ELSE ! jour+nuit
             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)
          mu0(i) = (zfrac1 * z1_mu+zfrac2 * z2_mu) / max(zfrac1+zfrac2, 1e-10)
       END IF
    END DO

  END SUBROUTINE zenang

end module zenang_m
1	guez	98	module zenang_m
2
3			IMPLICIT NONE
4
5			contains
6
7	guez	118	SUBROUTINE zenang(longi, gmtime, pdtrad, mu0, frac)
8	guez	98
9	guez	118	! Author: O. Boucher (LMD/CNRS), d'après les routines "zenith" et
10	guez	98	! "angle" de Z.X. Li
11
12	guez	118	! Date : première version le 13 décembre 1994, revu pour GCM le 30
13			! septembre 1996
14
15	guez	98	! Calcule les valeurs moyennes du cos de l'angle zénithal et
16	guez	125	! l'ensoleillement moyen entre "gmtime" et "gmtime + pdtrad"
17	guez	118	! connaissant la déclinaison, la latitude et la longitude.
18			! Différent de la routine "angle" en ce sens que "zenang" fournit
19	guez	125	! des moyennes de "mu0" et non des valeurs instantanées. Du coup
20			! "frac" prend toutes les valeurs entre 0 et 1. Cf. Capderou (2003
21			! 784, equation 9.11).
22	guez	98
23	guez	125	USE dimphy, ONLY: klon
24			USE yomcst, ONLY: r_incl
25			USE phyetat0_m, ONLY: rlat, rlon
26			use nr_util, only: assert, pi, twopi
27
28	guez	118	REAL, INTENT(IN):: longi
29			! longitude vraie de la terre dans son plan solaire à partir de
30			! l'équinoxe de printemps (in degrees)
31	guez	98
32	guez	118	REAL, INTENT(IN):: gmtime ! temps universel en fraction de jour
33			REAL, INTENT(IN):: pdtrad ! pas de temps du rayonnement (s)
34	guez	98
35	guez	118	REAL, INTENT(OUT):: mu0(:) ! (klon)
36			! cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad"
37	guez	98
38	guez	118	REAL, INTENT(OUT), OPTIONAL:: frac(:) ! (klon)
39			! ensoleillement moyen entre gmtime et gmtime+pdtrad
40	guez	98
41	guez	118	! Local:
42	guez	98
43			INTEGER i
44			REAL gmtime1, gmtime2
45	guez	118	REAL omega1, omega2 ! temps 1 et 2 exprimés en radians avec 0 à midi
46	guez	98
47	guez	118	REAL omega ! heure en rad du coucher de soleil
48			! - omega est donc l'heure en rad de lever du soleil
49
50	guez	98	REAL omegadeb, omegafin
51			REAL zfrac1, zfrac2, z1_mu, z2_mu
52			REAL lat_sun ! déclinaison en radians
53			REAL latr ! latitude du point de grille en radians
54
55			!----------------------------------------------------------------------
56
57	guez	118	if (present(frac)) call assert((/size(mu0), size(frac)/) == klon, "zenang")
58	guez	98
59			lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
60	guez	118	! Capderou (2003 784, equation 4.49)
61	guez	98
62	guez	206	gmtime1 = gmtime * 86400.
63			gmtime2 = gmtime * 86400. + pdtrad
64	guez	98
65			DO i = 1, klon
66	guez	206	latr = rlat(i) * pi / 180.
67			omega = 0. ! nuit polaire
68			IF (latr>=(pi / 2.-lat_sun) .OR. latr<=(-pi / 2.-lat_sun)) THEN
69	guez	118	omega = pi ! journée polaire
70	guez	98	END IF
71	guez	206	IF (latr<(pi / 2.+lat_sun) .AND. latr>(-pi / 2.+lat_sun) .AND. &
72			latr<(pi / 2.-lat_sun) .AND. latr>(-pi / 2.-lat_sun)) THEN
73			omega = -tan(latr) * tan(lat_sun)
74	guez	98	omega = acos(omega)
75			END IF
76
77	guez	206	omega1 = gmtime1 + rlon(i) * 86400. / 360.
78			omega1 = omega1 / 86400. * twopi
79	guez	125	omega1 = mod(omega1+twopi, twopi)
80	guez	98	omega1 = omega1 - pi
81
82	guez	206	omega2 = gmtime2 + rlon(i) * 86400. / 360.
83			omega2 = omega2 / 86400. * twopi
84	guez	125	omega2 = mod(omega2+twopi, twopi)
85	guez	98	omega2 = omega2 - pi
86
87	guez	118	IF (omega1<=omega2) THEN
88	guez	98	! on est dans la meme journee locale
89			IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN
90			! nuit
91	guez	206	IF (present(frac)) frac(i) = 0.
92			mu0(i) = 0.
93	guez	98	ELSE
94			! jour + nuit / jour
95			omegadeb = max(-omega, omega1)
96			omegafin = min(omega, omega2)
97	guez	206	IF (present(frac)) frac(i) = (omegafin-omegadeb) / (omega2-omega1)
98	guez	118	mu0(i) = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
99			* (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
100	guez	98	END IF
101	guez	118	ELSE
102			! omega1 > omega2, à cheval sur deux journées
103			! entre omega1 et pi
104			IF (omega1>=omega) THEN ! nuit
105	guez	206	zfrac1 = 0.
106			z1_mu = 0.
107	guez	118	ELSE ! jour+nuit
108	guez	98	omegadeb = max(-omega, omega1)
109			omegafin = omega
110			zfrac1 = omegafin - omegadeb
111	guez	125	z1_mu = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
112			* (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
113	guez	98	END IF
114	guez	118	! entre -pi et omega2
115			IF (omega2<=-omega) THEN ! nuit
116	guez	206	zfrac2 = 0.
117			z2_mu = 0.
118	guez	118	ELSE ! jour+nuit
119	guez	98	omegadeb = -omega
120			omegafin = min(omega, omega2)
121			zfrac2 = omegafin - omegadeb
122	guez	125	z2_mu = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
123			* (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
124	guez	98	END IF
125	guez	118	! moyenne
126	guez	206	IF (present(frac)) frac(i) = (zfrac1+zfrac2) / (omega2+twopi-omega1)
127			mu0(i) = (zfrac1 * z1_mu+zfrac2 * z2_mu) / max(zfrac1+zfrac2, 1e-10)
128	guez	118	END IF
129	guez	98	END DO
130
131			END SUBROUTINE zenang
132
133			end module zenang_m