| 1 | ! |
|---|
| 2 | ! $Header$ |
|---|
| 3 | ! |
|---|
| 4 | c====================================================================== |
|---|
| 5 | SUBROUTINE orbite(xjour,longi,dist) |
|---|
| 6 | IMPLICIT none |
|---|
| 7 | c====================================================================== |
|---|
| 8 | c Auteur(s): Z.X. Li (LMD/CNRS) (adapte du GCM du LMD) date: 19930818 |
|---|
| 9 | c Objet: pour un jour donne, calculer la longitude vraie de la terre |
|---|
| 10 | c (par rapport au point vernal-21 mars) dans son orbite solaire |
|---|
| 11 | c calculer aussi la distance terre-soleil (unite astronomique) |
|---|
| 12 | c====================================================================== |
|---|
| 13 | c Arguments: |
|---|
| 14 | c xjour--INPUT--R- jour de l'annee a compter du 1er janvier |
|---|
| 15 | c longi--OUTPUT-R- longitude vraie en degres par rapport au point |
|---|
| 16 | c vernal (21 mars) en degres |
|---|
| 17 | c dist---OUTPUT-R- distance terre-soleil (par rapport a la moyenne) |
|---|
| 18 | REAL xjour, longi, dist |
|---|
| 19 | c====================================================================== |
|---|
| 20 | #include "YOMCST.h" |
|---|
| 21 | C |
|---|
| 22 | C -- Variables dynamiques locales |
|---|
| 23 | REAL pir,xl,xllp,xee,xse,xlam,dlamm,anm,ranm,anv,ranv |
|---|
| 24 | C |
|---|
| 25 | pir = 4.0*ATAN(1.0) / 180.0 |
|---|
| 26 | xl=R_peri+180.0 |
|---|
| 27 | xllp=xl*pir |
|---|
| 28 | xee=R_ecc*R_ecc |
|---|
| 29 | xse=SQRT(1.0-xee) |
|---|
| 30 | xlam = (R_ecc/2.0+R_ecc*xee/8.0)*(1.0+xse)*SIN(xllp) |
|---|
| 31 | . - xee/4.0*(0.5+xse)*SIN(2.0*xllp) |
|---|
| 32 | . + R_ecc*xee/8.0*(1.0/3.0+xse)*SIN(3.0*xllp) |
|---|
| 33 | xlam=2.0*xlam/pir |
|---|
| 34 | dlamm=xlam+(xjour-81.0) |
|---|
| 35 | anm=dlamm-xl |
|---|
| 36 | ranm=anm*pir |
|---|
| 37 | xee=xee*R_ecc |
|---|
| 38 | ranv=ranm+(2.0*R_ecc-xee/4.0)*SIN(ranm) |
|---|
| 39 | . +5.0/4.0*R_ecc*R_ecc*SIN(2.0*ranm) |
|---|
| 40 | . +13.0/12.0*xee*SIN(3.0*ranm) |
|---|
| 41 | c |
|---|
| 42 | anv=ranv/pir |
|---|
| 43 | longi=anv+xl |
|---|
| 44 | C |
|---|
| 45 | dist = (1-R_ecc*R_ecc) |
|---|
| 46 | . /(1+R_ecc*COS(pir*(longi-(R_peri+180.0)))) |
|---|
| 47 | RETURN |
|---|
| 48 | END |
|---|
| 49 | c====================================================================== |
|---|
| 50 | SUBROUTINE angle(longi, lati, frac, muzero) |
|---|
| 51 | IMPLICIT none |
|---|
| 52 | c====================================================================== |
|---|
| 53 | c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818 |
|---|
| 54 | c Objet: Calculer la duree d'ensoleillement pour un jour et la hauteur |
|---|
| 55 | c du soleil (cosinus de l'angle zinithal) moyenne sur la journee |
|---|
| 56 | c====================================================================== |
|---|
| 57 | c Arguments: |
|---|
| 58 | c longi----INPUT-R- la longitude vraie de la terre dans son plan |
|---|
| 59 | c solaire a partir de l'equinoxe de printemps (degre) |
|---|
| 60 | c lati-----INPUT-R- la latitude d'un point sur la terre (degre) |
|---|
| 61 | c frac-----OUTPUT-R la duree d'ensoleillement dans la journee divisee |
|---|
| 62 | c par 24 heures (unite en fraction de 0 a 1) |
|---|
| 63 | c muzero---OUTPUT-R la moyenne du cosinus de l'angle zinithal sur |
|---|
| 64 | c la journee (0 a 1) |
|---|
| 65 | c====================================================================== |
|---|
| 66 | #include "dimensions.h" |
|---|
| 67 | #include "dimphy.h" |
|---|
| 68 | REAL longi |
|---|
| 69 | REAL lati(klon), frac(klon), muzero(klon) |
|---|
| 70 | #include "YOMCST.h" |
|---|
| 71 | REAL lat, omega, lon_sun, lat_sun |
|---|
| 72 | REAL pi_local, incl |
|---|
| 73 | INTEGER i |
|---|
| 74 | c |
|---|
| 75 | pi_local = 4.0 * ATAN(1.0) |
|---|
| 76 | incl=R_incl * pi_local / 180. |
|---|
| 77 | c |
|---|
| 78 | lon_sun = longi * pi_local / 180.0 |
|---|
| 79 | lat_sun = ASIN (sin(lon_sun)*SIN(incl) ) |
|---|
| 80 | c |
|---|
| 81 | DO i = 1, klon |
|---|
| 82 | lat = lati(i) * pi_local / 180.0 |
|---|
| 83 | c |
|---|
| 84 | IF ( lat .GE. (pi_local/2.+lat_sun) |
|---|
| 85 | . .OR. lat.LE.(-pi_local/2.+lat_sun)) THEN |
|---|
| 86 | omega = 0.0 ! nuit polaire |
|---|
| 87 | ELSE IF ( lat.GE.(pi_local/2.-lat_sun) |
|---|
| 88 | . .OR. lat.LE.(-pi_local/2.-lat_sun)) THEN |
|---|
| 89 | omega = pi_local ! journee polaire |
|---|
| 90 | ELSE |
|---|
| 91 | omega = -TAN(lat)*TAN(lat_sun) |
|---|
| 92 | omega = ACOS (omega) |
|---|
| 93 | ENDIF |
|---|
| 94 | c |
|---|
| 95 | frac(i) = omega / pi_local |
|---|
| 96 | c |
|---|
| 97 | IF (omega .GT. 0.0) THEN |
|---|
| 98 | muzero(i) = SIN(lat)*SIN(lat_sun) |
|---|
| 99 | . + COS(lat)*COS(lat_sun)*SIN(omega) / omega |
|---|
| 100 | ELSE |
|---|
| 101 | muzero(i) = 0.0 |
|---|
| 102 | ENDIF |
|---|
| 103 | ENDDO |
|---|
| 104 | c |
|---|
| 105 | RETURN |
|---|
| 106 | END |
|---|
| 107 | c==================================================================== |
|---|
| 108 | SUBROUTINE zenang(longi,gmtime,pdtrad,lat,long, |
|---|
| 109 | s pmu0,frac) |
|---|
| 110 | IMPLICIT none |
|---|
| 111 | c============================================================= |
|---|
| 112 | c Auteur : O. Boucher (LMD/CNRS) |
|---|
| 113 | c d'apres les routines zenith et angle de Z.X. Li |
|---|
| 114 | c Objet : calculer les valeurs moyennes du cos de l'angle zenithal |
|---|
| 115 | c et l'ensoleillement moyen entre gmtime1 et gmtime2 |
|---|
| 116 | c connaissant la declinaison, la latitude et la longitude. |
|---|
| 117 | c Rque : Different de la routine angle en ce sens que zenang |
|---|
| 118 | c fournit des moyennes de pmu0 et non des valeurs |
|---|
| 119 | c instantanees, du coup frac prend toutes les valeurs |
|---|
| 120 | c entre 0 et 1. |
|---|
| 121 | c Date : premiere version le 13 decembre 1994 |
|---|
| 122 | c revu pour GCM le 30 septembre 1996 |
|---|
| 123 | c=============================================================== |
|---|
| 124 | c longi----INPUT : la longitude vraie de la terre dans son plan |
|---|
| 125 | c solaire a partir de l'equinoxe de printemps (degre) |
|---|
| 126 | c gmtime---INPUT : temps universel en fraction de jour |
|---|
| 127 | c pdtrad---INPUT : pas de temps du rayonnement (secondes) |
|---|
| 128 | c lat------INPUT : latitude en degres |
|---|
| 129 | c long-----INPUT : longitude en degres |
|---|
| 130 | c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad |
|---|
| 131 | c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad |
|---|
| 132 | c================================================================ |
|---|
| 133 | #include "dimensions.h" |
|---|
| 134 | #include "dimphy.h" |
|---|
| 135 | #include "YOMCST.h" |
|---|
| 136 | c================================================================ |
|---|
| 137 | real longi, gmtime, pdtrad |
|---|
| 138 | real lat(klon), long(klon), pmu0(klon), frac(klon) |
|---|
| 139 | c================================================================ |
|---|
| 140 | integer i |
|---|
| 141 | real gmtime1, gmtime2 |
|---|
| 142 | real pi_local, deux_pi_local, incl |
|---|
| 143 | real omega1, omega2, omega |
|---|
| 144 | c omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi. |
|---|
| 145 | c omega : heure en radian du coucher de soleil |
|---|
| 146 | c -omega est donc l'heure en radian de lever du soleil |
|---|
| 147 | real omegadeb, omegafin |
|---|
| 148 | real zfrac1, zfrac2, z1_mu, z2_mu |
|---|
| 149 | real lat_sun ! declinaison en radian |
|---|
| 150 | real lon_sun ! longitude solaire en radian |
|---|
| 151 | real latr ! latitude du pt de grille en radian |
|---|
| 152 | c================================================================ |
|---|
| 153 | c |
|---|
| 154 | pi_local = 4.0 * ATAN(1.0) |
|---|
| 155 | deux_pi_local = 2.0 * pi_local |
|---|
| 156 | incl=R_incl * pi_local / 180. |
|---|
| 157 | c |
|---|
| 158 | lon_sun = longi * pi_local / 180.0 |
|---|
| 159 | lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) |
|---|
| 160 | c |
|---|
| 161 | gmtime1=gmtime*86400. |
|---|
| 162 | gmtime2=gmtime*86400.+pdtrad |
|---|
| 163 | c |
|---|
| 164 | DO i = 1, klon |
|---|
| 165 | c |
|---|
| 166 | latr = lat(i) * pi_local / 180. |
|---|
| 167 | c |
|---|
| 168 | c--pose probleme quand lat=+/-90 degres |
|---|
| 169 | c |
|---|
| 170 | c omega = -TAN(latr)*TAN(lat_sun) |
|---|
| 171 | c omega = ACOS(omega) |
|---|
| 172 | c IF (latr.GE.(pi_local/2.+lat_sun) |
|---|
| 173 | c . .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN |
|---|
| 174 | c omega = 0.0 ! nuit polaire |
|---|
| 175 | c ENDIF |
|---|
| 176 | c IF (latr.GE.(pi_local/2.-lat_sun) |
|---|
| 177 | c . .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN |
|---|
| 178 | c omega = pi_local ! journee polaire |
|---|
| 179 | c ENDIF |
|---|
| 180 | c |
|---|
| 181 | c--remplace par cela (le cas par defaut est different) |
|---|
| 182 | c |
|---|
| 183 | omega=0.0 !--nuit polaire |
|---|
| 184 | IF (latr.GE.(pi_local/2.-lat_sun) |
|---|
| 185 | . .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN |
|---|
| 186 | omega = pi_local ! journee polaire |
|---|
| 187 | ENDIF |
|---|
| 188 | IF (latr.LT.(pi_local/2.+lat_sun).AND. |
|---|
| 189 | . latr.GT.(-pi_local/2.+lat_sun).AND. |
|---|
| 190 | . latr.LT.(pi_local/2.-lat_sun).AND. |
|---|
| 191 | . latr.GT.(-pi_local/2.-lat_sun)) THEN |
|---|
| 192 | omega = -TAN(latr)*TAN(lat_sun) |
|---|
| 193 | omega = ACOS(omega) |
|---|
| 194 | ENDIF |
|---|
| 195 | c |
|---|
| 196 | omega1 = gmtime1 + long(i)*86400.0/360.0 |
|---|
| 197 | omega1 = omega1 / 86400.0*deux_pi_local |
|---|
| 198 | omega1 = MOD (omega1+deux_pi_local, deux_pi_local) |
|---|
| 199 | omega1 = omega1 - pi_local |
|---|
| 200 | c |
|---|
| 201 | omega2 = gmtime2 + long(i)*86400.0/360.0 |
|---|
| 202 | omega2 = omega2 / 86400.0*deux_pi_local |
|---|
| 203 | omega2 = MOD (omega2+deux_pi_local, deux_pi_local) |
|---|
| 204 | omega2 = omega2 - pi_local |
|---|
| 205 | c |
|---|
| 206 | IF (omega1.LE.omega2) THEN !--on est dans la meme journee locale |
|---|
| 207 | c |
|---|
| 208 | IF (omega2.LE.-omega .OR. omega1.GE.omega |
|---|
| 209 | . .OR. omega.LT.1e-5) THEN !--nuit |
|---|
| 210 | frac(i)=0.0 |
|---|
| 211 | pmu0(i)=0.0 |
|---|
| 212 | ELSE !--jour+nuit/jour |
|---|
| 213 | omegadeb=MAX(-omega,omega1) |
|---|
| 214 | omegafin=MIN(omega,omega2) |
|---|
| 215 | 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 | c |
|---|
| 222 | ELSE !---omega1 GT omega2 -- a cheval sur deux journees |
|---|
| 223 | c |
|---|
| 224 | c-------------------entre omega1 et pi |
|---|
| 225 | IF (omega1.GE.omega) THEN !--nuit |
|---|
| 226 | zfrac1=0.0 |
|---|
| 227 | z1_mu =0.0 |
|---|
| 228 | ELSE !--jour+nuit |
|---|
| 229 | omegadeb=MAX(-omega,omega1) |
|---|
| 230 | omegafin=omega |
|---|
| 231 | zfrac1=omegafin-omegadeb |
|---|
| 232 | z1_mu =SIN(latr)*SIN(lat_sun) + |
|---|
| 233 | . COS(latr)*COS(lat_sun)* |
|---|
| 234 | . (SIN(omegafin)-SIN(omegadeb))/ |
|---|
| 235 | . (omegafin-omegadeb) |
|---|
| 236 | ENDIF |
|---|
| 237 | c---------------------entre -pi et omega2 |
|---|
| 238 | IF (omega2.LE.-omega) THEN !--nuit |
|---|
| 239 | zfrac2=0.0 |
|---|
| 240 | z2_mu =0.0 |
|---|
| 241 | ELSE !--jour+nuit |
|---|
| 242 | omegadeb=-omega |
|---|
| 243 | omegafin=MIN(omega,omega2) |
|---|
| 244 | zfrac2=omegafin-omegadeb |
|---|
| 245 | z2_mu =SIN(latr)*SIN(lat_sun) + |
|---|
| 246 | . COS(latr)*COS(lat_sun)* |
|---|
| 247 | . (SIN(omegafin)-SIN(omegadeb))/ |
|---|
| 248 | . (omegafin-omegadeb) |
|---|
| 249 | c |
|---|
| 250 | ENDIF |
|---|
| 251 | c-----------------------moyenne |
|---|
| 252 | frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1) |
|---|
| 253 | pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10) |
|---|
| 254 | c |
|---|
| 255 | ENDIF !---comparaison omega1 et omega2 |
|---|
| 256 | c |
|---|
| 257 | ENDDO |
|---|
| 258 | c |
|---|
| 259 | END |
|---|
| 260 | c=================================================================== |
|---|
| 261 | SUBROUTINE zenith (longi, gmtime, lat, long, |
|---|
| 262 | s pmu0, fract) |
|---|
| 263 | IMPLICIT none |
|---|
| 264 | c |
|---|
| 265 | c Auteur(s): Z.X. Li (LMD/ENS) |
|---|
| 266 | c |
|---|
| 267 | c Objet: calculer le cosinus de l'angle zenithal du soleil en |
|---|
| 268 | c connaissant la declinaison du soleil, la latitude et la |
|---|
| 269 | c longitude du point sur la terre, et le temps universel |
|---|
| 270 | c |
|---|
| 271 | c Arguments d'entree: |
|---|
| 272 | c longi : declinaison du soleil (en degres) |
|---|
| 273 | c gmtime : temps universel en second qui varie entre 0 et 86400 |
|---|
| 274 | c lat : latitude en degres |
|---|
| 275 | c long : longitude en degres |
|---|
| 276 | c Arguments de sortie: |
|---|
| 277 | c pmu0 : cosinus de l'angle zenithal |
|---|
| 278 | c |
|---|
| 279 | c==================================================================== |
|---|
| 280 | #include "dimensions.h" |
|---|
| 281 | #include "dimphy.h" |
|---|
| 282 | #include "YOMCST.h" |
|---|
| 283 | c==================================================================== |
|---|
| 284 | REAL longi, gmtime |
|---|
| 285 | REAL lat(klon), long(klon), pmu0(klon), fract(klon) |
|---|
| 286 | c===================================================================== |
|---|
| 287 | INTEGER n |
|---|
| 288 | REAL zpi, zpir, omega, zgmtime |
|---|
| 289 | REAL incl, lat_sun, lon_sun |
|---|
| 290 | c---------------------------------------------------------------------- |
|---|
| 291 | zpi = 4.0*ATAN(1.0) |
|---|
| 292 | zpir = zpi / 180.0 |
|---|
| 293 | zgmtime=gmtime*86400. |
|---|
| 294 | c |
|---|
| 295 | incl=R_incl * zpir |
|---|
| 296 | c |
|---|
| 297 | lon_sun = longi * zpir |
|---|
| 298 | lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) |
|---|
| 299 | c |
|---|
| 300 | c--initialisation a la nuit |
|---|
| 301 | c |
|---|
| 302 | DO n =1, klon |
|---|
| 303 | pmu0(n)=0. |
|---|
| 304 | fract(n)=0.0 |
|---|
| 305 | ENDDO |
|---|
| 306 | c |
|---|
| 307 | c 1 degre en longitude = 240 secondes en temps |
|---|
| 308 | c |
|---|
| 309 | DO n = 1, klon |
|---|
| 310 | omega = zgmtime + long(n)*86400.0/360.0 |
|---|
| 311 | omega = omega / 86400.0 * 2.0 * zpi |
|---|
| 312 | omega = MOD(omega + 2.0 * zpi, 2.0 * zpi) |
|---|
| 313 | omega = omega - zpi |
|---|
| 314 | pmu0(n) = sin(lat(n)*zpir) * sin(lat_sun) |
|---|
| 315 | . + cos(lat(n)*zpir) * cos(lat_sun) |
|---|
| 316 | . * cos(omega) |
|---|
| 317 | pmu0(n) = MAX (pmu0(n), 0.0) |
|---|
| 318 | IF (pmu0(n).GT.1.E-6) fract(n)=1.0 |
|---|
| 319 | ENDDO |
|---|
| 320 | c |
|---|
| 321 | RETURN |
|---|
| 322 | END |
|---|
