trunk/phylmd/zenang.f

module zenang_m

  IMPLICIT NONE

contains

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

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

    ! 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.

    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 deux_pi

    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")
    
    deux_pi = 2*pi

    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.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.0/360.0
       omega1 = omega1/86400.0*deux_pi
       omega1 = mod(omega1+deux_pi, deux_pi)
       omega1 = omega1 - pi

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