trunk/phylmd/orbite.f

MODULE orbite_m

  ! From phylmd/orbite.F, v 1.1.1.1 2004/05/19 12:53:08

  IMPLICIT NONE

CONTAINS

  SUBROUTINE orbite(xjour, longi, dist)

    USE yomcst, ONLY : r_ecc, r_peri
    use nr_util, only: pi

    ! Auteur(s): Z.X. Li (LMD/CNRS)
    ! Date: 1993/08/18
    ! Pour un jour donné, calcule la longitude vraie de la Terre (par
    ! rapport au point vernal, 21 mars) dans son orbite solaire. Calcule aussi
    ! la distance Terre-Soleil, c'est-à-dire l'unité astronomique.

    REAL, INTENT (IN):: xjour ! jour de l'année à compter du premier janvier

    REAL, INTENT (OUT):: longi
    ! longitude vraie de la Terre dans son orbite solaire, par
    ! rapport au point vernal (21 mars), en degrés

    REAL, INTENT (OUT), OPTIONAL:: dist
    ! distance terre-soleil (par rapport a la moyenne)

    ! Variables locales
    REAL pir, xl, xllp, xee, xse, xlam, anm, ranm, ranv

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

    pir = pi / 180.
    xl = r_peri + 180.
    xllp = xl * pir
    xee = r_ecc * r_ecc
    xse = sqrt(1. - xee)
    xlam = (r_ecc / 2 + r_ecc * xee / 8.) * (1. + xse) * sin(xllp) &
         - xee / 4. * (0.5 + xse) * sin(2.*xllp) &
         + r_ecc * xee / 8. * (1. / 3. + xse) * sin(3.*xllp)
    xlam = 2. * xlam / pir
    anm = xlam + (xjour - 81.) - xl
    ranm = anm * pir
    xee = xee * r_ecc
    ranv = ranm + (2. * r_ecc - xee / 4.) * sin(ranm) + &
         5. / 4. * r_ecc * r_ecc * sin(2 * ranm) &
         + 13. / 12. * xee * sin(3.*ranm)

    longi = ranv / pir + xl

    IF (present(dist)) then
       dist = (1 - r_ecc*r_ecc) &
            / (1 + r_ecc*cos(pir*(longi - (r_peri + 180.))))
    end IF

  END SUBROUTINE orbite

  !*****************************************************************

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

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

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

    ! Calcule les valeurs moyennes du cos de l'angle zénithal et
    ! l'ensoleillement moyen entre "gmtime1" et "gmtime2" connaissant la
    ! déclinaison, la latitude et la longitude.
    ! Différent de la routine "angle" en ce sens que "zenang" fournit des
    ! moyennes de "pmu0" et non des valeurs instantanées.
    ! Du coup "frac" prend toutes les valeurs entre 0 et 1.

    ! Date : première version le 13 decembre 1994
    !        revu pour  GCM  le 30 septembre 1996

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

    REAL, INTENT (IN):: gmtime ! temps universel en fraction de jour
    REAL, INTENT (IN):: pdtrad ! pas de temps du rayonnement (secondes)

    REAL, INTENT (OUT):: pmu0(:) ! (klon)
    ! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")

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

    ! Variables local to the procedure:

    INTEGER i
    REAL gmtime1, gmtime2
    REAL deux_pi

    REAL omega1, omega2, omega
    ! omega1, omega2 : temps 1 et 2 exprimés en radians avec 0 à midi.
    ! omega : heure en radians du coucher de soleil 
    ! -omega est donc l'heure en radians 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(pmu0), size(frac)/) == klon, &
         "zenang")
    
    deux_pi = 2*pi

    lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
    ! Capderou (2003 #784, équation 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 ! journee 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

       TEST_OMEGA12: 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
             pmu0(i) = 0.0
          ELSE
             ! jour + nuit / jour
             omegadeb = max(-omega, omega1)
             omegafin = min(omega, omega2)
             IF (present(frac)) frac(i) = (omegafin-omegadeb)/(omega2-omega1)
             pmu0(i) = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &
                  omegafin)-sin(omegadeb))/(omegafin-omegadeb)
          END IF
       ELSE TEST_OMEGA12
          !---omega1 GT omega2 -- a cheval sur deux journees
          !-------------------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)
          pmu0(i) = (zfrac1*z1_mu+zfrac2*z2_mu)/max(zfrac1+zfrac2, 1.E-10)
       END IF TEST_OMEGA12
    END DO

  END SUBROUTINE zenang

