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 nr_util, only: assert, pi

    ! 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	guez	39	use nr_util, only: assert, pi
67	guez	3
68			! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et
69			! "angle" de Z.X. Li
70
71			! Calcule les valeurs moyennes du cos de l'angle zénithal et
72			! l'ensoleillement moyen entre "gmtime1" et "gmtime2" connaissant la
73			! déclinaison, la latitude et la longitude.
74			! Différent de la routine "angle" en ce sens que "zenang" fournit des
75			! moyennes de "pmu0" et non des valeurs instantanées.
76			! Du coup "frac" prend toutes les valeurs entre 0 et 1.
77
78	guez	20	! Date : première version le 13 decembre 1994
79			! revu pour GCM le 30 septembre 1996
80	guez	3
81	guez	22	REAL, INTENT (IN):: longi
82	guez	3	! (longitude vraie de la terre dans son plan solaire a partir de
83	guez	20	! l'equinoxe de printemps) (in degrees)
84	guez	3
85	guez	22	REAL, INTENT (IN):: gmtime ! temps universel en fraction de jour
86			REAL, INTENT (IN):: pdtrad ! pas de temps du rayonnement (secondes)
87	guez	3
88	guez	22	REAL, INTENT (OUT):: pmu0(:) ! (klon)
89	guez	20	! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")
90	guez	3
91	guez	22	REAL, INTENT (OUT), OPTIONAL:: frac(:) ! (klon)
92	guez	3	! (ensoleillement moyen entre gmtime et gmtime+pdtrad)
93
94			! Variables local to the procedure:
95
96	guez	22	INTEGER i
97			REAL gmtime1, gmtime2
98			REAL deux_pi
99	guez	20
100	guez	22	REAL omega1, omega2, omega
101	guez	20	! omega1, omega2 : temps 1 et 2 exprimés en radians avec 0 à midi.
102			! omega : heure en radians du coucher de soleil
103			! -omega est donc l'heure en radians de lever du soleil
104
105	guez	22	REAL omegadeb, omegafin
106			REAL zfrac1, zfrac2, z1_mu, z2_mu
107			REAL lat_sun ! déclinaison en radians
108			REAL latr ! latitude du point de grille en radians
109	guez	3
110			!----------------------------------------------------------------------
111
112	guez	23	if (present(frac)) call assert((/size(pmu0), size(frac)/) == klon, &
113			"zenang")
114
115	guez	22	deux_pi = 2*pi
116	guez	3
117	guez	22	lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
118			! Capderou (2003 #784, équation 4.49)
119	guez	3
120	guez	22	gmtime1 = gmtime*86400.
121			gmtime2 = gmtime*86400. + pdtrad
122
123	guez	3	DO i = 1, klon
124	guez	22	latr = rlat(i)*pi/180.
125			omega = 0.0 !--nuit polaire
126			IF (latr>=(pi/2.-lat_sun) .OR. latr<=(-pi/2.-lat_sun)) THEN
127			omega = pi ! journee polaire
128			END IF
129			IF (latr<(pi/2.+lat_sun) .AND. latr>(-pi/2.+lat_sun) .AND. &
130			latr<(pi/2.-lat_sun) .AND. latr>(-pi/2.-lat_sun)) THEN
131			omega = -tan(latr)*tan(lat_sun)
132			omega = acos(omega)
133			END IF
134	guez	3
135			omega1 = gmtime1 + rlon(i)*86400.0/360.0
136	guez	22	omega1 = omega1/86400.0*deux_pi
137			omega1 = mod(omega1+deux_pi, deux_pi)
138	guez	20	omega1 = omega1 - pi
139	guez	3
140			omega2 = gmtime2 + rlon(i)*86400.0/360.0
141	guez	22	omega2 = omega2/86400.0*deux_pi
142			omega2 = mod(omega2+deux_pi, deux_pi)
143	guez	20	omega2 = omega2 - pi
144	guez	3
145	guez	22	TEST_OMEGA12: IF (omega1<=omega2) THEN
146	guez	3	! on est dans la meme journee locale
147	guez	22	IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN
148	guez	3	! nuit
149	guez	22	IF (present(frac)) frac(i) = 0.0
150			pmu0(i) = 0.0
151	guez	3	ELSE
152			! jour + nuit / jour
153	guez	22	omegadeb = max(-omega, omega1)
154			omegafin = min(omega, omega2)
155			IF (present(frac)) frac(i) = (omegafin-omegadeb)/(omega2-omega1)
156			pmu0(i) = sin(latr)sin(lat_sun) + cos(latr)cos(lat_sun)*(sin( &
157			omegafin)-sin(omegadeb))/(omegafin-omegadeb)
158			END IF
159			ELSE TEST_OMEGA12
160	guez	3	!---omega1 GT omega2 -- a cheval sur deux journees
161			!-------------------entre omega1 et pi
162	guez	22	IF (omega1>=omega) THEN !--nuit
163			zfrac1 = 0.0
164			z1_mu = 0.0
165			ELSE !--jour+nuit
166			omegadeb = max(-omega, omega1)
167			omegafin = omega
168			zfrac1 = omegafin - omegadeb
169			z1_mu = sin(latr)sin(lat_sun) + cos(latr)cos(lat_sun)*(sin( &
170			omegafin)-sin(omegadeb))/(omegafin-omegadeb)
171			END IF
172	guez	3	!---------------------entre -pi et omega2
173	guez	22	IF (omega2<=-omega) THEN !--nuit
174			zfrac2 = 0.0
175			z2_mu = 0.0
176			ELSE !--jour+nuit
177			omegadeb = -omega
178			omegafin = min(omega, omega2)
179			zfrac2 = omegafin - omegadeb
180			z2_mu = sin(latr)sin(lat_sun) + cos(latr)cos(lat_sun)*(sin( &
181			omegafin)-sin(omegadeb))/(omegafin-omegadeb)
182	guez	3
183	guez	22	END IF
184	guez	3	!-----------------------moyenne
185	guez	22	IF (present(frac)) frac(i) = (zfrac1+zfrac2)/ &
186			(omega2+deux_pi-omega1)
187			pmu0(i) = (zfrac1z1_mu+zfrac2z2_mu)/max(zfrac1+zfrac2, 1.E-10)
188			END IF TEST_OMEGA12
189			END DO
190	guez	3
191			END SUBROUTINE zenang
192
193	guez	22	END MODULE orbite_m