source: LMDZ6/trunk/libf/phylmd/ecrad.v1.5.1/radiation_cloud_generator.F90 @ 5444

! radiation_cloud_generator.F90 - Generate water-content or optical-depth scalings for McICA
!
! (C) Copyright 2015- 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 clouds for McICA using a method modified from Raisanen et
! al. (2002)
!
! Modifications
!   2018-02-22  R. Hogan  Call masked version of PDF sampler for speed
!   2020-03-31  R. Hogan  More vectorizable version of Exp-Ran

module radiation_cloud_generator

  public

contains

  !---------------------------------------------------------------------
  ! Generate scaling factors for the cloud optical depth to represent
  ! cloud overlap, the overlap of internal cloud inhomogeneities and
  ! the fractional standard deviation of these inhomogeneities, for
  ! use in a Monte Carlo Independent Column Approximation radiation
  ! scheme. All returned profiles contain cloud, and the total cloud
  ! cover is also returned, so the calling function can then do a
  ! weighted average of clear and cloudy skies; this is a way to
  ! reduce the Monte Carlo noise in profiles with low cloud cover.
  subroutine cloud_generator(ng, nlev, i_overlap_scheme, &
       &  iseed, frac_threshold, &
       &  frac, overlap_param, decorrelation_scaling, &
       &  fractional_std, pdf_sampler, &
       &  od_scaling, total_cloud_cover, &
       &  use_beta_overlap, use_vectorizable_generator)

    use parkind1, only           : jprb
    use yomhook,  only           : lhook, dr_hook
    use radiation_io,   only     : nulerr, radiation_abort
    use random_numbers_mix, only : randomnumberstream, &
         initialize_random_numbers, uniform_distribution
    use radiation_pdf_sampler, only : pdf_sampler_type
    use radiation_cloud_cover, only : cum_cloud_cover_exp_ran, &
         &       cum_cloud_cover_max_ran, cum_cloud_cover_exp_exp, &
         &       IOverlapMaximumRandom, IOverlapExponentialRandom, &
         &       IOverlapExponential

    implicit none

    ! Inputs
    integer, intent(in)     :: ng    ! number of g points
    integer, intent(in)     :: nlev  ! number of model levels
    integer, intent(in)     :: i_overlap_scheme
    integer, intent(in)     :: iseed ! seed for random number generator

    ! Only cloud fractions above this threshold are considered to be
    ! clouds
    real(jprb), intent(in)  :: frac_threshold

    ! 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)

    ! Overlap parameter for internal inhomogeneities
    real(jprb), intent(in)  :: decorrelation_scaling

    ! Fractional standard deviation at each layer
    real(jprb), intent(in)  :: fractional_std(nlev)

    ! Object for sampling from a lognormal or gamma distribution
    type(pdf_sampler_type), intent(in) :: pdf_sampler

    ! 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 :: use_beta_overlap

    ! Do we use the more vectorizable cloud generator, at the expense
    ! of more random numbers being needed?
    logical, intent(in), optional :: use_vectorizable_generator

    ! Outputs

    ! Cloud optical depth scaling factor, with 0 indicating clear sky
    real(jprb), intent(out) :: od_scaling(ng,nlev)

    ! Total cloud cover using cloud fraction and overlap parameter
    real(jprb), intent(out) :: total_cloud_cover

    ! Local variables

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

    ! Scaled random number for finding cloud
    real(jprb) :: trigger

    ! Uniform deviates between 0 and 1
    real(jprb) :: rand_top(ng)

    ! Overlap parameter of inhomogeneities
    real(jprb) :: overlap_param_inhom(nlev-1)

    ! Seed for random number generator and stream for producing random
    ! numbers
    type(randomnumberstream) :: random_stream
    
    ! First and last cloudy layers
    integer :: ibegin, iend

    integer :: itrigger

    ! Loop index for model level and g-point
    integer :: jlev, jg

    ! Cloud cover of a pair of layers, and amount by which cloud at
    ! next level increases total cloud cover as seen from above
    real(jprb), dimension(nlev-1) :: pair_cloud_cover, overhang

    logical :: use_vec_gen

    real(jprb) :: hook_handle

    if (lhook) call dr_hook('radiation_cloud_generator:cloud_generator',0,hook_handle)

    if (i_overlap_scheme == IOverlapExponentialRandom) then
      call cum_cloud_cover_exp_ran(nlev, frac, overlap_param, &
           &   cum_cloud_cover, pair_cloud_cover, use_beta_overlap)
    else if (i_overlap_scheme == IOverlapMaximumRandom) then
      call cum_cloud_cover_max_ran(nlev, frac, &
           &   cum_cloud_cover, pair_cloud_cover)
    else if (i_overlap_scheme == IOverlapExponential) then
      call cum_cloud_cover_exp_exp(nlev, frac, overlap_param, &
           &   cum_cloud_cover, pair_cloud_cover, use_beta_overlap)
    else
      write(nulerr,'(a)') '*** Error: cloud overlap scheme not recognised'
      call radiation_abort()
    end if

    total_cloud_cover = cum_cloud_cover(nlev);
    do jlev = 1,nlev-1
      overhang(jlev) = cum_cloud_cover(jlev+1)-cum_cloud_cover(jlev)
    end do

    if (total_cloud_cover < frac_threshold) then
      ! Treat column as clear sky: calling function therefore will not
      ! use od_scaling so we don't need to calculate it
      total_cloud_cover = 0.0_jprb

    else
      ! Cloud is present: need to calculate od_scaling

      ! Find range of cloudy layers
      jlev = 1
      do while (frac(jlev) <= 0.0_jprb) 
        jlev = jlev + 1
      end do
      ibegin = jlev
      iend = jlev
      do jlev = jlev+1,nlev
        if (frac(jlev) > 0.0_jprb) then
          iend = jlev
        end if
      end do

      ! Set overlap parameter of inhomogeneities
      overlap_param_inhom = overlap_param

      do jlev = ibegin,iend-1
        if (overlap_param(jlev) > 0.0_jprb) then
          overlap_param_inhom(jlev) &
               &  = overlap_param(jlev)**(1.0_jprb/decorrelation_scaling)
        end if
      end do

      ! Reset optical depth scaling to clear skies
      od_scaling = 0.0_jprb

      if (present(use_vectorizable_generator)) then
        use_vec_gen = use_vectorizable_generator
      else
        use_vec_gen = .false.
      end if

      if (.not. use_vec_gen) then
        ! Original generator that minimizes the number of random
        ! numbers used, but is not vectorizable

        ! Expensive operation: initialize random number generator for
        ! this column
        call initialize_random_numbers(iseed, random_stream)

        ! Compute ng random numbers to use to locate cloud top
        call uniform_distribution(rand_top, random_stream)
        
        ! Loop over ng columns
        do jg = 1,ng
          ! Find the cloud top height corresponding to the current
          ! random number, and store in itrigger
          trigger = rand_top(jg) * total_cloud_cover
          jlev = ibegin
          do while (trigger > cum_cloud_cover(jlev) .and. jlev < iend)
            jlev = jlev + 1
          end do
          itrigger = jlev
          
          if (i_overlap_scheme /= IOverlapExponential) then
            call generate_column_exp_ran(ng, nlev, jg, random_stream, pdf_sampler, &
                 &  frac, pair_cloud_cover, &
                 &  cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, &
                 &  itrigger, iend, od_scaling)
          else
            call generate_column_exp_exp(ng, nlev, jg, random_stream, pdf_sampler, &
                 &  frac, pair_cloud_cover, &
                 &  cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, &
                 &  itrigger, iend, od_scaling)
          end if
          
        end do

      else
        ! Alternative generator (only for Exp-Ran overlap so far) that
        ! should be vectorizable but generates more random numbers,
        ! some of which are not used

        if (i_overlap_scheme == IOverlapExponential) then
          write(nulerr,'(a)') '*** Error: vectorizable cloud generator is not available with Exp-Exp overlap'
          call radiation_abort()
        end if

        call generate_columns_exp_ran(ng, nlev, iseed, pdf_sampler, &
             &  total_cloud_cover, frac_threshold, frac, pair_cloud_cover, &
             &  cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, &
             &  ibegin, iend, od_scaling)

      end if

    end if

    if (lhook) call dr_hook('radiation_cloud_generator:cloud_generator',1,hook_handle)

  end subroutine cloud_generator


  !---------------------------------------------------------------------
  ! Generate a column of optical depth scalings using
  ! exponential-random overlap (which includes maximum-random overlap
  ! as a limiting case)
  subroutine generate_column_exp_ran(ng, nlev, ig, random_stream, pdf_sampler, &
       &  frac, pair_cloud_cover, &
       &  cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, &
       &  itrigger, iend, od_scaling)

    use parkind1,              only : jprb
    use radiation_pdf_sampler, only : pdf_sampler_type
    use random_numbers_mix,    only : randomnumberstream, &
         initialize_random_numbers, uniform_distribution


    implicit none

    ! Number of g points / columns, and number of current column
    integer, intent(in) :: ng, ig

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

    ! Stream for producing random numbers
    type(randomnumberstream), intent(inout) :: random_stream

    ! Object for sampling from a lognormal or gamma distribution
    type(pdf_sampler_type), intent(in) :: pdf_sampler

    ! Cloud fraction, cumulative cloud cover and fractional standard
    ! deviation in each layer
    real(jprb), intent(in), dimension(nlev) :: frac, cum_cloud_cover, fractional_std

    ! Cloud cover of a pair of layers, and amount by which cloud at
    ! next level increases total cloud cover as seen from above
    real(jprb), intent(in), dimension(nlev-1) :: pair_cloud_cover, overhang

    ! Overlap parameter of inhomogeneities
    real(jprb), intent(in), dimension(nlev-1) :: overlap_param_inhom

    ! Top of highest cloudy layer (in this subcolumn) and base of
    ! lowest
    integer, intent(in) :: itrigger, iend

    ! Optical depth scaling to output
    real(jprb), intent(inout), dimension(ng,nlev) :: od_scaling

    ! Height indices
    integer :: jlev, jcloud

    ! Number of contiguous cloudy layers for which to compute optical
    ! depth scaling
    integer :: n_layers_to_scale

    integer :: iy

    ! Is it time to fill the od_scaling variable?
    logical :: do_fill_od_scaling

    real(jprb) :: rand_cloud(nlev)
    real(jprb) :: rand_inhom1(nlev), rand_inhom2(nlev)

    ! So far our vertically contiguous cloud contains only one layer
    n_layers_to_scale = 1
    iy = 0

    ! Locate the clouds below this layer: first generate some more
    ! random numbers
    call uniform_distribution(rand_cloud(1:(iend+1-itrigger)),random_stream)

    ! Loop from the layer below the local cloud top down to the
    ! bottom-most cloudy layer
    do jlev = itrigger+1,iend+1
      do_fill_od_scaling = .false.
      if (jlev <= iend) then
        iy = iy+1
        if (n_layers_to_scale > 0) then
          ! There is a cloud above, in which case the probability
          ! of cloud in the layer below is as follows
          if (rand_cloud(iy)*frac(jlev-1) &
               &  < frac(jlev) + frac(jlev-1) - pair_cloud_cover(jlev-1)) then
            ! Add another cloudy layer
            n_layers_to_scale = n_layers_to_scale + 1
          else 
            ! Reached the end of a contiguous set of cloudy layers and
            ! will compute the optical depth scaling immediately.
            do_fill_od_scaling = .true.
          end if
        else
          ! There is clear-sky above, in which case the
          ! probability of cloud in the layer below is as follows
          if (rand_cloud(iy)*(cum_cloud_cover(jlev-1) - frac(jlev-1)) &
               &  < pair_cloud_cover(jlev-1) - overhang(jlev-1) - frac(jlev-1)) then
            ! A new cloud top
            n_layers_to_scale = 1
          end if
        end if
      else
        ! We are at the bottom of the cloudy layers in the model,
        ! so in a moment need to populate the od_scaling array
        do_fill_od_scaling = .true.
      end if

      if (do_fill_od_scaling) then
        ! We have a contiguous range of layers for which we
        ! compute the od_scaling elements using some random
        ! numbers
        call uniform_distribution(rand_inhom1(1:n_layers_to_scale),random_stream)
        call uniform_distribution(rand_inhom2(1:n_layers_to_scale),random_stream)

        ! Loop through the sequence of cloudy layers
        do jcloud = 2,n_layers_to_scale
          ! Use second random number, and inhomogeneity overlap
          ! parameter, to decide whether the first random number
          ! should be repeated (corresponding to maximum overlap)
          ! or not (corresponding to random overlap)
          if (rand_inhom2(jcloud) &
               &  < overlap_param_inhom(jlev-n_layers_to_scale+jcloud-2)) then
            rand_inhom1(jcloud) = rand_inhom1(jcloud-1)
          end if
        end do
        
        ! Sample from a lognormal or gamma distribution to obtain
        ! the optical depth scalings
        call pdf_sampler%sample(fractional_std(jlev-n_layers_to_scale:jlev-1), &
             & rand_inhom1(1:n_layers_to_scale), od_scaling(ig,jlev-n_layers_to_scale:jlev-1))

        n_layers_to_scale = 0
      end if
          
    end do

  end subroutine generate_column_exp_ran



  !---------------------------------------------------------------------
  ! Generate a column of optical depth scalings using
  ! exponential-exponential overlap
  subroutine generate_column_exp_exp(ng, nlev, ig, random_stream, pdf_sampler, &
       &  frac, pair_cloud_cover, &
       &  cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, &
       &  itrigger, iend, od_scaling)

    use parkind1,              only : jprb
    use radiation_pdf_sampler, only : pdf_sampler_type
    use random_numbers_mix,    only : randomnumberstream, &
         initialize_random_numbers, uniform_distribution

    implicit none

    ! Number of g points / columns, and number of current column
    integer, intent(in) :: ng, ig

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

    ! Stream for producing random numbers
    type(randomnumberstream), intent(inout) :: random_stream

    ! Object for sampling from a lognormal or gamma distribution
    type(pdf_sampler_type), intent(in) :: pdf_sampler

    ! Cloud fraction, cumulative cloud cover and fractional standard
    ! deviation in each layer
    real(jprb), intent(in), dimension(nlev) :: frac, cum_cloud_cover, fractional_std

    ! Cloud cover of a pair of layers, and amount by which cloud at
    ! next level increases total cloud cover as seen from above
    real(jprb), intent(in), dimension(nlev-1) :: pair_cloud_cover, overhang

    ! Overlap parameter of inhomogeneities
    real(jprb), intent(in), dimension(nlev-1) :: overlap_param_inhom

    ! Top of highest cloudy layer (in this subcolumn) and base of
    ! lowest
    integer, intent(in) :: itrigger, iend

    ! Optical depth scaling to output
    real(jprb), intent(inout), dimension(ng,nlev) :: od_scaling

    ! Height indices
    integer :: jlev, jcloud

    integer :: iy

    real(jprb) :: rand_cloud(nlev)
    real(jprb) :: rand_inhom1(nlev), rand_inhom2(nlev)

    ! For each column analysed, this vector locates the clouds. It is
    ! only actually used for Exp-Exp overlap
    logical :: is_cloudy(nlev)

    ! Number of contiguous cloudy layers for which to compute optical
    ! depth scaling
    integer :: n_layers_to_scale

    iy = 0

    is_cloudy = .false.
    is_cloudy(itrigger) = .true.

    ! Locate the clouds below this layer: first generate some more
    ! random numbers
    call uniform_distribution(rand_cloud(1:(iend+1-itrigger)),random_stream)

    ! Loop from the layer below the local cloud top down to the
    ! bottom-most cloudy layer
    do jlev = itrigger+1,iend
      iy = iy+1
      if (is_cloudy(jlev-1)) then
        ! There is a cloud above, in which case the probability
        ! of cloud in the layer below is as follows
        if (rand_cloud(iy)*frac(jlev-1) &
             &  < frac(jlev) + frac(jlev-1) - pair_cloud_cover(jlev-1)) then
          ! Add another cloudy layer
          is_cloudy(jlev) = .true.
        end if
      else
        ! There is clear-sky above, in which case the
        ! probability of cloud in the layer below is as follows
        if (rand_cloud(iy)*(cum_cloud_cover(jlev-1) - frac(jlev-1)) &
             &  < pair_cloud_cover(jlev-1) - overhang(jlev-1) - frac(jlev-1)) then
            ! A new cloud top
          is_cloudy(jlev) = .true.
        end if
      end if
    end do

    ! We have a contiguous range of layers for which we compute the
    ! od_scaling elements using some random numbers

    ! In the Exp-Exp overlap scheme we do all layers at once
    n_layers_to_scale = iend+1 - itrigger
        
    call uniform_distribution(rand_inhom1(1:n_layers_to_scale),random_stream)
    call uniform_distribution(rand_inhom2(1:n_layers_to_scale),random_stream)
        
    ! Loop through the sequence of cloudy layers
    do jcloud = 2,n_layers_to_scale
      ! Use second random number, and inhomogeneity overlap
      ! parameter, to decide whether the first random number
      ! should be repeated (corresponding to maximum overlap)
      ! or not (corresponding to random overlap)
      if (rand_inhom2(jcloud) &
           &  < overlap_param_inhom(iend-n_layers_to_scale+jcloud-1)) then
        rand_inhom1(jcloud) = rand_inhom1(jcloud-1)
      end if
    end do
        
    ! Sample from a lognormal or gamma distribution to obtain the
    ! optical depth scalings

    ! Masked version assuming values outside the range itrigger:iend
    ! are already zero:
    call pdf_sampler%masked_sample(n_layers_to_scale, &
         &  fractional_std(itrigger:iend), &
         &  rand_inhom1(1:n_layers_to_scale), od_scaling(ig,itrigger:iend), &
         &  is_cloudy(itrigger:iend))
        
    ! ! IFS version:
    ! !$omp simd 
    ! do jlev=itrigger,iend
    !    if (.not. is_cloudy(jlev)) then
    !       od_scaling(ig,jlev) = 0.0_jprb
    !    else
    !       call sample_from_pdf_simd(&
    !            pdf_sampler,fractional_std(jlev),&
    !            rand_inhom1(jlev-itrigger+1), &
    !            od_scaling(ig,jlev))
    !    end if
    ! end do

  end subroutine generate_column_exp_exp


  !---------------------------------------------------------------------
  ! Extract the value of a lognormal distribution with fractional
  ! standard deviation "fsd" corresponding to the cumulative
  ! distribution function value "cdf", and return it in x. Since this
  ! is an elemental subroutine, fsd, cdf and x may be arrays. SIMD version.
  subroutine sample_from_pdf_simd(this, fsd, cdf, x)
    use parkind1,              only : jprb
    use radiation_pdf_sampler, only : pdf_sampler_type
    implicit none
    type(pdf_sampler_type), intent(in)  :: this

    ! Fractional standard deviation (0 to 4) and cumulative
    ! distribution function (0 to 1)
    real(jprb),              intent(in)  :: fsd, cdf

    ! Sample from distribution
    real(jprb),              intent(out) :: x

    ! Index to look-up table
    integer    :: ifsd, icdf

    ! Weights in bilinear interpolation
    real(jprb) :: wfsd, wcdf

    ! Bilinear interpolation with bounds
    wcdf = cdf * (this%ncdf-1) + 1.0_jprb
    icdf = max(1, min(int(wcdf), this%ncdf-1))
    wcdf = max(0.0_jprb, min(wcdf - icdf, 1.0_jprb))

    wfsd = (fsd-this%fsd1) * this%inv_fsd_interval + 1.0_jprb
    ifsd = max(1, min(int(wfsd), this%nfsd-1))
    wfsd = max(0.0_jprb, min(wfsd - ifsd, 1.0_jprb))

    x =      (1.0_jprb-wcdf)*(1.0_jprb-wfsd) * this%val(icdf  ,ifsd)   &
         & + (1.0_jprb-wcdf)*          wfsd  * this%val(icdf  ,ifsd+1) &
         & +           wcdf *(1.0_jprb-wfsd) * this%val(icdf+1,ifsd)   &
         & +           wcdf *          wfsd  * this%val(icdf+1,ifsd+1)

  end subroutine sample_from_pdf_simd


  !---------------------------------------------------------------------
  ! Generate columns of optical depth scalings using
  ! exponential-random overlap (which includes maximum-random overlap
  ! as a limiting case).  This version is intended to work better on
  ! hardware with long vector lengths.  As with all calculations in
  ! this file, we zoom into the fraction of the column with cloud at
  ! any height, so that all spectral intervals see a cloud somewhere.
  ! In the McICA solver, this is combined appropriately with the
  ! clear-sky calculation.
  subroutine generate_columns_exp_ran(ng, nlev, iseed, pdf_sampler, &
       &  total_cloud_cover, frac_threshold, frac, pair_cloud_cover, &
       &  cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, &
       &  ibegin, iend, od_scaling)

    use parkind1,              only : jprb
    use radiation_pdf_sampler, only : pdf_sampler_type
    use radiation_random_numbers, only : rng_type, IRngMinstdVector, IRngNative

    implicit none

    ! Number of g points / columns
    integer, intent(in) :: ng

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

    integer, intent(in) :: iseed ! seed for random number generator

    ! Stream for producing random numbers
    !type(randomnumberstream) :: random_stream
    type(rng_type) :: random_number_generator

    ! Object for sampling from a lognormal or gamma distribution
    type(pdf_sampler_type), intent(in) :: pdf_sampler

    ! Total cloud cover using cloud fraction and overlap parameter
    real(jprb), intent(in) :: total_cloud_cover

    real(jprb), intent(in) :: frac_threshold

    ! Cloud fraction, cumulative cloud cover and fractional standard
    ! deviation in each layer
    real(jprb), intent(in), dimension(nlev) :: frac, cum_cloud_cover, fractional_std

    ! Cloud cover of a pair of layers, and amount by which cloud at
    ! next level increases total cloud cover as seen from above
    real(jprb), intent(in), dimension(nlev-1) :: pair_cloud_cover, overhang

    ! Overlap parameter of inhomogeneities
    real(jprb), intent(in), dimension(nlev-1) :: overlap_param_inhom

    ! Top of highest cloudy layer and base of lowest
    integer, intent(inout) :: ibegin, iend

    ! Optical depth scaling to output
    real(jprb), intent(inout), dimension(ng,nlev) :: od_scaling

    ! Loop indices
    integer :: jlev, jg

    real(jprb) :: rand_cloud(ng,ibegin:iend)
    real(jprb) :: rand_inhom(ng,ibegin-1:iend), rand_inhom2(ng,ibegin:iend)

    ! Is the cloud fraction above the minimum threshold at each level
    logical :: is_any_cloud(ibegin:iend)

    ! Scaled random number for finding cloud
    real(jprb) :: trigger(ng)

    logical :: is_cloud(ng)    ! Is there cloud at this level and spectral interval?
    logical :: prev_cloud(ng)  ! Was there cloud at level above?
    logical :: first_cloud(ng) ! At level of first cloud counting down from top?
    logical :: found_cloud(ng) ! Cloud found in this column counting down from top?

    is_any_cloud = (frac(ibegin:iend) >= frac_threshold)

    ! Initialize random number generator for this column, and state
    ! that random numbers will be requested in blocks of length the
    ! number of spectral intervals ng.
    call random_number_generator%initialize(IRngMinstdVector, iseed=iseed, &
         &                                  nmaxstreams=ng)

    ! Random numbers to use to locate cloud top
    call random_number_generator%uniform_distribution(trigger)

    ! Random numbers to work out whether to transition vertically from
    ! clear to cloudy, cloudy to clear, clear to clear or cloudy to
    ! cloudy
    call random_number_generator%uniform_distribution(rand_cloud, is_any_cloud)

    ! Random numbers to generate sub-grid cloud structure
    call random_number_generator%uniform_distribution(rand_inhom)
    call random_number_generator%uniform_distribution(rand_inhom2, is_any_cloud)

    trigger = trigger * total_cloud_cover

    ! Initialize logicals for clear-sky above first cloudy layer
    found_cloud = .false.
    is_cloud    = .false.
    first_cloud = .false.

    ! Loop down through layers starting at the first cloudy layer
    do jlev = ibegin,iend

      if (is_any_cloud(jlev)) then