END MODULE orbite_m
1	guez	22	MODULE orbite_m
2	guez	3
3			! From phylmd/orbite.F, v 1.1.1.1 2004/05/19 12:53:08
4
5	guez	22	IMPLICIT NONE
6	guez	3
7	guez	22	CONTAINS
8	guez	3
9			SUBROUTINE orbite(xjour, longi, dist)
10
11	guez	22	USE yomcst, ONLY : r_ecc, r_peri
12	guez	36	use nr_util, only: pi
13	guez	3
14	guez	22	! Auteur(s): Z.X. Li (LMD/CNRS)
15			! Date: 1993/08/18
16			! Pour un jour donné, calcule la longitude vraie de la Terre (par
17			! rapport au point vernal, 21 mars) dans son orbite solaire. Calcule aussi
18	guez	3	! la distance Terre-Soleil, c'est-à-dire l'unité astronomique.
19
20	guez	22	REAL, INTENT (IN):: xjour ! jour de l'année à compter du premier janvier
21	guez	3
22	guez	22	REAL, INTENT (OUT):: longi
23			! longitude vraie de la Terre dans son orbite solaire, par
24			! rapport au point vernal (21 mars), en degrés
25	guez	3
26	guez	22	REAL, INTENT (OUT), OPTIONAL:: dist
27			! distance terre-soleil (par rapport a la moyenne)
28	guez	3
29			! Variables locales
30	guez	22	REAL pir, xl, xllp, xee, xse, xlam, anm, ranm, ranv
31	guez	3
32			!----------------------------------------------------------------------
33
34	guez	22	pir = pi / 180.
35			xl = r_peri + 180.
36			xllp = xl * pir
37			xee = r_ecc * r_ecc
38			xse = sqrt(1. - xee)
39			xlam = (r_ecc / 2 + r_ecc * xee / 8.) * (1. + xse) * sin(xllp) &
40			- xee / 4. * (0.5 + xse) * sin(2.*xllp) &
41			+ r_ecc * xee / 8. * (1. / 3. + xse) * sin(3.*xllp)
42			xlam = 2. * xlam / pir
43			anm = xlam + (xjour - 81.) - xl
44			ranm = anm * pir
45			xee = xee * r_ecc
46			ranv = ranm + (2. * r_ecc - xee / 4.) * sin(ranm) + &
47			5. / 4. * r_ecc * r_ecc * sin(2 * ranm) &
48			+ 13. / 12. * xee * sin(3.*ranm)
49	guez	3
50	guez	22	longi = ranv / pir + xl
51	guez	3
52	guez	22	IF (present(dist)) then
53			dist = (1 - r_ecc*r_ecc) &
54			/ (1 + r_ecccos(pir(longi - (r_peri + 180.))))
55			end IF
56	guez	3
57			END SUBROUTINE orbite
58
59			!*****************************************************************
60
61			SUBROUTINE zenang(longi, gmtime, pdtrad, pmu0, frac)
62
63	guez	22	USE dimphy, ONLY : klon
64			USE yomcst, ONLY : r_incl
65			USE phyetat0_m, ONLY : rlat, rlon
66			USE comconst, ONLY : pi
67	guez	36	use nr_util, only: assert
68	guez	3
69			! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et
70			! "angle" de Z.X. Li
71
72			! Calcule les valeurs moyennes du cos de l'angle zénithal et
73			! l'ensoleillement moyen entre "gmtime1" et "gmtime2" connaissant la
74			! déclinaison, la latitude et la longitude.
75			! Différent de la routine "angle" en ce sens que "zenang" fournit des
76			! moyennes de "pmu0" et non des valeurs instantanées.
77			! Du coup "frac" prend toutes les valeurs entre 0 et 1.
78
79	guez	20	! Date : première version le 13 decembre 1994
80			! revu pour GCM le 30 septembre 1996
81	guez	3
82	guez	22	REAL, INTENT (IN):: longi
83	guez	3	! (longitude vraie de la terre dans son plan solaire a partir de
84	guez	20	! l'equinoxe de printemps) (in degrees)
85	guez	3
86	guez	22	REAL, INTENT (IN):: gmtime ! temps universel en fraction de jour
87			REAL, INTENT (IN):: pdtrad ! pas de temps du rayonnement (secondes)
88	guez	3
89	guez	22	REAL, INTENT (OUT):: pmu0(:) ! (klon)
90	guez	20	! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")
91	guez	3
92	guez	22	REAL, INTENT (OUT), OPTIONAL:: frac(:) ! (klon)
93	guez	3	! (ensoleillement moyen entre gmtime et gmtime+pdtrad)
94
95			! Variables local to the procedure:
96
97	guez	22	INTEGER i
98			REAL gmtime1, gmtime2
99			REAL deux_pi
100	guez	20
101	guez	22	REAL omega1, omega2, omega
102	guez	20	! omega1, omega2 : temps 1 et 2 exprimés en radians avec 0 à midi.
103			! omega : heure en radians du coucher de soleil
104			! -omega est donc l'heure en radians de lever du soleil
105
106	guez	22	REAL omegadeb, omegafin
107			REAL zfrac1, zfrac2, z1_mu, z2_mu
108			REAL lat_sun ! déclinaison en radians
109			REAL latr ! latitude du point de grille en radians
110	guez	3
111			!----------------------------------------------------------------------
112
113	guez	23	if (present(frac)) call assert((/size(pmu0), size(frac)/) == klon, &
114			"zenang")
115
116	guez	22	deux_pi = 2*pi
117	guez	3
118	guez	22	lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
119			! Capderou (2003 #784, équation 4.49)
120	guez	3
121	guez	22	gmtime1 = gmtime*86400.
122			gmtime2 = gmtime*86400. + pdtrad
123
124	guez	3	DO i = 1, klon
125	guez	22	latr = rlat(i)*pi/180.
126			omega = 0.0 !--nuit polaire
127			IF (latr>=(pi/2.-lat_sun) .OR. latr<=(-pi/2.-lat_sun)) THEN
128			omega = pi ! journee polaire
129			END IF
130			IF (latr<(pi/2.+lat_sun) .AND. latr>(-pi/2.+lat_sun) .AND. &
131			latr<(pi/2.-lat_sun) .AND. latr>(-pi/2.-lat_sun)) THEN
132			omega = -tan(latr)*tan(lat_sun)
133			omega = acos(omega)
134			END IF
135	guez	3
136			omega1 = gmtime1 + rlon(i)*86400.0/360.0
137	guez	22	omega1 = omega1/86400.0*deux_pi
138			omega1 = mod(omega1+deux_pi, deux_pi)
139	guez	20	omega1 = omega1 - pi
140	guez	3
141			omega2 = gmtime2 + rlon(i)*86400.0/360.0
142	guez	22	omega2 = omega2/86400.0*deux_pi
143			omega2 = mod(omega2+deux_pi, deux_pi)
144	guez	20	omega2 = omega2 - pi
145	guez	3
146	guez	22	TEST_OMEGA12: IF (omega1<=omega2) THEN
147	guez	3	! on est dans la meme journee locale
148	guez	22	IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN
149	guez	3	! nuit
150	guez	22	IF (present(frac)) frac(i) = 0.0
151			pmu0(i) = 0.0
152	guez	3	ELSE
153			! jour + nuit / jour
154	guez	22	omegadeb = max(-omega, omega1)
155			omegafin = min(omega, omega2)
156			IF (present(frac)) frac(i) = (omegafin-omegadeb)/(omega2-omega1)
157			pmu0(i) = sin(latr)sin(lat_sun) + cos(latr)cos(lat_sun)*(sin( &
158			omegafin)-sin(omegadeb))/(omegafin-omegadeb)
159			END IF
160			ELSE TEST_OMEGA12
161	guez	3	!---omega1 GT omega2 -- a cheval sur deux journees
162			!-------------------entre omega1 et pi
163	guez	22	IF (omega1>=omega) THEN !--nuit
164			zfrac1 = 0.0
165			z1_mu = 0.0
166			ELSE !--jour+nuit
167			omegadeb = max(-omega, omega1)
168			omegafin = omega
169			zfrac1 = omegafin - omegadeb
170			z1_mu = sin(latr)sin(lat_sun) + cos(latr)cos(lat_sun)*(sin( &
171			omegafin)-sin(omegadeb))/(omegafin-omegadeb)
172			END IF
173	guez	3	!---------------------entre -pi et omega2
174	guez	22	IF (omega2<=-omega) THEN !--nuit
175			zfrac2 = 0.0
176			z2_mu = 0.0
177			ELSE !--jour+nuit
178			omegadeb = -omega
179			omegafin = min(omega, omega2)
180			zfrac2 = omegafin - omegadeb
181			z2_mu = sin(latr)sin(lat_sun) + cos(latr)cos(lat_sun)*(sin( &
182			omegafin)-sin(omegadeb))/(omegafin-omegadeb)
183	guez	3
184	guez	22	END IF
185	guez	3	!-----------------------moyenne
186	guez	22	IF (present(frac)) frac(i) = (zfrac1+zfrac2)/ &
187			(omega2+deux_pi-omega1)
188			pmu0(i) = (zfrac1z1_mu+zfrac2z2_mu)/max(zfrac1+zfrac2, 1.E-10)
189			END IF TEST_OMEGA12
190			END DO
191	guez	3
192			END SUBROUTINE zenang
193
194	guez	22	END MODULE orbite_m