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	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	USE dimphy, ONLY: klon
10			USE yomcst, ONLY: r_incl
11			USE phyetat0_m, ONLY: rlat, rlon
12	guez	98	use nr_util, only: assert, pi
13
14	guez	118	! Author: O. Boucher (LMD/CNRS), d'après les routines "zenith" et
15	guez	98	! "angle" de Z.X. Li
16
17	guez	118	! Date : première version le 13 décembre 1994, revu pour GCM le 30
18			! septembre 1996
19
20	guez	98	! Calcule les valeurs moyennes du cos de l'angle zénithal et
21	guez	118	! 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	guez	98
27	guez	118	REAL, INTENT(IN):: longi
28			! longitude vraie de la terre dans son plan solaire à partir de
29			! l'équinoxe de printemps (in degrees)
30	guez	98
31	guez	118	REAL, INTENT(IN):: gmtime ! temps universel en fraction de jour
32			REAL, INTENT(IN):: pdtrad ! pas de temps du rayonnement (s)
33	guez	98
34	guez	118	REAL, INTENT(OUT):: mu0(:) ! (klon)
35			! cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad"
36	guez	98
37	guez	118	REAL, INTENT(OUT), OPTIONAL:: frac(:) ! (klon)
38			! ensoleillement moyen entre gmtime et gmtime+pdtrad
39	guez	98
40	guez	118	! Local:
41	guez	98
42			INTEGER i
43			REAL gmtime1, gmtime2
44			REAL deux_pi
45
46	guez	118	REAL omega1, omega2 ! temps 1 et 2 exprimés en radians avec 0 à midi
47	guez	98
48	guez	118	REAL omega ! heure en rad du coucher de soleil
49			! - omega est donc l'heure en rad de lever du soleil
50
51	guez	98	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	guez	118	if (present(frac)) call assert((/size(mu0), size(frac)/) == klon, "zenang")
59	guez	98
60			deux_pi = 2*pi
61
62			lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
63	guez	118	! Capderou (2003 784, equation 4.49)
64	guez	98
65			gmtime1 = gmtime*86400.
66			gmtime2 = gmtime*86400. + pdtrad
67
68			DO i = 1, klon
69			latr = rlat(i)*pi/180.
70	guez	118	omega = 0.0 ! nuit polaire
71	guez	98	IF (latr>=(pi/2.-lat_sun) .OR. latr<=(-pi/2.-lat_sun)) THEN
72	guez	118	omega = pi ! journée polaire
73	guez	98	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	guez	118	IF (omega1<=omega2) THEN
91	guez	98	! 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	guez	118	mu0(i) = 0.0
96	guez	98	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	guez	118	mu0(i) = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
102			* (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
103	guez	98	END IF
104	guez	118	ELSE
105			! omega1 > omega2, à cheval sur deux journées
106			! entre omega1 et pi
107			IF (omega1>=omega) THEN ! nuit
108	guez	98	zfrac1 = 0.0
109			z1_mu = 0.0
110	guez	118	ELSE ! jour+nuit
111	guez	98	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	guez	118	! entre -pi et omega2
118			IF (omega2<=-omega) THEN ! nuit
119	guez	98	zfrac2 = 0.0
120			z2_mu = 0.0
121	guez	118	ELSE ! jour+nuit
122	guez	98	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	guez	118	! 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	guez	98	END DO
134
135			END SUBROUTINE zenang
136
137			end module zenang_m