rrtmgp_lw_main.F90

 
 
module rrtmgp_lw_main
  use mpi_f08
  use machine,                only: kind_phys, kind_dbl_prec
  use mo_optical_props,       only: ty_optical_props_1scl, ty_optical_props_2str
  use mo_rte_lw,              only: rte_lw
  use mo_gas_optics_rrtmgp,   only: ty_gas_optics_rrtmgp
  use mo_gas_concentrations,  only: ty_gas_concs
  use mo_fluxes_byband,       only: ty_fluxes_byband
  use mo_source_functions,    only: ty_source_func_lw
  use radiation_tools,        only: check_error_msg
  use rrtmgp_lw_gas_optics,   only: lw_gas_props,rrtmgp_lw_gas_optics_init
  use rrtmgp_lw_cloud_optics, only: lw_cloud_props, rrtmgp_lw_cloud_optics_init, abssnow0,   &
                                    abssnow1, absrain
  use gfs_rrtmgp_pre,         only: istr_h2o, istr_co2, istr_o3, istr_n2o, istr_ch4,         &
                                    istr_o2, istr_ccl4, istr_cfc11, istr_cfc12, istr_cfc22,  &
                                    eps, oneminus, ftiny
  use mersenne_twister,       only: random_setseed, random_number, random_stat 
  use rrtmgp_sampling,        only: sampled_mask, draw_samples
  use mo_rte_kind,            only: rte_wp => wp
  implicit none
 
  public rrtmgp_lw_main_init, rrtmgp_lw_main_run
