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