source: LMDZ6/branches/blowing_snow/libf/phylmd/ecrad/radiation_cloud_cover.F90 @ 5450

! radiation_cloud_cover.F90 - Compute cumulative cloud cover for McICA
!
! (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
!
! Generate profiles of the cumulative cloud cover as seen from TOA,
! used in the McICA cloud generator.
!
! Modifications
!   2020-10-07  R. Hogan  Ensure iobj1 initialized in case of alpha_obj==0

module radiation_cloud_cover

  use parkind1, only           : jprb

  public

  ! Three overlap schemes.  Note that "Exponential" means that
  ! clear-sky regions have no special significance for computing the
  ! cumulative cloud cover: non-contiguous clouds are exponentially
  ! rather than randomly overlapped. This is the situaition in the
  ! McRad radiation scheme at ECMWF.
  enum, bind(c)
    enumerator IOverlapMaximumRandom, IOverlapExponentialRandom, &
         &     IOverlapExponential
  end enum
  character(len=*), parameter :: OverlapName(0:2) = (/ 'Max-Ran', &
       &                                               'Exp-Ran', &
       &                                               'Exp-Exp' /)

  ! Maximum cloud fraction distinguishable from 1
  real(jprb), parameter :: MaxCloudFrac = 1.0_jprb-epsilon(1.0_jprb)*10.0_jprb


contains

  !---------------------------------------------------------------------
  ! Convert "beta" overlap parameter of Shonk et al. (2010) to "alpha"
  ! overlap parameter of Hogan and Illingworth (2000)
  elemental function beta2alpha(beta, frac1, frac2)

    implicit none

    ! Beta overlap parameter and the cloud fractions in the upper and
    ! lower layers
    real(jprb), intent(in) :: beta, frac1, frac2

    real(jprb)             :: beta2alpha

    ! Absolute difference in cloud fraction
    real(jprb)             :: frac_diff

    if (beta < 1.0_jprb) then
      frac_diff = abs(frac1-frac2)
      beta2alpha = beta &
           &  + (1.0_jprb-beta)*frac_diff &
           &  / (frac_diff + 1.0_jprb/beta - 1.0_jprb)
    else
      beta2alpha = 1.0_jprb
    end if

  end function beta2alpha


  !---------------------------------------------------------------------
  ! Compute total cloud cover according to the specified overlap
  ! rule. This can be used to compute the high, mid and low cloud
  ! cover by passing in subsets of the cloud fraction array
  function cloud_cover(nlev, i_overlap_scheme, frac, overlap_param, &
       &               is_beta_overlap)

    implicit none
    
    ! Number of levels and the overlap scheme to be applied
    integer, intent(in)    :: nlev, i_overlap_scheme

    ! Cloud fraction and the overlap parameter between adjacent pairs
    ! of levels
    real(jprb), intent(in) :: frac(nlev), overlap_param(nlev-1)

    ! Do we use the "beta" overlap scheme of Shonk et al. (2010)?
    ! Default is false.
    logical, intent(in), optional :: is_beta_overlap

    ! Return cloud cover
    real(jprb)             :: cloud_cover

    ! Cumulative cloud cover from TOA to the base of each layer
    real(jprb) :: cum_cloud_cover(nlev)

    ! Cloud cover of a pair of layers
    real(jprb) :: pair_cloud_cover(nlev-1)

    if (i_overlap_scheme == IOverlapExponentialRandom) then
      call cum_cloud_cover_exp_ran(nlev, frac, overlap_param, &
           &   cum_cloud_cover, pair_cloud_cover, is_beta_overlap)
    else if (i_overlap_scheme == IOverlapExponential) then
      call cum_cloud_cover_exp_exp(nlev, frac, overlap_param, &
           &   cum_cloud_cover, pair_cloud_cover, is_beta_overlap)
    else
      call cum_cloud_cover_max_ran(nlev, frac, cum_cloud_cover, &
           &                       pair_cloud_cover)
    end if

    cloud_cover = cum_cloud_cover(nlev)

  end function cloud_cover


  !---------------------------------------------------------------------
  ! Maximum-random overlap: Geleyn & Hollingsworth formula
  subroutine cum_cloud_cover_max_ran(nlev, frac, &
       & cum_cloud_cover, pair_cloud_cover)

    use yomhook,  only           : lhook, dr_hook

    implicit none

    ! Inputs
    integer, intent(in)     :: nlev  ! number of model levels

    ! Cloud fraction on full levels
    real(jprb), intent(in)  :: frac(nlev)

    ! Outputs

    ! Cumulative cloud cover from TOA to the base of each layer
    real(jprb), intent(out) :: cum_cloud_cover(nlev)

    ! Cloud cover of a pair of layers
    real(jprb), intent(out) :: pair_cloud_cover(nlev-1)

    ! Local variables

    ! Cumulative product needed in computation of total_cloud_cover
    real(jprb) :: cum_product

    ! Loop index for model level
    integer :: jlev

    real(jprb) :: hook_handle

    if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_max_ran',0,hook_handle)

    ! Loop to compute total cloud cover and the cumulative cloud cover

    ! down to the base of each layer
    cum_product = 1.0_jprb - frac(1)
    cum_cloud_cover(1) = frac(1)
    do jlev = 1,nlev-1
      ! Compute the combined cloud cover of layers jlev and jlev+1
      pair_cloud_cover(jlev) = max(frac(jlev),frac(jlev+1))

      if (frac(jlev) >= MaxCloudFrac) then
        ! Cloud cover has reached one
        cum_product = 0.0_jprb
      else
        cum_product = cum_product * (1.0_jprb - pair_cloud_cover(jlev)) &
             &  / (1.0_jprb - frac(jlev))
      end if
      cum_cloud_cover(jlev+1) = 1.0_jprb - cum_product
    end do

    if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_max_ran',1,hook_handle)

  end subroutine cum_cloud_cover_max_ran
  

  !---------------------------------------------------------------------
  ! Exponential-random overlap: exponential overlap for contiguous
  ! clouds, random overlap for non-contiguous clouds
  subroutine cum_cloud_cover_exp_ran(nlev, frac, overlap_param, &
       & cum_cloud_cover, pair_cloud_cover, is_beta_overlap)

    use yomhook,  only           : lhook, dr_hook

    implicit none

    ! Inputs
    integer, intent(in)     :: nlev  ! number of model levels

    ! Cloud fraction on full levels
    real(jprb), intent(in)  :: frac(nlev)

    ! Cloud overlap parameter for interfaces between model layers,
    ! where 0 indicates random overlap and 1 indicates maximum-random
    ! overlap
    real(jprb), intent(in)  :: overlap_param(nlev-1)

    ! This routine has been coded using the "alpha" overlap parameter
    ! of Hogan and Illingworth (2000). If the following logical is
    ! present and true then the input is interpretted to be the "beta"
    ! overlap parameter of Shonk et al. (2010), and needs to be
    ! converted to alpha.
    logical, intent(in), optional :: is_beta_overlap

    ! Outputs

    ! Cumulative cloud cover from TOA to the base of each layer
    real(jprb), intent(out) :: cum_cloud_cover(nlev)

    ! Cloud cover of a pair of layers
    real(jprb), intent(out) :: pair_cloud_cover(nlev-1)

    ! Local variables

    ! Cumulative product needed in computation of total_cloud_cover
    real(jprb) :: cum_product

    ! "Alpha" overlap parameter
    real(jprb) :: overlap_alpha
    logical    :: do_overlap_conversion

    ! Loop index for model level
    integer :: jlev

    real(jprb) :: hook_handle

    if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_exp_ran',0,hook_handle)

    
    if (present(is_beta_overlap)) then
      do_overlap_conversion = is_beta_overlap
    else
      do_overlap_conversion = .false.
    end if

    ! Loop to compute total cloud cover and the cumulative cloud cover

    ! down to the base of each layer
    cum_product = 1.0_jprb - frac(1)
    cum_cloud_cover(1) = frac(1)
    do jlev = 1,nlev-1
      ! Convert to "alpha" overlap parameter if necessary
      if (do_overlap_conversion) then
        overlap_alpha = beta2alpha(overlap_param(jlev), &
             &                     frac(jlev), frac(jlev+1))
      else
        overlap_alpha = overlap_param(jlev)
      end if

      ! Compute the combined cloud cover of layers jlev and jlev+1
      pair_cloud_cover(jlev) = overlap_alpha*max(frac(jlev),frac(jlev+1)) &
           &  + (1.0_jprb - overlap_alpha) &
           &  * (frac(jlev)+frac(jlev+1)-frac(jlev)*frac(jlev+1))

      if (frac(jlev) >= MaxCloudFrac) then
        ! Cloud cover has reached one
        cum_product = 0.0_jprb
      else
        cum_product = cum_product * (1.0_jprb - pair_cloud_cover(jlev)) &
             &  / (1.0_jprb - frac(jlev))
      end if
      cum_cloud_cover(jlev+1) = 1.0_jprb - cum_product
    end do


    if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_exp_ran',1,hook_handle)

  end subroutine cum_cloud_cover_exp_ran



  !---------------------------------------------------------------------
  ! Exponential-exponential overlap: exponential overlap for both
  ! contiguous and non-contiguous clouds. This is the result of the
  ! simple Raisanen cloud generator, but unfortunately it has no
  ! (known) analytic formula for the total cloud cover, or the
  ! cumulative cloud cover.  In partially cloudy columns, The McICA
  ! scheme needs this info in order to devote all the cloudy g-points
  ! to columns containing cloud, which reduces McICA noise. The
  ! following routine provides an approximate estimate of cumulative
  ! cloud cover consistent with the exponential-exponential scheme.
  subroutine cum_cloud_cover_exp_exp(nlev, frac, overlap_param, &
       & cum_cloud_cover, pair_cloud_cover, is_beta_overlap)

    use yomhook,  only           : lhook, dr_hook

    implicit none

    ! Inputs
    integer, intent(in)     :: nlev  ! number of model levels

    ! Cloud fraction on full levels
    real(jprb), intent(in)  :: frac(nlev)

    ! Cloud overlap parameter for interfaces between model layers,
    ! where 0 indicates random overlap and 1 indicates maximum-random
    ! overlap
    real(jprb), intent(in)  :: overlap_param(nlev-1)

    ! This routine has been coded using the "alpha" overlap parameter
    ! of Hogan and Illingworth (2000). If the following logical is
    ! present and true then the input is interpretted to be the "beta"
    ! overlap parameter of Shonk et al. (2010), and needs to be
    ! converted to alpha.
    logical, intent(in), optional :: is_beta_overlap

    ! Outputs

    ! Cumulative cloud cover from TOA to the base of each layer
    real(jprb), intent(out) :: cum_cloud_cover(nlev)

    ! Cloud cover of a pair of layers
    real(jprb), intent(out) :: pair_cloud_cover(nlev-1)

    ! Local variables

    ! If this routine is called from the radiation_interface tree then
    ! very low cloud fractions have already been set to zero, but if
    ! it is called as a cloud cover diagnostic then this can't be
    ! guaranteed so a small non-zero numbers is required
    real(jprb), parameter :: min_frac = 1.0e-6_jprb

    ! "Alpha" overlap parameter
    real(jprb) :: overlap_alpha(nlev-1)
    logical    :: do_overlap_conversion

    ! Variables describing "concave cloud objects", i.e. contiguous
    ! layers where the cloud fraction monotonically increases then
    ! monotonically decreases

    ! Number of objects
    integer    :: nobj

    ! Indices to the location of the top, maximum cloud fraction, and
    ! base, of each concave cloud object
    integer    :: i_top_obj(nlev)
    integer    :: i_max_obj(nlev)
    integer    :: i_base_obj(nlev)
    ! Poor-man's linked list to allow for deletion of objects: this
    ! variable points to the index of the next active object
    integer    :: i_next_obj(nlev)

    ! Cloud cover of object
    real(jprb) :: cc_obj(nlev)

    ! Overlap parameter between objects
    real(jprb) :: alpha_obj(nlev)

    ! Do (while) loop index for model level
    integer :: jlev

    ! Do loop index for object
    integer :: jobj

    ! Maximum correlation between adjacent objects
    real(jprb) :: alpha_max

    ! Indices to pair of objects
    integer    :: iobj1, iobj2

    ! Combined cloud cover of pair of objects, and scaling to modify
    ! cumulative cloud cover of lower layer
    real(jprb) :: cc_pair, scaling

    real(jprb) :: hook_handle

    if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_exp_exp',0,hook_handle)

    
    if (present(is_beta_overlap)) then
      do_overlap_conversion = is_beta_overlap
    else
      do_overlap_conversion = .false.
    end if

    ! Loop down through atmosphere to locate objects and compute their
    ! basic properties
    jlev = 1
    nobj = 0
    do while (jlev <= nlev)
      if (frac(jlev) > min_frac) then
        ! Starting a new object: store its top
        nobj = nobj + 1
        i_top_obj(nobj) = jlev;
        ! Find its maximum cloud fraction
        jlev = jlev + 1
        do while (jlev <= nlev) 
          if (frac(jlev) < frac(jlev-1)) then
            exit
          end if
          jlev = jlev + 1
        end do
        i_max_obj(nobj) = jlev - 1
        ! Find its base
        do while (jlev <= nlev)
          if (frac(jlev) > frac(jlev-1) .or. frac(jlev) <= min_frac) then
            exit
          end if
          jlev = jlev + 1
        end do
        ! In the case of cloud fraction profile starting from the top
        ! like this: 0.1 0.2 0.1 0.2 0.1, we may want the object grouping
        ! to be (0.1 0.2) (0.1 0.2 0.1), not (0.1 0.2 0.1) (0.2 0.1),
        ! in which case the following should be uncommented
        !        if (jlev < nlev) then
        !          if (frac(jlev) > frac(jlev-1)) then
        !            jlev = jlev - 1
        !          end if
        !        end if
        i_base_obj(nobj) = jlev - 1
        ! Index to the next object
        i_next_obj(nobj) = nobj+1
      else
        jlev = jlev + 1
      end if
    end do

    ! Array assignments
    cum_cloud_cover = 0.0_jprb
    pair_cloud_cover = 0.0_jprb

    if (nobj > 0) then
      ! Only do any more work if there is cloud present
      
      ! To minimize the potential calls to beta2alpha, we do all the
      ! computations related to overlap parameter here
      if (.not. do_overlap_conversion) then

        ! Compute the combined cloud cover of pairs of layers
        do jlev = 1,nlev-1
          pair_cloud_cover(jlev) &
               &  = overlap_param(jlev)*max(frac(jlev),frac(jlev+1)) &
               &  + (1.0_jprb - overlap_param(jlev)) &
               &  * (frac(jlev)+frac(jlev+1)-frac(jlev)*frac(jlev+1))
        end do
        ! Estimate the effective overlap parameter "alpha_obj" between
        ! adjacent objects as the product of the layerwise overlap
        ! parameters between their layers of maximum cloud fraction
        do jobj = 1,nobj-1
          alpha_obj(jobj) &
               &  = product(overlap_param(i_max_obj(jobj):i_max_obj(jobj+1)-1))
        end do

      else

        ! Convert Shonk et al overlap parameter to Hogan and
        ! Illingworth definition
        overlap_alpha = beta2alpha(overlap_param, &
             &                     frac(1:nlev-1), frac(2:nlev))
        ! Compute the combined cloud cover of pairs of layers
        do jlev = 1,nlev-1
          pair_cloud_cover(jlev) &
               &  = overlap_alpha(jlev)*max(frac(jlev),frac(jlev+1)) &
               &  + (1.0_jprb - overlap_alpha(jlev)) &
               &  * (frac(jlev)+frac(jlev+1)-frac(jlev)*frac(jlev+1))          
        end do
        ! Estimate the effective overlap parameter "alpha_obj" between
        ! adjacent objects as the product of the layerwise overlap
        ! parameters between their layers of maximum cloud fraction
        do jobj = 1,nobj-1
          alpha_obj(jobj) &
               &  = product(overlap_alpha(i_max_obj(jobj):i_max_obj(jobj+1)-1))
        end do

      end if


      ! Compute the cumulative cloud cover working down from the top
      ! of each object: this will later be converted to the cumulative
      ! cloud cover working down from TOA
      do jobj = 1,nobj
        cum_cloud_cover(i_top_obj(jobj)) = frac(i_top_obj(jobj))
        do jlev = i_top_obj(jobj), i_base_obj(jobj)-1
          if (frac(jlev) >= MaxCloudFrac) then
            ! Cloud cover has reached one
            cum_cloud_cover(jlev+1) = 1.0_jprb
          else
            cum_cloud_cover(jlev+1) = 1.0_jprb &
                 &  - (1.0_jprb - cum_cloud_cover(jlev)) &
                 &  * (1.0_jprb - pair_cloud_cover(jlev))  &
                 &  / (1.0_jprb - frac(jlev))
          end if
        end do
        cc_obj(jobj) = cum_cloud_cover(i_base_obj(jobj))
      end do

      iobj1 = 1

      ! Sequentially combine objects until there is only one left
      ! covering the full vertical extent of clouds in the profile
      do while (nobj > 1)
        ! Find the most correlated adjacent pair of objects
        alpha_max = 0.0_jprb

        ! Need to re-initialize iobj1 here in case alpha_obj(:)==0.0,
        ! which would mean that the "if" statement in the following
        ! loop would never get triggered
        iobj1 = 1

        jobj = 1
        do while (jobj < nobj)
          if (alpha_obj(jobj) > alpha_max) then
            alpha_max = alpha_obj(jobj)
            iobj1 = jobj
          end if
          jobj = i_next_obj(jobj)
        end do

        ! iobj1 is the index to the first object in the pair, set
        ! iobj2 to the second
        iobj2 = i_next_obj(iobj1)

        ! Set the cumulative cloud cover in the clear-sky gap between
        ! the objects to the value at the base of the upper object
        cum_cloud_cover(i_base_obj(iobj1)+1:i_top_obj(iobj2)-1) &
             &  = cum_cloud_cover(i_base_obj(iobj1))

        ! Calculate the combined cloud cover of the pair of objects
        cc_pair = alpha_obj(iobj1)*max(cc_obj(iobj1), cc_obj(iobj2)) &
             &  + (1.0_jprb - alpha_obj(iobj1)) &
             &  * (cc_obj(iobj1) + cc_obj(iobj2) - cc_obj(iobj1)*cc_obj(iobj2))
        scaling = min(max((cc_pair-cc_obj(iobj1)) / max(min_frac, cc_obj(iobj2)), &
             &            0.0_jprb), &
             &        1.0_jprb)
        
        ! Scale the combined cloud cover of the lower object to
        ! account for its overlap with the upper object
        do jlev = i_top_obj(iobj2),i_base_obj(iobj2)
          cum_cloud_cover(jlev) = cum_cloud_cover(i_base_obj(iobj1)) &
               +  cum_cloud_cover(jlev) * scaling
        end do
        
        ! Merge the objects by setting the properties of the upper
        ! object to the combined properties of both.  Note that
        ! i_max_obj is not modified because it is no longer needed.
        cc_obj(iobj1) = cc_pair
        i_base_obj(iobj1) = i_base_obj(iobj2)
        i_next_obj(iobj1) = i_next_obj(iobj2)
        alpha_obj(iobj1)  = alpha_obj(iobj2)
        nobj = nobj - 1
      end do

      ! Finish off the total cloud cover below cloud
      cum_cloud_cover(i_base_obj(iobj1)+1:nlev) &
           &  = cum_cloud_cover(i_base_obj(iobj1)) 

      ! Ensure that the combined cloud cover of pairs of layers is
      ! consistent with the overhang
      do jlev = 1,nlev-1
        pair_cloud_cover(jlev) = max(pair_cloud_cover(jlev), &
             &     frac(jlev)+cum_cloud_cover(jlev+1)-cum_cloud_cover(jlev))
      end do

      ! Sometimes round-off error can lead to cloud cover just above
      ! one, which in turn can lead to direct shortwave fluxes just
      ! below zero
      cum_cloud_cover = min(cum_cloud_cover, 1.0_jprb)

    end if ! cloud is present in profile

    if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_exp_exp',1,hook_handle)

  end subroutine cum_cloud_cover_exp_exp

end module radiation_cloud_cover