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