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

    ! 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	guez	3	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	guez	20	use comconst, only: pi
125	guez	3
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	guez	20	! Date : première version le 13 decembre 1994
137			! revu pour GCM le 30 septembre 1996
138	guez	3
139			real, intent(in):: longi
140			! (longitude vraie de la terre dans son plan solaire a partir de
141	guez	20	! l'equinoxe de printemps) (in degrees)
142	guez	3
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	guez	20	! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")
148	guez	3
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	guez	20	real deux_pi, incl
157
158	guez	3	real omega1, omega2, omega
159	guez	20	! 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	guez	3	real omegadeb, omegafin
164			real zfrac1, zfrac2, z1_mu, z2_mu
165	guez	20	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	guez	3
169			!----------------------------------------------------------------------
170
171	guez	20	deux_pi = 2 * pi
172			incl=R_incl * pi / 180.
173	guez	3
174	guez	20	lon_sun = longi * pi / 180.
175			lat_sun = ASIN(SIN(lon_sun)*SIN(incl))
176	guez	3
177			gmtime1=gmtime*86400.
178			gmtime2=gmtime*86400.+pdtrad
179
180			DO i = 1, klon
181	guez	20	latr = rlat(i) * pi / 180.
182	guez	3	omega=0.0 !--nuit polaire
183	guez	20	IF (latr >= (pi/2.-lat_sun) &
184			.OR. latr <= (-pi/2.-lat_sun)) THEN
185			omega = pi ! journee polaire
186	guez	3	ENDIF
187	guez	20	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	guez	3	omega = -TAN(latr)*TAN(lat_sun)
192			omega = ACOS(omega)
193			ENDIF
194
195			omega1 = gmtime1 + rlon(i)*86400.0/360.0
196	guez	20	omega1 = omega1 / 86400.0*deux_pi
197			omega1 = MOD (omega1+deux_pi, deux_pi)
198			omega1 = omega1 - pi
199	guez	3
200			omega2 = gmtime2 + rlon(i)*86400.0/360.0
201	guez	20	omega2 = omega2 / 86400.0*deux_pi
202			omega2 = MOD (omega2+deux_pi, deux_pi)
203			omega2 = omega2 - pi
204	guez	3
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	guez	20	frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi-omega1)
253	guez	3	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