source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_regions.F90 @ 4062

! radiation_regions.F90 -- Properties of horizontal regions in Tripleclouds & SPARTACUS
!
! (C) Copyright 2016- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!
! Author:  Robin Hogan
! Email:   r.j.hogan@ecmwf.int
!
! Modifications
!   2017-07-14  R. Hogan  Incorporate gamma distribution option
!   2018-10-06  R. Hogan  Merged from radiation_optical_depth_scaling.h and radiation_overlap.F90

module radiation_regions

  implicit none

  public :: calc_region_properties

contains

  !---------------------------------------------------------------------
  ! Compute the optical depth scalings for the optically "thick" and
  ! "thin" regions of a Tripleclouds representation of a sub-grid PDF
  ! of cloud optical depth. Following Shonk and Hogan (2008), the 16th
  ! percentile is used for the thin region, and the formulas estimate
  ! this for both lognormal and gamma distributions. However, an
  ! adjustment is needed for the gamma distribution at large
  ! fractional standard deviations.
  subroutine calc_region_properties(nlev, nreg, istartcol, iendcol, do_gamma, &
       &  cloud_fraction, frac_std, reg_fracs, od_scaling, cloud_fraction_threshold)

    use parkind1,     only : jprb
    use yomhook,      only : lhook, dr_hook
    use radiation_io, only : nulerr, radiation_abort

    ! Minimum od_scaling in the case of a Gamma distribution
    real(jprb), parameter :: MinGammaODScaling = 0.025_jprb

    ! At large fractional standard deviations (FSDs), we cannot
    ! capture the behaviour of a gamma distribution with two equally
    ! weighted points; we need to weight the first ("lower") point
    ! more.  The weight of the first point is normally 0.5, but for
    ! FSDs between 1.5 and 3.725 it rises linearly to 0.9, and for
    ! higher FSD it is capped at this value.  The weight of the second
    ! point is one minus the first point.
    real(jprb), parameter :: MinLowerFrac      = 0.5_jprb
    real(jprb), parameter :: MaxLowerFrac      = 0.9_jprb
    real(jprb), parameter :: FSDAtMinLowerFrac = 1.5_jprb
    real(jprb), parameter :: FSDAtMaxLowerFrac = 3.725_jprb
    ! Between FSDAtMinLowerFrac and FSDAtMaxLowerFrac,
    ! LowerFrac=LowerFracFSDIntercept+FSD*LowerFracFSDGradient
    real(jprb), parameter :: LowerFracFSDGradient &
         &  = (MaxLowerFrac-MinLowerFrac) / (FSDAtMaxLowerFrac-FSDAtMinLowerFrac)
    real(jprb), parameter :: LowerFracFSDIntercept &
         &  = MinLowerFrac - FSDAtMinLowerFrac*LowerFracFSDGradient

    ! Number of levels and regions
    integer, intent(in) :: nlev, nreg

    ! Range of columns to process
    integer, intent(in) :: istartcol, iendcol

    ! Do we do a lognormal or gamma distribution?
    logical, intent(in) :: do_gamma

    ! Cloud fraction, i.e. the fraction of the gridbox assigned to all
    ! regions numbered 2 and above (region 1 is clear sky)
    real(jprb), intent(in), dimension(:,:)  :: cloud_fraction ! (ncol,nlev)

    ! Fractional standard deviation of in-cloud water content
    real(jprb), intent(in), dimension(:,:)  :: frac_std       ! (ncol,nlev)

    ! Fractional area coverage of each region
    real(jprb), intent(out) :: reg_fracs(1:nreg,nlev,istartcol:iendcol)

    ! Optical depth scaling for the cloudy regions
    real(jprb), intent(out) :: od_scaling(2:nreg,nlev,istartcol:iendcol)

    ! Regions smaller than this are ignored
    real(jprb), intent(in), optional :: cloud_fraction_threshold

    ! In case the user doesn't supply cloud_fraction_threshold we use
    ! a default value
    real(jprb) :: frac_threshold

    ! Loop indices
    integer :: jcol, jlev

    real(jprb) :: hook_handle

    if (lhook) call dr_hook('radiation_region_properties:calc_region_properties',0,hook_handle)

    if (present(cloud_fraction_threshold)) then
      frac_threshold = cloud_fraction_threshold
    else
      frac_threshold = 1.0e-20_jprb
    end if

    if (nreg == 2) then
      ! Only one clear-sky and one cloudy region: cloudy region is
      ! homogeneous
      reg_fracs(2,1:nlev,istartcol:iendcol)  = transpose(cloud_fraction(istartcol:iendcol,1:nlev))
      reg_fracs(1,1:nlev,istartcol:iendcol)  = 1.0_jprb - reg_fracs(2,1:nlev,istartcol:iendcol)
      od_scaling(2,1:nlev,istartcol:iendcol) = 1.0_jprb

    else if (nreg == 3) then
      ! Two cloudy regions with optical depth scaled by 1-x and
      ! 1+x.
      ! Simple version which fails when fractional_std >= 1:
      !od_scaling(2) = 1.0_jprb-cloud%fractional_std(jcol,jlev)
      ! According to Shonk and Hogan (2008), 1-FSD should correspond to
      ! the 16th percentile. 
      if (.not. do_gamma) then
        ! If we treat the distribution as a lognormal such that the
        ! equivalent Normal has a mean mu and standard deviation
        ! sigma, then the 16th percentile of the lognormal is very
        ! close to exp(mu-sigma).
        do jcol = istartcol,iendcol
          do jlev = 1,nlev
            if (cloud_fraction(jcol,jlev) < frac_threshold) then
              reg_fracs(1,jlev,jcol)       = 1.0_jprb
              reg_fracs(2:3,jlev,jcol)  = 0.0_jprb
              od_scaling(2:3,jlev,jcol) = 1.0_jprb
            else
              reg_fracs(1,jlev,jcol) = 1.0_jprb - cloud_fraction(jcol,jlev)
              reg_fracs(2:3,jlev,jcol) = cloud_fraction(jcol,jlev)*0.5_jprb
              od_scaling(2,jlev,jcol) &
                   &  = exp(-sqrt(log(frac_std(jcol,jlev)**2+1))) &
                   &  / sqrt(frac_std(jcol,jlev)**2+1)
              od_scaling(3,jlev,jcol) = 2.0_jprb - od_scaling(2,jlev,jcol)
            end if
          end do
        end do
      else
        ! If we treat the distribution as a gamma then the 16th
        ! percentile is close to the following.  Note that because it
        ! becomes vanishingly small for FSD >~ 2, we have a lower
        ! limit of 1/40, and for higher FSDs reduce the fractional
        ! cover of the denser region - see region_fractions routine
        ! below
        do jcol = istartcol,iendcol
          do jlev = 1,nlev
            if (cloud_fraction(jcol,jlev) < frac_threshold) then
              reg_fracs(1,jlev,jcol)    = 1.0_jprb
              reg_fracs(2:3,jlev,jcol)  = 0.0_jprb
              od_scaling(2:3,jlev,jcol) = 1.0_jprb
            else
              ! Fraction of the clear-sky region
              reg_fracs(1,jlev,jcol) = 1.0_jprb - cloud_fraction(jcol,jlev)
