trunk/phylmd/zenang.f

module zenang_m

  IMPLICIT NONE

contains

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

    ! 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 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" parce que "zenang" fournit des
    ! moyennes de "mu0" et non des valeurs instantanées. Du coup
    ! "fract" 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:: fract(:) ! (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(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 = gmtime1 + pdtrad

    DO i = 1, klon
       latr = rlat(i) * pi / 180.
       
       IF (latr >= pi / 2. - lat_sun .OR. latr <= - pi / 2. - lat_sun) then
          omega = pi ! journée polaire
       else IF (latr < pi / 2. + lat_sun .AND. latr > - pi / 2. + lat_sun) THEN
          omega = acos(- tan(latr) * tan(lat_sun))
       else
          omega = 0. ! nuit polaire
       END IF

       omega1 = mod((gmtime1 + rlon(i) * 86400. / 360.) / 86400. * twopi &
            + twopi, twopi) - pi
       omega2 = mod((gmtime2 + rlon(i) * 86400. / 360.) / 86400. * twopi &
            + twopi, twopi) - 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(fract)) fract(i) = 0.
             mu0(i) = 0.
          ELSE
             ! jour + nuit / jour
             omegadeb = max(- omega, omega1)
             omegafin = min(omega, omega2)
             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
             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(fract)) fract(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	213	SUBROUTINE zenang(longi, gmtime, pdtrad, mu0, fract)
8	guez	98
9	guez	118	! Author: O. Boucher (LMD/CNRS), d'après les routines "zenith" et
10	guez	210	! "angle" de Z. X. Li
11	guez	98
12	guez	118	! Date : première version le 13 décembre 1994, revu pour GCM le 30
13			! septembre 1996
14
15	guez	213	! Calcule les valeurs moyennes du cosinus 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	guez	213	! Différent de la routine "angle" parce que "zenang" fournit des
19			! moyennes de "mu0" et non des valeurs instantanées. Du coup
20			! "fract" prend toutes les valeurs entre 0 et 1. Cf. Capderou (2003
21	guez	125	! 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	guez	213	! cosine of mean zenith angle between "gmtime" and "gmtime + pdtrad"
37	guez	98
38	guez	213	REAL, INTENT(OUT), OPTIONAL:: fract(:) ! (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	guez	213	! "- omega" est donc l'heure en rad de lever du soleil.
49	guez	118
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	213	if (present(fract)) call assert((/size(mu0), size(fract)/) == klon, &
58			"zenang")
59	guez	98
60			lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
61	guez	118	! Capderou (2003 784, equation 4.49)
62	guez	98
63	guez	206	gmtime1 = gmtime * 86400.
64	guez	213	gmtime2 = gmtime1 + pdtrad
65	guez	98
66			DO i = 1, klon
67	guez	206	latr = rlat(i) * pi / 180.
68	guez	213
69			IF (latr >= pi / 2. - lat_sun .OR. latr <= - pi / 2. - lat_sun) then
70	guez	118	omega = pi ! journée polaire
71	guez	213	else IF (latr < pi / 2. + lat_sun .AND. latr > - pi / 2. + lat_sun) THEN
72			omega = acos(- tan(latr) * tan(lat_sun))
73			else
74			omega = 0. ! nuit polaire
75	guez	98	END IF
76
77	guez	213	omega1 = mod((gmtime1 + rlon(i) * 86400. / 360.) / 86400. * twopi &
78			+ twopi, twopi) - pi
79			omega2 = mod((gmtime2 + rlon(i) * 86400. / 360.) / 86400. * twopi &
80			+ twopi, twopi) - pi
81	guez	98
82	guez	213	IF (omega1 <= omega2) THEN
83	guez	98	! on est dans la meme journee locale
84	guez	213	IF (omega2 <= - omega .OR. omega1 >= omega .OR. omega < 1E-5) THEN
85	guez	98	! nuit
86	guez	213	IF (present(fract)) fract(i) = 0.
87	guez	206	mu0(i) = 0.
88	guez	98	ELSE
89			! jour + nuit / jour
90	guez	213	omegadeb = max(- omega, omega1)
91	guez	98	omegafin = min(omega, omega2)
92	guez	213	IF (present(fract)) fract(i) = (omegafin - omegadeb) &
93			/ (omega2 - omega1)
94	guez	118	mu0(i) = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
95			* (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
96	guez	98	END IF
97	guez	118	ELSE
98			! omega1 > omega2, à cheval sur deux journées
99			! entre omega1 et pi
100	guez	213	IF (omega1 >= omega) THEN ! nuit
101	guez	206	zfrac1 = 0.
102			z1_mu = 0.
103	guez	213	ELSE ! jour + nuit
104			omegadeb = max(- omega, omega1)
105	guez	98	omegafin = omega
106			zfrac1 = omegafin - omegadeb
107	guez	125	z1_mu = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
108			* (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
109	guez	98	END IF
110	guez	213	! entre - pi et omega2
111			IF (omega2 <= - omega) THEN ! nuit
112	guez	206	zfrac2 = 0.
113			z2_mu = 0.
114	guez	213	ELSE ! jour + nuit
115			omegadeb = - omega
116	guez	98	omegafin = min(omega, omega2)
117			zfrac2 = omegafin - omegadeb
118	guez	125	z2_mu = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
119			* (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
120	guez	98	END IF
121	guez	118	! moyenne
122	guez	213	IF (present(fract)) fract(i) = (zfrac1 + zfrac2) &
123			/ (omega2 + twopi - omega1)
124			mu0(i) = (zfrac1 * z1_mu + zfrac2 * z2_mu) &
125			/ max(zfrac1 + zfrac2, 1e-10)
126	guez	118	END IF
127	guez	98	END DO
128
129			END SUBROUTINE zenang
130
131			end module zenang_m