| 1 | MODULE solar |
|---|
| 2 | |
|---|
| 3 | #include "use_logging.h" |
|---|
| 4 | |
|---|
| 5 | IMPLICIT NONE |
|---|
| 6 | PRIVATE |
|---|
| 7 | |
|---|
| 8 | REAL, PARAMETER :: pi_local = 4.0 * ATAN(1.0), deux_pi_local = 2.0 * pi_local |
|---|
| 9 | |
|---|
| 10 | PUBLIC :: solang, zenang, mucorr |
|---|
| 11 | |
|---|
| 12 | CONTAINS |
|---|
| 13 | |
|---|
| 14 | PURE SUBROUTINE solang ( kgrid,psilon,pcolon,psilat,pcolat, & |
|---|
| 15 | ptim1,ptim2,ptim3,pmu0,pfract ) |
|---|
| 16 | |
|---|
| 17 | !----------------------------------------------------------------------- |
|---|
| 18 | ! CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID |
|---|
| 19 | ! |
|---|
| 20 | ! ==== INPUTS === |
|---|
| 21 | ! |
|---|
| 22 | ! PSILON(KGRID) : SINUS OF THE LONGITUDE |
|---|
| 23 | ! PCOLON(KGRID) : COSINUS OF THE LONGITUDE |
|---|
| 24 | ! PSILAT(KGRID) : SINUS OF THE LATITUDE |
|---|
| 25 | ! PCOLAT(KGRID) : COSINUS OF THE LATITUDE |
|---|
| 26 | ! PTIM1 : SIN(DECLI) |
|---|
| 27 | ! PTIM2 : COS(DECLI)*COS(TIME) |
|---|
| 28 | ! PTIM3 : SIN(DECLI)*SIN(TIME) |
|---|
| 29 | ! |
|---|
| 30 | ! ==== OUTPUTS === |
|---|
| 31 | ! |
|---|
| 32 | ! PMU0 (KGRID) : SOLAR ANGLE |
|---|
| 33 | ! PFRACT(KGRID) : DAY FRACTION OF THE TIME INTERVAL |
|---|
| 34 | ! |
|---|
| 35 | ! REFERENCE. |
|---|
| 36 | ! ---------- |
|---|
| 37 | ! |
|---|
| 38 | ! RADIATIVE PROCESSES IN METEOROLOGIE AND CLIMATOLOGIE |
|---|
| 39 | ! PALTRIDGE AND PLATT |
|---|
| 40 | ! |
|---|
| 41 | ! AUTHOR. |
|---|
| 42 | ! ------- |
|---|
| 43 | ! FREDERIC HOURDIN |
|---|
| 44 | ! |
|---|
| 45 | ! MODIFICATIONS. |
|---|
| 46 | ! -------------- |
|---|
| 47 | ! ORIGINAL :90-01-14 |
|---|
| 48 | ! 92-02-14 CALCULATIONS DONE THE ENTIER GRID (J.Polcher) |
|---|
| 49 | !----------------------------------------------------------------------- |
|---|
| 50 | |
|---|
| 51 | INTEGER, INTENT(IN) :: kgrid |
|---|
| 52 | REAL, INTENT(IN) :: ptim1,ptim2,ptim3 |
|---|
| 53 | REAL, INTENT(IN) :: psilon(kgrid), pcolon(kgrid) |
|---|
| 54 | REAL, INTENT(IN) :: psilat(kgrid), pcolat(kgrid) |
|---|
| 55 | REAL, INTENT(OUT) :: pmu0(kgrid), pfract(kgrid) |
|---|
| 56 | |
|---|
| 57 | INTEGER jl |
|---|
| 58 | REAL ztim1,ztim2,ztim3 |
|---|
| 59 | !------------------------------------------------------------------------ |
|---|
| 60 | !------------------------------------------------------------------------ |
|---|
| 61 | !------------------------------------------------------------------------ |
|---|
| 62 | ! |
|---|
| 63 | !------------------------------------------------------------------------ |
|---|
| 64 | ! |
|---|
| 65 | !* 1. INITIALISATION |
|---|
| 66 | ! -------------- |
|---|
| 67 | ! |
|---|
| 68 | ! |
|---|
| 69 | DO jl=1,kgrid |
|---|
| 70 | pmu0(jl)=0. |
|---|
| 71 | pfract(jl)=0. |
|---|
| 72 | ENDDO |
|---|
| 73 | ! |
|---|
| 74 | !* 1.1 COMPUTATION OF THE SOLAR ANGLE |
|---|
| 75 | ! ------------------------------ |
|---|
| 76 | ! |
|---|
| 77 | DO jl=1,kgrid |
|---|
| 78 | ztim1=psilat(jl)*ptim1 |
|---|
| 79 | ztim2=pcolat(jl)*ptim2 |
|---|
| 80 | ztim3=pcolat(jl)*ptim3 |
|---|
| 81 | pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl) |
|---|
| 82 | ENDDO |
|---|
| 83 | ! |
|---|
| 84 | !* 1.2 DISTINCTION BETWEEN DAY AND NIGHT |
|---|
| 85 | ! --------------------------------- |
|---|
| 86 | ! |
|---|
| 87 | DO jl=1,kgrid |
|---|
| 88 | IF (pmu0(jl).gt.0.) THEN |
|---|
| 89 | pfract(jl)=1. |
|---|
| 90 | ELSE |
|---|
| 91 | pmu0(jl)=0. |
|---|
| 92 | pfract(jl)=0. |
|---|
| 93 | ENDIF |
|---|
| 94 | ENDDO |
|---|
| 95 | |
|---|
| 96 | END SUBROUTINE solang |
|---|
| 97 | |
|---|
| 98 | SUBROUTINE zenang(klon,longi,gmtime,pdtrad,lat,long, & |
|---|
| 99 | & pmu0,frac) |
|---|
| 100 | USE astronomy |
|---|
| 101 | |
|---|
| 102 | !============================================================= |
|---|
| 103 | ! Auteur : O. Boucher (LMD/CNRS) |
|---|
| 104 | ! d'apres les routines zenith et angle de Z.X. Li |
|---|
| 105 | ! Objet : calculer les valeurs moyennes du cos de l'angle zenithal |
|---|
| 106 | ! et l'ensoleillement moyen entre gmtime1 et gmtime2 |
|---|
| 107 | ! connaissant la declinaison, la latitude et la longitude. |
|---|
| 108 | ! Rque : Different de la routine angle en ce sens que zenang |
|---|
| 109 | ! fournit des moyennes de pmu0 et non des valeurs |
|---|
| 110 | ! instantanees, du coup frac prend toutes les valeurs |
|---|
| 111 | ! entre 0 et 1. |
|---|
| 112 | ! Date : premiere version le 13 decembre 1994 |
|---|
| 113 | ! revu pour GCM le 30 septembre 1996 |
|---|
| 114 | !=============================================================== |
|---|
| 115 | ! longi----INPUT : la longitude vraie de la terre dans son plan |
|---|
| 116 | ! solaire a partir de l'equinoxe de printemps (degre) |
|---|
| 117 | ! gmtime---INPUT : temps universel en fraction de jour |
|---|
| 118 | ! pdtrad---INPUT : pas de temps du rayonnement (secondes) |
|---|
| 119 | ! lat------INPUT : latitude en degres |
|---|
| 120 | ! long-----INPUT : longitude en degres |
|---|
| 121 | ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad |
|---|
| 122 | ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad |
|---|
| 123 | !================================================================ |
|---|
| 124 | integer klon |
|---|
| 125 | !================================================================ |
|---|
| 126 | real longi, gmtime, pdtrad |
|---|
| 127 | real lat(klon), long(klon), pmu0(klon), frac(klon) |
|---|
| 128 | !================================================================ |
|---|
| 129 | integer i |
|---|
| 130 | real gmtime1, gmtime2 |
|---|
| 131 | real incl |
|---|
| 132 | real omega1, omega2, omega |
|---|
| 133 | ! omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi. |
|---|
| 134 | ! omega : heure en radian du coucher de soleil |
|---|
| 135 | ! -omega est donc l'heure en radian de lever du soleil |
|---|
| 136 | real omegadeb, omegafin |
|---|
| 137 | real zfrac1, zfrac2, z1_mu, z2_mu |
|---|
| 138 | ! declinaison en radian |
|---|
| 139 | real lat_sun |
|---|
| 140 | ! longitude solaire en radian |
|---|
| 141 | real lon_sun |
|---|
| 142 | ! latitude du pt de grille en radian |
|---|
| 143 | real latr |
|---|
| 144 | !================================================================ |
|---|
| 145 | ! |
|---|
| 146 | ! incl=R_incl * pi_local / 180. |
|---|
| 147 | WRITELOG(*,*) 'Obliquite =' ,obliquit |
|---|
| 148 | LOG_INFO('solar') |
|---|
| 149 | |
|---|
| 150 | incl=obliquit * pi_local / 180. |
|---|
| 151 | ! |
|---|
| 152 | ! lon_sun = longi * pi_local / 180.0 |
|---|
| 153 | lon_sun = longi |
|---|
| 154 | lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) |
|---|
| 155 | ! |
|---|
| 156 | gmtime1=gmtime*86400. |
|---|
| 157 | gmtime2=gmtime*86400.+pdtrad |
|---|
| 158 | ! |
|---|
| 159 | DO i = 1, klon |
|---|
| 160 | ! |
|---|
| 161 | ! latr = lat(i) * pi_local / 180. |
|---|
| 162 | latr = lat(i) |
|---|
| 163 | ! |
|---|
| 164 | !--pose probleme quand lat=+/-90 degres |
|---|
| 165 | ! |
|---|
| 166 | ! omega = -TAN(latr)*TAN(lat_sun) |
|---|
| 167 | ! omega = ACOS(omega) |
|---|
| 168 | ! IF (latr.GE.(pi_local/2.+lat_sun) |
|---|
| 169 | ! . .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN |
|---|
| 170 | ! omega = 0.0 ! nuit polaire |
|---|
| 171 | ! ENDIF |
|---|
| 172 | ! IF (latr.GE.(pi_local/2.-lat_sun) |
|---|
| 173 | ! . .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN |
|---|
| 174 | ! omega = pi_local ! journee polaire |
|---|
| 175 | ! ENDIF |
|---|
| 176 | ! |
|---|
| 177 | !--remplace par cela (le cas par defaut est different) |
|---|
| 178 | ! |
|---|
| 179 | !--nuit polaire |
|---|
| 180 | omega=0.0 |
|---|
| 181 | IF (latr.GE.(pi_local/2.-lat_sun) & |
|---|
| 182 | & .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN |
|---|
| 183 | ! journee polaire |
|---|
| 184 | omega = pi_local |
|---|
| 185 | ENDIF |
|---|
| 186 | IF (latr.LT.(pi_local/2.+lat_sun).AND. & |
|---|
| 187 | & latr.GT.(-pi_local/2.+lat_sun).AND. & |
|---|
| 188 | & latr.LT.(pi_local/2.-lat_sun).AND. & |
|---|
| 189 | & latr.GT.(-pi_local/2.-lat_sun)) THEN |
|---|
| 190 | omega = -TAN(latr)*TAN(lat_sun) |
|---|
| 191 | omega = ACOS(omega) |
|---|
| 192 | ENDIF |
|---|
| 193 | ! |
|---|
| 194 | omega1 = gmtime1 + long(i)*86400.0/360.0 |
|---|
| 195 | omega1 = omega1 / 86400.0*deux_pi_local |
|---|
| 196 | omega1 = MOD (omega1+deux_pi_local, deux_pi_local) |
|---|
| 197 | omega1 = omega1 - pi_local |
|---|
| 198 | ! |
|---|
| 199 | omega2 = gmtime2 + long(i)*86400.0/360.0 |
|---|
| 200 | omega2 = omega2 / 86400.0*deux_pi_local |
|---|
| 201 | omega2 = MOD (omega2+deux_pi_local, deux_pi_local) |
|---|
| 202 | omega2 = omega2 - pi_local |
|---|
| 203 | ! |
|---|
| 204 | !--on est dans la meme journee locale |
|---|
| 205 | IF (omega1.LE.omega2) THEN |
|---|
| 206 | ! |
|---|
| 207 | IF (omega2.LE.-omega .OR. omega1.GE.omega .OR. omega.LT.1e-5) & |
|---|
| 208 | THEN |
|---|
| 209 | !--nuit |
|---|
| 210 | frac(i)=0.0 |
|---|
| 211 | pmu0(i)=0.0 |
|---|
| 212 | !--jour+nuit/jou |
|---|
| 213 | ELSE |
|---|
| 214 | omegadeb=MAX(-omega,omega1) |
|---|
| 215 | omegafin=MIN(omega,omega2) |
|---|
| 216 | frac(i)=(omegafin-omegadeb)/(omega2-omega1) |
|---|
| 217 | pmu0(i)=SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)* & |
|---|
| 218 | (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb) |
|---|
| 219 | ENDIF |
|---|
| 220 | ! |
|---|
| 221 | !---omega1 GT omega2 -- a cheval sur deux journees |
|---|
| 222 | ELSE |
|---|
| 223 | ! |
|---|
| 224 | !-------------------entre omega1 et pi |
|---|
| 225 | !--nuit |
|---|
| 226 | IF (omega1.GE.omega) THEN |
|---|
| 227 | zfrac1=0.0 |
|---|
| 228 | z1_mu =0.0 |
|---|
| 229 | !--jour+nuit |
|---|
| 230 | ELSE |
|---|
| 231 | omegadeb=MAX(-omega,omega1) |
|---|
| 232 | omegafin=omega |
|---|
| 233 | zfrac1=omegafin-omegadeb |
|---|
| 234 | z1_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)* & |
|---|
| 235 | (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb) |
|---|
| 236 | ENDIF |
|---|
| 237 | !---------------------entre -pi et omega2 |
|---|
| 238 | !--nuit |
|---|
| 239 | IF (omega2.LE.-omega) THEN |
|---|
| 240 | zfrac2=0.0 |
|---|
| 241 | z2_mu =0.0 |
|---|
| 242 | !--jour+nuit |
|---|
| 243 | ELSE |
|---|
| 244 | omegadeb=-omega |
|---|
| 245 | omegafin=MIN(omega,omega2) |
|---|
| 246 | zfrac2=omegafin-omegadeb |
|---|
| 247 | z2_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)* & |
|---|
| 248 | (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb) |
|---|
| 249 | ! |
|---|
| 250 | ENDIF |
|---|
| 251 | !-----------------------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 | ! |
|---|
| 255 | !---comparaison omega1 et omega2 |
|---|
| 256 | ENDIF |
|---|
| 257 | ! |
|---|
| 258 | ENDDO |
|---|
| 259 | ! |
|---|
| 260 | END SUBROUTINE zenang |
|---|
| 261 | |
|---|
| 262 | PURE SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad) |
|---|
| 263 | |
|---|
| 264 | !======================================================================= |
|---|
| 265 | ! |
|---|
| 266 | ! Calcul of equivalent solar angle and and fraction of day whithout |
|---|
| 267 | ! diurnal cycle. |
|---|
| 268 | ! |
|---|
| 269 | ! parmeters : |
|---|
| 270 | ! ----------- |
|---|
| 271 | ! |
|---|
| 272 | ! Input : |
|---|
| 273 | ! ------- |
|---|
| 274 | ! npts number of points |
|---|
| 275 | ! pdeclin solar declinaison |
|---|
| 276 | ! plat(npts) latitude |
|---|
| 277 | ! phaut hauteur typique de l'atmosphere |
|---|
| 278 | ! prad rayon planetaire |
|---|
| 279 | ! |
|---|
| 280 | ! Output : |
|---|
| 281 | ! -------- |
|---|
| 282 | ! pmu(npts) equivalent cosinus of the solar angle |
|---|
| 283 | ! pfract(npts) fractionnal day |
|---|
| 284 | ! |
|---|
| 285 | !======================================================================= |
|---|
| 286 | |
|---|
| 287 | !----------------------------------------------------------------------- |
|---|
| 288 | ! |
|---|
| 289 | ! 0. Declarations : |
|---|
| 290 | ! ----------------- |
|---|
| 291 | |
|---|
| 292 | ! Arguments : |
|---|
| 293 | ! ----------- |
|---|
| 294 | INTEGER, INTENT(IN) :: npts |
|---|
| 295 | REAL, INTENT(IN) :: phaut, prad, pdeclin, plat(npts) |
|---|
| 296 | REAL, INTENT(OUT) :: pmu(npts), pfract(npts) |
|---|
| 297 | ! |
|---|
| 298 | ! Local variables : |
|---|
| 299 | ! ----------------- |
|---|
| 300 | INTEGER j |
|---|
| 301 | REAL z,cz,sz,tz,phi,cphi,sphi,tphi |
|---|
| 302 | REAL ap,a,t,b |
|---|
| 303 | REAL alph |
|---|
| 304 | |
|---|
| 305 | REAL, PARAMETER :: pi=2.*asin(1.) |
|---|
| 306 | |
|---|
| 307 | !----------------------------------------------------------------------- |
|---|
| 308 | |
|---|
| 309 | z = pdeclin |
|---|
| 310 | cz = cos (z) |
|---|
| 311 | sz = sin (z) |
|---|
| 312 | |
|---|
| 313 | DO j = 1, npts |
|---|
| 314 | |
|---|
| 315 | phi = plat(j) |
|---|
| 316 | cphi = cos(phi) |
|---|
| 317 | if (cphi.le.1.e-9) cphi=1.e-9 |
|---|
| 318 | sphi = sin(phi) |
|---|
| 319 | tphi = sphi / cphi |
|---|
| 320 | b = cphi * cz |
|---|
| 321 | t = -tphi * sz / cz |
|---|
| 322 | a = 1.0 - t*t |
|---|
| 323 | ap = a |
|---|
| 324 | |
|---|
| 325 | IF(t.eq.0.) then |
|---|
| 326 | t=0.5*pi |
|---|
| 327 | ELSE |
|---|
| 328 | IF (a.lt.0.) a = 0. |
|---|
| 329 | t = sqrt(a) / t |
|---|
| 330 | IF (t.lt.0.) then |
|---|
| 331 | t = -atan (-t) + pi |
|---|
| 332 | ELSE |
|---|
| 333 | t = atan(t) |
|---|
| 334 | ENDIF |
|---|
| 335 | ENDIF |
|---|
| 336 | |
|---|
| 337 | pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi |
|---|
| 338 | pfract(j) = t / pi |
|---|
| 339 | IF (ap .lt.0.) then |
|---|
| 340 | pmu(j) = sphi * sz |
|---|
| 341 | pfract(j) = 1.0 |
|---|
| 342 | ENDIF |
|---|
| 343 | IF (pmu(j).le.0.0) pmu(j) = 0.0 |
|---|
| 344 | pmu(j) = pmu(j) / pfract(j) |
|---|
| 345 | IF (pmu(j).eq.0.) pfract(j) = 0. |
|---|
| 346 | |
|---|
| 347 | END DO |
|---|
| 348 | |
|---|
| 349 | !----------------------------------------------------------------------- |
|---|
| 350 | ! correction de rotondite: |
|---|
| 351 | ! ------------------------ |
|---|
| 352 | |
|---|
| 353 | alph=phaut/prad |
|---|
| 354 | DO j=1,npts |
|---|
| 355 | ! !!!!!! |
|---|
| 356 | pmu(j)=sqrt(1224.*pmu(j)*pmu(j)+1.)/35. |
|---|
| 357 | ! $ (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j)) |
|---|
| 358 | END DO |
|---|
| 359 | |
|---|
| 360 | END SUBROUTINE mucorr |
|---|
| 361 | |
|---|
| 362 | END MODULE solar |
|---|