! Added for DWD (2020)
!NEC$ shortloop
        do jg = 1,ng
          ! The intention is that all these operations are vectorizable,
          ! since all are vector operations on vectors of length ng...

          ! Copy the cloud mask between levels
          prev_cloud(jg) = is_cloud(jg)

          ! For each spectral interval, has the first cloud appeared at this level?
          first_cloud(jg) = (trigger(jg) <= cum_cloud_cover(jlev) .and. .not. found_cloud(jg))

          ! ...if so, add to found_cloud
          found_cloud(jg) = found_cloud(jg) .or. first_cloud(jg)

          ! There is cloud at this level either if a first cloud has
          ! appeared, or using separate probability calculations
          ! depending on whether there is a cloud above (given by
          ! prev_cloud)
          is_cloud(jg) = first_cloud(jg) &
               &  .or. found_cloud(jg) .and. merge(rand_cloud(jg,jlev)*frac(jlev-1) &
               &               < frac(jlev)+frac(jlev-1)-pair_cloud_cover(jlev-1), &
               &             rand_cloud(jg,jlev)*(cum_cloud_cover(jlev-1) - frac(jlev-1)) &
               &               < pair_cloud_cover(jlev-1) - overhang(jlev-1) - frac(jlev-1), &
               &             prev_cloud(jg))
          ! The random number determining cloud structure decorrelates
          ! with the one above it according to the overlap parameter,
          ! but always decorrelates if there is clear-sky above.  If
          ! there is clear-sky in the present level, the random number
          ! is set to zero to ensure that the optical depth scaling is
          ! also zero.
          rand_inhom(jg,jlev) = merge(merge(rand_inhom(jg,jlev-1), rand_inhom(jg,jlev), &
               &                           rand_inhom2(jg,jlev) < overlap_param_inhom(jlev-1) &
               &                           .and. prev_cloud(jg)), &
               &                     0.0_jprb, is_cloud(jg))
        end do
      else
        ! No cloud at this level
        is_cloud = .false.
      end if
    end do
       
    ! Sample from a lognormal or gamma distribution to obtain the
    ! optical depth scalings, calling the faster masked version and
    ! assuming values outside the range ibegin:iend are already zero
    call pdf_sampler%masked_block_sample(iend-ibegin+1, ng, &
         &  fractional_std(ibegin:iend), &
         &  rand_inhom(:,ibegin:iend), od_scaling(:,ibegin:iend), &
         &  is_any_cloud)

  end subroutine generate_columns_exp_ran

end module radiation_cloud_generator