| 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_ecc*xee/8.0)*(1.0+xse)*SIN(xllp) & |
| 38 |
|
|
- xee/4.0*(0.5+xse)*SIN(2.0*xllp) & |
| 39 |
|
|
+ R_ecc*xee/8.0*(1.0/3.0+xse)*SIN(3.0*xllp) |
| 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.0*R_ecc-xee/4.0)*SIN(ranm) & |
| 46 |
|
|
+5.0/4.0*R_ecc*R_ecc*SIN(2.0*ranm) & |
| 47 |
|
|
+13.0/12.0*xee*SIN(3.0*ranm) |
| 48 |
|
|
|
| 49 |
|
|
anv=ranv/pir |
| 50 |
|
|
longi=anv+xl |
| 51 |
|
|
|
| 52 |
|
|
if (present(dist)) dist & |
| 53 |
|
|
= (1-R_ecc*R_ecc) /(1+R_ecc*COS(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 |
|
|
|
| 125 |
|
|
! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et |
| 126 |
|
|
! "angle" de Z.X. Li |
| 127 |
|
|
|
| 128 |
|
|
! Calcule les valeurs moyennes du cos de l'angle zénithal et |
| 129 |
|
|
! l'ensoleillement moyen entre "gmtime1" et "gmtime2" connaissant la |
| 130 |
|
|
! déclinaison, la latitude et la longitude. |
| 131 |
|
|
! Différent de la routine "angle" en ce sens que "zenang" fournit des |
| 132 |
|
|
! moyennes de "pmu0" et non des valeurs instantanées. |
| 133 |
|
|
! Du coup "frac" prend toutes les valeurs entre 0 et 1. |
| 134 |
|
|
|
| 135 |
|
|
! Date : premiere version le 13 decembre 1994 |
| 136 |
|
|
! revu pour GCM le 30 septembre 1996 |
| 137 |
|
|
|
| 138 |
|
|
real, intent(in):: longi |
| 139 |
|
|
! (longitude vraie de la terre dans son plan solaire a partir de |
| 140 |
|
|
! l'equinoxe de printemps (degre)) |
| 141 |
|
|
|
| 142 |
|
|
real, intent(in):: gmtime ! temps universel en fraction de jour |
| 143 |
|
|
real, intent(in):: pdtrad ! pas de temps du rayonnement (secondes) |
| 144 |
|
|
|
| 145 |
|
|
real, intent(out):: pmu0(klon) |
| 146 |
|
|
! (cosinus de l'angle zenithal moyen entre gmtime et gmtime+pdtrad) |
| 147 |
|
|
|
| 148 |
|
|
real, intent(out), optional:: frac(klon) |
| 149 |
|
|
! (ensoleillement moyen entre gmtime et gmtime+pdtrad) |
| 150 |
|
|
|
| 151 |
|
|
! Variables local to the procedure: |
| 152 |
|
|
|
| 153 |
|
|
integer i |
| 154 |
|
|
real gmtime1, gmtime2 |
| 155 |
|
|
real pi_local, deux_pi_local, incl |
| 156 |
|
|
real omega1, omega2, omega |
| 157 |
|
|
! omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi. |
| 158 |
|
|
! omega : heure en radian du coucher de soleil |
| 159 |
|
|
! -omega est donc l'heure en radian de lever du soleil |
| 160 |
|
|
real omegadeb, omegafin |
| 161 |
|
|
real zfrac1, zfrac2, z1_mu, z2_mu |
| 162 |
|
|
real lat_sun ! declinaison en radian |
| 163 |
|
|
real lon_sun ! longitude solaire en radian |
| 164 |
|
|
real latr ! latitude du pt de grille en radian |
| 165 |
|
|
|
| 166 |
|
|
!---------------------------------------------------------------------- |
| 167 |
|
|
|
| 168 |
|
|
pi_local = 4.0 * ATAN(1.0) |
| 169 |
|
|
deux_pi_local = 2.0 * pi_local |
| 170 |
|
|
incl=R_incl * pi_local / 180. |
| 171 |
|
|
|
| 172 |
|
|
lon_sun = longi * pi_local / 180.0 |
| 173 |
|
|
lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) |
| 174 |
|
|
|
| 175 |
|
|
gmtime1=gmtime*86400. |
| 176 |
|
|
gmtime2=gmtime*86400.+pdtrad |
| 177 |
|
|
|
| 178 |
|
|
DO i = 1, klon |
| 179 |
|
|
latr = rlat(i) * pi_local / 180. |
| 180 |
|
|
omega=0.0 !--nuit polaire |
| 181 |
|
|
IF (latr >= (pi_local/2.-lat_sun) & |
| 182 |
|
|
.OR. latr <= (-pi_local/2.-lat_sun)) THEN |
| 183 |
|
|
omega = pi_local ! journee polaire |
| 184 |
|
|
ENDIF |
| 185 |
|
|
IF (latr < (pi_local/2.+lat_sun).AND. & |
| 186 |
|
|
latr.GT.(-pi_local/2.+lat_sun).AND. & |
| 187 |
|
|
latr < (pi_local/2.-lat_sun).AND. & |
| 188 |
|
|
latr.GT.(-pi_local/2.-lat_sun)) THEN |
| 189 |
|
|
omega = -TAN(latr)*TAN(lat_sun) |
| 190 |
|
|
omega = ACOS(omega) |
| 191 |
|
|
ENDIF |
| 192 |
|
|
|
| 193 |
|
|
omega1 = gmtime1 + rlon(i)*86400.0/360.0 |
| 194 |
|
|
omega1 = omega1 / 86400.0*deux_pi_local |
| 195 |
|
|
omega1 = MOD (omega1+deux_pi_local, deux_pi_local) |
| 196 |
|
|
omega1 = omega1 - pi_local |
| 197 |
|
|
|
| 198 |
|
|
omega2 = gmtime2 + rlon(i)*86400.0/360.0 |
| 199 |
|
|
omega2 = omega2 / 86400.0*deux_pi_local |
| 200 |
|
|
omega2 = MOD (omega2+deux_pi_local, deux_pi_local) |
| 201 |
|
|
omega2 = omega2 - pi_local |
| 202 |
|
|
|
| 203 |
|
|
test_omega12: IF (omega1 <= omega2) THEN |
| 204 |
|
|
! on est dans la meme journee locale |
| 205 |
|
|
IF (omega2 <= -omega .OR. omega1 >= omega .OR. omega < 1e-5) THEN |
| 206 |
|
|
! nuit |
| 207 |
|
|
if (present(frac)) frac(i)=0.0 |
| 208 |
|
|
pmu0(i)=0.0 |
| 209 |
|
|
ELSE |
| 210 |
|
|
! jour + nuit / jour |
| 211 |
|
|
omegadeb=MAX(-omega, omega1) |
| 212 |
|
|
omegafin=MIN(omega, omega2) |
| 213 |
|
|
if (present(frac)) frac(i)=(omegafin-omegadeb)/(omega2-omega1) |
| 214 |
|
|
pmu0(i)=SIN(latr)*SIN(lat_sun) + & |
| 215 |
|
|
COS(latr)*COS(lat_sun)* & |
| 216 |
|
|
(SIN(omegafin)-SIN(omegadeb))/ & |
| 217 |
|
|
(omegafin-omegadeb) |
| 218 |
|
|
ENDIF |
| 219 |
|
|
ELSE test_omega12 |
| 220 |
|
|
!---omega1 GT omega2 -- a cheval sur deux journees |
| 221 |
|
|
!-------------------entre omega1 et pi |
| 222 |
|
|
IF (omega1 >= omega) THEN !--nuit |
| 223 |
|
|
zfrac1=0.0 |
| 224 |
|
|
z1_mu =0.0 |
| 225 |
|
|
ELSE !--jour+nuit |
| 226 |
|
|
omegadeb=MAX(-omega, omega1) |
| 227 |
|
|
omegafin=omega |
| 228 |
|
|
zfrac1=omegafin-omegadeb |
| 229 |
|
|
z1_mu =SIN(latr)*SIN(lat_sun) + & |
| 230 |
|
|
COS(latr)*COS(lat_sun)* & |
| 231 |
|
|
(SIN(omegafin)-SIN(omegadeb))/ & |
| 232 |
|
|
(omegafin-omegadeb) |
| 233 |
|
|
ENDIF |
| 234 |
|
|
!---------------------entre -pi et omega2 |
| 235 |
|
|
IF (omega2 <= -omega) THEN !--nuit |
| 236 |
|
|
zfrac2=0.0 |
| 237 |
|
|
z2_mu =0.0 |
| 238 |
|
|
ELSE !--jour+nuit |
| 239 |
|
|
omegadeb=-omega |
| 240 |
|
|
omegafin=MIN(omega, omega2) |
| 241 |
|
|
zfrac2=omegafin-omegadeb |
| 242 |
|
|
z2_mu =SIN(latr)*SIN(lat_sun) + & |
| 243 |
|
|
COS(latr)*COS(lat_sun)* & |
| 244 |
|
|
(SIN(omegafin)-SIN(omegadeb))/ & |
| 245 |
|
|
(omegafin-omegadeb) |
| 246 |
|
|
|
| 247 |
|
|
ENDIF |
| 248 |
|
|
!-----------------------moyenne |
| 249 |
|
|
if (present(frac)) & |
| 250 |
|
|
frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1) |
| 251 |
|
|
pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2, 1.E-10) |
| 252 |
|
|
ENDIF test_omega12 |
| 253 |
|
|
ENDDO |
| 254 |
|
|
|
| 255 |
|
|
END SUBROUTINE zenang |
| 256 |
|
|
|
| 257 |
|
|
!***************************************************************** |
| 258 |
|
|
|
| 259 |
|
|
SUBROUTINE zenith (longi, gmtime, lat, long, pmu0, fract) |
| 260 |
|
|
|
| 261 |
|
|
use dimphy, only: klon |
| 262 |
|
|
use YOMCST, only: R_incl |
| 263 |
|
|
use comconst, only: pi |
| 264 |
|
|
|
| 265 |
|
|
! Author: Z.X. Li (LMD/ENS) |
| 266 |
|
|
|
| 267 |
|
|
! Calcule le cosinus de l'angle zénithal du Soleil en connaissant |
| 268 |
|
|
! la déclinaison du Soleil, la latitude et la longitude du point |
| 269 |
|
|
! sur la Terre, et le temps universel. |
| 270 |
|
|
|
| 271 |
|
|
REAL, intent(in):: longi ! déclinaison du soleil (en degrés) |
| 272 |
|
|
|
| 273 |
|
|
REAL, intent(in):: gmtime |
| 274 |
|
|
! (temps universel en fraction de jour, entre 0 et 1) |
| 275 |
|
|
|
| 276 |
|
|
REAL, intent(in):: lat(klon), long(klon) ! latitude et longitude en degres |
| 277 |
|
|
REAL, intent(out):: pmu0(klon) ! cosinus de l'angle zenithal |
| 278 |
|
|
REAL, intent(out):: fract(klon) |
| 279 |
|
|
|
| 280 |
|
|
! Variables local to the procedure: |
| 281 |
|
|
|
| 282 |
|
|
INTEGER n |
| 283 |
|
|
REAL zpir, omega |
| 284 |
|
|
REAL lat_sun |
| 285 |
|
|
|
| 286 |
|
|
!---------------------------------------------------------------------- |
| 287 |
|
|
|
| 288 |
|
|
zpir = pi / 180.0 |
| 289 |
|
|
lat_sun = ASIN (SIN(longi * zpir)*SIN(R_incl * zpir) ) |
| 290 |
|
|
|
| 291 |
|
|
! Initialisation à la nuit: |
| 292 |
|
|
pmu0 = 0. |
| 293 |
|
|
fract = 0. |
| 294 |
|
|
|
| 295 |
|
|
! 1 degree in longitude = 240 s in time |
| 296 |
|
|
|
| 297 |
|
|
DO n = 1, klon |
| 298 |
|
|
omega = (gmtime + long(n) / 360.) * 2 * pi |
| 299 |
|
|
omega = MOD(omega + 2 * pi, 2 * pi) |
| 300 |
|
|
omega = omega - pi |
| 301 |
|
|
pmu0(n) = sin(lat(n)*zpir) * sin(lat_sun) & |
| 302 |
|
|
+ cos(lat(n)*zpir) * cos(lat_sun) & |
| 303 |
|
|
* cos(omega) |
| 304 |
|
|
pmu0(n) = MAX (pmu0(n), 0.) |
| 305 |
|
|
IF (pmu0(n) > 1.E-6) fract(n)=1.0 |
| 306 |
|
|
ENDDO |
| 307 |
|
|
|
| 308 |
|
|
END SUBROUTINE zenith |
| 309 |
|
|
|
| 310 |
|
|
end module orbite_m |
Parent Directory
| 
Revision Log