rrtmgp_sw_main.F90

 
 
module rrtmgp_sw_main
  use mpi_f08
  use machine,                only: kind_phys, kind_dbl_prec
  use mo_optical_props,       only: ty_optical_props_2str
  use module_radsw_parameters, only: cmpfsw_type
  use mo_rte_sw,              only: rte_sw
  use mo_gas_concentrations,  only: ty_gas_concs
  use mo_fluxes_byband,       only: ty_fluxes_byband
  use radiation_tools,        only: check_error_msg
  use rrtmgp_sw_gas_optics,   only: sw_gas_props,rrtmgp_sw_gas_optics_init
  use rrtmgp_sw_cloud_optics, only: sw_cloud_props, rrtmgp_sw_cloud_optics_init, a0r, a0s,  &
                                    a1s, b0r, b0s, b1s, c0r, c0s
  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_sw_main_init, rrtmgp_sw_main_run
 
contains
 
  subroutine rrtmgp_sw_main_init(rrtmgp_root_dir, rrtmgp_sw_file_gas, rrtmgp_sw_file_clouds,&
       active_gases_array, nrghice, mpicomm, mpirank, mpiroot, nLay, rrtmgp_phys_blksz,     &
       errmsg, errflg)
 
    ! Inputs
    character(len=128),intent(in) :: &
         rrtmgp_root_dir,       & !< RTE-RRTMGP root directory
         rrtmgp_sw_file_clouds, & !< RRTMGP file containing K-distribution data
         rrtmgp_sw_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
    character(len=*), intent(out) :: &
         errmsg
    integer,          intent(out) :: &
         errflg
 
    ! Initialize CCPP error handling variables
    errmsg = ''
    errflg = 0
 
    ! RRTMGP shortwave gas-optics (k-distribution) initialization
    call rrtmgp_sw_gas_optics_init(rrtmgp_root_dir, rrtmgp_sw_file_gas, active_gases_array,&
         mpicomm, mpirank, mpiroot, errmsg, errflg)
 
    ! RRTMGP shortwave cloud-optics initialization
    call rrtmgp_sw_cloud_optics_init(rrtmgp_root_dir, rrtmgp_sw_file_clouds,               &
         nrghice, mpicomm, mpirank, mpiroot, errmsg, errflg)
 
  end subroutine rrtmgp_sw_main_init
 
  subroutine rrtmgp_sw_main_run(doSWrad, doSWclrsky, top_at_1,   &
       nCol, nDay, nLay, nGases, rrtmgp_phys_blksz, idx, icseed_sw, iovr, iovr_convcld,     &
       iovr_max, iovr_maxrand, iovr_rand, iovr_dcorr, iovr_exp, iovr_exprand, isubc_sw,     &
       iSFC, sfc_alb_nir_dir, sfc_alb_nir_dif, sfc_alb_uvvis_dir, sfc_alb_uvvis_dif, coszen,&
       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, aersw_tau, aersw_ssa, aersw_g, solcon, scmpsw,                   &
       fluxswUP_allsky, fluxswDOWN_allsky, fluxswUP_clrsky, fluxswDOWN_clrsky, cldtausw,    &
       errmsg, errflg)
 
    ! Inputs
    logical, intent(in) :: &
         doswrad,             & ! Flag to perform shortwave calculation
         doswclrsky,          & ! Flag to compute clear-sky fluxes
         top_at_1               ! Flag for vertical ordering convention
    integer,intent(in) :: &
         ncol,                & ! Number of horizontal points
         nday,                & ! Number of daytime points
         nlay,                & ! Number of vertical grid points.
         ngases,              & ! Number of active gases
         rrtmgp_phys_blksz,   & ! Number of horizontal points to process at once.
         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_sw,            & !
         isfc
    integer,intent(in),dimension(:) :: &
         idx                    ! Index array for daytime points
    integer,intent(in),dimension(:) :: &
         icseed_sw              ! Seed for random number generation for shortwave radiation
    real(kind_phys), dimension(:), intent(in) :: &
         sfc_alb_nir_dir,     & ! Surface albedo (direct)
         sfc_alb_nir_dif,     & ! Surface albedo (diffuse)
         sfc_alb_uvvis_dir,   & ! Surface albedo (direct)
         sfc_alb_uvvis_dif,   & ! Surface albedo (diffuse)
         coszen                 ! Cosize of SZA
    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
         cloud_overlap_param    !
    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) :: &
         aersw_tau,           & ! Aerosol optical depth
         aersw_ssa,           & ! Aerosol single scattering albedo
         aersw_g                ! Aerosol asymmetry paramter
    character(len=*), dimension(:), intent(in) :: &
         active_gases_array     ! List of active gases from namelist as array
    real(kind_phys), intent(in) :: &
         solcon                 ! Solar constant
 
    ! Outputs
    character(len=*), intent(out) :: &
         errmsg                ! CCPP error message
    integer, intent(out) :: &
         errflg                ! CCPP error flag
    real(kind_phys), dimension(:,:), intent(inout) :: &
         cldtausw              ! Approx 10.mu band layer cloud optical depth  
    real(kind_phys), dimension(:,:), intent(inout) :: &
         fluxswup_allsky,    & ! RRTMGP upward all-sky flux profiles (W/m2)
         fluxswdown_allsky,  & ! RRTMGP downward all-sky flux profiles (W/m2)
         fluxswup_clrsky,    & ! RRTMGP upward clear-sky flux profiles (W/m2)
         fluxswdown_clrsky     ! RRTMGP downward clear-sky flux profiles (W/m2)
    type(cmpfsw_type), dimension(:), intent(inout) :: &
         scmpsw                ! 2D surface fluxes, components:
                               ! uvbfc - total sky downward uv-b flux (W/m2)
                               ! uvbf0 - clear sky downward uv-b flux (W/m2)
                               ! nirbm - downward nir direct beam flux (W/m2)
                               ! nirdf - downward nir diffused flux (W/m2)
                               ! visbm - downward uv+vis direct beam flux (W/m2)
                               ! visdf - downward uv+vis diffused flux (W/m2)
 
    ! Local variables
    type(cmpfsw_type), dimension(rrtmgp_phys_blksz) :: scmpsw_clrsky, scmpsw_allsky
    type(ty_fluxes_byband)      :: flux_allsky, flux_clrsky
    real(kind_phys) :: tau_rain, tau_snow, ssa_rain, ssa_snow, asy_rain, asy_snow, &
         tau_prec, asy_prec, ssa_prec, asyw, ssaw, za1, za2, flux_dir, flux_dif
    real(kind_phys), dimension(rrtmgp_phys_blksz) :: zcf0, zcf1
    real(kind_dbl_prec), dimension(sw_gas_props%get_ngpt()) :: rng1d
    real(rte_wp), dimension(sw_gas_props%get_ngpt(),nLay,rrtmgp_phys_blksz) :: rng3d,rng3d2
    real(kind_dbl_prec), dimension(sw_gas_props%get_ngpt()*nLay) :: rng2d
    logical, dimension(rrtmgp_phys_blksz,nLay,sw_gas_props%get_ngpt()) :: maskmcica
    logical :: cloudy_column, clear_column, dogp_sgs_pbl, dogp_sgs_cnv
    real(rte_wp), dimension(sw_gas_props%get_nband(),rrtmgp_phys_blksz) :: &
         sfc_alb_dir, sfc_alb_dif
    real(rte_wp), dimension(rrtmgp_phys_blksz,nLay+1,sw_gas_props%get_nband()),target :: &
         fluxsw_up_allsky, fluxsw_up_clrsky, fluxsw_dn_dir_clrsky, fluxsw_dn_allsky, &
         fluxsw_dn_clrsky, fluxsw_dn_dir_allsky
    integer :: iband, ibd, ibd_uv, icol, igas, ilay, ix, ix2, iblck
    integer, dimension(rrtmgp_phys_blksz) :: ipseed_sw, icols
    type(random_stat) :: rng_stat
    real(kind_phys), dimension(2,sw_gas_props%get_nband()) :: bandlimits
    real(kind_phys), dimension(2), parameter :: &
         nir_uvvis_bnd = (/12850,16000/), &
         uvb_bnd       = (/29000,38000/)
    real(rte_wp), dimension(rrtmgp_phys_blksz,sw_gas_props%get_ngpt()) :: toa_src_sw
 
    type(ty_gas_concs)          :: gas_concs
    type(ty_optical_props_2str) :: sw_optical_props_accum, sw_optical_props_aerosol_local,    &
         sw_optical_props_cloudsbyband, sw_optical_props_cnvcloudsbyband,                     &
         sw_optical_props_pblcloudsbyband, sw_optical_props_precipbyband,                     &
         sw_optical_props_clouds
 
    ! Initialize CCPP error handling variables
    errmsg = ''
    errflg = 0
 
    if (.not. doswrad) 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
    
    ! ty_gas_concs
    call check_error_msg('rrtmgp_sw_main_gas_concs_init',gas_concs%init(active_gases_array))
 
    ! ty_optical_props
    call check_error_msg('rrtmgp_sw_main_accumulated_optics_init',&
         sw_optical_props_accum%alloc_2str(rrtmgp_phys_blksz, nlay, sw_gas_props))
    call check_error_msg('rrtmgp_sw_main_cloud_optics_init',&
         sw_optical_props_cloudsbyband%alloc_2str(rrtmgp_phys_blksz, nlay, sw_gas_props%get_band_lims_wavenumber()))
    call check_error_msg('rrtmgp_sw_main_precip_optics_init',&
         sw_optical_props_precipbyband%alloc_2str(rrtmgp_phys_blksz, nlay, sw_gas_props%get_band_lims_wavenumber()))
    call check_error_msg('rrtmgp_sw_mian_cloud_sampling_init', &
         sw_optical_props_clouds%alloc_2str(rrtmgp_phys_blksz, nlay, sw_gas_props))
    call check_error_msg('rrtmgp_sw_main_aerosol_optics_init',&
         sw_optical_props_aerosol_local%alloc_2str(rrtmgp_phys_blksz, nlay, sw_gas_props%get_band_lims_wavenumber()))
    if (dogp_sgs_cnv) then
       call check_error_msg('rrtmgp_sw_main_cnv_cloud_optics_init',&
            sw_optical_props_cnvcloudsbyband%alloc_2str(rrtmgp_phys_blksz, nlay, sw_gas_props%get_band_lims_wavenumber()))
    endif
    if (dogp_sgs_pbl) then
       call check_error_msg('rrtmgp_sw_main_pbl_cloud_optics_init',&
            sw_optical_props_pblcloudsbyband%alloc_2str(rrtmgp_phys_blksz, nlay, sw_gas_props%get_band_lims_wavenumber()))
    endif
 
    if (nday .gt. 0) then
 
       bandlimits = sw_gas_props%get_band_lims_wavenumber()
       ! ######################################################################################
       !
       ! Loop over all (daylit) columns...
       !
       ! ######################################################################################
       do icol=1,nday,rrtmgp_phys_blksz
          !ix  = idx(iCol)
          !ix2 = idx(iCol + rrtmgp_phys_blksz - 1)
          icols = idx(icol:icol + rrtmgp_phys_blksz - 1)
 
          ! 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(icols(iblck),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
          cloudy_column = any(zcf1 .gt. eps)
          clear_column  = .true.
          if (cloudy_column) clear_column = .false.
 
          ! ###################################################################################
          !
          ! Initialize/reset
          !
          ! ###################################################################################
          sw_optical_props_clouds%tau             = 0._kind_phys
          sw_optical_props_clouds%ssa             = 0._kind_phys
          sw_optical_props_clouds%g               = 0._kind_phys
          sw_optical_props_accum%tau              = 0._kind_phys
          sw_optical_props_accum%ssa              = 0._kind_phys
          sw_optical_props_accum%g                = 0._kind_phys
          sw_optical_props_cloudsbyband%tau       = 0._kind_phys
          sw_optical_props_cloudsbyband%ssa       = 0._kind_phys
          sw_optical_props_cloudsbyband%g         = 0._kind_phys
          sw_optical_props_precipbyband%tau       = 0._kind_phys
          sw_optical_props_precipbyband%ssa       = 0._kind_phys
          sw_optical_props_precipbyband%g         = 0._kind_phys
          if (dogp_sgs_cnv) then
             sw_optical_props_cnvcloudsbyband%tau = 0._kind_phys
             sw_optical_props_cnvcloudsbyband%ssa = 0._kind_phys
             sw_optical_props_cnvcloudsbyband%g   = 0._kind_phys
          endif
          if (dogp_sgs_pbl) then
             sw_optical_props_pblcloudsbyband%tau = 0._kind_phys
             sw_optical_props_pblcloudsbyband%ssa = 0._kind_phys
             sw_optical_props_pblcloudsbyband%g   = 0._kind_phys
          endif
          scmpsw_clrsky= cmpfsw_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 )
          scmpsw_allsky= cmpfsw_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 )
          cldtausw = 0._kind_phys
 
          ! ty_fluxes_byband
          fluxsw_up_allsky            = 0._rte_wp
          fluxsw_dn_allsky            = 0._rte_wp
          fluxsw_dn_dir_allsky        = 0._rte_wp
          fluxsw_up_clrsky            = 0._rte_wp
          fluxsw_dn_clrsky            = 0._rte_wp
          flux_allsky%bnd_flux_up     => fluxsw_up_allsky
          flux_allsky%bnd_flux_dn     => fluxsw_dn_allsky
          flux_allsky%bnd_flux_dn_dir => fluxsw_dn_dir_allsky
          flux_clrsky%bnd_flux_up     => fluxsw_up_clrsky
          flux_clrsky%bnd_flux_dn     => fluxsw_dn_clrsky
 
          ! ###################################################################################
          !
          ! Set gas-concentrations
          !
          ! ###################################################################################
          call check_error_msg('rrtmgp_sw_main_set_vmr_o2',  &
               gas_concs%set_vmr(trim(active_gases_array(istr_o2)), real(vmr_o2(icols,:),  kind=rte_wp)))
          call check_error_msg('rrtmgp_sw_main_set_vmr_co2', &
               gas_concs%set_vmr(trim(active_gases_array(istr_co2)),real(vmr_co2(icols,:), kind=rte_wp)))
          call check_error_msg('rrtmgp_sw_main_set_vmr_ch4', &
               gas_concs%set_vmr(trim(active_gases_array(istr_ch4)),real(vmr_ch4(icols,:), kind=rte_wp)))
          call check_error_msg('rrtmgp_sw_main_set_vmr_n2o', &
               gas_concs%set_vmr(trim(active_gases_array(istr_n2o)),real(vmr_n2o(icols,:), kind=rte_wp)))
          call check_error_msg('rrtmgp_sw_main_set_vmr_h2o', &
               gas_concs%set_vmr(trim(active_gases_array(istr_h2o)),real(vmr_h2o(icols,:), kind=rte_wp)))
          call check_error_msg('rrtmgp_sw_main_set_vmr_o3',  &
               gas_concs%set_vmr(trim(active_gases_array(istr_o3)), real(vmr_o3(icols,:),  kind=rte_wp)))
 
          ! ###################################################################################
          !
          ! Compute gas-optics
          !
          ! ###################################################################################
 
          call check_error_msg('rrtmgp_sw_main_gas_optics',sw_gas_props%gas_optics(&
               real(p_lay(icols,:), kind=rte_wp),          & ! IN  - Pressure @ layer-centers (Pa)
               real(p_lev(icols,:), kind=rte_wp),          & ! IN  - Pressure @ layer-interfaces (Pa)
               real(t_lay(icols,:), kind=rte_wp),          & ! IN  - Temperature @ layer-centers (K)
               gas_concs,                                  & ! IN  - RRTMGP DDT: trace gas volumne mixing-ratios
               sw_optical_props_accum,                     & ! OUT - RRTMGP DDT: Shortwave optical properties, by
                                                             !                   spectral point (tau,ssa,g)
               toa_src_sw))                                  ! OUT - TOA incident shortwave radiation (spectral)
          ! Scale incident flux
          do iblck = 1, rrtmgp_phys_blksz
             toa_src_sw(iblck,:) = toa_src_sw(iblck,:)*solcon / sum(toa_src_sw(iblck,:))
          enddo
 
          ! ###################################################################################
          !
          ! Set surface albedo
          !
          ! Use near-IR albedo for bands with wavenumbers extending to 12850cm-1
          ! Use uv-vis albedo for bands with wavenumbers greater than 16000cm-1
          ! For overlapping band, average near-IR and us-vis albedos.
          !
          ! ###################################################################################
          do iblck = 1, rrtmgp_phys_blksz
             do iband=1,sw_gas_props%get_nband()
                if (bandlimits(1,iband) .lt. nir_uvvis_bnd(1)) then
                   sfc_alb_dir(iband,iblck) = sfc_alb_nir_dir(icols(iblck))
                   sfc_alb_dif(iband,iblck) = sfc_alb_nir_dif(icols(iblck))
                endif
                if (bandlimits(1,iband) .eq. nir_uvvis_bnd(1)) then
                   sfc_alb_dir(iband,iblck) = 0.5_kind_phys*(sfc_alb_nir_dir(icols(iblck)) +    &
                                                             sfc_alb_uvvis_dir(icols(iblck)))
                   sfc_alb_dif(iband,iblck) = 0.5_kind_phys*(sfc_alb_nir_dif(icols(iblck)) +    &
                                                             sfc_alb_uvvis_dif(icols(iblck)))
                   ibd = iband
                endif
                if (bandlimits(1,iband) .ge. nir_uvvis_bnd(2)) then
                   sfc_alb_dir(iband,iblck) = sfc_alb_uvvis_dir(icols(iblck))
                   sfc_alb_dif(iband,iblck) = sfc_alb_uvvis_dif(icols(iblck))
                endif
                if (bandlimits(1,iband) .eq. uvb_bnd(1)) ibd_uv = iband
             enddo
          enddo
 
          ! ###################################################################################
          !
          ! Compute optics for cloud(s) and precipitation, sample clouds...
          !
          ! ###################################################################################
          if (cloudy_column) then
             ! Gridmean/mp-clouds
             call check_error_msg('rrtmgp_sw_main_cloud_optics',sw_cloud_props%cloud_optics(&
                  real(cld_lwp(icols,:),  kind=rte_wp),         & ! IN  - Cloud liquid water path
                  real(cld_iwp(icols,:),  kind=rte_wp),         & ! IN  - Cloud ice water path
                  real(cld_reliq(icols,:),kind=rte_wp),         & ! IN  - Cloud liquid effective radius
                  real(cld_reice(icols,:),kind=rte_wp),         & ! IN  - Cloud ice effective radius
                  sw_optical_props_cloudsbyband))                 ! OUT - RRTMGP DDT: Shortwave optical properties, 
                                                                  !       in each band (tau,ssa,g)
             cldtausw(icols,:) = sw_optical_props_cloudsbyband%tau(:,:,11)
          
             ! Include convective clouds?
             if (dogp_sgs_cnv) then
                ! Compute
                call check_error_msg('rrtmgp_sw_main_cnv_cloud_optics',sw_cloud_props%cloud_optics(&
                     real(cld_cnv_lwp(icols,:),  kind=rte_wp),  & ! IN  - Convective cloud liquid water path (g/m2)
                     real(cld_cnv_iwp(icols,:),  kind=rte_wp),  & ! IN  - Convective cloud ice water path (g/m2)
                     real(cld_cnv_reliq(icols,:),kind=rte_wp),  & ! IN  - Convective cloud liquid effective radius (microns)
                     real(cld_cnv_reice(icols,:),kind=rte_wp),  & ! IN  - Convective cloud ice effective radius (microns)
                     sw_optical_props_cnvcloudsbyband))           ! OUT - RRTMGP DDT containing convective cloud radiative properties
                                                                  !       in each band
                ! Increment
                call check_error_msg('rrtmgp_sw_main_increment_cnvclouds_to_clouds',&
                     sw_optical_props_cnvcloudsbyband%increment(sw_optical_props_cloudsbyband))
             endif
 
             ! Include PBL clouds?
             if (dogp_sgs_pbl) then
                ! Compute
                call check_error_msg('rrtmgp_sw_main_pbl_cloud_optics',sw_cloud_props%cloud_optics(&
                     real(cld_pbl_lwp(icols,:),   kind=rte_wp), & ! IN  - PBL cloud liquid water path (g/m2)
                     real(cld_pbl_iwp(icols,:),   kind=rte_wp), & ! IN  - PBL cloud ice water path (g/m2)
                     real(cld_pbl_reliq(icols,:), kind=rte_wp), & ! IN  - PBL cloud liquid effective radius (microns)
                     real(cld_pbl_reice(icols,:), kind=rte_wp), & ! IN  - PBL cloud ice effective radius (microns)
                     sw_optical_props_pblcloudsbyband))           ! OUT - RRTMGP DDT containing PBL cloud radiative properties
                                                                  !       in each band
                ! Increment
                call check_error_msg('rrtmgp_sw_main_increment_pblclouds_to_clouds',&
                     sw_optical_props_pblcloudsbyband%increment(sw_optical_props_cloudsbyband))
             endif
          
             ! Cloud precipitation optics: rain and snow(+groupel)
             do iblck = 1, rrtmgp_phys_blksz
                do ilay=1,nlay
                   if (cld_frac(icols(iblck),ilay) .gt. ftiny) then
                      ! Rain/Snow optical depth (No band dependence)
                      tau_rain = cld_rwp(icols(iblck),ilay)*a0r
                      if (cld_swp(icols(iblck),ilay) .gt. 0. .and. cld_resnow(icols(iblck),ilay) .gt. 10._kind_phys) then
                         tau_snow = cld_swp(icols(iblck),ilay)*1.09087*(a0s + a1s/(1.0315*cld_resnow(icols(iblck),ilay)))     ! fu's formula 
                      else
                         tau_snow = 0._kind_phys
                      endif
                      
                      ! Rain/Snow single-scattering albedo and asymmetry (Band dependent)
                      do iband=1,sw_gas_props%get_nband()
                         ! By species
                         ssa_rain = tau_rain*(1.-b0r(iband))
                         asy_rain = ssa_rain*c0r(iband)
                         ssa_snow = tau_snow*(1.-(b0s(iband)+b1s(iband)*1.0315*cld_resnow(icols(iblck),ilay)))
                         asy_snow = ssa_snow*c0s(iband)
                         ! Combine
                         tau_prec = max(1.e-12_kind_phys, tau_rain + tau_snow)
                         ssa_prec = max(1.e-12_kind_phys, ssa_rain + ssa_snow)
                         asy_prec = max(1.e-12_kind_phys, asy_rain + asy_snow)
                         asyw     = asy_prec/max(1.e-12_kind_phys, ssa_prec)
                         ssaw     = min(1._kind_phys-0.000001, ssa_prec/tau_prec)
                         za1      = asyw * asyw
                         za2      = ssaw * za1                      
                         sw_optical_props_precipbyband%tau(iblck,ilay,iband) = (1._kind_phys - za2) * tau_prec
                         sw_optical_props_precipbyband%ssa(iblck,ilay,iband) = (ssaw - za2) / (1._kind_phys - za2)
                         sw_optical_props_precipbyband%g(iblck,ilay,iband)   = asyw/(1+asyw)
                      enddo
                   endif
                enddo
             enddo
             ! Increment 
             call check_error_msg('rrtmgp_sw_main_increment_precip_to_clouds',&
                  sw_optical_props_precipbyband%increment(sw_optical_props_cloudsbyband))
          
             ! ###################################################################################
             !
             ! Cloud-sampling
             !
             ! ###################################################################################
             ! Change random number seed value for each radiation invocation (isubc_sw =1 or 2).
             if(isubc_sw == 1) then      ! advance prescribed permutation seed
                do iblck = 1, rrtmgp_phys_blksz
                   ipseed_sw(iblck) = sw_gas_props%get_ngpt() + icols(iblck)
                enddo
             elseif (isubc_sw == 2) then ! use input array of permutation seeds
                do iblck = 1, rrtmgp_phys_blksz
                   ipseed_sw(iblck) = icseed_sw(icols(iblck))
                enddo
             endif
 
             ! Call RNG
             do iblck = 1, rrtmgp_phys_blksz
                call random_setseed(ipseed_sw(iblck),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,iblck) = rng1d
                   enddo
                else
                   do ilay=1,nlay
                      call random_number(rng1d,rng_stat)
                      rng3d(:,ilay,iblck) = 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(icols,:), kind=rte_wp), maskmcica)
             endif
             ! Exponential decorrelation length overlap
             if (iovr == iovr_dcorr) then
                do iblck = 1, rrtmgp_phys_blksz
                   ! Generate second RNG
                   call random_setseed(ipseed_sw(iblck),rng_stat)
                   call random_number(rng2d,rng_stat)
                   rng3d2(:,:,iblck) = reshape(source = rng2d,shape=[sw_gas_props%get_ngpt(),nlay])
                enddo
                !
                call sampled_mask(rng3d, real(cld_frac(icols,:), kind=rte_wp), &
                     maskmcica, overlap_param = real(cloud_overlap_param(icols,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(icols,:), kind=rte_wp), &
                     maskmcica, overlap_param = real(cloud_overlap_param(icols,1:nlay-1), kind=rte_wp))
             endif
             ! Sampling. Map band optical depth to each g-point using McICA
             call check_error_msg('rrtmgp_sw_main_cloud_sampling',&
                  draw_samples(maskmcica, .true., &
                  sw_optical_props_cloudsbyband, sw_optical_props_clouds))
          endif ! cloudy_column
 
          ! ###################################################################################
          !
          ! Compute clear-sky fluxes (gaseous+aerosol)
          !
          ! ###################################################################################
          ! Increment optics (always)
          sw_optical_props_aerosol_local%tau = aersw_tau(icols,:,:)
          sw_optical_props_aerosol_local%ssa = aersw_ssa(icols,:,:)
          sw_optical_props_aerosol_local%g   = aersw_g(icols,:,:)
          call check_error_msg('rrtmgp_sw_main_increment_aerosol_to_clrsky', & 
               sw_optical_props_aerosol_local%increment(sw_optical_props_accum))
 
          ! Compute clear-sky fluxes (Yes for no-clouds. Optional for cloudy scenes)
          if (clear_column .or. doswclrsky) then
             call check_error_msg('rrtmgp_sw_main_rte_sw_clrsky',rte_sw(     &
                  sw_optical_props_accum,           & ! IN  - optical-properties
                  top_at_1,                         & ! IN  - veritcal ordering flag
                  real(coszen(icols), kind=rte_wp), & ! IN  - Cosine of solar zenith angle
                  toa_src_sw,                       & ! IN  - incident solar flux at TOA
                  sfc_alb_dir,                      & ! IN  - Shortwave surface albedo (direct)
                  sfc_alb_dif,                      & ! IN  - Shortwave surface albedo (diffuse)
                  flux_clrsky))                       ! OUT - Fluxes, clear-sky, 3D (1,nLay,nBand) 
             
             ! Store fluxes
             fluxswup_clrsky(icols,:)   = sum(flux_clrsky%bnd_flux_up, dim=3)
             fluxswdown_clrsky(icols,:) = sum(flux_clrsky%bnd_flux_dn, dim=3)
 
             ! Compute surface downward beam/diffused flux components
             do iblck = 1, rrtmgp_phys_blksz
                do iband=1,sw_gas_props%get_nband()
                   flux_dir = flux_clrsky%bnd_flux_dn(iblck,isfc,iband)
                   flux_dif = 0._kind_phys
                   ! Near-IR bands
                   if (iband < ibd) then
                      scmpsw_clrsky(iblck)%nirbm = scmpsw_clrsky(iblck)%nirbm + flux_dir
                      scmpsw_clrsky(iblck)%nirdf = scmpsw_clrsky(iblck)%nirdf + flux_dif
                   endif
                   ! Transition band
                   if (iband == ibd) then
                      scmpsw_clrsky(iblck)%nirbm = scmpsw_clrsky(iblck)%nirbm + flux_dir*0.5_kind_phys
                      scmpsw_clrsky(iblck)%nirdf = scmpsw_clrsky(iblck)%nirdf + flux_dif*0.5_kind_phys
                      scmpsw_clrsky(iblck)%visbm = scmpsw_clrsky(iblck)%visbm + flux_dir*0.5_kind_phys
                      scmpsw_clrsky(iblck)%visdf = scmpsw_clrsky(iblck)%visdf + flux_dif*0.5_kind_phys
                   endif
                   ! UV-VIS bands
                   if (iband > ibd) then
                      scmpsw_clrsky(iblck)%visbm = scmpsw_clrsky(iblck)%visbm + flux_dir
                      scmpsw_clrsky(iblck)%visdf = scmpsw_clrsky(iblck)%visdf + flux_dif
                   endif
                   ! uv-b surface downward flux
                   scmpsw_clrsky(iblck)%uvbfc    = flux_clrsky%bnd_flux_dn(iblck,isfc,ibd_uv)
                enddo
             enddo
          else
             fluxswup_clrsky(icols,:)   = 0._kind_phys
             fluxswdown_clrsky(icols,:) = 0._kind_phys
             scmpsw                     = cmpfsw_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 )
          endif
 
          ! ###################################################################################
          !
          ! All-sky fluxes (clear-sky + clouds + precipitation)
          !
          ! ###################################################################################
          if (cloudy_column) then
             ! Delta scale
             !call check_error_msg('rrtmgp_sw_main_delta_scale',sw_optical_props_clouds%delta_scale())
 
             ! Increment
             call check_error_msg('rrtmgp_sw_main_increment_clouds_to_clrsky', & 
                  sw_optical_props_clouds%increment(sw_optical_props_accum))
 
             ! Compute fluxes
             call check_error_msg('rrtmgp_sw_main_rte_sw_allsky',rte_sw(     &
                  sw_optical_props_accum,           & ! IN  - optical-properties
                  top_at_1,                         & ! IN  - veritcal ordering flag
                  real(coszen(icols), kind=rte_wp), & ! IN  - Cosine of solar zenith angle
                  toa_src_sw,                       & ! IN  - incident solar flux at TOA
                  sfc_alb_dir,                      & ! IN  - Shortwave surface albedo (direct)
                  sfc_alb_dif,                      & ! IN  - Shortwave surface albedo (diffuse)
                  flux_allsky))                       ! OUT - Fluxes, clear-sky, 3D (1,nLay,nBand)
             
             ! Store fluxes
             fluxswup_allsky(icols,:)   = sum(flux_allsky%bnd_flux_up, dim=3)
             fluxswdown_allsky(icols,:) = sum(flux_allsky%bnd_flux_dn, dim=3)
             
             ! Compute and store downward beam/diffused flux components
             do iblck = 1, rrtmgp_phys_blksz
                ! Loop over bands, sum fluxes...
                do iband=1,sw_gas_props%get_nband()
                   flux_dir = flux_allsky%bnd_flux_dn_dir(iblck,isfc,iband) 
                   flux_dif = flux_allsky%bnd_flux_dn(iblck,isfc,iband) - flux_allsky%bnd_flux_dn_dir(iblck,isfc,iband)
                   ! Near-IR bands
                   if (iband < ibd) then
                      scmpsw_allsky(iblck)%nirbm = scmpsw_allsky(iblck)%nirbm + flux_dir
                      scmpsw_allsky(iblck)%nirdf = scmpsw_allsky(iblck)%nirdf + flux_dif
                   endif
                   ! Transition band
                   if (iband == ibd) then
                      scmpsw_allsky(iblck)%nirbm = scmpsw_allsky(iblck)%nirbm + flux_dir*0.5_kind_phys
                      scmpsw_allsky(iblck)%nirdf = scmpsw_allsky(iblck)%nirdf + flux_dif*0.5_kind_phys
                      scmpsw_allsky(iblck)%visbm = scmpsw_allsky(iblck)%visbm + flux_dir*0.5_kind_phys
                      scmpsw_allsky(iblck)%visdf = scmpsw_allsky(iblck)%visdf + flux_dif*0.5_kind_phys
                   endif
                   ! UV-VIS bands
                   if (iband > ibd) then
                      scmpsw_allsky(iblck)%visbm = scmpsw_allsky(iblck)%visbm + flux_dir
                      scmpsw_allsky(iblck)%visdf = scmpsw_allsky(iblck)%visdf + flux_dif
                   endif
                   ! uv-b surface downward flux 
                   scmpsw_allsky(iblck)%uvbfc    = flux_allsky%bnd_flux_dn(iblck,isfc,ibd_uv)
                enddo
                ! Store surface downward beam/diffused flux components
                if (zcf1(iblck) .gt. eps) then
                   scmpsw(icols(iblck))%nirbm = scmpsw_allsky(iblck)%nirbm
                   scmpsw(icols(iblck))%nirdf = scmpsw_allsky(iblck)%nirdf
                   scmpsw(icols(iblck))%visbm = scmpsw_allsky(iblck)%visbm
                   scmpsw(icols(iblck))%visdf = scmpsw_allsky(iblck)%visdf
                   scmpsw(icols(iblck))%uvbfc = flux_allsky%bnd_flux_dn(iblck,isfc,ibd_uv)
                else
                   scmpsw(icols(iblck))%nirbm = scmpsw_clrsky(iblck)%nirbm
                   scmpsw(icols(iblck))%nirdf = scmpsw_clrsky(iblck)%nirdf
                   scmpsw(icols(iblck))%visbm = scmpsw_clrsky(iblck)%visbm
                   scmpsw(icols(iblck))%visdf = scmpsw_clrsky(iblck)%visdf
                   scmpsw(icols(iblck))%uvbfc = flux_clrsky%bnd_flux_dn(iblck,isfc,ibd_uv)
                endif
                scmpsw(icols(iblck))%uvbf0    = flux_clrsky%bnd_flux_dn(iblck,isfc,ibd_uv)
             enddo
          else ! No clouds
             fluxswup_allsky(icols,:)   = sum(flux_clrsky%bnd_flux_up, dim=3)
             fluxswdown_allsky(icols,:) = sum(flux_clrsky%bnd_flux_dn, dim=3)
             do iblck = 1, rrtmgp_phys_blksz
                scmpsw(icols(iblck))%nirbm = scmpsw_clrsky(iblck)%nirbm
                scmpsw(icols(iblck))%nirdf = scmpsw_clrsky(iblck)%nirdf
                scmpsw(icols(iblck))%visbm = scmpsw_clrsky(iblck)%visbm
                scmpsw(icols(iblck))%visdf = scmpsw_clrsky(iblck)%visdf
                scmpsw(icols(iblck))%uvbfc = flux_clrsky%bnd_flux_dn(iblck,isfc,ibd_uv)
                scmpsw(icols(iblck))%uvbf0 = flux_clrsky%bnd_flux_dn(iblck,isfc,ibd_uv)
             enddo
          endif
          !
       enddo ! nday
    else
       fluxswup_allsky(:,:)   = 0._kind_phys
       fluxswdown_allsky(:,:) = 0._kind_phys
       fluxswup_clrsky(:,:)   = 0._kind_phys
       fluxswdown_clrsky(:,:) = 0._kind_phys
       scmpsw                 = cmpfsw_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 )
    endif
  end subroutine rrtmgp_sw_main_run
end module rrtmgp_sw_main