| 1 | ! ================================================================================================================================= |
|---|
| 2 | ! MODULE : solar |
|---|
| 3 | ! |
|---|
| 4 | ! CONTACT : orchidee-help _at_ listes.ipsl.fr |
|---|
| 5 | ! |
|---|
| 6 | ! LICENCE : IPSL (2011) |
|---|
| 7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
|---|
| 8 | ! |
|---|
| 9 | !>\BRIEF This module define solar |
|---|
| 10 | !! |
|---|
| 11 | !!\n DESCRIPTION: |
|---|
| 12 | !! |
|---|
| 13 | !! RECENT CHANGE(S): None |
|---|
| 14 | !! |
|---|
| 15 | !! REFERENCE(S) : None |
|---|
| 16 | !! |
|---|
| 17 | !! SVN : |
|---|
| 18 | !! $HeadURL$ |
|---|
| 19 | !! $Date$ |
|---|
| 20 | !! $Revision$ |
|---|
| 21 | !! \n |
|---|
| 22 | !_ ================================================================================================================================ |
|---|
| 23 | |
|---|
| 24 | MODULE solar |
|---|
| 25 | |
|---|
| 26 | USE defprec |
|---|
| 27 | USE constantes |
|---|
| 28 | USE ioipsl_para |
|---|
| 29 | |
|---|
| 30 | IMPLICIT NONE |
|---|
| 31 | |
|---|
| 32 | |
|---|
| 33 | |
|---|
| 34 | CONTAINS |
|---|
| 35 | |
|---|
| 36 | |
|---|
| 37 | !! ================================================================================================================================ |
|---|
| 38 | !! SUBROUTINE : solarang |
|---|
| 39 | !! |
|---|
| 40 | !>\BRIEF This subroutine computes the solar angle according to the method used by GSWP and developed by J.C. Morill. |
|---|
| 41 | !! |
|---|
| 42 | !! DESCRIPTION : None |
|---|
| 43 | !! |
|---|
| 44 | !! RECENT CHANGE(S): None |
|---|
| 45 | !! |
|---|
| 46 | !! MAIN OUTPUT VARIABLE(S): ::csang |
|---|
| 47 | !! |
|---|
| 48 | !! REFERENCE(S) : See for further details : |
|---|
| 49 | !! http://www.atmo.arizona.edu/~morrill/swrad.html |
|---|
| 50 | !! |
|---|
| 51 | !! FLOWCHART : None |
|---|
| 52 | !! \n |
|---|
| 53 | !_ ================================================================================================================================ |
|---|
| 54 | |
|---|
| 55 | SUBROUTINE solarang (julian, julian0, iim, jjm, lon, lat, csang) |
|---|
| 56 | |
|---|
| 57 | USE calendar |
|---|
| 58 | USE time, ONLY : one_year, calendar_str |
|---|
| 59 | |
|---|
| 60 | IMPLICIT NONE |
|---|
| 61 | |
|---|
| 62 | !! 0. Parameters and variables declaration |
|---|
| 63 | |
|---|
| 64 | !! 0.1 Input variables |
|---|
| 65 | |
|---|
| 66 | INTEGER, INTENT(in) :: iim, jjm !! |
|---|
| 67 | REAL, INTENT(in) :: julian !! |
|---|
| 68 | REAL, INTENT(in) :: julian0 !! |
|---|
| 69 | REAL, DIMENSION(iim,jjm), INTENT(in) :: lon, lat !! |
|---|
| 70 | |
|---|
| 71 | !! 0.2 Output variables |
|---|
| 72 | |
|---|
| 73 | REAL, DIMENSION(iim,jjm), INTENT(out) :: csang !! |
|---|
| 74 | |
|---|
| 75 | !! 0.4 Local variables |
|---|
| 76 | |
|---|
| 77 | REAL :: gamma !! |
|---|
| 78 | REAL :: dec !! |
|---|
| 79 | REAL :: decd !! |
|---|
| 80 | REAL :: et !! |
|---|
| 81 | REAL :: gmt !! |
|---|
| 82 | REAL :: le !! |
|---|
| 83 | REAL :: ls !! |
|---|
| 84 | REAL :: lcorr !! |
|---|
| 85 | REAL :: latime !! |
|---|
| 86 | REAL :: omega !! |
|---|
| 87 | REAL :: omegad !! |
|---|
| 88 | REAL :: llatd, llat !! |
|---|
| 89 | INTEGER :: igmt !! |
|---|
| 90 | INTEGER :: ilon, ilat !! |
|---|
| 91 | INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: zone !! |
|---|
| 92 | !$OMP THREADPRIVATE(zone) |
|---|
| 93 | REAL, SAVE, ALLOCATABLE, DIMENSION(:) :: lhour !! |
|---|
| 94 | !$OMP THREADPRIVATE(lhour) |
|---|
| 95 | |
|---|
| 96 | !_ ================================================================================================================================ |
|---|
| 97 | |
|---|
| 98 | IF (printlev >= 4) WRITE(numout,*) 'Calendar information', julian, julian0 |
|---|
| 99 | !- |
|---|
| 100 | !- 1) Day angle gamma |
|---|
| 101 | !- |
|---|
| 102 | ! gamma = 2.*pi*MOD(julian,one_year)/one_year |
|---|
| 103 | gamma = 2.*pi*(julian-julian0)/one_year |
|---|
| 104 | !- |
|---|
| 105 | !- 2) Solar declination (assumed constant for a 24 hour period) page 7 |
|---|
| 106 | !- in radians |
|---|
| 107 | !- |
|---|
| 108 | IF (printlev >= 4) WRITE(numout,*) 'Solar declination', gamma, one_year |
|---|
| 109 | ! |
|---|
| 110 | dec = ( 0.006918-0.399912*COS(gamma)+0.070257*SIN(gamma) & |
|---|
| 111 | & -0.006758*COS(2*gamma)+0.000907*SIN(2*gamma) & |
|---|
| 112 | & -0.002697*COS(3*gamma)+0.00148*SIN(3*gamma)) |
|---|
| 113 | decd = dec*(180/pi) |
|---|
| 114 | !- |
|---|
| 115 | !- maximum error 0.0006 rad (<3'), leads to error |
|---|
| 116 | !- of less than 1/2 degree in zenith angle |
|---|
| 117 | !- |
|---|
| 118 | IF (printlev >= 4) WRITE(numout,*) 'Equation of time', decd |
|---|
| 119 | !- 3) Equation of time page 11 |
|---|
| 120 | !- |
|---|
| 121 | et = ( 0.000075+0.001868*COS(gamma)-0.032077*SIN(gamma)& |
|---|
| 122 | & -0.014615*COS(2*gamma)-0.04089*SIN(2*gamma))*229.18 |
|---|
| 123 | !- |
|---|
| 124 | !- Get the time zones for the current time |
|---|
| 125 | !- |
|---|
| 126 | gmt = 24.*(julian-INT(julian)) |
|---|
| 127 | IF (.NOT.ALLOCATED(zone)) ALLOCATE(zone(iim)) |
|---|
| 128 | IF (.NOT.ALLOCATED(lhour)) ALLOCATE(lhour(iim)) |
|---|
| 129 | !- |
|---|
| 130 | !igmt = NINT(gmt) |
|---|
| 131 | IF (printlev >= 4) WRITE(numout,*) 'Get time zone', et, gmt |
|---|
| 132 | CALL time_zone(gmt, iim, jjm, lon, zone, lhour) |
|---|
| 133 | !- |
|---|
| 134 | !- Start a loop through the grid |
|---|
| 135 | !- |
|---|
| 136 | IF (printlev >= 4) WRITE(numout,*) 'Start a loop through the grid' |
|---|
| 137 | DO ilon=1,iim |
|---|
| 138 | !--- |
|---|
| 139 | !--- 4) Local apparent time page 13 |
|---|
| 140 | !--- |
|---|
| 141 | !--- ls standard longitude (nearest 15 degree meridian) |
|---|
| 142 | !--- le local longtitude |
|---|
| 143 | !--- lhour local standard time |
|---|
| 144 | !--- latime local apparent time |
|---|
| 145 | !--- lcorr longitudunal correction (minutes) |
|---|
| 146 | !--- |
|---|
| 147 | le = lon(ilon,1) |
|---|
| 148 | ls = ((zone(ilon)-1)*15)-180. |
|---|
| 149 | lcorr = 4.*(ls-le)*(-1) |
|---|
| 150 | latime = lhour(ilon)+lcorr/60.+et/60. |
|---|
| 151 | IF (latime < 0.) latime = latime+24 |
|---|
| 152 | IF (latime > 24.) latime = latime-24 |
|---|
| 153 | !--- |
|---|
| 154 | !--- 5) Hour angle omega page 15 |
|---|
| 155 | !--- |
|---|
| 156 | !--- hour angle is zero at noon, positive in the morning |
|---|
| 157 | !--- It ranges from 180 to -180 |
|---|
| 158 | !--- |
|---|
| 159 | !--- omegad is omega in degrees, omega is in radians |
|---|
| 160 | !--- |
|---|
| 161 | omegad = (latime-12.)*(-15.) |
|---|
| 162 | omega = omegad*pi/180. |
|---|
| 163 | !--- |
|---|
| 164 | DO ilat=1,jjm |
|---|
| 165 | !----- |
|---|
| 166 | !----- 6) Zenith angle page 15 |
|---|
| 167 | !----- |
|---|
| 168 | !----- csang cosine of zenith angle (radians) |
|---|
| 169 | !----- llatd = local latitude in degrees |
|---|
| 170 | !----- llat = local latitude in radians |
|---|
| 171 | !----- |
|---|
| 172 | llatd = lat(1,ilat) |
|---|
| 173 | llat = llatd*pi/180. |
|---|
| 174 | csang(ilon,ilat) = & |
|---|
| 175 | & MAX(zero,SIN(dec)*SIN(llat)+COS(dec)*COS(llat)*COS(omega)) |
|---|
| 176 | ENDDO |
|---|
| 177 | ENDDO |
|---|
| 178 | !---------------------- |
|---|
| 179 | END SUBROUTINE solarang |
|---|
| 180 | |
|---|
| 181 | |
|---|
| 182 | !! ================================================================================================================================ |
|---|
| 183 | !! SUBROUTINE : time_zone |
|---|
| 184 | !! |
|---|
| 185 | !>\BRIEF |
|---|
| 186 | !! |
|---|
| 187 | !! DESCRIPTION : None |
|---|
| 188 | !! |
|---|
| 189 | !! RECENT CHANGE(S): None |
|---|
| 190 | !! |
|---|
| 191 | !! MAIN OUTPUT VARIABLE(S): ::zone, ::lhour |
|---|
| 192 | !! |
|---|
| 193 | !! REFERENCE(S) : None |
|---|
| 194 | !! |
|---|
| 195 | !! FLOWCHART : None |
|---|
| 196 | !! \n |
|---|
| 197 | !_ ================================================================================================================================ |
|---|
| 198 | |
|---|
| 199 | SUBROUTINE time_zone (gmt, iim, jjm, lon, zone, lhour) |
|---|
| 200 | |
|---|
| 201 | IMPLICIT NONE |
|---|
| 202 | |
|---|
| 203 | !! 0. Parameters and variables declaration |
|---|
| 204 | |
|---|
| 205 | !! 0.1 Input variables |
|---|
| 206 | |
|---|
| 207 | INTEGER, INTENT(in) :: iim, jjm !! |
|---|
| 208 | REAL, DIMENSION(iim,jjm), INTENT(in) :: lon !! |
|---|
| 209 | REAL, INTENT(in) :: gmt !! |
|---|
| 210 | |
|---|
| 211 | !! 0.2 Output variables |
|---|
| 212 | |
|---|
| 213 | INTEGER, DIMENSION(iim), INTENT(out) :: zone !! |
|---|
| 214 | REAL, DIMENSION(iim), INTENT(out) :: lhour !! |
|---|
| 215 | |
|---|
| 216 | !! 0.4 Local variables |
|---|
| 217 | |
|---|
| 218 | INTEGER :: deg !! |
|---|
| 219 | !!?? REAL :: deg |
|---|
| 220 | INTEGER :: i, ilon, change !! Indices |
|---|
| 221 | |
|---|
| 222 | !_ ================================================================================================================================ |
|---|
| 223 | |
|---|
| 224 | DO ilon=1,iim |
|---|
| 225 | !--- |
|---|
| 226 | !--- Convert longitude index to longtitude in degrees |
|---|
| 227 | !--- |
|---|
| 228 | deg = lon(ilon,1) |
|---|
| 229 | !--- |
|---|
| 230 | !--- Determine into which time zone (15 degree interval) the |
|---|
| 231 | !--- longitude falls. |
|---|
| 232 | !--- |
|---|
| 233 | DO i=1,25 |
|---|
| 234 | IF (deg < (-187.5+(15*i))) THEN |
|---|
| 235 | zone(ilon) = i |
|---|
| 236 | IF (zone(ilon) == 25) zone(ilon) = 1 |
|---|
| 237 | EXIT |
|---|
| 238 | ENDIF |
|---|
| 239 | ENDDO |
|---|
| 240 | !--- |
|---|
| 241 | !--- Calculate change (in number of hours) from GMT time to |
|---|
| 242 | !--- local hour. Change will be negative for zones < 13 and |
|---|
| 243 | !--- positive for zones > 13. |
|---|
| 244 | !--- |
|---|
| 245 | !--- There is also a correction for lhour < 0 and lhour > 23 |
|---|
| 246 | !--- to lhour between 0 and 23. |
|---|
| 247 | !--- |
|---|
| 248 | IF (zone(ilon) < 13) THEN |
|---|
| 249 | change = zone(ilon)-13 |
|---|
| 250 | lhour(ilon) = gmt+change |
|---|
| 251 | ENDIF |
|---|
| 252 | !--- |
|---|
| 253 | IF (zone(ilon) == 13) THEN |
|---|
| 254 | lhour(ilon) = gmt |
|---|
| 255 | ENDIF |
|---|
| 256 | !--- |
|---|
| 257 | IF (zone(ilon) > 13) THEN |
|---|
| 258 | change = zone(ilon)-13 |
|---|
| 259 | lhour(ilon) = gmt+change |
|---|
| 260 | ENDIF |
|---|
| 261 | IF (lhour(ilon) < 0) lhour(ilon) = lhour(ilon)+24 |
|---|
| 262 | IF (lhour(ilon) >= 24) lhour(ilon) = lhour(ilon)-24 |
|---|
| 263 | !--- |
|---|
| 264 | ENDDO |
|---|
| 265 | !----------------------- |
|---|
| 266 | END SUBROUTINE time_zone |
|---|
| 267 | |
|---|
| 268 | |
|---|
| 269 | |
|---|
| 270 | !! ================================================================================================================================ |
|---|
| 271 | !! SUBROUTINE : downward_solar_flux |
|---|
| 272 | !! |
|---|
| 273 | !>\BRIEF |
|---|
| 274 | !! |
|---|
| 275 | !! DESCRIPTION : None |
|---|
| 276 | !! |
|---|
| 277 | !! RECENT CHANGE(S): None |
|---|
| 278 | !! |
|---|
| 279 | !! MAIN OUTPUT VARIABLE(S): ::solad, ::solai |
|---|
| 280 | !! |
|---|
| 281 | !! REFERENCE(S) :See for further details : |
|---|
| 282 | !! http://www.atmo.arizona.edu/~morrill/swrad.html |
|---|
| 283 | !! |
|---|
| 284 | !! FLOWCHART : None |
|---|
| 285 | !! \n |
|---|
| 286 | !_ ================================================================================================================================ |
|---|
| 287 | |
|---|
| 288 | SUBROUTINE downward_solar_flux (npoi,latitude,jday,rtime,cloud,nband,solad,solai) |
|---|
| 289 | |
|---|
| 290 | USE time, ONLY : one_year, calendar_str |
|---|
| 291 | |
|---|
| 292 | IMPLICIT NONE |
|---|
| 293 | |
|---|
| 294 | !! 0. Parameter and variables declaration |
|---|
| 295 | |
|---|
| 296 | !! 0.1 Input variables |
|---|
| 297 | |
|---|
| 298 | INTEGER, INTENT(in) :: npoi !! number of grid points (unitless) |
|---|
| 299 | REAL, DIMENSION(npoi), INTENT(in) :: latitude !! latitude in degrees |
|---|
| 300 | REAL, INTENT(in) :: jday |
|---|
| 301 | REAL,INTENT(in) :: rtime |
|---|
| 302 | REAL, DIMENSION(npoi), INTENT(in) :: cloud !! cloud fraction (0-1, unitless) |
|---|
| 303 | INTEGER,INTENT(in) :: nband !! number of visible bands (unitless) |
|---|
| 304 | |
|---|
| 305 | !! 0.2 Output variables |
|---|
| 306 | |
|---|
| 307 | REAL, DIMENSION(npoi,nband), INTENT(out) :: solad !! direct downward solar flux |
|---|
| 308 | !! @tex $(W.m^{-2})$ @endtex |
|---|
| 309 | REAL, DIMENSION(npoi,nband), INTENT(out) :: solai !! diffuse downward solar flux |
|---|
| 310 | !! @tex $(W.m^{-2})$ @endtex |
|---|
| 311 | |
|---|
| 312 | !! 0.4 Local variables |
|---|
| 313 | |
|---|
| 314 | ! Parametres orbitaux: |
|---|
| 315 | |
|---|
| 316 | REAL, SAVE :: ecc !! Eccentricity |
|---|
| 317 | !$OMP THREADPRIVATE(ecc) |
|---|
| 318 | REAL, SAVE :: perh !! Longitude of perihelie |
|---|
| 319 | !$OMP THREADPRIVATE(perh) |
|---|
| 320 | REAL, SAVE :: xob !! obliquity |
|---|
| 321 | !$OMP THREADPRIVATE(xob) |
|---|
| 322 | REAL, PARAMETER :: pir = pi/180. !! |
|---|
| 323 | REAL, SAVE :: step !! |
|---|
| 324 | !$OMP THREADPRIVATE(step) |
|---|
| 325 | REAL :: xl !! |
|---|
| 326 | REAL :: so !! |
|---|
| 327 | REAL :: xllp !! |
|---|
| 328 | REAL :: xee !! |
|---|
| 329 | REAL :: xse !! |
|---|
| 330 | REAL :: xlam !! |
|---|
| 331 | REAL :: dlamm !! |
|---|
| 332 | REAL :: anm !! |
|---|
| 333 | REAL :: ranm !! |
|---|
| 334 | REAL :: ranv !! |
|---|
| 335 | REAL :: anv !! |
|---|
| 336 | REAL :: tls !! |
|---|
| 337 | REAL :: rlam !! |
|---|
| 338 | REAL :: sd !! |
|---|
| 339 | REAL :: cd !! |
|---|
| 340 | REAL :: deltar !! |
|---|
| 341 | REAL :: delta !! |
|---|
| 342 | REAL :: Dis_ST !! |
|---|
| 343 | REAL :: ddt !! |
|---|
| 344 | REAL :: orbit !! |
|---|
| 345 | REAL :: angle !! |
|---|
| 346 | REAL :: xdecl !! |
|---|
| 347 | REAL :: xlat !! |
|---|
| 348 | REAL, DIMENSION(npoi) :: trans !! |
|---|
| 349 | REAL, DIMENSION(npoi) :: fdiffuse !! |
|---|
| 350 | REAL, DIMENSION(npoi) :: coszen !! cosine of solar zenith angle |
|---|
| 351 | REAL :: sw !! |
|---|
| 352 | REAL :: frac !! |
|---|
| 353 | INTEGER :: i,ib !! |
|---|
| 354 | LOGICAL, SAVE :: lstep_init = .TRUE. !! |
|---|
| 355 | !$OMP THREADPRIVATE(lstep_init) |
|---|
| 356 | |
|---|
| 357 | !_ ================================================================================================================================ |
|---|
| 358 | |
|---|
| 359 | IF (lstep_init) THEN |
|---|
| 360 | IF ( TRIM(calendar_str) .EQ. 'gregorian' ) THEN |
|---|
| 361 | step = 1.0 |
|---|
| 362 | ELSE |
|---|
| 363 | step = one_year/365.2425 |
|---|
| 364 | ENDIF |
|---|
| 365 | !- |
|---|
| 366 | ! Read Orbital Parameters |
|---|
| 367 | !- |
|---|
| 368 | |
|---|
| 369 | !Config Key = ECCENTRICITY |
|---|
| 370 | !Config Desc = Use prescribed values |
|---|
| 371 | !Config If = ALLOW_WEATHERGEN |
|---|
| 372 | !Config Def = 0.016724 |
|---|
| 373 | !Config Help = |
|---|
| 374 | !Config Units = [-] |
|---|
| 375 | ecc = 0.016724 |
|---|
| 376 | CALL getin_p ('ECCENTRICITY',ecc) |
|---|
| 377 | IF (printlev >= 1) WRITE(numout,*) 'ECCENTRICITY: ',ecc |
|---|
| 378 | ! |
|---|
| 379 | !Config Key = PERIHELIE |
|---|
| 380 | !Config Desc = Use prescribed values |
|---|
| 381 | !Config If = ALLOW_WEATHERGEN |
|---|
| 382 | !Config Def = 102.04 |
|---|
| 383 | !Config Help = |
|---|
| 384 | !Config Units = [-] |
|---|
| 385 | perh = 102.04 |
|---|
| 386 | CALL getin_p ('PERIHELIE',perh) |
|---|
| 387 | IF (printlev >= 1) WRITE(numout,*) 'PERIHELIE: ',perh |
|---|
| 388 | ! |
|---|
| 389 | !Config Key = OBLIQUITY |
|---|
| 390 | !Config Desc = Use prescribed values |
|---|
| 391 | !Config If = ALLOW_WEATHERGEN |
|---|
| 392 | !Config Def = 23.446 |
|---|
| 393 | !Config Help = |
|---|
| 394 | !Config Units = [Degrees] |
|---|
| 395 | xob = 23.446 |
|---|
| 396 | CALL getin_p ('OBLIQUITY',xob) |
|---|
| 397 | IF (printlev >= 1) WRITE(numout,*) 'OBLIQUITY: ',xob |
|---|
| 398 | |
|---|
| 399 | lstep_init = .FALSE. |
|---|
| 400 | ENDIF |
|---|
| 401 | |
|---|
| 402 | !- |
|---|
| 403 | ! calendar and orbital calculations |
|---|
| 404 | !- |
|---|
| 405 | |
|---|
| 406 | !- |
|---|
| 407 | ! calculate the earth's orbital angle (around the sun) in radians |
|---|
| 408 | orbit = 2.0*pi*jday/365.2425 |
|---|
| 409 | !- |
|---|
| 410 | ! calculate the solar hour angle in radians |
|---|
| 411 | angle = 2.0*pi*(rtime-12.0)/24.0 |
|---|
| 412 | !- |
|---|
| 413 | ! calculate the current solar declination angle |
|---|
| 414 | ! ref: global physical climatology, hartmann, appendix a |
|---|
| 415 | ! |
|---|
| 416 | ! xdecl = 0.006918 & |
|---|
| 417 | ! -0.399912*cos(orbit) & |
|---|
| 418 | ! +0.070257*sin(orbit) & |
|---|
| 419 | ! -0.006758*cos(2.0*orbit) & |
|---|
| 420 | ! +0.000907*sin(2.0*orbit) & |
|---|
| 421 | ! -0.002697*cos(3.0*orbit) & |
|---|
| 422 | ! +0.001480*sin(3.0*orbit) |
|---|
| 423 | ! |
|---|
| 424 | ! calculate the effective solar constant, |
|---|
| 425 | ! including effects of eccentricity |
|---|
| 426 | ! ref: global physical climatology, hartmann, appendix a |
|---|
| 427 | ! |
|---|
| 428 | ! sw = 1370.*( 1.000110 & |
|---|
| 429 | ! +0.034221*cos(orbit) & |
|---|
| 430 | ! +0.001280*sin(orbit) & |
|---|
| 431 | ! +0.000719*cos(2.0*orbit) & |
|---|
| 432 | ! +0.000077*sin(2.0*orbit)) |
|---|
| 433 | ! |
|---|
| 434 | ! correction Marie-France Loutre |
|---|
| 435 | ! |
|---|
| 436 | ! orbital parameters |
|---|
| 437 | ! |
|---|
| 438 | ! ecc = 0.016724 |
|---|
| 439 | ! perh = 102.04 |
|---|
| 440 | ! xob = 23.446 |
|---|
| 441 | !- |
|---|
| 442 | xl = perh+180.0 |
|---|
| 443 | ! so : sinus of obliquity |
|---|
| 444 | so = sin(xob*pir) |
|---|
| 445 | !- |
|---|
| 446 | xllp = xl*pir |
|---|
| 447 | xee = ecc*ecc |
|---|
| 448 | xse = sqrt(1.0d0-xee) |
|---|
| 449 | ! xlam : true long. sun for mean long. = 0 |
|---|
| 450 | xlam = (ecc/2.0+ecc*xee/8.0d0)*(1.0+xse)*sin(xllp)-xee/4.0 & |
|---|
| 451 | & *(0.5+xse)*sin(2.0*xllp)+ecc*xee/8.0*(1.0/3.0+xse) & |
|---|
| 452 | & *sin(3.0*xllp) |
|---|
| 453 | xlam =2.0d0*xlam/pir |
|---|
| 454 | ! dlamm : mean long. sun for ma-ja |
|---|
| 455 | dlamm =xlam+(INT(jday)-79)*step |
|---|
| 456 | anm = dlamm-xl |
|---|
| 457 | ranm = anm*pir |
|---|
| 458 | xee = xee*ecc |
|---|
| 459 | ! ranv : anomalie vraie (radian) |
|---|
| 460 | ranv = ranm+(2.0*ecc-xee/4.0)*sin(ranm)+5.0/4.0*ecc*ecc & |
|---|
| 461 | & *sin(2.0*ranm)+13.0/12.0*xee*sin(3.0*ranm) |
|---|
| 462 | ! anv : anomalie vraie (degrees) |
|---|
| 463 | anv = ranv/pir |
|---|
| 464 | ! tls : longitude vraie (degrees) |
|---|
| 465 | tls = anv+xl |
|---|
| 466 | ! rlam : longitude vraie (radian) |
|---|
| 467 | rlam = tls*pir |
|---|
| 468 | ! sd and cd: cosinus and sinus of solar declination angle (delta) |
|---|
| 469 | ! sinus delta = sin (obl)*sin(lambda) with lambda = real longitude |
|---|
| 470 | ! (Phd. thesis of Marie-France Loutre, ASTR-UCL, Belgium, 1993) |
|---|
| 471 | sd = so*sin(rlam) |
|---|
| 472 | cd = sqrt(1.0d0-sd*sd) |
|---|
| 473 | ! delta : Solar Declination (degrees, angle sun at equator) |
|---|
| 474 | deltar = atan(sd/cd) |
|---|
| 475 | delta = deltar/pir |
|---|
| 476 | !- |
|---|
| 477 | ! Eccentricity Effect |
|---|
| 478 | !- |
|---|
| 479 | Dis_ST =(1.0-ecc*ecc)/(1.0+ecc*cos(ranv)) |
|---|
| 480 | ! ddt : 1 / normalized earth's sun distance |
|---|
| 481 | ddt = 1.0/Dis_ST |
|---|
| 482 | !- |
|---|
| 483 | ! Insolation normal to the atmosphere (W/m2) |
|---|
| 484 | !- |
|---|
| 485 | sw = ddt *ddt *1370.d0 |
|---|
| 486 | !- |
|---|
| 487 | xdecl = deltar |
|---|
| 488 | !- |
|---|
| 489 | ! solar calculations |
|---|
| 490 | !- |
|---|
| 491 | do i=1,npoi |
|---|
| 492 | !--- |
|---|
| 493 | !-- calculate the latitude in radians |
|---|
| 494 | !--- |
|---|
| 495 | xlat = latitude(i)*pi/180.0 |
|---|
| 496 | !--- |
|---|
| 497 | !-- calculate the cosine of the solar zenith angle |
|---|
| 498 | !--- |
|---|
| 499 | coszen(i) = MAX(0.0, (sin(xlat)*sin(xdecl) & |
|---|
| 500 | & + cos(xlat)*cos(xdecl)*cos(angle))) |
|---|
| 501 | !--- |
|---|
| 502 | !-- calculate the solar transmission through the atmosphere |
|---|
| 503 | !-- using simple linear function of tranmission and cloud cover |
|---|
| 504 | !--- |
|---|
| 505 | !-- note that the 'cloud cover' data is typically obtained from |
|---|
| 506 | !-- sunshine hours -- not direct cloud observations |
|---|
| 507 | !--- |
|---|
| 508 | !-- where, cloud cover = 1 - sunshine fraction |
|---|
| 509 | !--- |
|---|
| 510 | !-- different authors present different values for the slope and |
|---|
| 511 | !-- intercept terms of this equation |
|---|
| 512 | !--- |
|---|
| 513 | !-- Friend, A: Parameterization of a global daily weather generator |
|---|
| 514 | !-- for terrestrial ecosystem and biogeochemical modelling, |
|---|
| 515 | !-- Ecological Modelling |
|---|
| 516 | !--- |
|---|
| 517 | !-- Spitters et al., 1986: Separating the diffuse and direct component |
|---|
| 518 | !-- of global radiation and its implications for modeling canopy |
|---|
| 519 | !-- photosynthesis, Part I: Components of incoming radiation, |
|---|
| 520 | !-- Agricultural and Forest Meteorology, 38, 217-229. |
|---|
| 521 | !--- |
|---|
| 522 | !-- A. Friend : trans = 0.251+0.509*(1.0-cloud(i)) |
|---|
| 523 | !-- Spitters et al. : trans = 0.200+0.560*(1.0-cloud(i)) |
|---|
| 524 | !--- |
|---|
| 525 | !-- we are using the values from A. Friend |
|---|
| 526 | !--- |
|---|
| 527 | trans(i) = 0.251+0.509*(1.0-cloud(i)) |
|---|
| 528 | !--- |
|---|
| 529 | !-- calculate the fraction of indirect (diffuse) solar radiation |
|---|
| 530 | !-- based upon the cloud cover |
|---|
| 531 | !--- |
|---|
| 532 | !-- note that these relationships typically are measured for either |
|---|
| 533 | !-- monthly or daily timescales, and may not be exactly appropriate |
|---|
| 534 | !-- for hourly calculations -- however, in ibis, cloud cover is fixed |
|---|
| 535 | !-- through the entire day so it may not make much difference |
|---|
| 536 | !--- |
|---|
| 537 | !-- method i -- |
|---|
| 538 | !--- |
|---|
| 539 | !-- we use a simple empirical relationships |
|---|
| 540 | !-- from Nikolov and Zeller (1992) |
|---|
| 541 | !--- |
|---|
| 542 | !-- Nikolov, N. and K.F. Zeller, 1992: A solar radiation algorithm |
|---|
| 543 | !-- for ecosystem dynamics models, Ecological Modelling, 61, 149-168. |
|---|
| 544 | !--- |
|---|
| 545 | fdiffuse(i) = 1.0045+trans(i)*( 0.0435+trans(i) & |
|---|
| 546 | & *(-3.5227+trans(i)*2.6313)) |
|---|
| 547 | !--- |
|---|
| 548 | IF (trans(i) > 0.75) fdiffuse(i) = 0.166 |
|---|
| 549 | !--- |
|---|
| 550 | !-- method ii -- |
|---|
| 551 | !--- |
|---|
| 552 | !-- another method was suggested by Spitters et al. (1986) based on |
|---|
| 553 | !-- long-term data from the Netherlands |
|---|
| 554 | !-- |
|---|
| 555 | !-- Spitters et al., 1986: Separating the diffuse and direct component |
|---|
| 556 | !-- of global radiation and its implications for modeling canopy |
|---|
| 557 | !-- photosynthesis, Part I: Components of incoming radiation, |
|---|
| 558 | !-- Agricultural and Forest Meteorology, 38, 217-229. |
|---|
| 559 | !--- |
|---|
| 560 | !-- if ((trans == 0.00).and.(trans < 0.07)) then |
|---|
| 561 | !-- fdiffuse = 1.0 |
|---|
| 562 | !-- else if ((trans >= 0.07).and.(trans < 0.35)) then |
|---|
| 563 | !-- fdiffuse = 1.0-2.3*(trans-0.07)**2 |
|---|
| 564 | !-- else if ((trans >= 0.35).and.(trans < 0.75)) then |
|---|
| 565 | !-- fdiffuse = 1.33-1.46*trans |
|---|
| 566 | !-- ELSE |
|---|
| 567 | !-- fdiffuse = 0.23 |
|---|
| 568 | !-- endif |
|---|
| 569 | !--- |
|---|
| 570 | ENDDO |
|---|
| 571 | !----------------------- |
|---|
| 572 | !- |
|---|
| 573 | ! do for each waveband |
|---|
| 574 | !- |
|---|
| 575 | DO ib=1,nband |
|---|
| 576 | !--- |
|---|
| 577 | !-- calculate the fraction in each waveband |
|---|
| 578 | !--- |
|---|
| 579 | !-- GK010200 |
|---|
| 580 | IF (nband == 2) then |
|---|
| 581 | frac = 0.46+0.08*REAL(ib-1) |
|---|
| 582 | ELSE |
|---|
| 583 | frac = 1./REAL(nband) |
|---|
| 584 | ENDIF |
|---|
| 585 | !--- |
|---|
| 586 | DO i=1,npoi |
|---|
| 587 | !----- |
|---|
| 588 | !---- calculate the direct and indirect solar radiation |
|---|
| 589 | !----- |
|---|
| 590 | solad(i,ib) = sw*coszen(i)*frac*trans(i)*(1.-fdiffuse(i)) |
|---|
| 591 | solai(i,ib) = sw*coszen(i)*frac*trans(i)*fdiffuse(i) |
|---|
| 592 | ENDDO |
|---|
| 593 | ENDDO |
|---|
| 594 | |
|---|
| 595 | |
|---|
| 596 | END SUBROUTINE downward_solar_flux |
|---|
| 597 | |
|---|
| 598 | |
|---|
| 599 | END MODULE solar |
|---|
