phys/module_ra_farms.F

   1   module module_ra_farms
   2
   3     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   4     !                                                 !
   5     ! Purpose: Couples the FARMS model to WRF         !
   6     !                                                 !
   7     ! Author: Yu Xie coded the FARMS model            !
   8     !         Pedro A. Jimenez coupled FARMS to WRF   !
   9     !                                                 !
  10     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  11
  12     use module_model_constants, only : G
  13
  14     implicit none
  15
  16     private
  17
  18     public :: Farms_driver
  19
  20     real, parameter :: THREE_OVER_TWO = 3.0 / 2.0
  21     integer, parameter :: TAU_ICE_METHOD = 1
  22     logical, parameter :: USE_REST2 = .false.
  23
  24     real, parameter :: PTAU_MIN = 0.000001
  25     real, parameter :: DELTA_MIN = 0.000001
  26     real, parameter :: DE_ICE_MIN = 5.0, DE_ICE_MAX = 140.0
  27     real, parameter :: DE_CLOUD_MIN = 5.0, DE_CLOUD_MAX = 120.0
  28     real, parameter :: TAU_MIN = 0.0001, TAU_MAX = 300.0
  29     real, parameter :: RE_CLOUD_CLIM = 8.E-6, RE_ICE_CLIM = 24.E-6, RE_SNOW_CLIM = 24.E-6
  30
  31   contains
  32
  33     subroutine Farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, kte,  &
  34            p8w, rho, dz8w, albedo, aer_opt, aerssa2d, aerasy2d, aod5502d, angexp2d,  &
  35            coszen_loc, qv, qi, qs, qc, re_cloud, re_ice, re_snow,    &
  36            julian, swdown2, swddir2, swddni2, swddif2, swdownc2, swddnic2, &
  37            has_reqc, has_reqi, has_reqs, cldfra)
  38
  39       implicit None
  40
  41       integer, intent(in) :: ims, ime, jms, jme, its, ite, jts, jte, kms, kme, &
  42           kts, kte
  43
  44       integer, intent(in) :: aer_opt, has_reqc, has_reqi, has_reqs
  45       real,    intent(in) :: julian
  46
  47       real, dimension(ims:ime, jms:jme), intent(in) :: albedo, coszen_loc
  48
  49       real, dimension(ims:ime, kms:kme, jms:jme ), intent(in) :: qv, qi, qs, qc, &
  50           p8w, rho, dz8w, cldfra
  51       real, dimension(ims:ime, kms:kme, jms:jme ), intent(in) :: re_cloud, re_ice, re_snow
  52
  53       real, dimension(ims:ime, jms:jme), intent(inout) :: aerssa2d, aerasy2d, aod5502d, angexp2d
  54       real, dimension(ims:ime,jms:jme), intent(inout) :: swddir2, swdown2, &
  55                                                          swddif2, swddni2, &
  56                                                          swdownc2, swddnic2
  57
  58         ! Local
  59       integer :: i, j, k
  60       real    :: tau_qv, tau_qi, tau_qs, pmw, swp, iwp, lwp, beta
  61       real    :: re_cloud_path, re_ice_path, re_snow_path, q_aux, cldfra_h
  62       real, dimension(kms:kme) :: rhodz, re_cloud_k, re_ice_k, re_snow_k
  63
  64
  65       j_loop: do j = jts, jte
  66          i_loop: do i = its, ite
  67            daytime_if: if (coszen_loc(i, j) <= 0.0 ) then
  68              swdown2(i, j) = 0.0
  69              swddni2(i, j) = 0.0
  70              swddir2(i, j) = 0.0
  71              swddif2(i, j) = 0.0
  72              swdownc2(i, j) = 0.0
  73              swddnic2(i, j) = 0.0
  74            else
  75              rhodz(:) = rho(i, :, j) * dz8w(i, :, j) / (1. + qv(i, :, j))
  76              re_cloud_k(:) = re_cloud(i, :, j)
  77              re_ice_k(:) = re_ice(i, :, j)
  78              re_snow_k(:) = re_snow(i, :, j)
  79
  80              if (has_reqc == 1) then
  81                do k = kts, kte
  82                  if (CLDFRA (i, k, j) > 0.0 .and. re_cloud_k(k) < 2.5E-6) re_cloud_k(k) = RE_CLOUD_CLIM
  83                end do
  84              else
  85                re_cloud_k(:) = RE_CLOUD_CLIM
  86              end if
  87
  88              if (has_reqi == 1) then
  89                do k = kts, kte
  90                  if (cldfra(i, k, j) > 0.0 .and. re_ice_k(k) < 5.0E-6) re_ice_k(k) = RE_ICE_CLIM
  91                end do
  92              else
  93                re_ice_k(:) = RE_ICE_CLIM
  94              end if
  95
  96              if (has_reqs == 1) then
  97                do k = kts, kte
  98                  if (cldfra(i, k, j) > 0.0 .and. re_snow_k(k) < 10.0E-6) re_snow_k(k) = RE_SNOW_CLIM
  99                end do
 100              else
 101                re_snow_k(:) = RE_SNOW_CLIM
 102              end if
 103
 104                ! PMW
 105              pmw = integrate_1var (rhodz, qv(i, :, j), kms, kme, kts, kte)
 106
 107                ! Calc effective radius water
 108              q_aux = integrate_1var (rhodz, qc(i, :, j), kms, kme, kts, kte)
 109              lwp = q_aux
 110
 111              if (q_aux > 0.0) then
 112                re_cloud_path = integrate_2var (rhodz, qc(i, :, j), &
 113                    re_cloud_k, kms, kme, kts, kte)
 114                re_cloud_path = re_cloud_path / q_aux
 115              else
 116                re_cloud_path = 0.0
 117              end if
 118
 119                ! Calc effective radius ice
 120              q_aux = integrate_1var (rhodz, qi(i, :, j), kms, kme, kts, kte)
 121              iwp = q_aux
 122
 123              if (q_aux > 0.0) then
 124                re_ice_path = integrate_2var (rhodz, qi(i, :, j), &
 125                    re_ice_k, kms, kme, kts, kte)
 126                re_ice_path = re_ice_path / q_aux
 127              else
 128                re_ice_path = 0.0
 129              end if
 130
 131                ! Calc effective radius snow
 132              q_aux = integrate_1var (rhodz, qs(i, :, j), kms, kme, kts, kte)
 133              swp = q_aux
 134
 135              if (q_aux > 0.0) then
 136                re_snow_path = integrate_2var (rhodz, qs(i, :, j), &
 137                    re_snow_k, kms, kme, kts, kte)
 138                re_snow_path = re_snow_path / q_aux
 139              else
 140                re_snow_path = 0.0
 141              end if
 142
 143              ! Calculate horizontal cloud fraction
 144              q_aux = integrate_1var (rhodz, qc(i, :, j) + qi(i, :, j) + qs(i, :, j), kms, kme, kts, kte)
 145              if (q_aux > 0.0) then
 146                cldfra_h = integrate_2var (rhodz, qc(i, :, j) + qi(i, :, j) + qs(i, :, j), &
 147                    cldfra(i, :, j), kms, kme, kts, kte)
 148                cldfra_h = cldfra_h / q_aux
 149              else
 150                cldfra_h = 0.0
 151              end if
 152
 153                ! optical thickness  water
 154              if (re_cloud_path > 0.0) then
 155                tau_qv = THREE_OVER_TWO * lwp / re_cloud_path / 1000.0
 156              else
 157                tau_qv = 0.0
 158              end if
 159
 160                ! Optical thickness ice
 161              if (re_ice_path > 0.0) then
 162                if (TAU_ICE_METHOD == 1) then
 163                    ! Eq 10 in Matrosov et al. (2002)
 164                  tau_qi = iwp * 1000.0 * (0.02 + 4.2 / (2.0 * re_ice_path * 1.0e+6))
 165                else
 166                  tau_qi = iwp * 1000.0 * (-0.006656 + 3.686 / (2.0 * re_ice_path * 1.0e+6))
 167                end if
 168              else
 169                tau_qi = 0.0
 170              end if
 171
 172                ! Optical thickness snow
 173              if (re_snow_path > 0.0) then
 174                if (TAU_ICE_METHOD == 1) then
 175                    ! Eq 10 in Matrosov et al. (2002)
 176                  tau_qs = swp * 1000.0 * (0.02 + 4.2 / (2.0 * re_snow_path * 1.0e+6))
 177                else
 178                  tau_qs = swp * 1000.0 * (-0.006656 + 3.686 / (2.0 * re_snow_path * 1.0e+6))
 179                end if
 180              else
 181                tau_qs = 0.0
 182              end if
 183
 184              beta = aod5502d(i, j) * (1000.0/ 550.0) ** (- angexp2d(i, j))
 185
 186              Call Farms (p8w(i, 1, j), albedo(i, j), aerssa2d(i, j),         &
 187                  aerasy2d(i, j), coszen_loc(i, j), beta,                     &
 188                  angexp2d(i, j), pmw, tau_qv, tau_qi, tau_qs, cldfra_h,      &
 189                  re_cloud_path, re_ice_path, re_snow_path, int(julian),      &
 190                  swdown2(i, j), swddni2(i, j), swddif2(i, j), swddir2(i, j), &
 191                  swdownc2(i, j), swddnic2(i, j))
 192
 193           end if daytime_if
 194         end do i_loop
 195       end do j_loop
 196
 197     end subroutine Farms_driver
 198
 199
 200     function Integrate_1var(rhodz, var1_1d, kms, kme, kts, kte) &
 201         result (return_value)
 202
 203       implicit none
 204
 205       integer, intent(in) :: kts, kte, kms, kme
 206       real, dimension(kms:kme), intent(in) :: var1_1d, rhodz
 207
 208         ! Local
 209       real :: return_value
 210       integer :: k
 211
 212       return_value = 0.0
 213       do k = kts, kte - 1
 214         return_value = return_value + var1_1d(k) * rhodz(k)
 215       end do
 216
 217     end function Integrate_1var
 218
 219
 220     function Integrate_2var(rhodz, var1_1d, var2_1d, kms, kme, kts, kte) &
 221         result (return_value)
 222
 223       implicit none
 224
 225       integer, intent(in) :: kts, kte, kms, kme
 226       real, dimension(kms:kme), intent(in) :: var1_1d, var2_1d, rhodz
 227
 228         ! Local
 229       real :: return_value
 230       integer :: k
 231
 232       return_value = 0.0
 233       do k = kts, kte - 1
 234         return_value = return_value + var1_1d(k) * var2_1d(k) * rhodz(k)
 235       end do
 236
 237     end function Integrate_2var
 238
 239
 240     subroutine FARMS (p_pa, albdo, ssa, g, solarangle, beta, alpha, w_mm, &
 241         tau_qv, tau_qi, tau_qs, cldfra_h, re_cloud_path_m, re_ice_path_m, re_snow_path_m, &
 242         juday, ghi, dni, dif, dir, ghi_clear, dni_clear)
 243
 244         !!!!!! This Fast All-sky Radiation Model for Solar applications (FARMS) was developed by
 245         !!!!!! Yu Xie (Yu.Xie@nrel.gov). Please contact him for more
 246         !!!!!! information.
 247
 248         ! Adapted by PAJ to couple it with WRF
 249
 250         ! Input values
 251         ! p_pa: surface air pressure (Pa)
 252         ! albdo: surface albedo
 253         ! SSA: single-scattering albedo for aerosol
 254         ! g: asymmetric factor of aerosol
 255         ! Z: solar zenith angle
 256         ! beta: Angstrom turbidity coeff., i.e. AOD at 1000 nm.
 257         ! alpha: Angstrom wavelength exponent
 258         !        For a cloudy sky, the suggested values are beta=0.1 and alpha=1.2,
 259         !        unless you have more accurate measurements.
 260         ! w: PWV (cm)
 261         ! tau_qv: cloud optical thickness (liquid water)
 262         ! tau_qi: cloud optical thickness (ice)
 263         ! tau_qs: cloud optical thickness (snow)
 264         ! re_cloud_path_m: effective radious liquid
 265         ! re_ice_path_m: effective radious ice
 266         ! re_snow_path_m: effective radious snow
 267         ! juday: day of year (1-366)
 268         ! ozone: ozone amount (1000DU, i.e. 0.4=400DU)
 269         ! phase: cloud thermodynamic phase (1=water, 2=ice)
 270         ! De: Cloud effective particle size (micron)
 271         !     (For a water cloud, De=2*effective radius)
 272         !     (For an ice cloud, the defination of De follows Fu. 1996.
 273         !     Another useful defination is given by "Determination of ice cloud models using MODIS and MISR data" )
 274
 275
 276       implicit none
 277
 278       real, intent(in) :: p_pa, albdo, ssa, g, solarangle, beta, alpha, w_mm, &
 279           tau_qv, tau_qi, tau_qs, re_cloud_path_m, re_ice_path_m, re_snow_path_m, cldfra_h
 280       integer, intent(in) :: juday
 281       real, intent(out) :: ghi, dni, dir, dif, ghi_clear, dni_clear
 282
 283         ! Local vars
 284       real :: de_cloud, de_cloud2, de_ice, de_ice2, de_snow, de_snow2, f0, f1, ftotal, ozone, &
 285           p, radius, Ruucld, pi, Ruucld_water, Ruucld_ice, Ruucld_snow, Ruuclr, tau_qv2, tau_qi2, &
 286           tau_qs2, tau, Tddcld, Tddclr, Tducld, Tducld_ice, tducld_snow, tducld_water, Tduclr, Tuuclr, &
 287           w, Z, tau_tot
 288
 289
 290       PI = acos(-1.0)
 291
 292       p = p_pa / 100.0
 293       ozone = 0.265
 294       w = w_mm / 10.0
 295
 296       de_cloud = 2.0 * re_cloud_path_m / 1.0e-6
 297       de_ice = 2.0 * re_ice_path_m / 1.0e-6
 298       de_snow = 2.0 * re_snow_path_m / 1.0e-6
 299
 300       de_cloud = 2.0 * re_cloud_path_m / 1.0e-6
 301       de_cloud2 = Max (de_cloud, DE_CLOUD_MIN)
 302       de_cloud2 = Min (de_cloud2, DE_CLOUD_MAX)
 303
 304       de_ice = 2.0 * re_ice_path_m / 1.0e-6
 305       de_ice2 = Max (de_ice, DE_ICE_MIN)
 306       de_ice2 = Min (de_ice2, DE_ICE_MAX)
 307
 308       de_snow = 2.0 * re_snow_path_m / 1.0e-6
 309       de_snow2 = Max (de_snow, DE_ICE_MIN)
 310       de_snow2 = Min (de_snow2, DE_ICE_MAX)
 311
 312       tau_qv2 = Max (tau_qv, TAU_MIN)
 313       tau_qv2 = Min (tau_qv2, TAU_MAX)
 314       tau_qi2 = Max (tau_qi, TAU_MIN)
 315       tau_qi2 = Min (tau_qi2, TAU_MAX)
 316       tau_qs2 = Max (tau_qs, TAU_MIN)
 317       tau_qs2 = Min (tau_qs2, TAU_MAX)
 318
 319       tau_tot = 0.0
 320       if (tau_qv > TAU_MIN) tau_tot = tau_tot + tau_qv2
 321       if (tau_qi > TAU_MIN) tau_tot = tau_tot + tau_qi2
 322       if (tau_qs > TAU_MIN) tau_tot = tau_tot + tau_qs2
 323       tau_tot = Min (tau_tot, TAU_MAX)
 324
 325       Z = acos(solarangle) * 180.0 / PI
 326
 327       call SUNEARTH( juday, Radius )
 328       F0 = 1361.2/(Radius*Radius)
 329
 330       call CLEARSKYALL(p, albdo, SSA, g, Z, Radius, beta,&
 331           alpha, ozone, w, Tddclr, Tduclr, Ruuclr, Tuuclr)
 332
 333         ! Liquid hydrometeors
 334       if (tau_qv > TAU_MIN) then
 335         Call Watermodel (tau_tot, de_cloud2, solarangle, Tducld_water, Ruucld_water)
 336       else
 337         Tducld_water = 0.0
 338         Ruucld_water = 0.0
 339         tau_qv2 = 0.0
 340       end if
 341
 342         ! Ice hydrometeors
 343       if (tau_qi > TAU_MIN) then
 344         call Icemodel (tau_tot, de_ice2, solarangle, Tducld_ice, Ruucld_ice)
 345       else
 346         Tducld_ice = 0.0
 347         Ruucld_ice = 0.0
 348         tau_qi2 = 0.0
 349       end if
 350
 351         ! Snow hydrometeors
 352       if (tau_qs > TAU_MIN) then
 353         call Icemodel (tau_tot, de_snow2, solarangle, Tducld_snow, Ruucld_snow)
 354       else
 355         Tducld_snow = 0.0
 356         Ruucld_snow = 0.0
 357         tau_qs2 = 0.0
 358       end if
 359
 360       if (tau_tot > 0.0) then
 361         Tducld =  (tau_qv2 * Tducld_water + tau_qi2 * Tducld_ice + tau_qs2 * Tducld_snow) / tau_tot
 362         Tducld = min(Tducld, 1.0)
 363         Tducld = max(Tducld, 0.0)
 364
 365         Ruucld = (tau_qv2 * Ruucld_water + tau_qi2 * Ruucld_ice + tau_qs2 * Ruucld_snow)  / tau_tot
 366         Ruucld = min(Ruucld, 1.0)
 367         Ruucld = max(Ruucld, 0.0)
 368       else
 369         Tducld = 0.0
 370         Ruucld = 0.0
 371       end if
 372
 373       tau = tau_qv2 + tau_qi2 + tau_qs2
 374       tau = Min (tau, TAU_MAX)
 375       tau = Max (tau, TAU_MIN)
 376
 377       Tddcld = exp(-tau/solarangle)
 378
 379       dni_clear = F0 * Tddclr
 380       dni = dni_clear * Tddcld
 381       dir = dni * solarangle
 382
 383       F1 = solarangle*F0*( Tddcld*(Tddclr+Tduclr) + Tducld*Tuuclr )
 384
 385       Ftotal = F1/(1.0-albdo*(Ruuclr + Ruucld*Tuuclr*Tuuclr))
 386       ghi = Ftotal
 387       ghi_clear =  solarangle * F0 * ((Tddclr + Tduclr) / (1.0 - albdo * Ruuclr))
 388
 389       ghi = cldfra_h * ghi + (1.0 - cldfra_h) * ghi_clear
 390       dni = cldfra_h * dni + (1.0 - cldfra_h) * dni_clear
 391
 392       dif = ghi - dir
 393
 394     end subroutine farms
 395
 396
 397     subroutine SUNEARTH( juday, R )
 398
 399       implicit none
 400
 401       integer, intent(in) :: juday
 402       real, intent(out) :: R
 403
 404       real :: pi, b, R1
 405
 406
 407       PI = acos(-1.0)
 408
 409       b = 2.0*PI*juday/365.0
 410       R1 = 1.00011 + 0.034221*cos(b) + 0.001280*sin(b) + &
 411            0.000719*cos(2.0*b) +0.000077*sin(2.0*b)
 412
 413       R = R1**(-0.5)
 414
 415     end subroutine SUNEARTH
 416
 417
 418     subroutine CLEARSKYALL( p, albdo, SSA, g, Z, Radius, beta,&
 419         alpha, ozone, w, Tddclr, Tduclr, Ruuclr, Tuuclr )
 420
 421         implicit none
 422
 423         real, intent(in) :: p, albdo, SSA, g, Z, Radius, beta, alpha, ozone, w
 424         real, intent(out) :: Tddclr, Tduclr, Ruuclr, Tuuclr
 425
 426           ! Local vars
 427         integer, parameter :: nangle = 10
 428         integer :: i
 429         real :: mu(nangle), angle(nangle), aa(nangle)
 430         real :: PI
 431
 432
 433         PI = acos(-1.0)
 434
 435         do i=1, nangle
 436           mu(i) = (i-1.0)*0.1+0.1
 437           angle(i) = acos(mu(i))*180.0/PI
 438         end do
 439
 440         do i=1, nangle
 441             ! PAJ: can use both REST2 and BIRD
 442           if (USE_REST2) then
 443             call CLEARSKY (p, albdo, SSA, g, angle(i), Radius, beta,&
 444             alpha, ozone, w, Tddclr, Tduclr, Ruuclr)
 445           else
 446             call BIRD(p, albdo, SSA, g, angle(i), Radius, beta,&
 447             alpha, ozone, w, Tddclr, Tduclr, Ruuclr)
 448           end if
 449           aa(i) = Tddclr
 450         end do
 451
 452         Tuuclr = 0.0
 453         do i=1, nangle
 454           Tuuclr = Tuuclr + 2.0*mu(i)*aa(i)*0.1
 455         end do
 456
 457           ! PAJ: can use both REST2 and BIRD
 458         if (USE_REST2) then
 459           call CLEARSKY(p, albdo, SSA, g, Z, Radius, beta,&
 460           alpha, ozone, w, Tddclr, Tduclr, Ruuclr)
 461         else
 462           call BIRD(p, albdo, SSA, g, Z, Radius, beta,&
 463           alpha, ozone, w, Tddclr, Tduclr, Ruuclr)
 464         end if
 465
 466     end subroutine CLEARSKYALL
 467
 468
 469     subroutine CLEARSKY(p, albdo, SSA, g, Z, Radius, beta,&
 470         alpha, ozone, w, Tddclr, Tduclr, Ruuclr)
 471
 472       implicit none
 473
 474       real, intent(in) :: p, albdo, SSA, g, Z, Radius, beta, alpha, ozone, w
 475       real, intent(out) :: Tddclr, Tduclr, Ruuclr
 476
 477
 478       return
 479
 480     end subroutine CLEARSKY
 481
 482
 483     subroutine BIRD(p, albdo, SSA, g, Z, Radius, beta,&
 484         alpha, ozone, w, Tddclr, Tduclr, Ruuclr)
 485
 486         ! This clear-sky model follows the equations given by Bird (1981)
 487         ! This subroutine for the all-sky fast model is given by
 488         ! Yu Xie (Yu.Xie@nrel.gov)
 489
 490       implicit none
 491
 492       real, intent (in) :: p, albdo, SSA, g, Z, Radius, beta, alpha, ozone, w
 493       real, intent (out) :: Tddclr, Tduclr, Ruuclr
 494
 495         ! Local vars
 496       real :: degrad, airmass, airmassp, T_rayleigh, x0, T_o, T_gases, xw, T_water, &
 497           tau038, tau050, taua, T_aerosol, T_AA, T_AS, F0, Fdif, Ftotal, Fddclr
 498
 499
 500       degrad=.017453293d+00
 501       airmass = 1/(cos(Z*degrad)+0.15*(93.885-Z)**(-1.25))
 502       airmassp = p*airmass/1013.0
 503
 504       T_rayleigh = exp(-0.0903*(airmassp**0.84)*(1.0 + airmassp -&
 505           airmassp**1.01))
 506
 507       x0 = ozone*airmass
 508       T_o = 1.0 - 0.1611*x0*(1.0+139.48*x0)**(-0.3035) -&
 509           0.002715*x0/(1.0+0.044*x0 + 0.0003*x0**2.0)
 510       T_gases = exp( -0.0127*(airmassp**0.26) )
 511       xw = w*airmass
 512       T_water = 1.0 - 2.4959*xw/( (1.0+79.034*xw)**0.6828 + 6.385*xw )
 513
 514       tau038 = beta*(0.38**(-alpha))
 515       tau050 = beta*(0.50**(-alpha))
 516       taua = 0.2758*tau038 + 0.35*tau050
 517
 518       T_aerosol = exp( -(taua**0.873)*(1.0+taua-taua**0.7088)&
 519           *(airmass**0.9108) )
 520       T_AA = 1.0 - 0.1*(1.0 - airmass + airmass**1.06)*&
 521           (1.0 - T_aerosol)
 522       T_AS = T_aerosol/T_AA
 523
 524       Ruuclr = 0.0685 + (1.0-g)*(1.0-T_AS)
 525
 526       F0 = 1.0
 527       Fddclr = F0*(cos(Z*degrad))*0.9662*T_rayleigh*T_o*T_gases*T_water*T_aerosol
 528
 529       Fdif = F0*( cos(Z*degrad) )*0.79*T_o*T_gases*T_water*T_AA*&
 530           (0.5*(1.0-T_rayleigh) + g*(1.0-T_AS) )&
 531          /(1.0 - airmass + airmass**1.02 )
 532
 533       Ftotal = ( Fddclr + Fdif )/(1.0 - albdo*Ruuclr)
 534
 535       Tddclr = Fddclr/( F0*cos(Z*degrad) )
 536       Tduclr = (Ftotal - Fddclr)/(F0*cos(Z*degrad))
 537
 538     end subroutine BIRD
 539
 540
 541     subroutine WATERMODEL( tau, De, solarangle, Tducld, Ruucld )
 542
 543       implicit none
 544
 545       real, intent (in) :: tau, De, solarangle
 546       real, intent (out) :: Tducld, Ruucld
 547
 548         ! Local vars
 549       real :: solarconst, Ptau, PDHI, delta, y, PPDHI
 550
 551
 552       solarconst = 1385.72180
 553       Ptau = (2.8850+0.002*(De-60.0))*solarangle-0.007347
 554       Ptau = max(Ptau, PTAU_MIN)
 555       PDHI = (1.0+(De-60.0)*0.0002)*1087.24*solarangle**1.1605
 556
 557       delta = -0.644531*solarangle+1.20117+0.129807/solarangle &
 558           -0.00121096/(solarangle*solarangle) + &
 559           1.52587e-07/(solarangle*solarangle*solarangle)
 560       delta = Max (delta, DELTA_MIN)
 561
 562       y = 0.012*(tau-Ptau)*solarangle
 563       PPDHI = (1.0+SINH(y))*PDHI*&
 564           exp(-( (log10(tau)-log10(Ptau))**2.0 )/delta)
 565
 566       Tducld = PPDHI/(solarconst*solarangle)
 567
 568       if (tau .lt. 1.0) then
 569         Ruucld = 0.107359*tau
 570       endif
 571
 572       if (tau .ge. 1.0) then
 573         Ruucld = 1.03 - exp(-(0.5+log10(tau))*(0.5+log10(tau))/3.105 )
 574       endif
 575
 576     end subroutine WATERMODEL
 577
 578
 579     subroutine ICEMODEL( tau, De, solarangle, Tducld, Ruucld )
 580
 581       implicit none
 582
 583       real, intent (in) :: tau, De, solarangle
 584       real, intent (out) :: Tducld, Ruucld
 585
 586         ! Local vars
 587
 588       real :: solarconst, Ptau, PDHI, delta, y, PPDHI
 589
 590       solarconst = 1385.72180
 591
 592       if ( De .le. 26.0 ) then
 593         Ptau = 2.8487*solarangle- 0.0029
 594       endif
 595       if ( De .gt. 26.0 ) then
 596         Ptau =  (2.8355 + (100.0-De)*0.006)*solarangle - 0.00612
 597       endif
 598       Ptau = max(Ptau, PTAU_MIN)
 599
 600       PDHI = 1047.6367*solarangle**1.0883
 601
 602       delta = -0.0549531*solarangle+0.617632+0.17876/(solarangle) &
 603           -0.002174/(solarangle*solarangle)
 604       delta = Max (delta, DELTA_MIN)
 605
 606
 607       y = 0.01*(tau-Ptau)*solarangle
 608       PPDHI = (1.0+SINH(y))*PDHI*&
 609           exp(-( (log10(tau)-log10(Ptau))**2.0 )/delta)
 610
 611       Tducld = PPDHI/(solarconst*solarangle)
 612
 613       if (tau .lt. 1.0) then
 614         Ruucld = 0.094039*tau
 615       endif
 616
 617       if (tau .ge. 1.0) then
 618         Ruucld = 1.02 - exp(-(0.5+log10(tau))*(0.5+log10(tau))/3.25 )
 619       endif
 620
 621     end subroutine ICEMODEL
 622
 623   end module module_ra_farms
 624
 625
 626