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