[7541] | 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: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_global/solar.f90 $ |
---|
| 19 | !! $Date: 2017-10-18 11:15:06 +0200 (Wed, 18 Oct 2017) $ |
---|
| 20 | !! $Revision: 4693 $ |
---|
| 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 |
---|