[3] | 1 | c==================================================================== |
---|
| 2 | SUBROUTINE zenang(longi,gmtime,pdtrad,lat,long, |
---|
| 3 | s pmu0,frac) |
---|
| 4 | c============================================================= |
---|
| 5 | c Auteur : O. Boucher (LMD/CNRS) |
---|
| 6 | c d'apres les routines zenith et angle de Z.X. Li |
---|
| 7 | c Objet : calculer les valeurs moyennes du cos de l'angle zenithal |
---|
| 8 | c et l'ensoleillement moyen entre gmtime1 et gmtime2 |
---|
| 9 | c connaissant la declinaison, la latitude et la longitude. |
---|
| 10 | c Rque : Different de la routine angle en ce sens que zenang |
---|
| 11 | c fournit des moyennes de pmu0 et non des valeurs |
---|
| 12 | c instantanees, du coup frac prend toutes les valeurs |
---|
| 13 | c entre 0 et 1. |
---|
| 14 | c Date : premiere version le 13 decembre 1994 |
---|
| 15 | c revu pour GCM le 30 septembre 1996 |
---|
| 16 | c=============================================================== |
---|
| 17 | c longi----INPUT : la longitude vraie de la terre dans son plan |
---|
| 18 | c solaire a partir de l'equinoxe de printemps (degre) |
---|
| 19 | c gmtime---INPUT : temps universel en fraction de jour |
---|
| 20 | c pdtrad---INPUT : pas de temps du rayonnement (secondes) |
---|
| 21 | c lat------INPUT : latitude en degres |
---|
| 22 | c long-----INPUT : longitude en degres |
---|
| 23 | c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad |
---|
| 24 | c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad |
---|
| 25 | c================================================================ |
---|
[101] | 26 | use dimphy |
---|
| 27 | IMPLICIT none |
---|
[3] | 28 | #include "dimensions.h" |
---|
| 29 | #include "YOMCST.h" |
---|
| 30 | c================================================================ |
---|
| 31 | real longi, gmtime, pdtrad |
---|
| 32 | real lat(klon), long(klon), pmu0(klon), frac(klon) |
---|
| 33 | c================================================================ |
---|
| 34 | integer i |
---|
| 35 | real gmtime1, gmtime2 |
---|
[892] | 36 | real incl |
---|
[3] | 37 | real omega1, omega2, omega |
---|
| 38 | c omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi. |
---|
| 39 | c omega : heure en radian du coucher de soleil |
---|
| 40 | c -omega est donc l'heure en radian de lever du soleil |
---|
| 41 | real omegadeb, omegafin |
---|
| 42 | real zfrac1, zfrac2, z1_mu, z2_mu |
---|
| 43 | real lat_sun ! declinaison en radian |
---|
| 44 | real lon_sun ! longitude solaire en radian |
---|
| 45 | real latr ! latitude du pt de grille en radian |
---|
| 46 | c================================================================ |
---|
| 47 | c |
---|
[892] | 48 | incl=R_incl * RPI / 180. |
---|
[3] | 49 | c |
---|
[892] | 50 | lon_sun = longi * RPI / 180.0 |
---|
[3] | 51 | lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) |
---|
| 52 | c |
---|
| 53 | gmtime1=gmtime*RDAY |
---|
| 54 | gmtime2=gmtime*RDAY+pdtrad |
---|
| 55 | c |
---|
| 56 | DO i = 1, klon |
---|
| 57 | c |
---|
[892] | 58 | latr = lat(i) * RPI / 180. |
---|
[3] | 59 | c |
---|
| 60 | c--pose probleme quand lat=+/-90 degres |
---|
| 61 | c |
---|
| 62 | c omega = -TAN(latr)*TAN(lat_sun) |
---|
| 63 | c omega = ACOS(omega) |
---|
[892] | 64 | c IF (latr.GE.(RPI/2.+lat_sun) |
---|
| 65 | c . .OR. latr.LE.(-RPI/2.+lat_sun)) THEN |
---|
[3] | 66 | c omega = 0.0 ! nuit polaire |
---|
| 67 | c ENDIF |
---|
[892] | 68 | c IF (latr.GE.(RPI/2.-lat_sun) |
---|
| 69 | c . .OR. latr.LE.(-RPI/2.-lat_sun)) THEN |
---|
| 70 | c omega = RPI ! journee polaire |
---|
[3] | 71 | c ENDIF |
---|
| 72 | c |
---|
| 73 | c--remplace par cela (le cas par defaut est different) |
---|
| 74 | c |
---|
| 75 | omega=0.0 !--nuit polaire |
---|
[892] | 76 | IF (latr.GE.(RPI/2.-lat_sun) |
---|
| 77 | . .OR. latr.LE.(-RPI/2.-lat_sun)) THEN |
---|
| 78 | omega = RPI ! journee polaire |
---|
[3] | 79 | ENDIF |
---|
[892] | 80 | IF (latr.LT.(RPI/2.+lat_sun).AND. |
---|
| 81 | . latr.GT.(-RPI/2.+lat_sun).AND. |
---|
| 82 | . latr.LT.(RPI/2.-lat_sun).AND. |
---|
| 83 | . latr.GT.(-RPI/2.-lat_sun)) THEN |
---|
[3] | 84 | omega = -TAN(latr)*TAN(lat_sun) |
---|
| 85 | omega = ACOS(omega) |
---|
| 86 | ENDIF |
---|
| 87 | c |
---|
| 88 | omega1 = gmtime1 + long(i)*RDAY/360.0 |
---|
[892] | 89 | omega1 = omega1*2*RPI / RDAY |
---|
| 90 | omega1 = MOD (omega1+2*RPI, 2*RPI) |
---|
| 91 | omega1 = omega1 - RPI |
---|
[3] | 92 | c |
---|
| 93 | omega2 = gmtime2 + long(i)*RDAY/360.0 |
---|
[892] | 94 | omega2 = omega2*2*RPI / RDAY |
---|
| 95 | omega2 = MOD (omega2+2*RPI, 2*RPI) |
---|
| 96 | omega2 = omega2 - RPI |
---|
[3] | 97 | c |
---|
| 98 | IF (omega1.LE.omega2) THEN !--on est dans la meme journee locale |
---|
| 99 | c |
---|
| 100 | IF (omega2.LE.-omega .OR. omega1.GE.omega |
---|
| 101 | . .OR. omega.LT.1e-5) THEN !--nuit |
---|
| 102 | frac(i)=0.0 |
---|
| 103 | pmu0(i)=0.0 |
---|
| 104 | ELSE !--jour+nuit/jour |
---|
| 105 | omegadeb=MAX(-omega,omega1) |
---|
| 106 | omegafin=MIN(omega,omega2) |
---|
| 107 | frac(i)=(omegafin-omegadeb)/(omega2-omega1) |
---|
| 108 | pmu0(i)=SIN(latr)*SIN(lat_sun) + |
---|
| 109 | . COS(latr)*COS(lat_sun)* |
---|
| 110 | . (SIN(omegafin)-SIN(omegadeb))/ |
---|
| 111 | . (omegafin-omegadeb) |
---|
| 112 | ENDIF |
---|
| 113 | c |
---|
| 114 | ELSE !---omega1 GT omega2 -- a cheval sur deux journees |
---|
| 115 | c |
---|
| 116 | c-------------------entre omega1 et pi |
---|
| 117 | IF (omega1.GE.omega) THEN !--nuit |
---|
| 118 | zfrac1=0.0 |
---|
| 119 | z1_mu =0.0 |
---|
| 120 | ELSE !--jour+nuit |
---|
| 121 | omegadeb=MAX(-omega,omega1) |
---|
| 122 | omegafin=omega |
---|
| 123 | zfrac1=omegafin-omegadeb |
---|
| 124 | z1_mu =SIN(latr)*SIN(lat_sun) + |
---|
| 125 | . COS(latr)*COS(lat_sun)* |
---|
| 126 | . (SIN(omegafin)-SIN(omegadeb))/ |
---|
| 127 | . (omegafin-omegadeb) |
---|
| 128 | ENDIF |
---|
| 129 | c---------------------entre -pi et omega2 |
---|
| 130 | IF (omega2.LE.-omega) THEN !--nuit |
---|
| 131 | zfrac2=0.0 |
---|
| 132 | z2_mu =0.0 |
---|
| 133 | ELSE !--jour+nuit |
---|
| 134 | omegadeb=-omega |
---|
| 135 | omegafin=MIN(omega,omega2) |
---|
| 136 | zfrac2=omegafin-omegadeb |
---|
| 137 | z2_mu =SIN(latr)*SIN(lat_sun) + |
---|
| 138 | . COS(latr)*COS(lat_sun)* |
---|
| 139 | . (SIN(omegafin)-SIN(omegadeb))/ |
---|
| 140 | . (omegafin-omegadeb) |
---|
| 141 | c |
---|
| 142 | ENDIF |
---|
| 143 | c-----------------------moyenne |
---|
[892] | 144 | frac(i)=(zfrac1+zfrac2)/(omega2+2*RPI-omega1) |
---|
[3] | 145 | pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10) |
---|
| 146 | c |
---|
| 147 | ENDIF !---comparaison omega1 et omega2 |
---|
| 148 | c |
---|
| 149 | ENDDO |
---|
| 150 | c |
---|
| 151 | END |
---|