libf/phylmd/orbite.f90

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 numer_rec, 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 numer_rec, 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	MODULE orbite_m
2
3	! From phylmd/orbite.F, v 1.1.1.1 2004/05/19 12:53:08
4
5	IMPLICIT NONE
6
7	CONTAINS
8
9	SUBROUTINE orbite(xjour, longi, dist)
10
11	USE yomcst, ONLY : r_ecc, r_peri
12	use numer_rec, only: pi
13
14	! 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	! la distance Terre-Soleil, c'est-à-dire l'unité astronomique.
19
20	REAL, INTENT (IN):: xjour ! jour de l'année à compter du premier janvier
21
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
26	REAL, INTENT (OUT), OPTIONAL:: dist
27	! distance terre-soleil (par rapport a la moyenne)
28
29	! Variables locales
30	REAL pir, xl, xllp, xee, xse, xlam, anm, ranm, ranv
31
32	!----------------------------------------------------------------------
33
34	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
50	longi = ranv / pir + xl
51
52	IF (present(dist)) then
53	dist = (1 - r_ecc*r_ecc) &
54	/ (1 + r_ecccos(pir(longi - (r_peri + 180.))))
55	end IF
56
57	END SUBROUTINE orbite
58
59	!*****************************************************************
60
61	SUBROUTINE zenang(longi, gmtime, pdtrad, pmu0, frac)
62
63	USE dimphy, ONLY : klon
64	USE yomcst, ONLY : r_incl
65	USE phyetat0_m, ONLY : rlat, rlon
66	USE comconst, ONLY : pi
67	use numer_rec, only: assert
68
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	! Date : première version le 13 decembre 1994
80	! revu pour GCM le 30 septembre 1996
81
82	REAL, INTENT (IN):: longi
83	! (longitude vraie de la terre dans son plan solaire a partir de
84	! l'equinoxe de printemps) (in degrees)
85
86	REAL, INTENT (IN):: gmtime ! temps universel en fraction de jour
87	REAL, INTENT (IN):: pdtrad ! pas de temps du rayonnement (secondes)
88
89	REAL, INTENT (OUT):: pmu0(:) ! (klon)
90	! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")
91
92	REAL, INTENT (OUT), OPTIONAL:: frac(:) ! (klon)
93	! (ensoleillement moyen entre gmtime et gmtime+pdtrad)
94
95	! Variables local to the procedure:
96
97	INTEGER i
98	REAL gmtime1, gmtime2
99	REAL deux_pi
100
101	REAL omega1, omega2, omega
102	! 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	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
111	!----------------------------------------------------------------------
112
113	if (present(frac)) call assert((/size(pmu0), size(frac)/) == klon, &
114	"zenang")
115
116	deux_pi = 2*pi
117
118	lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
119	! Capderou (2003 #784, équation 4.49)
120
121	gmtime1 = gmtime*86400.
122	gmtime2 = gmtime*86400. + pdtrad
123
124	DO i = 1, klon
125	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
136	omega1 = gmtime1 + rlon(i)*86400.0/360.0
137	omega1 = omega1/86400.0*deux_pi
138	omega1 = mod(omega1+deux_pi, deux_pi)
139	omega1 = omega1 - pi
140
141	omega2 = gmtime2 + rlon(i)*86400.0/360.0
142	omega2 = omega2/86400.0*deux_pi
143	omega2 = mod(omega2+deux_pi, deux_pi)
144	omega2 = omega2 - pi
145
146	TEST_OMEGA12: IF (omega1<=omega2) THEN
147	! on est dans la meme journee locale
148	IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN
149	! nuit
150	IF (present(frac)) frac(i) = 0.0
151	pmu0(i) = 0.0
152	ELSE
153	! jour + nuit / jour
154	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	!---omega1 GT omega2 -- a cheval sur deux journees
162	!-------------------entre omega1 et pi
163	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	!---------------------entre -pi et omega2
174	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
184	END IF
185	!-----------------------moyenne
186	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
192	END SUBROUTINE zenang
193
194	END MODULE orbite_m