!#define OLD_GAMMA_REGION_BEHAVIOUR 1
#ifdef OLD_GAMMA_REGION_BEHAVIOUR
              ! Use previous behaviour (ecRad version 1.1.5 and
              ! earlier): cloudy fractions are the same and there is
              ! no minimum optical depth scaling; this tends to lead
              ! to an overprediction of the reflection from scenes
              ! with a large fractional standard deviation of optical
              ! depth.
              ! Fraction and optical-depth scaling of the lower of the
              ! two cloudy regions
              reg_fracs(2,jlev,jcol) = cloud_fraction(jcol,jlev) * 0.5_jprb
              od_scaling(2,jlev,jcol) = &
                   &  exp(-frac_std(jcol,jlev)*(1.0_jprb + 0.5_jprb*frac_std(jcol,jlev) &
                   &                     *(1.0_jprb+0.5_jprb*frac_std(jcol,jlev))))

#else
              ! Improved behaviour.
              ! Fraction and optical-depth scaling of the lower of the
              ! two cloudy regions
              reg_fracs(2,jlev,jcol) = cloud_fraction(jcol,jlev) &
                   &  * max(MinLowerFrac, min(MaxLowerFrac, &
                   &  LowerFracFSDIntercept + frac_std(jcol,jlev)*LowerFracFSDGradient))
              od_scaling(2,jlev,jcol) = MinGammaODScaling &
                   &  + (1.0_jprb - MinGammaODScaling) &
                   &    * exp(-frac_std(jcol,jlev)*(1.0_jprb + 0.5_jprb*frac_std(jcol,jlev) &
                   &                     *(1.0_jprb+0.5_jprb*frac_std(jcol,jlev))))
#endif
              ! Fraction of the upper of the two cloudy regions
              reg_fracs(3,jlev,jcol) = 1.0_jprb-reg_fracs(1,jlev,jcol)-reg_fracs(2,jlev,jcol)
              ! Ensure conservation of the mean optical depth
              od_scaling(3,jlev,jcol) = (cloud_fraction(jcol,jlev) &
                   &  -reg_fracs(2,jlev,jcol)*od_scaling(2,jlev,jcol)) / reg_fracs(3,jlev,jcol)
            end if
          end do
        end do
      end if
    else ! nreg > 3
      write(nulerr,'(a)') '*** Error: only 2 or 3 regions may be specified'
      call radiation_abort()
    end if

    if (lhook) call dr_hook('radiation_region_properties:calc_region_properties',1,hook_handle)

  end subroutine calc_region_properties

end module radiation_regions