contains
 
  subroutine rrtmgp_lw_main_init(rrtmgp_root_dir, rrtmgp_lw_file_gas, rrtmgp_lw_file_clouds,&
       active_gases_array, nrghice, mpicomm, mpirank, mpiroot, nLay, rrtmgp_phys_blksz,     &
       is_init_gas_optics, is_init_cloud_optics, errmsg, errflg)
 
    ! Inputs
    character(len=128),intent(in) :: &
         rrtmgp_root_dir,       & !< RTE-RRTMGP root directory
         rrtmgp_lw_file_clouds, & !< RRTMGP file containing coefficients used to compute
 
         rrtmgp_lw_file_gas       
    character(len=*), dimension(:), intent(in), optional :: &
         active_gases_array
    integer, intent(inout) :: &
         nrghice
    type(mpi_comm),intent(in) :: &
         mpicomm
    integer,intent(in) :: &
         mpirank,               & !< Current MPI rank
         mpiroot,               & !< Master MPI rank
         rrtmgp_phys_blksz,     & !< Number of horizontal points to process at once.
         nlay
 
    ! Outputs
    logical,          intent(inout) :: &
         is_init_gas_optics,    & !< Initialization flag
         is_init_cloud_optics
    character(len=*), intent(  out) :: &
         errmsg
    integer,          intent(  out) :: &
         errflg
 
    ! Initialize CCPP error handling variables 
    errmsg = ''
    errflg = 0
 
    ! RRTMGP longwave gas-optics (k-distribution) initialization
    call rrtmgp_lw_gas_optics_init(rrtmgp_root_dir, rrtmgp_lw_file_gas,                  &
         active_gases_array, mpicomm, mpirank, mpiroot, is_init_gas_optics, errmsg, errflg)
 
    ! RRTMGP longwave cloud-optics initialization
    call rrtmgp_lw_cloud_optics_init(rrtmgp_root_dir, rrtmgp_lw_file_clouds,             &
         nrghice, mpicomm, mpirank, mpiroot, is_init_cloud_optics, errmsg, errflg)
 
  end subroutine rrtmgp_lw_main_init
 
  subroutine rrtmgp_lw_main_run(doLWrad, doLWclrsky, top_at_1, doGP_lwscat,              &
       nCol, nLay, nGases,rrtmgp_phys_blksz,                 &
       nGauss_angles, icseed_lw, iovr, iovr_convcld, iovr_max, iovr_maxrand, iovr_rand,  &
       iovr_dcorr, iovr_exp, iovr_exprand, isubc_lw, semis, tsfg, p_lay, p_lev, t_lay,   &
       t_lev,  vmr_o2, vmr_h2o, vmr_o3, vmr_ch4, vmr_n2o, vmr_co2,                       &
       cld_frac, cld_lwp, cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow,            &
       cld_rwp, cld_rerain, precip_frac, cld_cnv_lwp, cld_cnv_reliq, cld_cnv_iwp,        &
       cld_cnv_reice, cld_pbl_lwp, cld_pbl_reliq, cld_pbl_iwp, cld_pbl_reice,            &
       cloud_overlap_param, active_gases_array, aerlw_tau, aerlw_ssa, aerlw_g,           &
       fluxlwUP_allsky, fluxlwDOWN_allsky, fluxlwUP_clrsky, fluxlwDOWN_clrsky,           &
       fluxlwUP_radtime, fluxlwDOWN_radtime, fluxlwUP_jac, errmsg, errflg)
 
    ! Inputs
    logical, intent(in) :: &
         dolwrad,            & ! Flag to perform longwave calculation
         dolwclrsky,         & ! Flag to compute clear-sky fluxes
         top_at_1,           & ! Flag for vertical ordering convention
         dogp_lwscat           ! Flag to include scattering in clouds
    integer,intent(in) :: &
         ncol,               & ! Number of horizontal points
         nlay,               & ! Number of vertical grid points.
         ngases,             & ! Number of active gases
         rrtmgp_phys_blksz,  & ! Number of horizontal points to process at once.
         ngauss_angles,      & ! Number of gaussian quadrature angles used
         iovr,               & ! Choice of cloud-overlap method
         iovr_convcld,       & ! Choice of convective cloud-overlap
         iovr_max,           & ! Flag for maximum cloud overlap method
         iovr_maxrand,       & ! Flag for maximum-random cloud overlap method
         iovr_rand,          & ! Flag for random cloud overlap method
         iovr_dcorr,         & ! Flag for decorrelation-length cloud overlap method
         iovr_exp,           & ! Flag for exponential cloud overlap method
         iovr_exprand,       & ! Flag for exponential-random cloud overlap method
         isubc_lw              ! Flag for cloud-seeding (rng) for cloud-sampling
    integer,intent(in),dimension(:) :: &
         icseed_lw             ! Seed for random number generation for longwave radiation
    real(kind_phys), dimension(:), intent(in) :: &
         semis,              & ! Surface-emissivity (1)
         tsfg                  ! Skin temperature (K)
    real(kind_phys), dimension(:,:), intent(in) :: &
         p_lay,               & ! Pressure @ model layer-centers (Pa)
         t_lay,               & ! Temperature (K)
         p_lev,               & ! Pressure @ model layer-interfaces (Pa)
         t_lev,               & ! Temperature @ model levels (K)
         vmr_o2,              & ! Molar-mixing ratio oxygen
         vmr_h2o,             & ! Molar-mixing ratio water vapor
         vmr_o3,              & ! Molar-mixing ratio ozone
         vmr_ch4,             & ! Molar-mixing ratio methane
         vmr_n2o,             & ! Molar-mixing ratio nitrous oxide
         vmr_co2                ! Molar-mixing ratio carbon dioxide
    real(kind_phys), dimension(:,:), intent(in) :: &    
         cld_frac,            & ! Cloud-fraction for   stratiform   clouds
         cld_lwp,             & ! Water path for       stratiform   liquid cloud-particles
         cld_reliq,           & ! Effective radius for stratiform   liquid cloud-particles
         cld_iwp,             & ! Water path for       stratiform   ice    cloud-particles
         cld_reice,           & ! Effective radius for stratiform   ice    cloud-particles
         cld_swp,             & ! Water path for                    snow   hydrometeors
         cld_resnow,          & ! Effective radius for              snow   hydrometeors
         cld_rwp,             & ! Water path for                    rain   hydrometeors
         cld_rerain,          & ! Effective radius for              rain   hydrometeors
         precip_frac,         & ! Precipitation fraction (not active, currently precipitation optics uses cloud-fraction)
         cloud_overlap_param    ! Cloud overlap parameter
    real(kind_phys), dimension(:,:), intent(in), optional :: &         
         cld_cnv_lwp,         & ! Water path for       convective   liquid cloud-particles
         cld_cnv_reliq,       & ! Effective radius for convective   liquid cloud-particles
         cld_cnv_iwp,         & ! Water path for       convective   ice    cloud-particles
         cld_cnv_reice,       & ! Effective radius for convective   ice    cloud-particles
         cld_pbl_lwp,         & ! Water path for       PBL          liquid cloud-particles
         cld_pbl_reliq,       & ! Effective radius for PBL          liquid cloud-particles
         cld_pbl_iwp,         & ! Water path for       PBL          ice    cloud-particles
         cld_pbl_reice          ! Effective radius for PBL          ice    cloud-particles
    real(kind_phys), dimension(:,:,:), intent(in) :: &
         aerlw_tau,           & ! Aerosol optical depth
         aerlw_ssa,           & ! Aerosol single scattering albedo
         aerlw_g                ! Aerosol asymmetry paramter
    character(len=*), dimension(:), intent(in) :: &
         active_gases_array     ! List of active gases from namelist as array
 
    ! Outputs
    real(kind_phys), dimension(:,:), intent(inout), optional :: &
         fluxlwup_jac           ! Jacobian of upwelling LW surface radiation (W/m2/K)
    real(kind_phys), dimension(:,:), intent(inout) :: &
         fluxlwup_allsky,     & ! All-sky flux (W/m2)
         fluxlwdown_allsky,   & ! All-sky flux (W/m2)
         fluxlwup_clrsky,     & ! Clear-sky flux (W/m2)
         fluxlwdown_clrsky,   & ! All-sky flux (W/m2)
         fluxlwup_radtime,    & ! Copy of fluxes (Used for coupling)
         fluxlwdown_radtime     !
    character(len=*), intent(out) :: & 
         errmsg                 ! CCPP error message
    integer, intent(out) :: & 
         errflg                 ! CCPP error flag
 
    ! Local variables
    type(ty_fluxes_byband) :: flux_allsky, flux_clrsky
    integer :: icol, ilay, igas, iband, icol2, ix, iblck
    integer, dimension(rrtmgp_phys_blksz) :: ipseed_lw
    type(random_stat) :: rng_stat
    real(kind_phys), dimension(rrtmgp_phys_blksz) :: zcf0, zcf1
    logical, dimension(rrtmgp_phys_blksz,nLay,lw_gas_props%get_ngpt()) :: maskmcica
    real(kind_phys), dimension(rrtmgp_phys_blksz) :: tau_rain, tau_snow
    real(kind_dbl_prec), dimension(lw_gas_props%get_ngpt()) :: rng1d
    real(rte_wp), dimension(lw_gas_props%get_ngpt(),nLay,rrtmgp_phys_blksz) :: rng3d,rng3d2
    real(kind_dbl_prec), dimension(lw_gas_props%get_ngpt()*nLay) :: rng2d
    real(rte_wp), dimension(rrtmgp_phys_blksz,nLay+1,lw_gas_props%get_nband()),target :: &
         fluxlw_up_allsky, fluxlw_up_clrsky, fluxlw_dn_allsky, fluxlw_dn_clrsky
    real(rte_wp), dimension(rrtmgp_phys_blksz,lw_gas_props%get_ngpt()) :: lw_ds
    real(rte_wp), dimension(lw_gas_props%get_nband(),rrtmgp_phys_blksz) :: sfc_emiss_byband
    real(rte_wp), dimension(rrtmgp_phys_blksz,nLay+1) :: fluxlw_up_jac
    logical :: dogp_sgs_cnv, dogp_sgs_pbl
    ! Local RRTMGP DDTs.
    type(ty_gas_concs)          :: gas_concs
    type(ty_optical_props_1scl) :: lw_optical_props_clrsky, lw_optical_props_aerosol_local
    type(ty_optical_props_2str) :: lw_optical_props_clouds, lw_optical_props_cloudsbyband,    &
         lw_optical_props_cnvcloudsbyband, lw_optical_props_pblcloudsbyband,                  &
         lw_optical_props_precipbyband
    type(ty_source_func_lw)     :: sources 
 
    ! Initialize CCPP error handling variables
    errmsg = ''
    errflg = 0
 
    if (.not. dolwrad) return
 
    ! Do we have convective cloud properties?
    dogp_sgs_cnv = .false.
    if (present(cld_cnv_lwp) .and. present(cld_cnv_reliq) .and. &
         present(cld_cnv_iwp) .and. present(cld_cnv_reice)) then
       dogp_sgs_cnv = .true.
    endif
    
    ! Do we have pbl cloud prperties?
    dogp_sgs_pbl = .false.
    if (present(cld_pbl_lwp) .and. present(cld_pbl_reliq) .and. &
         present(cld_pbl_iwp) .and. present(cld_pbl_reice)) then
       dogp_sgs_pbl = .true.
    endif
    
    !
    ! Initialize RRTMGP DDTs (local)
    !
 
    ! ty_gas_concs
    call check_error_msg('rrtmgp_lw_main_gas_concs_run',gas_concs%init(active_gases_array))
 
    ! ty_optical_props
    call check_error_msg('rrtmgp_lw_main_gas_optics_run',&
         lw_optical_props_clrsky%alloc_1scl(rrtmgp_phys_blksz, nlay, lw_gas_props))
    call check_error_msg('rrtmgp_lw_main_sources_run',&
         sources%alloc(rrtmgp_phys_blksz, nlay, lw_gas_props))
    call check_error_msg('rrtmgp_lw_main_cloud_optics_run',&
         lw_optical_props_cloudsbyband%alloc_2str(rrtmgp_phys_blksz, nlay, lw_gas_props%get_band_lims_wavenumber()))
    call check_error_msg('rrtmgp_lw_main_precip_optics_run',&
         lw_optical_props_precipbyband%alloc_2str(rrtmgp_phys_blksz, nlay, lw_gas_props%get_band_lims_wavenumber()))
    call check_error_msg('rrtmgp_lw_mian_cloud_sampling_run', &
         lw_optical_props_clouds%alloc_2str(rrtmgp_phys_blksz, nlay, lw_gas_props))
    call check_error_msg('rrtmgp_lw_main_aerosol_optics_run',&
         lw_optical_props_aerosol_local%alloc_1scl(rrtmgp_phys_blksz, nlay, lw_gas_props%get_band_lims_wavenumber()))
    if (dogp_sgs_cnv) then
       call check_error_msg('rrtmgp_lw_main_cnv_cloud_optics_run',&
            lw_optical_props_cnvcloudsbyband%alloc_2str(rrtmgp_phys_blksz, nlay, lw_gas_props%get_band_lims_wavenumber()))
    endif
    if (dogp_sgs_pbl) then
       call check_error_msg('rrtmgp_lw_main_pbl_cloud_optics_run',&
            lw_optical_props_pblcloudsbyband%alloc_2str(rrtmgp_phys_blksz, nlay, lw_gas_props%get_band_lims_wavenumber()))
    endif
 
    ! ######################################################################################
    !
    ! Loop over all columns...
    !
    ! ###################################################################################### 
    do icol=1,ncol,rrtmgp_phys_blksz
       icol2 = icol + rrtmgp_phys_blksz - 1
 
       ! Initialize/reset
 
       ! ty_optical_props
       lw_optical_props_clrsky%tau       = 0._kind_phys
       lw_optical_props_precipbyband%tau = 0._kind_phys
       lw_optical_props_precipbyband%ssa = 0._kind_phys
       lw_optical_props_precipbyband%g   = 0._kind_phys
       lw_optical_props_cloudsbyband%tau = 0._kind_phys
       lw_optical_props_cloudsbyband%ssa = 0._kind_phys
       lw_optical_props_cloudsbyband%g   = 0._kind_phys
       lw_optical_props_clouds%tau       = 0._kind_phys
       lw_optical_props_clouds%ssa       = 0._kind_phys
       lw_optical_props_clouds%g         = 0._kind_phys
       sources%sfc_source                = 0._kind_phys
       sources%lay_source                = 0._kind_phys
       sources%sfc_source_Jac            = 0._kind_phys
       fluxlw_up_allsky                  = 0._rte_wp
       fluxlw_dn_allsky                  = 0._rte_wp
       fluxlw_up_clrsky                  = 0._rte_wp
       fluxlw_dn_clrsky                  = 0._rte_wp
       if (dogp_sgs_cnv) lw_optical_props_cnvcloudsbyband%tau = 0._kind_phys
       if (dogp_sgs_pbl) lw_optical_props_pblcloudsbyband%tau = 0._kind_phys
 
       ! ty_fluxes_byband
       fluxlw_up_allsky        = 0._rte_wp
       fluxlw_dn_allsky        = 0._rte_wp
       fluxlw_up_clrsky        = 0._rte_wp
       fluxlw_dn_clrsky        = 0._rte_wp
       flux_allsky%bnd_flux_up => fluxlw_up_allsky
       flux_allsky%bnd_flux_dn => fluxlw_dn_allsky
       flux_clrsky%bnd_flux_up => fluxlw_up_clrsky
       flux_clrsky%bnd_flux_dn => fluxlw_dn_clrsky
 
       ! ###################################################################################
       !
       ! Set gas-concentrations
       !
       ! ###################################################################################
       call check_error_msg('rrtmgp_lw_main_set_vmr_o2',  &
            gas_concs%set_vmr(trim(active_gases_array(istr_o2)), real(vmr_o2(icol:icol2,:),kind=rte_wp)))
       call check_error_msg('rrtmgp_lw_main_set_vmr_co2', &
            gas_concs%set_vmr(trim(active_gases_array(istr_co2)),real(vmr_co2(icol:icol2,:),kind=rte_wp)))
       call check_error_msg('rrtmgp_lw_main_set_vmr_ch4', &
            gas_concs%set_vmr(trim(active_gases_array(istr_ch4)),real(vmr_ch4(icol:icol2,:),kind=rte_wp)))
       call check_error_msg('rrtmgp_lw_main_set_vmr_n2o', &
            gas_concs%set_vmr(trim(active_gases_array(istr_n2o)),real(vmr_n2o(icol:icol2,:),kind=rte_wp)))
       call check_error_msg('rrtmgp_lw_main_set_vmr_h2o', &
            gas_concs%set_vmr(trim(active_gases_array(istr_h2o)),real(vmr_h2o(icol:icol2,:),kind=rte_wp)))
       call check_error_msg('rrtmgp_lw_main_set_vmr_o3',  &
            gas_concs%set_vmr(trim(active_gases_array(istr_o3)), real(vmr_o3(icol:icol2,:),kind=rte_wp)))
 
       ! ###################################################################################
       !
       ! Surface emissity in each band
       !
       ! ###################################################################################
       ! Assign same emissivity to all band
       do iblck=1,rrtmgp_phys_blksz
          if (semis(icol+iblck-1) > eps .and. semis(icol+iblck-1) <= 1._kind_phys) then
             do iband=1,lw_gas_props%get_nband()
                sfc_emiss_byband(iband,iblck) = semis(icol+iblck-1)
             enddo
          else
             sfc_emiss_byband(1:lw_gas_props%get_nband(),iblck) = 1.0
          endif
       enddo
 
       ! ###################################################################################
       !
       ! Compute gas-optics...
       !
       ! ###################################################################################
       call check_error_msg('rrtmgp_lw_main_gas_optics',lw_gas_props%gas_optics(&
            real(p_lay(icol:icol2,:),kind=rte_wp),              & ! IN  - Pressure @ layer-centers (Pa)
            real(p_lev(icol:icol2,:),kind=rte_wp),              & ! IN  - Pressure @ layer-interfaces (Pa)
            real(t_lay(icol:icol2,:),kind=rte_wp),              & ! IN  - Temperature @ layer-centers (K)
            real(tsfg(icol:icol2),kind=rte_wp),                 & ! IN  - Skin-temperature (K)
            gas_concs,                        & ! IN  - RRTMGP DDT: trace gas volumne mixing-ratios
            lw_optical_props_clrsky,          & ! OUT - RRTMGP DDT: longwave optical properties
            sources,                          & ! OUT - RRTMGP DDT: source functions
            tlev=real(t_lev(icol:icol2,:),kind=rte_wp)))          ! IN  - Temperature @ layer-interfaces (K) (optional)
 
       ! ###################################################################################
       !
       ! Compute cloud-optics...
       !
       ! ###################################################################################
       ! Create clear/cloudy indicator
       zcf0(:) = 1._kind_phys
       zcf1(:) = 1._kind_phys
       do iblck = 1, rrtmgp_phys_blksz
          do ilay=1,nlay
             zcf0(iblck) = min(zcf0(iblck), 1._kind_phys - cld_frac(icol+iblck-1,ilay))
          enddo
          if (zcf0(iblck) <= ftiny)   zcf0(iblck) = 0._kind_phys
          if (zcf0(iblck) > oneminus) zcf0(iblck) = 1._kind_phys
          zcf1(iblck) = 1._kind_phys - zcf0(iblck)
       enddo
 
       if (any(zcf1 .gt. eps)) then
          ! Microphysical (gridmean) cloud optics
          call check_error_msg('rrtmgp_lw_main_cloud_optics',lw_cloud_props%cloud_optics(&
               real(cld_lwp(icol:icol2,:),kind=rte_wp),                & ! IN  - Cloud liquid water path (g/m2)
               real(cld_iwp(icol:icol2,:),kind=rte_wp),                & ! IN  - Cloud ice water path (g/m2)
               real(cld_reliq(icol:icol2,:),kind=rte_wp),              & ! IN  - Cloud liquid effective radius (microns)
               real(cld_reice(icol:icol2,:),kind=rte_wp),              & ! IN  - Cloud ice effective radius (microns)
               lw_optical_props_cloudsbyband))         ! OUT - RRTMGP DDT containing cloud radiative properties
                                                       !       in each band
          ! Include convective (subgrid scale) clouds?
          if (dogp_sgs_cnv) then
             ! Compute
             call check_error_msg('rrtmgp_lw_main_cnv_cloud_optics',lw_cloud_props%cloud_optics(&
                  real(cld_cnv_lwp(icol:icol2,:),kind=rte_wp),         & ! IN  - Convective cloud liquid water path (g/m2)
                  real(cld_cnv_iwp(icol:icol2,:),kind=rte_wp),         & ! IN  - Convective cloud ice water path (g/m2)
                  real(cld_cnv_reliq(icol:icol2,:),kind=rte_wp),       & ! IN  - Convective cloud liquid effective radius (microns)
                  real(cld_cnv_reice(icol:icol2,:),kind=rte_wp),       & ! IN  - Convective cloud ice effective radius (microns)
                  lw_optical_props_cnvcloudsbyband))   ! OUT - RRTMGP DDT containing convective cloud radiative properties
                                                       !       in each band
             ! Increment
             call check_error_msg('rrtmgp_lw_main_increment_cnvclouds_to_clouds',&
                  lw_optical_props_cnvcloudsbyband%increment(lw_optical_props_cloudsbyband))
          endif
 
          ! Include PBL (subgrid scale) clouds?
          if (dogp_sgs_pbl) then
             ! Compute
             call check_error_msg('rrtmgp_lw_main_pbl_cloud_optics',lw_cloud_props%cloud_optics(&
                  real(cld_pbl_lwp(icol:icol2,:),kind=rte_wp),         & ! IN  - PBL cloud liquid water path (g/m2)
                  real(cld_pbl_iwp(icol:icol2,:),kind=rte_wp),         & ! IN  - PBL cloud ice water path (g/m2)
                  real(cld_pbl_reliq(icol:icol2,:),kind=rte_wp),       & ! IN  - PBL cloud liquid effective radius (microns)
                  real(cld_pbl_reice(icol:icol2,:),kind=rte_wp),       & ! IN  - PBL cloud ice effective radius (microns)
                  lw_optical_props_pblcloudsbyband))   ! OUT - RRTMGP DDT containing PBL cloud radiative properties
                                                       !       in each band
             ! Increment
             call check_error_msg('rrtmgp_lw_main_increment_pblclouds_to_clouds',&
                  lw_optical_props_pblcloudsbyband%increment(lw_optical_props_cloudsbyband))
          endif
       endif
 
       ! ###################################################################################
       !
       ! Cloud precipitation optics: rain and snow(+groupel)
       !
       ! ###################################################################################
       tau_rain(:) = 0._kind_phys
       tau_snow(:) = 0._kind_phys
       do ix=1,rrtmgp_phys_blksz
          do ilay=1,nlay
             if (cld_frac(icol+ix-1,ilay) .gt. eps) then
                ! Rain optical-depth (No band dependence)
                tau_rain(ix) = absrain*cld_rwp(icol+ix-1,ilay)
                
                ! Snow (+groupel) optical-depth (No band dependence)
                if (cld_swp(icol+ix-1,ilay) .gt. 0. .and. cld_resnow(icol+ix-1,ilay) .gt. 10._kind_phys) then
                   tau_snow(ix) = abssnow0*1.05756*cld_swp(icol+ix-1,ilay)/cld_resnow(icol+ix-1,ilay)
                else
                   tau_snow(ix) = 0.0
                endif
                do iband=1,lw_gas_props%get_nband()
                   lw_optical_props_precipbyband%tau(ix,ilay,iband) = tau_rain(ix) + tau_snow(ix)
                enddo
             endif
          enddo
       enddo
       ! Increment
       call check_error_msg('rrtmgp_lw_main_increment_precip_to_clouds',&
            lw_optical_props_precipbyband%increment(lw_optical_props_cloudsbyband))
 
       ! ###################################################################################
       !
       ! Cloud-sampling
       ! *Note* All of the included cloud-types are sampled together, not independently.
       !
       ! ###################################################################################
       if (any(zcf1 .gt. eps)) then
          ! Change random number seed value for each radiation invocation (isubc_lw =1 or 2).
          if(isubc_lw == 1) then      ! advance prescribed permutation seed
             do ix=1,rrtmgp_phys_blksz
                ipseed_lw(ix) = lw_gas_props%get_ngpt() + icol + ix - 1
             enddo
          elseif (isubc_lw == 2) then ! use input array of permutation seeds
             do ix=1,rrtmgp_phys_blksz
                ipseed_lw(ix) = icseed_lw(icol+ix-1)
             enddo
          endif
 
          ! Call RNG
          do ix=1,rrtmgp_phys_blksz
             call random_setseed(ipseed_lw(ix),rng_stat)
             ! Use same rng for each layer
             if (iovr == iovr_max) then
                call random_number(rng1d,rng_stat)
                do ilay=1,nlay
                   rng3d(:,ilay,ix) = rng1d
                enddo
             else
                do ilay=1,nlay
                   call random_number(rng1d,rng_stat)
                   rng3d(:,ilay,ix) = rng1d
                enddo
             endif
          enddo
 
          ! Cloud-overlap.
          ! Maximum-random, random or maximum.
          if (iovr == iovr_maxrand .or. iovr == iovr_rand .or. iovr == iovr_max) then
             call sampled_mask(rng3d, real(cld_frac(icol:icol2,:),kind=rte_wp), maskmcica)
          endif
          ! Exponential decorrelation length overlap
          if (iovr == iovr_dcorr) then
             do ix=1,rrtmgp_phys_blksz
                ! Generate second RNG
                call random_setseed(ipseed_lw(ix),rng_stat)
                call random_number(rng2d,rng_stat)
                rng3d2(:,:,ix) = reshape(source = rng2d,shape=[lw_gas_props%get_ngpt(),nlay])
             enddo
             !
             call sampled_mask(rng3d, real(cld_frac(icol:icol2,:),kind=rte_wp), &
                  maskmcica, overlap_param = real(cloud_overlap_param(icol:icol2,1:nlay-1),kind=rte_wp), &
                  randoms2 = rng3d2)
          endif
          ! Exponential or Exponential-random
          if (iovr == iovr_exp .or. iovr == iovr_exprand) then
             call sampled_mask(rng3d, real(cld_frac(icol:icol2,:),kind=rte_wp), &
                  maskmcica, overlap_param = real(cloud_overlap_param(icol:icol2,1:nlay-1),kind=rte_wp))
          endif
          ! Sampling. Map band optical depth to each g-point using McICA
          call check_error_msg('rrtmgp_lw_main_cloud_sampling',&
               draw_samples(maskmcica, .true., &
               lw_optical_props_cloudsbyband, lw_optical_props_clouds))
       endif
 
       ! ###################################################################################
       !
       ! Compute clear-sky fluxes (gaseous+aerosol) (optional)
       !
       ! ###################################################################################
       ! Increment
       lw_optical_props_aerosol_local%tau = aerlw_tau(icol:icol2,:,:)
       call check_error_msg('rrtmgp_lw_main_increment_aerosol_to_clrsky',&
            lw_optical_props_aerosol_local%increment(lw_optical_props_clrsky))
 
       ! Call RTE solver
       if (dolwclrsky) then
          call check_error_msg('rrtmgp_lw_main_opt_angle',&
               lw_gas_props%compute_optimal_angles(lw_optical_props_clrsky,lw_ds))
          if (ngauss_angles .gt. 1) then
             call check_error_msg('rrtmgp_lw_main_lw_rte_clrsky',rte_lw(           &
                  lw_optical_props_clrsky,         & ! IN  - optical-properties
                  top_at_1,                        & ! IN  - veritcal ordering flag
                  sources,                         & ! IN  - source function
                  sfc_emiss_byband,                & ! IN  - surface emissivity in each LW band
                  flux_clrsky,                     & ! OUT - Fluxes
                  n_gauss_angles = ngauss_angles))   ! IN  - Number of angles in Gaussian quadrature
          else
             call check_error_msg('rrtmgp_lw_main_lw_rte_clrsky',rte_lw(           &
                  lw_optical_props_clrsky,         & ! IN  - optical-properties
                  top_at_1,                        & ! IN  - veritcal ordering flag
                  sources,                         & ! IN  - source function
                  sfc_emiss_byband,                & ! IN  - surface emissivity in each LW band
                  flux_clrsky,                     & ! OUT - Fluxes
                  lw_ds = lw_ds))
          endif
          
          ! Store fluxes
          fluxlwup_clrsky(icol:icol2,:)   = sum(flux_clrsky%bnd_flux_up, dim=3)
          fluxlwdown_clrsky(icol:icol2,:) = sum(flux_clrsky%bnd_flux_dn, dim=3)
       else
          fluxlwup_clrsky(icol:icol2,:)   = 0.0
          fluxlwdown_clrsky(icol:icol2,:) = 0.0   
       endif
 
       ! ###################################################################################
       !
       ! All-sky fluxes (clear-sky + clouds + precipitation)
       ! *Note* CCPP does not allow for polymorphic types, they are ambiguous to the CCPP
       ! framework. rte-rrtmgp uses polymorphic types extensively, for example, querying the
       ! type to determine physics configuration/pathway/etc...
       !
       ! The logic in the code below is to satisfy the polymorphishm in the rte-rrtmgp code.
       ! The rte-rrtmgp "increment" procedures are utilized to provide the correct type to the
       ! rte solver (rte_lw). Rte_lw quieries the type determine if scattering is to be 
       ! included in the calculation. The increment procedures are called so that the correct
       ! optical properties are inherited. ugh...
       ! 
       ! ###################################################################################
 
       ! Include LW cloud-scattering?
       if (dogp_lwscat) then 
          ! Increment
          call check_error_msg('rrtmgp_lw_main_increment_clrsky_to_clouds',&
               lw_optical_props_clrsky%increment(lw_optical_props_clouds))
          
          if (present(fluxlwup_jac)) then
             ! Compute LW Jacobians
             call check_error_msg('rrtmgp_lw_main_lw_rte_allsky',rte_lw(           &
                  lw_optical_props_clouds,         & ! IN  - optical-properties
                  top_at_1,                        & ! IN  - veritcal ordering flag
                  sources,                         & ! IN  - source function
                  sfc_emiss_byband,                & ! IN  - surface emissivity in each LW band
                  flux_allsky,                     & ! OUT - Flxues 
                  n_gauss_angles = ngauss_angles,  & ! IN  - Number of angles in Gaussian quadrature
                  flux_up_jac    = fluxlw_up_jac))   ! OUT - surface temperature flux (upward) Jacobian (W/m2/K)
          else
             call check_error_msg('rrtmgp_lw_main_lw_rte_allsky',rte_lw(           &
                  lw_optical_props_clouds,         & ! IN  - optical-properties
                  top_at_1,                        & ! IN  - veritcal ordering flag
                  sources,                         & ! IN  - source function
                  sfc_emiss_byband,                & ! IN  - surface emissivity in each LW band
                  flux_allsky,                     & ! OUT - Flxues 
                  n_gauss_angles = ngauss_angles))   ! IN  - Number of angles in Gaussian quadrature    
          end if
       ! No scattering in LW clouds.   
       else
          ! Increment
          call check_error_msg('rrtmgp_lw_main_increment_clouds_to_clrsky', &
               lw_optical_props_clouds%increment(lw_optical_props_clrsky))
          
          if (present(fluxlwup_jac)) then
             ! Compute LW Jacobians
             call check_error_msg('rrtmgp_lw_rte_run',rte_lw(           &
                  lw_optical_props_clrsky,         & ! IN  - optical-properties
                  top_at_1,                        & ! IN  - veritcal ordering flag
                  sources,                         & ! IN  - source function
                  sfc_emiss_byband,                & ! IN  - surface emissivity in each LW band
                  flux_allsky,                     & ! OUT - Flxues 
                  n_gauss_angles = ngauss_angles,  & ! IN  - Number of angles in Gaussian quadrature
                  flux_up_jac    = fluxlw_up_jac))   ! OUT - surface temperature flux (upward) Jacobian (W/m2/K)
          else
             call check_error_msg('rrtmgp_lw_rte_run',rte_lw(           &
                  lw_optical_props_clrsky,         & ! IN  - optical-properties
                  top_at_1,                        & ! IN  - veritcal ordering flag
                  sources,                         & ! IN  - source function
                  sfc_emiss_byband,                & ! IN  - surface emissivity in each LW band
                  flux_allsky,                     & ! OUT - Flxues 
                  n_gauss_angles = ngauss_angles))   ! IN  - Number of angles in Gaussian quadrature    
          end if
       endif
       
       ! Store fluxes
       fluxlwup_allsky(icol:icol2,:)   = sum(flux_allsky%bnd_flux_up, dim=3)
       fluxlwdown_allsky(icol:icol2,:) = sum(flux_allsky%bnd_flux_dn, dim=3)
 
       ! Save fluxes for coupling
       if (present(fluxlwup_jac)) then
          fluxlwup_jac(icol:icol2,:)    = fluxlw_up_jac
       endif
       fluxlwup_radtime(icol:icol2,:)   = fluxlwup_allsky(icol:icol2,:)
       fluxlwdown_radtime(icol:icol2,:) = fluxlwdown_allsky(icol:icol2,:)
 
    enddo
 
  end subroutine rrtmgp_lw_main_run
end module rrtmgp_lw_main