Sources/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

    ! 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'annee a compter du 1er 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, dlamm, anm, ranm, anv, ranv

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

    pir = 4.0*ATAN(1.0) / 180.0
    xl=R_peri+180.0
    xllp=xl*pir
    xee=R_ecc*R_ecc
    xse=SQRT(1.0-xee)
    xlam = (R_ecc/2.0+R_ecc*xee/8.0)*(1.0+xse)*SIN(xllp) &
         - xee/4.0*(0.5+xse)*SIN(2.0*xllp) &
         + R_ecc*xee/8.0*(1.0/3.0+xse)*SIN(3.0*xllp)
    xlam=2.0*xlam/pir
    dlamm=xlam+(xjour-81.0)
    anm=dlamm-xl
    ranm=anm*pir
    xee=xee*R_ecc
    ranv=ranm+(2.0*R_ecc-xee/4.0)*SIN(ranm) &
         +5.0/4.0*R_ecc*R_ecc*SIN(2.0*ranm) &
         +13.0/12.0*xee*SIN(3.0*ranm)

    anv=ranv/pir
    longi=anv+xl

    if (present(dist)) dist &
         = (1-R_ecc*R_ecc) /(1+R_ecc*COS(pir*(longi-(R_peri+180.0))))

  END SUBROUTINE orbite

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

  SUBROUTINE angle(longi, lati, frac, muzero)

    use dimphy, only: klon
    use YOMCST, only: R_incl

    ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18
    ! Calcule la durée d'ensoleillement pour un jour et la hauteur du
    ! soleil (cosinus de l'angle zénithal) moyennée sur la journee.

    ! Arguments:
    ! longi----INPUT-R- la longitude vraie de la terre dans son plan 
    !                   solaire a partir de l'equinoxe de printemps (degre)
    ! lati-----INPUT-R- la latitude d'un point sur la terre (degre)
    ! frac-----OUTPUT-R la duree d'ensoleillement dans la journee divisee
    !                   par 24 heures (unite en fraction de 0 a 1)
    ! muzero---OUTPUT-R la moyenne du cosinus de l'angle zinithal sur
    !                   la journee (0 a 1)

    REAL longi
    REAL lati(klon), frac(klon), muzero(klon)
    REAL lat, omega, lon_sun, lat_sun
    REAL pi_local, incl
    INTEGER i

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

    pi_local = 4.0 * ATAN(1.0)
    incl=R_incl * pi_local / 180.

    lon_sun = longi * pi_local / 180.0
    lat_sun = ASIN (sin(lon_sun)*SIN(incl) )

    DO i = 1, klon
       lat = lati(i) * pi_local / 180.0

       IF ( lat .GE. (pi_local/2.+lat_sun) &
            .OR. lat <= (-pi_local/2.+lat_sun)) THEN
          omega = 0.0   ! nuit polaire
       ELSE IF ( lat.GE.(pi_local/2.-lat_sun) &
            .OR. lat <= (-pi_local/2.-lat_sun)) THEN
          omega = pi_local   ! journee polaire
       ELSE
          omega = -TAN(lat)*TAN(lat_sun)
          omega = ACOS (omega)
       ENDIF

       frac(i) = omega / pi_local

       IF (omega .GT. 0.0) THEN
          muzero(i) = SIN(lat)*SIN(lat_sun) &
               + COS(lat)*COS(lat_sun)*SIN(omega) / omega
       ELSE
          muzero(i) = 0.0
       ENDIF
    ENDDO

  END SUBROUTINE angle

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

  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

    ! 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, incl

    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 lon_sun          ! longitude solaire en radians
    real latr             ! latitude du point de grille en radians

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

    deux_pi = 2 * pi
    incl=R_incl * pi / 180.

    lon_sun = longi * pi / 180.
    lat_sun = ASIN(SIN(lon_sun)*SIN(incl))

    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
       ENDIF
       IF (latr < (pi/2.+lat_sun).AND. &
            latr.GT.(-pi/2.+lat_sun).AND. &
            latr < (pi/2.-lat_sun).AND. &
            latr.GT.(-pi/2.-lat_sun)) THEN
          omega = -TAN(latr)*TAN(lat_sun)
          omega = ACOS(omega)
       ENDIF

       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)        
          ENDIF
       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)
          ENDIF
          !---------------------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)

          ENDIF
          !-----------------------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)
       ENDIF test_omega12
    ENDDO

  END SUBROUTINE zenang

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

  SUBROUTINE zenith (longi, gmtime, lat, long, pmu0, fract)

    use dimphy, only: klon
    use YOMCST, only: R_incl
    use comconst, only: pi

    ! Author: Z.X. Li (LMD/ENS)

    ! Calcule le cosinus de l'angle zénithal du Soleil en connaissant
    ! la déclinaison du Soleil, la latitude et la longitude du point
    ! sur la Terre, et le temps universel.

    REAL, intent(in):: longi ! déclinaison du soleil (en degrés)

    REAL, intent(in):: gmtime
    ! (temps universel en fraction de jour, entre 0 et 1)

    REAL, intent(in):: lat(klon), long(klon) ! latitude et longitude en degres
    REAL, intent(out):: pmu0(klon) ! cosinus de l'angle zenithal
    REAL, intent(out):: fract(klon)

    ! Variables local to the procedure:

    INTEGER n
    REAL zpir, omega
    REAL lat_sun

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

    zpir = pi / 180.0
    lat_sun = ASIN (SIN(longi * zpir)*SIN(R_incl * zpir) )

    ! Initialisation à la nuit:
    pmu0 = 0.
    fract = 0.

    ! 1 degree in longitude = 240 s in time

    DO n = 1, klon
       omega = (gmtime + long(n) / 360.) * 2 * pi
       omega = MOD(omega + 2 * pi, 2 * pi)
       omega = omega - pi
       pmu0(n) = sin(lat(n)*zpir) * sin(lat_sun) &
            + cos(lat(n)*zpir) * cos(lat_sun) &
            * cos(omega)
       pmu0(n) = MAX (pmu0(n), 0.)
       IF (pmu0(n) > 1.E-6) fract(n)=1.0
    ENDDO

  END SUBROUTINE zenith

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
12
13	! Auteur(s): Z.X. Li (LMD/CNRS) Date: 1993/08/18 Pour un jour
14	! donné, calcule la longitude vraie de la Terre (par rapport au
15	! point vernal, 21 mars) dans son orbite solaire. Calcule aussi
16	! la distance Terre-Soleil, c'est-à-dire l'unité astronomique.
17
18	REAL, intent(in):: xjour ! jour de l'annee a compter du 1er janvier
19
20	real, intent(out):: longi
21	! (longitude vraie de la Terre dans son orbite solaire, par
22	! rapport au point vernal (21 mars), en degrés)
23
24	real, intent(out), optional:: dist
25	! (distance terre-soleil (par rapport a la moyenne))
26
27	! Variables locales
28	REAL pir, xl, xllp, xee, xse, xlam, dlamm, anm, ranm, anv, ranv
29
30	!----------------------------------------------------------------------
31
32	pir = 4.0*ATAN(1.0) / 180.0
33	xl=R_peri+180.0
34	xllp=xl*pir
35	xee=R_ecc*R_ecc
36	xse=SQRT(1.0-xee)
37	xlam = (R_ecc/2.0+R_eccxee/8.0)(1.0+xse)*SIN(xllp) &
38	- xee/4.0(0.5+xse)SIN(2.0*xllp) &
39	+ R_eccxee/8.0(1.0/3.0+xse)SIN(3.0xllp)
40	xlam=2.0*xlam/pir
41	dlamm=xlam+(xjour-81.0)
42	anm=dlamm-xl
43	ranm=anm*pir
44	xee=xee*R_ecc
45	ranv=ranm+(2.0R_ecc-xee/4.0)SIN(ranm) &
46	+5.0/4.0R_eccR_eccSIN(2.0ranm) &
47	+13.0/12.0xeeSIN(3.0*ranm)
48
49	anv=ranv/pir
50	longi=anv+xl
51
52	if (present(dist)) dist &
53	= (1-R_eccR_ecc) /(1+R_eccCOS(pir*(longi-(R_peri+180.0))))
54
55	END SUBROUTINE orbite
56
57	!*****************************************************************
58
59	SUBROUTINE angle(longi, lati, frac, muzero)
60
61	use dimphy, only: klon
62	use YOMCST, only: R_incl
63
64	! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18
65	! Calcule la durée d'ensoleillement pour un jour et la hauteur du
66	! soleil (cosinus de l'angle zénithal) moyennée sur la journee.
67
68	! Arguments:
69	! longi----INPUT-R- la longitude vraie de la terre dans son plan
70	! solaire a partir de l'equinoxe de printemps (degre)
71	! lati-----INPUT-R- la latitude d'un point sur la terre (degre)
72	! frac-----OUTPUT-R la duree d'ensoleillement dans la journee divisee
73	! par 24 heures (unite en fraction de 0 a 1)
74	! muzero---OUTPUT-R la moyenne du cosinus de l'angle zinithal sur
75	! la journee (0 a 1)
76
77	REAL longi
78	REAL lati(klon), frac(klon), muzero(klon)
79	REAL lat, omega, lon_sun, lat_sun
80	REAL pi_local, incl
81	INTEGER i
82
83	!----------------------------------------------------------------------
84
85	pi_local = 4.0 * ATAN(1.0)
86	incl=R_incl * pi_local / 180.
87
88	lon_sun = longi * pi_local / 180.0
89	lat_sun = ASIN (sin(lon_sun)*SIN(incl) )
90
91	DO i = 1, klon
92	lat = lati(i) * pi_local / 180.0
93
94	IF ( lat .GE. (pi_local/2.+lat_sun) &
95	.OR. lat <= (-pi_local/2.+lat_sun)) THEN
96	omega = 0.0 ! nuit polaire
97	ELSE IF ( lat.GE.(pi_local/2.-lat_sun) &
98	.OR. lat <= (-pi_local/2.-lat_sun)) THEN
99	omega = pi_local ! journee polaire
100	ELSE
101	omega = -TAN(lat)*TAN(lat_sun)
102	omega = ACOS (omega)
103	ENDIF
104
105	frac(i) = omega / pi_local
106
107	IF (omega .GT. 0.0) THEN
108	muzero(i) = SIN(lat)*SIN(lat_sun) &
109	+ COS(lat)COS(lat_sun)SIN(omega) / omega
110	ELSE
111	muzero(i) = 0.0
112	ENDIF
113	ENDDO
114
115	END SUBROUTINE angle
116
117	!*****************************************************************
118
119	SUBROUTINE zenang(longi, gmtime, pdtrad, pmu0, frac)
120
121	use dimphy, only: klon
122	use YOMCST, only: r_incl
123	use phyetat0_m, only: rlat, rlon
124	use comconst, only: pi
125
126	! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et
127	! "angle" de Z.X. Li
128
129	! Calcule les valeurs moyennes du cos de l'angle zénithal et
130	! l'ensoleillement moyen entre "gmtime1" et "gmtime2" connaissant la
131	! déclinaison, la latitude et la longitude.
132	! Différent de la routine "angle" en ce sens que "zenang" fournit des
133	! moyennes de "pmu0" et non des valeurs instantanées.
134	! Du coup "frac" prend toutes les valeurs entre 0 et 1.
135
136	! Date : première version le 13 decembre 1994
137	! revu pour GCM le 30 septembre 1996
138
139	real, intent(in):: longi
140	! (longitude vraie de la terre dans son plan solaire a partir de
141	! l'equinoxe de printemps) (in degrees)
142
143	real, intent(in):: gmtime ! temps universel en fraction de jour
144	real, intent(in):: pdtrad ! pas de temps du rayonnement (secondes)
145
146	real, intent(out):: pmu0(klon)
147	! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")
148
149	real, intent(out), optional:: frac(klon)
150	! (ensoleillement moyen entre gmtime et gmtime+pdtrad)
151
152	! Variables local to the procedure:
153
154	integer i
155	real gmtime1, gmtime2
156	real deux_pi, incl
157
158	real omega1, omega2, omega
159	! omega1, omega2 : temps 1 et 2 exprimés en radians avec 0 à midi.
160	! omega : heure en radians du coucher de soleil
161	! -omega est donc l'heure en radians de lever du soleil
162
163	real omegadeb, omegafin
164	real zfrac1, zfrac2, z1_mu, z2_mu
165	real lat_sun ! déclinaison en radians
166	real lon_sun ! longitude solaire en radians
167	real latr ! latitude du point de grille en radians
168
169	!----------------------------------------------------------------------
170
171	deux_pi = 2 * pi
172	incl=R_incl * pi / 180.
173
174	lon_sun = longi * pi / 180.
175	lat_sun = ASIN(SIN(lon_sun)*SIN(incl))
176
177	gmtime1=gmtime*86400.
178	gmtime2=gmtime*86400.+pdtrad
179
180	DO i = 1, klon
181	latr = rlat(i) * pi / 180.
182	omega=0.0 !--nuit polaire
183	IF (latr >= (pi/2.-lat_sun) &
184	.OR. latr <= (-pi/2.-lat_sun)) THEN
185	omega = pi ! journee polaire
186	ENDIF
187	IF (latr < (pi/2.+lat_sun).AND. &
188	latr.GT.(-pi/2.+lat_sun).AND. &
189	latr < (pi/2.-lat_sun).AND. &
190	latr.GT.(-pi/2.-lat_sun)) THEN
191	omega = -TAN(latr)*TAN(lat_sun)
192	omega = ACOS(omega)
193	ENDIF
194
195	omega1 = gmtime1 + rlon(i)*86400.0/360.0
196	omega1 = omega1 / 86400.0*deux_pi
197	omega1 = MOD (omega1+deux_pi, deux_pi)
198	omega1 = omega1 - pi
199
200	omega2 = gmtime2 + rlon(i)*86400.0/360.0
201	omega2 = omega2 / 86400.0*deux_pi
202	omega2 = MOD (omega2+deux_pi, deux_pi)
203	omega2 = omega2 - pi
204
205	test_omega12: IF (omega1 <= omega2) THEN
206	! on est dans la meme journee locale
207	IF (omega2 <= -omega .OR. omega1 >= omega .OR. omega < 1e-5) THEN
208	! nuit
209	if (present(frac)) frac(i)=0.0
210	pmu0(i)=0.0
211	ELSE
212	! jour + nuit / jour
213	omegadeb=MAX(-omega, omega1)
214	omegafin=MIN(omega, omega2)
215	if (present(frac)) frac(i)=(omegafin-omegadeb)/(omega2-omega1)
216	pmu0(i)=SIN(latr)*SIN(lat_sun) + &
217	COS(latr)COS(lat_sun) &
218	(SIN(omegafin)-SIN(omegadeb))/ &
219	(omegafin-omegadeb)
220	ENDIF
221	ELSE test_omega12
222	!---omega1 GT omega2 -- a cheval sur deux journees
223	!-------------------entre omega1 et pi
224	IF (omega1 >= omega) THEN !--nuit
225	zfrac1=0.0
226	z1_mu =0.0
227	ELSE !--jour+nuit
228	omegadeb=MAX(-omega, omega1)
229	omegafin=omega
230	zfrac1=omegafin-omegadeb
231	z1_mu =SIN(latr)*SIN(lat_sun) + &
232	COS(latr)COS(lat_sun) &
233	(SIN(omegafin)-SIN(omegadeb))/ &
234	(omegafin-omegadeb)
235	ENDIF
236	!---------------------entre -pi et omega2
237	IF (omega2 <= -omega) THEN !--nuit
238	zfrac2=0.0
239	z2_mu =0.0
240	ELSE !--jour+nuit
241	omegadeb=-omega
242	omegafin=MIN(omega, omega2)
243	zfrac2=omegafin-omegadeb
244	z2_mu =SIN(latr)*SIN(lat_sun) + &
245	COS(latr)COS(lat_sun) &
246	(SIN(omegafin)-SIN(omegadeb))/ &
247	(omegafin-omegadeb)
248
249	ENDIF
250	!-----------------------moyenne
251	if (present(frac)) &
252	frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi-omega1)
253	pmu0(i)=(zfrac1z1_mu+zfrac2z2_mu)/MAX(zfrac1+zfrac2, 1.E-10)
254	ENDIF test_omega12
255	ENDDO
256
257	END SUBROUTINE zenang
258
259	!*****************************************************************
260
261	SUBROUTINE zenith (longi, gmtime, lat, long, pmu0, fract)
262
263	use dimphy, only: klon
264	use YOMCST, only: R_incl
265	use comconst, only: pi
266
267	! Author: Z.X. Li (LMD/ENS)
268
269	! Calcule le cosinus de l'angle zénithal du Soleil en connaissant
270	! la déclinaison du Soleil, la latitude et la longitude du point
271	! sur la Terre, et le temps universel.
272
273	REAL, intent(in):: longi ! déclinaison du soleil (en degrés)
274
275	REAL, intent(in):: gmtime
276	! (temps universel en fraction de jour, entre 0 et 1)
277
278	REAL, intent(in):: lat(klon), long(klon) ! latitude et longitude en degres
279	REAL, intent(out):: pmu0(klon) ! cosinus de l'angle zenithal
280	REAL, intent(out):: fract(klon)
281
282	! Variables local to the procedure:
283
284	INTEGER n
285	REAL zpir, omega
286	REAL lat_sun
287
288	!----------------------------------------------------------------------
289
290	zpir = pi / 180.0
291	lat_sun = ASIN (SIN(longi * zpir)SIN(R_incl zpir) )
292
293	! Initialisation à la nuit:
294	pmu0 = 0.
295	fract = 0.
296
297	! 1 degree in longitude = 240 s in time
298
299	DO n = 1, klon
300	omega = (gmtime + long(n) / 360.) * 2 * pi
301	omega = MOD(omega + 2 * pi, 2 * pi)
302	omega = omega - pi
303	pmu0(n) = sin(lat(n)zpir) sin(lat_sun) &
304	+ cos(lat(n)zpir) cos(lat_sun) &
305	* cos(omega)
306	pmu0(n) = MAX (pmu0(n), 0.)
307	IF (pmu0(n) > 1.E-6) fract(n)=1.0
308	ENDDO
309
310	END SUBROUTINE zenith
311
312	end module orbite_m