LCOV - code coverage report
Current view: | top level - physics/rrtmgp - rrtmgp_inputs.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 252 | 299 | 84.3 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 11 | 11 | 100.0 %

              Line data    Source code
           1             : module rrtmgp_inputs
           2             :
           3             : !--------------------------------------------------------------------------------
           4             : ! Transform data for inputs from CAM's data structures to those used by
           5             : ! RRTMGP.  Subset the number of model levels if CAM's top exceeds RRTMGP's
           6             : ! valid domain.  Add an extra layer if CAM's top is below 1 Pa.
           7             : ! The vertical indexing increases from top to bottom of atmosphere in both
           8             : ! CAM and RRTMGP arrays.
           9             : !--------------------------------------------------------------------------------
          10             :
          11             : use shr_kind_mod,     only: r8=>shr_kind_r8
          12             : use ppgrid,           only: pcols, pver, pverp
          13             :
          14             : use physconst,        only: stebol, pi
          15             :
          16             : use physics_types,    only: physics_state
          17             : use physics_buffer,   only: physics_buffer_desc
          18             : use camsrfexch,       only: cam_in_t
          19             :
          20             : use radconstants,     only: nradgas, gaslist, nswbands, nlwbands, nswgpts, nlwgpts,   &
          21             :                             get_sw_spectral_boundaries, idx_sw_diag, idx_sw_cloudsim, &
          22             :                             idx_lw_cloudsim
          23             :
          24             : use rad_constituents, only: rad_cnst_get_gas
          25             :
          26             : use cloud_rad_props,  only: get_liquid_optics_sw, liquid_cloud_get_rad_props_lw, &
          27             :                             get_ice_optics_sw,    ice_cloud_get_rad_props_lw,    &
          28             :                             get_snow_optics_sw,   snow_cloud_get_rad_props_lw,   &
          29             :                             get_grau_optics_sw,   grau_cloud_get_rad_props_lw
          30             :
          31             : use mcica_subcol_gen, only: mcica_subcol_sw, mcica_subcol_lw
          32             :
          33             : use aer_rad_props,    only: aer_rad_props_sw, aer_rad_props_lw
          34             :
          35             : use mo_gas_concentrations, only: ty_gas_concs
          36             : use mo_gas_optics_rrtmgp,  only: ty_gas_optics_rrtmgp
          37             : use mo_optical_props,      only: ty_optical_props_2str, ty_optical_props_1scl
          38             :
          39             : use cam_history_support,   only: fillvalue
          40             : use cam_logfile,      only: iulog
          41             : use cam_abortutils,   only: endrun
          42             : use error_messages,   only: alloc_err
          43             :
          44             : implicit none
          45             : private
          46             : save
          47             :
          48             : public :: &
          49             :    rrtmgp_inputs_init,  &
          50             :    rrtmgp_set_state,    &
          51             :    rrtmgp_set_gases_lw, &
          52             :    rrtmgp_set_gases_sw, &
          53             :    rrtmgp_set_cloud_lw, &
          54             :    rrtmgp_set_cloud_sw, &
          55             :    rrtmgp_set_aer_lw,   &
          56             :    rrtmgp_set_aer_sw
          57             :
          58             :
          59             : ! This value is to match the arbitrary small value used in RRTMG to decide
          60             : ! when a quantity is effectively zero.
          61             : real(r8), parameter :: tiny = 1.0e-80_r8
          62             :
          63             : ! Indices for copying data between cam and rrtmgp arrays
          64             : integer :: ktopcam ! Index in CAM arrays of top level (layer or interface) at which
          65             :                    ! RRTMGP is active.
          66             : integer :: ktoprad ! Index in RRTMGP arrays of the layer or interface corresponding
          67             :                    ! to CAM's top layer or interface
          68             :
          69             : ! wavenumber (cm^-1) boundaries of shortwave bands
          70             : real(r8) :: sw_low_bounds(nswbands), sw_high_bounds(nswbands)
          71             :
          72             : ! Mapping from RRTMG shortwave bands to RRTMGP.  Currently needed to continue using
          73             : ! the SW optics datasets from RRTMG (even thought there is a slight mismatch in the
          74             : ! band boundaries of the 2 bands that overlap with the LW bands).
          75             : integer, parameter, dimension(14) :: rrtmg_to_rrtmgp_swbands = &
          76             :    [ 14, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 ]
          77             :
          78             : !==================================================================================================
          79             : contains
          80             : !==================================================================================================
          81             :
          82        1536 : subroutine rrtmgp_inputs_init(ktcam, ktrad)
          83             :
          84             :    ! Note that this routine must be called after the calls to set_wavenumber_bands which set
          85             :    ! the sw/lw band boundaries in the radconstants module.
          86             :
          87             :    integer, intent(in) :: ktcam
          88             :    integer, intent(in) :: ktrad
          89             :
          90        1536 :    ktopcam = ktcam
          91        1536 :    ktoprad = ktrad
          92             :
          93             :    ! Initialize the module data containing the SW band boundaries.
          94        1536 :    call get_sw_spectral_boundaries(sw_low_bounds, sw_high_bounds, 'cm^-1')
          95             :
          96        1536 : end subroutine rrtmgp_inputs_init
          97             :
          98             : !=========================================================================================
          99             :
         100      746136 : subroutine rrtmgp_set_state( &
         101             :    state, cam_in, ncol, nlay, nday,           &
         102     2238408 :    idxday, coszrs, kdist_sw, t_sfc, emis_sfc,  &
         103      746136 :    t_rad, pmid_rad, pint_rad, t_day, pmid_day, &
         104      746136 :    pint_day, coszrs_day, alb_dir, alb_dif)
         105             :
         106             :    ! arguments
         107             :    type(physics_state),         intent(in) :: state     ! CAM physics state
         108             :    type(cam_in_t),              intent(in) :: cam_in    ! CAM import state
         109             :    integer,                     intent(in) :: ncol      ! # cols in CAM chunk
         110             :    integer,                     intent(in) :: nlay      ! # layers in rrtmgp grid
         111             :    integer,                     intent(in) :: nday      ! # daylight columns
         112             :    integer,                     intent(in) :: idxday(:) ! chunk indicies of daylight columns
         113             :    real(r8),                    intent(in) :: coszrs(:) ! cosine of solar zenith angle
         114             :    class(ty_gas_optics_rrtmgp), intent(in) :: kdist_sw  ! spectral information
         115             :
         116             :    real(r8), intent(out) :: t_sfc(ncol)              ! surface temperature [K]
         117             :    real(r8), intent(out) :: emis_sfc(nlwbands,ncol)  ! emissivity at surface []
         118             :    real(r8), intent(out) :: t_rad(ncol,nlay)         ! layer midpoint temperatures [K]
         119             :    real(r8), intent(out) :: pmid_rad(ncol,nlay)      ! layer midpoint pressures [Pa]
         120             :    real(r8), intent(out) :: pint_rad(ncol,nlay+1)    ! layer interface pressures [Pa]
         121             :    real(r8), intent(out) :: t_day(nday,nlay)         ! layer midpoint temperatures [K]
         122             :    real(r8), intent(out) :: pmid_day(nday,nlay)      ! layer midpoint pressure [Pa]
         123             :    real(r8), intent(out) :: pint_day(nday,nlay+1)    ! layer interface pressures [Pa]
         124             :    real(r8), intent(out) :: coszrs_day(nday)         ! cosine of solar zenith angle
         125             :    real(r8), intent(out) :: alb_dir(nswbands,nday)   ! surface albedo, direct radiation
         126             :    real(r8), intent(out) :: alb_dif(nswbands,nday)   ! surface albedo, diffuse radiation
         127             :
         128             :    ! local variables
         129             :    integer :: i, k, iband
         130             :
         131             :    real(r8) :: tref_min, tref_max
         132             :
         133             :    character(len=*), parameter :: sub='rrtmgp_set_state'
         134             :    character(len=512) :: errmsg
         135             :    !--------------------------------------------------------------------------------
         136             :
         137    12458736 :    t_sfc = sqrt(sqrt(cam_in%lwup(:ncol)/stebol))  ! Surface temp set based on longwave up flux.
         138             :
         139             :    ! Set surface emissivity to 1.0.
         140             :    ! The land model *does* have its own surface emissivity, but is not spectrally resolved.
         141             :    ! The LW upward flux is calculated with that land emissivity, and the "radiative temperature"
         142             :    ! t_sfc is derived from that flux. We assume, therefore, that the emissivity is unity
         143             :    ! to be consistent with t_sfc.
         144   199860336 :    emis_sfc(:,:) = 1._r8
         145             :
         146             :    ! Level ordering is the same for both CAM and RRTMGP (top to bottom)
         147  1159408584 :    t_rad(:,ktoprad:) = state%t(:ncol,ktopcam:)
         148  1159408584 :    pmid_rad(:,ktoprad:) = state%pmid(:ncol,ktopcam:)
         149  1171867320 :    pint_rad(:,ktoprad:) = state%pint(:ncol,ktopcam:)
         150             :
         151             :    ! Add extra layer values if needed.
         152      746136 :    if (nlay == pverp) then
         153    12458736 :       t_rad(:,1) = state%t(:ncol,1)
         154             :       ! The top reference pressure from the RRTMGP coefficients datasets is 1.005183574463 Pa
         155             :       ! Set the top of the extra layer just below that.
         156    12458736 :       pint_rad(:,1) = 1.01_r8
         157             :
         158             :       ! next interface down in LT will always be > 1Pa
         159             :       ! but in MT we apply adjustment to have it be 1.02 Pa if it was too high
         160    12458736 :       where (pint_rad(:,2) <= pint_rad(:,1)) pint_rad(:,2) = pint_rad(:,1)+0.01_r8
         161             :
         162             :       ! set the highest pmid (in the "extra layer") to the midpoint (guarantees > 1Pa)
         163    12458736 :       pmid_rad(:,1)   = pint_rad(:,1) + 0.5_r8 * (pint_rad(:,2) - pint_rad(:,1))
         164             :
         165             :       ! For case of CAM MT, also ensure pint_rad(:,2) > pint_rad(:,1) & pmid_rad(:,2) > max(pmid_rad(:,1), min_pressure)
         166    12458736 :       where (pmid_rad(:,2) <= kdist_sw%get_press_min()) pmid_rad(:,2) = pint_rad(:,2) + 0.01_r8
         167             :    else
         168             :       ! nlay < pverp, thus the 1 Pa level is within a CAM layer.  Assuming the top interface of
         169             :       ! this layer is at a pressure < 1 Pa, we need to adjust the top of this layer so that it
         170             :       ! is within the valid pressure range of RRTMGP (otherwise RRTMGP issues an error).  Then
         171             :       ! set the midpoint pressure halfway between the interfaces.
         172           0 :       pint_rad(:,1) = 1.01_r8
         173           0 :       pmid_rad(:,1) = 0.5_r8 * (pint_rad(:,1) + pint_rad(:,2))
         174             :    end if
         175             :
         176             :    ! Limit temperatures to be within the limits of RRTMGP validity.
         177      746136 :    tref_min = kdist_sw%get_temp_min()
         178      746136 :    tref_max = kdist_sw%get_temp_max()
         179  1171867320 :    t_rad = merge(t_rad, tref_min, t_rad > tref_min)
         180  1171867320 :    t_rad = merge(t_rad, tref_max, t_rad < tref_max)
         181    12458736 :    t_sfc = merge(t_sfc, tref_min, t_sfc > tref_min)
         182    12458736 :    t_sfc = merge(t_sfc, tref_max, t_sfc < tref_max)
         183             :
         184             :    ! Construct arrays containing only daylight columns
         185     6602436 :    do i = 1, nday
         186   556348500 :       t_day(i,:)    = t_rad(idxday(i),:)
         187   556348500 :       pmid_day(i,:) = pmid_rad(idxday(i),:)
         188   562204800 :       pint_day(i,:) = pint_rad(idxday(i),:)
         189     6602436 :       coszrs_day(i) = coszrs(idxday(i))
         190             :    end do
         191             :
         192             :    ! Assign albedos to the daylight columns (from E3SM implementation)
         193             :    ! Albedos are imported from the surface models as broadband (visible, and near-IR),
         194             :    ! and we need to map these to appropriate narrower bands used in RRTMGP. Bands
         195             :    ! are categorized broadly as "visible/UV" or "infrared" based on wavenumber.
         196             :    ! Loop over bands, and determine for each band whether it is broadly in the
         197             :    ! visible or infrared part of the spectrum based on a dividing line of
         198             :    ! 0.7 micron, or 14286 cm^-1
         199    11192040 :    do iband = 1,nswbands
         200    10445904 :       if (is_visible(sw_low_bounds(iband)) .and. &
         201      746136 :          is_visible(sw_high_bounds(iband))) then
         202             :
         203             :          ! Entire band is in the visible
         204    26409744 :          do i = 1, nday
         205    23425200 :             alb_dir(iband,i) = cam_in%asdir(idxday(i))
         206    26409744 :             alb_dif(iband,i) = cam_in%asdif(idxday(i))
         207             :          end do
         208             :
         209     7461360 :       else if (.not.is_visible(sw_low_bounds(iband)) .and. &
         210             :                .not.is_visible(sw_high_bounds(iband))) then
         211             :          ! Entire band is in the longwave (near-infrared)
         212    59421924 :          do i = 1, nday
         213    52706700 :             alb_dir(iband,i) = cam_in%aldir(idxday(i))
         214    59421924 :             alb_dif(iband,i) = cam_in%aldif(idxday(i))
         215             :          end do
         216             :       else
         217             :          ! Band straddles the visible to near-infrared transition, so we take
         218             :          ! the albedo to be the average of the visible and near-infrared
         219             :          ! broadband albedos
         220     6602436 :          do i = 1, nday
         221     5856300 :             alb_dir(iband,i) = 0.5_r8 * (cam_in%aldir(idxday(i)) + cam_in%asdir(idxday(i)))
         222     6602436 :             alb_dif(iband,i) = 0.5_r8 * (cam_in%aldif(idxday(i)) + cam_in%asdif(idxday(i)))
         223             :          end do
         224             :       end if
         225             :    end do
         226             :
         227             :    ! Strictly enforce albedo bounds
         228    88590636 :    where (alb_dir < 0)
         229             :        alb_dir = 0.0_r8
         230             :    end where
         231    88590636 :    where (alb_dir > 1)
         232             :        alb_dir = 1.0_r8
         233             :    end where
         234    88590636 :    where (alb_dif < 0)
         235             :        alb_dif = 0.0_r8
         236             :    end where
         237    88590636 :    where (alb_dif > 1)
         238             :        alb_dif = 1.0_r8
         239             :    end where
         240             :
         241      746136 : end subroutine rrtmgp_set_state
         242             :
         243             : !=========================================================================================
         244             :
         245    20891808 : pure logical function is_visible(wavenumber)
         246             :
         247             :    ! Wavenumber is in the visible if it is above the visible threshold
         248             :    ! wavenumber, and in the infrared if it is below the threshold
         249             :    ! This function doesn't distinquish between visible and UV.
         250             :
         251             :    ! wavenumber in inverse cm (cm^-1)
         252             :    real(r8), intent(in) :: wavenumber
         253             :
         254             :    ! Set threshold between visible and infrared to 0.7 micron, or 14286 cm^-1
         255             :    real(r8), parameter :: visible_wavenumber_threshold = 14286._r8  ! cm^-1
         256             :
         257    20891808 :    if (wavenumber > visible_wavenumber_threshold) then
         258             :       is_visible = .true.
         259             :    else
         260    14176584 :       is_visible = .false.
         261             :    end if
         262             :
         263    20891808 : end function is_visible
         264             :
         265             : !=========================================================================================
         266             :
         267     9083192 : function get_molar_mass_ratio(gas_name) result(massratio)
         268             :
         269             :    ! return the molar mass ratio of dry air to gas based on gas_name
         270             :
         271             :    character(len=*),intent(in) :: gas_name
         272             :    real(r8)                    :: massratio
         273             :
         274             :    ! local variables
         275             :    real(r8), parameter :: amdw = 1.607793_r8    ! Molecular weight of dry air / water vapor
         276             :    real(r8), parameter :: amdc = 0.658114_r8    ! Molecular weight of dry air / carbon dioxide
         277             :    real(r8), parameter :: amdo = 0.603428_r8    ! Molecular weight of dry air / ozone
         278             :    real(r8), parameter :: amdm = 1.805423_r8    ! Molecular weight of dry air / methane
         279             :    real(r8), parameter :: amdn = 0.658090_r8    ! Molecular weight of dry air / nitrous oxide
         280             :    real(r8), parameter :: amdo2 = 0.905140_r8   ! Molecular weight of dry air / oxygen
         281             :    real(r8), parameter :: amdc1 = 0.210852_r8   ! Molecular weight of dry air / CFC11
         282             :    real(r8), parameter :: amdc2 = 0.239546_r8   ! Molecular weight of dry air / CFC12
         283             :
         284             :    character(len=*), parameter :: sub='get_molar_mass_ratio'
         285             :    !----------------------------------------------------------------------------
         286             :
         287    18166384 :    select case (trim(gas_name))
         288             :       case ('H2O')
         289     1135399 :          massratio = amdw
         290             :       case ('CO2')
         291     1135399 :          massratio = amdc
         292             :       case ('O3')
         293     1135399 :          massratio = amdo
         294             :       case ('CH4')
         295     1135399 :          massratio = amdm
         296             :       case ('N2O')
         297     1135399 :          massratio = amdn
         298             :       case ('O2')
         299     1135399 :          massratio = amdo2
         300             :       case ('CFC11')
         301     1135399 :          massratio = amdc1
         302             :       case ('CFC12')
         303     1135399 :          massratio = amdc2
         304             :       case default
         305    18166384 :          call endrun(sub//": Invalid gas: "//trim(gas_name))
         306             :    end select
         307             :
         308     9083192 : end function get_molar_mass_ratio
         309             :
         310             : !=========================================================================================
         311             :
         312     9083192 : subroutine rad_gas_get_vmr(icall, gas_name, state, pbuf, nlay, numactivecols, gas_concs, idxday)
         313             :
         314             :    ! Set volume mixing ratio in gas_concs object.
         315             :    ! The gas_concs%set_vmr method copies data into internally allocated storage.
         316             :
         317             :    integer,                     intent(in) :: icall      ! index of climate/diagnostic radiation call
         318             :    character(len=*),            intent(in) :: gas_name
         319             :    type(physics_state), target, intent(in) :: state
         320             :    type(physics_buffer_desc),   pointer    :: pbuf(:)
         321             :    integer,                     intent(in) :: nlay           ! number of layers in radiation calculation
         322             :    integer,                     intent(in) :: numactivecols  ! number of columns, ncol for LW, nday for SW
         323             :
         324             :    type(ty_gas_concs),       intent(inout) :: gas_concs  ! the result is VRM inside gas_concs
         325             :
         326             :    integer, optional,          intent(in) :: idxday(:)   ! indices of daylight columns in a chunk
         327             :
         328             :    ! Local variables
         329    18166384 :    integer :: i, idx(numactivecols)
         330             :    integer :: istat
         331     9083192 :    real(r8), pointer     :: gas_mmr(:,:)
         332     9083192 :    real(r8), allocatable :: gas_vmr(:,:)
         333     9083192 :    real(r8), allocatable :: mmr(:,:)
         334             :    real(r8) :: massratio
         335             :
         336             :    ! For ozone profile above model
         337             :    real(r8) :: P_top, P_int, P_mid, alpha, beta, a, b, chi_mid, chi_0, chi_eff
         338             :
         339             :    character(len=128)          :: errmsg
         340             :    character(len=*), parameter :: sub = 'rad_gas_get_vmr'
         341             :    !----------------------------------------------------------------------------
         342             :
         343             :    ! set the column indices; when idxday is provided (e.g. daylit columns) use them, otherwise just count.
         344   149634392 :    do i = 1, numactivecols
         345   149634392 :       if (present(idxday)) then
         346    46850400 :          idx(i) = idxday(i)
         347             :       else
         348    93700800 :          idx(i) = i
         349             :       end if
         350             :    end do
         351             :
         352             :    ! gas_mmr points to a "chunk" in either the state or pbuf objects.  That storage is
         353             :    ! dimensioned (pcols,pver).
         354     9083192 :    call rad_cnst_get_gas(icall, gas_name, state, pbuf, gas_mmr)
         355             :
         356             :    ! Copy into storage for RRTMGP
         357    36332768 :    allocate(mmr(numactivecols, nlay), stat=istat)
         358     9083192 :    call alloc_err(istat, sub, 'mmr', numactivecols*nlay)
         359    27249576 :    allocate(gas_vmr(numactivecols, nlay), stat=istat)
         360     9083192 :    call alloc_err(istat, sub, 'gas_vmr', numactivecols*nlay)
         361             :
         362   149634392 :    do i = 1, numactivecols
         363 13220895992 :       mmr(i,ktoprad:) = gas_mmr(idx(i),ktopcam:)
         364             :    end do
         365             :
         366             :    ! If an extra layer is being used, copy mmr from the top layer of CAM to the extra layer.
         367     9083192 :    if (nlay == pverp) then
         368   149634392 :       mmr(:,1) = mmr(:,2)
         369             :    end if
         370             :
         371             :    ! special case: H2O is specific humidity, not mixing ratio. Use r = q/(1-q):
         372     9083192 :    if (gas_name == 'H2O') then
         373  1759339505 :       mmr = mmr / (1._r8 - mmr)
         374             :    end if
         375             :
         376             :    ! convert MMR to VMR, multipy by ratio of dry air molar mas to gas molar mass.
         377     9083192 :    massratio = get_molar_mass_ratio(gas_name)
         378 14083799232 :    gas_vmr = mmr * massratio
         379             :
         380             :    ! special case: Setting O3 in the extra layer:
         381             :    !
         382             :    ! For the purpose of attenuating solar fluxes above the CAM model top, we assume that ozone
         383             :    ! mixing decreases linearly in each column from the value in the top layer of CAM to zero at
         384             :    ! the pressure level set by P_top. P_top has been set to 50 Pa (0.5 hPa) based on model tuning
         385             :    ! to produce temperatures at the top of CAM that are most consistent with WACCM at similar pressure levels.
         386             :
         387     9083192 :    if ((gas_name == 'O3') .and. (nlay == pverp)) then
         388             :       P_top = 50.0_r8
         389    18704299 :       do i = 1, numactivecols
         390    17568900 :             P_int = state%pint(idx(i),1) ! pressure (Pa) at upper interface of CAM
         391    17568900 :             P_mid = state%pmid(idx(i),1) ! pressure (Pa) at midpoint of top layer of CAM
         392    17568900 :             alpha = log(P_int/P_top)
         393    17568900 :             beta =  log(P_mid/P_int)/log(P_mid/P_top)
         394             :
         395    17568900 :             a =  ( (1._r8 + alpha) * exp(-alpha) - 1._r8 ) / alpha
         396    17568900 :             b =  1._r8 - exp(-alpha)
         397             :
         398    18704299 :             if (alpha .gt. 0) then             ! only apply where top level is below 80 km
         399           0 :                chi_mid = gas_vmr(i,1)          ! molar mixing ratio of O3 at midpoint of top layer
         400           0 :                chi_0 = chi_mid /  (1._r8 + beta)
         401           0 :                chi_eff = chi_0 * (a + b)
         402           0 :                gas_vmr(i,1) = chi_eff
         403             :             end if
         404             :       end do
         405             :    end if
         406             :
         407     9083192 :    errmsg = gas_concs%set_vmr(gas_name, gas_vmr)
         408     9083192 :    if (len_trim(errmsg) > 0) then
         409           0 :       call endrun(sub//': ERROR, gas_concs%set_vmr: '//trim(errmsg))
         410             :    end if
         411             :
         412     9083192 :    deallocate(gas_vmr)
         413     9083192 :    deallocate(mmr)
         414             :
         415    18166384 : end subroutine rad_gas_get_vmr
         416             :
         417             : !==================================================================================================
         418             :
         419      746136 : subroutine rrtmgp_set_gases_lw(icall, state, pbuf, nlay, gas_concs)
         420             :
         421             :    ! Set gas vmr for the gases in the radconstants module's gaslist.
         422             :
         423             :    ! The memory management for the gas_concs object is internal.  The arrays passed to it
         424             :    ! are copied to the internally allocated memory.  Each call to the set_vmr method checks
         425             :    ! whether the gas already has memory allocated, and if it does that memory is deallocated
         426             :    ! and new memory is allocated.
         427             :
         428             :    ! arguments
         429             :    integer,                     intent(in)    :: icall      ! index of climate/diagnostic radiation call
         430             :    type(physics_state), target, intent(in)    :: state
         431             :    type(physics_buffer_desc),   pointer       :: pbuf(:)
         432             :    integer,                     intent(in)    :: nlay
         433             :    type(ty_gas_concs),          intent(inout) :: gas_concs
         434             :
         435             :    ! local variables
         436             :    integer :: i, ncol
         437             :    character(len=*), parameter :: sub = 'rrtmgp_set_gases_lw'
         438             :    !--------------------------------------------------------------------------------
         439             :
         440      746136 :    ncol = state%ncol
         441     6715224 :    do i = 1, nradgas
         442     6715224 :       call rad_gas_get_vmr(icall, gaslist(i), state, pbuf, nlay, ncol, gas_concs)
         443             :    end do
         444      746136 : end subroutine rrtmgp_set_gases_lw
         445             :
         446             : !==================================================================================================
         447             :
         448      389263 : subroutine rrtmgp_set_gases_sw( &
         449             :    icall, state, pbuf, nlay, nday, &
         450      389263 :    idxday, gas_concs)
         451             :
         452             :    ! Return gas_concs with gas volume mixing ratio on DAYLIT columns.
         453             :    ! Set all gases in radconstants gaslist.
         454             :
         455             :    ! arguments
         456             :    integer,                     intent(in)    :: icall      ! index of climate/diagnostic radiation call
         457             :    type(physics_state), target, intent(in)    :: state
         458             :    type(physics_buffer_desc),   pointer       :: pbuf(:)
         459             :    integer,                     intent(in)    :: nlay
         460             :    integer,                     intent(in)    :: nday
         461             :    integer,                     intent(in)    :: idxday(:)
         462             :    type(ty_gas_concs),          intent(inout) :: gas_concs
         463             :
         464             :    ! local variables
         465             :    integer :: i
         466             :    character(len=*), parameter :: sub = 'rrtmgp_set_gases_sw'
         467             :    !----------------------------------------------------------------------------
         468             :
         469             :    ! use the optional argument idxday to specify which columns are sunlit
         470     3503367 :     do i = 1,nradgas
         471     3503367 :       call rad_gas_get_vmr(icall, gaslist(i), state, pbuf, nlay, nday, gas_concs, idxday=idxday)
         472             :    end do
         473             :
         474      389263 : end subroutine rrtmgp_set_gases_sw
         475             :
         476             : !==================================================================================================
         477             :
         478      746136 : subroutine rrtmgp_set_cloud_lw( &
         479             :    state, pbuf, ncol, nlay, nlaycam, &
         480             :    cld, cldfsnow, cldfgrau, cldfprime, graupel_in_rad, &
         481             :    kdist_lw, cloud_lw, cld_lw_abs_cloudsim, snow_lw_abs_cloudsim, grau_lw_abs_cloudsim)
         482             :
         483             :    ! Compute combined cloud optical properties.
         484             :    ! Create MCICA stochastic arrays for cloud LW optical properties.
         485             :    ! Initialize optical properties object (cloud_lw) and load with MCICA columns.
         486             :
         487             :    ! arguments
         488             :    type(physics_state),         intent(in)  :: state
         489             :    type(physics_buffer_desc),   pointer     :: pbuf(:)
         490             :    integer,  intent(in) :: ncol           ! number of columns in CAM chunk
         491             :    integer,  intent(in) :: nlay           ! number of layers in radiation calculation (may include "extra layer")
         492             :    integer,  intent(in) :: nlaycam        ! number of CAM layers in radiation calculation
         493             :    real(r8), pointer    :: cld(:,:)       ! cloud fraction (liq+ice)
         494             :    real(r8), pointer    :: cldfsnow(:,:)  ! cloud fraction of just "snow clouds"
         495             :    real(r8), pointer    :: cldfgrau(:,:)  ! cloud fraction of just "graupel clouds"
         496             :    real(r8), intent(in) :: cldfprime(pcols,pver) ! combined cloud fraction
         497             :
         498             :    logical,                     intent(in)  :: graupel_in_rad ! use graupel in radiation code
         499             :    class(ty_gas_optics_rrtmgp), intent(in)  :: kdist_lw
         500             :    type(ty_optical_props_1scl), intent(out) :: cloud_lw
         501             :
         502             :    ! Diagnostic outputs
         503             :    real(r8), intent(out) :: cld_lw_abs_cloudsim(pcols,pver) ! in-cloud liq+ice optical depth (for COSP)
         504             :    real(r8), intent(out) :: snow_lw_abs_cloudsim(pcols,pver)! in-cloud snow optical depth (for COSP)
         505             :    real(r8), intent(out) :: grau_lw_abs_cloudsim(pcols,pver)! in-cloud Graupel optical depth (for COSP)
         506             :
         507             :    ! Local variables
         508             :
         509             :    integer :: i, k
         510             :
         511             :    ! cloud radiative parameters are "in cloud" not "in cell"
         512             :    real(r8) :: liq_lw_abs(nlwbands,pcols,pver)   ! liquid absorption optics depth (LW)
         513             :    real(r8) :: ice_lw_abs(nlwbands,pcols,pver)   ! ice absorption optics depth (LW)
         514             :    real(r8) :: cld_lw_abs(nlwbands,pcols,pver)   ! cloud absorption optics depth (LW)
         515             :    real(r8) :: snow_lw_abs(nlwbands,pcols,pver)  ! snow absorption optics depth (LW)
         516             :    real(r8) :: grau_lw_abs(nlwbands,pcols,pver)  ! graupel absorption optics depth (LW)
         517             :    real(r8) :: c_cld_lw_abs(nlwbands,pcols,pver) ! combined cloud absorption optics depth (LW)
         518             :
         519             :    ! Arrays for converting from CAM chunks to RRTMGP inputs.
         520     1492272 :    real(r8) :: cldf(ncol,nlaycam)
         521     1492272 :    real(r8) :: tauc(nlwbands,ncol,nlaycam)
         522      746136 :    real(r8) :: taucmcl(nlwgpts,ncol,nlaycam)
         523             :
         524             :    character(len=128) :: errmsg
         525             :    character(len=*), parameter :: sub = 'rrtmgp_set_cloud_lw'
         526             :    !--------------------------------------------------------------------------------
         527             :
         528             :    ! Combine the cloud optical properties.  These calculations are done on CAM "chunks".
         529             :
         530             :    ! gammadist liquid optics
         531      746136 :    call liquid_cloud_get_rad_props_lw(state, pbuf, liq_lw_abs)
         532             :    ! Mitchell ice optics
         533      746136 :    call ice_cloud_get_rad_props_lw(state, pbuf, ice_lw_abs)
         534             :
         535 18587757384 :    cld_lw_abs(:,:ncol,:) = liq_lw_abs(:,:ncol,:) + ice_lw_abs(:,:ncol,:)
         536             :
         537      746136 :    if (associated(cldfsnow)) then
         538             :       ! add in snow
         539      746136 :       call snow_cloud_get_rad_props_lw(state, pbuf, snow_lw_abs)
         540    12458736 :       do i = 1, ncol
         541  1101730536 :          do k = 1, pver
         542  1100984400 :             if (cldfprime(i,k) > 0._r8) then
         543           0 :                c_cld_lw_abs(:,i,k) = ( cldfsnow(i,k)*snow_lw_abs(:,i,k) &
         544  3972169782 :                                             + cld(i,k)*cld_lw_abs(:,i,k) )/cldfprime(i,k)
         545             :             else
         546 14545450818 :                c_cld_lw_abs(:,i,k) = 0._r8
         547             :             end if
         548             :          end do
         549             :       end do
         550             :    else
         551           0 :       c_cld_lw_abs(:,:ncol,:) = cld_lw_abs(:,:ncol,:)
         552             :    end if
         553             :
         554             :    ! add in graupel
         555      746136 :    if (associated(cldfgrau) .and. graupel_in_rad) then
         556           0 :       call grau_cloud_get_rad_props_lw(state, pbuf, grau_lw_abs)
         557           0 :       do i = 1, ncol
         558           0 :          do k = 1, pver
         559           0 :             if (cldfprime(i,k) > 0._r8) then
         560           0 :                c_cld_lw_abs(:,i,k) = ( cldfgrau(i,k)*grau_lw_abs(:,i,k) &
         561           0 :                                             + cld(i,k)*c_cld_lw_abs(:,i,k) )/cldfprime(i,k)
         562             :             else
         563           0 :                c_cld_lw_abs(:,i,k) = 0._r8
         564             :             end if
         565             :          end do
         566             :       end do
         567             :    end if
         568             :
         569             :    ! Cloud optics for COSP
         570  1180387152 :    cld_lw_abs_cloudsim  = cld_lw_abs(idx_lw_cloudsim,:,:)
         571  1180387152 :    snow_lw_abs_cloudsim = snow_lw_abs(idx_lw_cloudsim,:,:)
         572  1180387152 :    grau_lw_abs_cloudsim = grau_lw_abs(idx_lw_cloudsim,:,:)
         573             :
         574             :    ! Extract just the layers of CAM where RRTMGP does calculations.
         575             :
         576             :    ! Subset "chunk" data so just the number of CAM layers in the
         577             :    ! radiation calculation are used by MCICA to produce subcolumns.
         578  1159408584 :    cldf = cldfprime(:ncol, ktopcam:)
         579 18587757384 :    tauc = c_cld_lw_abs(:, :ncol, ktopcam:)
         580             :
         581             :    ! Enforce tauc >= 0.
         582 18587757384 :    tauc = merge(tauc, 0.0_r8, tauc > 0.0_r8)
         583             :
         584             :    ! MCICA uses spectral data (on bands) to construct subcolumns (one per g-point)
         585             :    call mcica_subcol_lw( &
         586             :       kdist_lw, nlwbands, nlwgpts, ncol, nlaycam, &
         587      746136 :       nlwgpts, state%pmid, cldf, tauc, taucmcl    )
         588             :
         589      746136 :    errmsg =cloud_lw%alloc_1scl(ncol, nlay, kdist_lw)
         590      746136 :    if (len_trim(errmsg) > 0) then
         591           0 :       call endrun(sub//': ERROR: cloud_lw%alloc_1scalar: '//trim(errmsg))
         592             :    end if
         593             :
         594             :    ! If there is an extra layer in the radiation then this initialization
         595             :    ! will provide zero optical depths there.
         596 >14999*10^7 :    cloud_lw%tau = 0.0_r8
         597             :
         598             :    ! Set the properties on g-points.
         599    96251544 :    do i = 1, nlwgpts
         600 >14840*10^7 :       cloud_lw%tau(:,ktoprad:,i) = taucmcl(i,:,:)
         601             :    end do
         602             :
         603             :    ! validate checks that: tau > 0
         604      746136 :    errmsg = cloud_lw%validate()
         605      746136 :    if (len_trim(errmsg) > 0) then
         606           0 :       call endrun(sub//': ERROR: cloud_lw%validate: '//trim(errmsg))
         607             :    end if
         608             :
         609      746136 : end subroutine rrtmgp_set_cloud_lw
         610             :
         611             : !==================================================================================================
         612             :
         613      746136 : subroutine rrtmgp_set_cloud_sw( &
         614             :    state, pbuf, nlay, nday, idxday, &
         615      746136 :    nnite, idxnite, pmid, cld, cldfsnow, &
         616             :    cldfgrau, cldfprime, graupel_in_rad, kdist_sw, cloud_sw, &
         617             :    tot_cld_vistau, tot_icld_vistau, liq_icld_vistau, ice_icld_vistau, snow_icld_vistau, &
         618             :    grau_icld_vistau, cld_tau_cloudsim, snow_tau_cloudsim, grau_tau_cloudsim)
         619             :
         620             :    ! Compute combined cloud optical properties.
         621             :    ! Create MCICA stochastic arrays for cloud SW optical properties.
         622             :    ! Initialize optical properties object (cloud_sw) and load with MCICA columns.
         623             :
         624             :    ! arguments
         625             :    type(physics_state), target, intent(in) :: state
         626             :    type(physics_buffer_desc), pointer      :: pbuf(:)
         627             :    integer,  intent(in) :: nlay           ! number of layers in radiation calculation (may include "extra layer")
         628             :    integer,  intent(in) :: nday           ! number of daylight columns
         629             :    integer,  intent(in) :: idxday(pcols)  ! indices of daylight columns in the chunk
         630             :    integer,  intent(in) :: nnite          ! number of night columns
         631             :    integer,  intent(in) :: idxnite(pcols) ! indices of night columns in the chunk
         632             :
         633             :    real(r8), intent(in) :: pmid(nday,nlay)! pressure at layer midpoints (Pa) used to seed RNG.
         634             :
         635             :    real(r8), pointer    :: cld(:,:)      ! cloud fraction (liq+ice)
         636             :    real(r8), pointer    :: cldfsnow(:,:) ! cloud fraction of just "snow clouds"
         637             :    real(r8), pointer    :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds"
         638             :    real(r8), intent(in) :: cldfprime(pcols,pver) ! combined cloud fraction
         639             :
         640             :    logical,                     intent(in)  :: graupel_in_rad ! graupel in radiation code
         641             :    class(ty_gas_optics_rrtmgp), intent(in)  :: kdist_sw  ! shortwave gas optics object
         642             :    type(ty_optical_props_2str), intent(out) :: cloud_sw  ! SW cloud optical properties object
         643             :
         644             :    ! Diagnostic outputs
         645             :    real(r8), intent(out) :: tot_cld_vistau(pcols,pver)   ! gbx total cloud optical depth
         646             :    real(r8), intent(out) :: tot_icld_vistau(pcols,pver)  ! in-cld total cloud optical depth
         647             :    real(r8), intent(out) :: liq_icld_vistau(pcols,pver)  ! in-cld liq cloud optical depth
         648             :    real(r8), intent(out) :: ice_icld_vistau(pcols,pver)  ! in-cld ice cloud optical depth
         649             :    real(r8), intent(out) :: snow_icld_vistau(pcols,pver) ! snow in-cloud visible sw optical depth
         650             :    real(r8), intent(out) :: grau_icld_vistau(pcols,pver) ! Graupel in-cloud visible sw optical depth
         651             :    real(r8), intent(out) :: cld_tau_cloudsim(pcols,pver) ! in-cloud liq+ice optical depth (for COSP)
         652             :    real(r8), intent(out) :: snow_tau_cloudsim(pcols,pver)! in-cloud snow optical depth (for COSP)
         653             :    real(r8), intent(out) :: grau_tau_cloudsim(pcols,pver)! in-cloud Graupel optical depth (for COSP)
         654             :
         655             :    ! Local variables
         656             :
         657             :    integer :: i, k, ncol
         658             :    integer :: igpt, nver
         659             :    integer :: istat
         660             :    integer, parameter :: changeseed = 1
         661             :
         662             :    ! cloud radiative parameters are "in cloud" not "in cell"
         663             :    real(r8) :: ice_tau    (nswbands,pcols,pver) ! ice extinction optical depth
         664             :    real(r8) :: ice_tau_w  (nswbands,pcols,pver) ! ice single scattering albedo * tau
         665             :    real(r8) :: ice_tau_w_g(nswbands,pcols,pver) ! ice asymmetry parameter * tau * w
         666             :    real(r8) :: liq_tau    (nswbands,pcols,pver) ! liquid extinction optical depth
         667             :    real(r8) :: liq_tau_w  (nswbands,pcols,pver) ! liquid single scattering albedo * tau
         668             :    real(r8) :: liq_tau_w_g(nswbands,pcols,pver) ! liquid asymmetry parameter * tau * w
         669             :    real(r8) :: cld_tau    (nswbands,pcols,pver) ! cloud extinction optical depth
         670             :    real(r8) :: cld_tau_w  (nswbands,pcols,pver) ! cloud single scattering albedo * tau
         671             :    real(r8) :: cld_tau_w_g(nswbands,pcols,pver) ! cloud asymmetry parameter * w * tau
         672             :    real(r8) :: snow_tau    (nswbands,pcols,pver) ! snow extinction optical depth
         673             :    real(r8) :: snow_tau_w  (nswbands,pcols,pver) ! snow single scattering albedo * tau
         674             :    real(r8) :: snow_tau_w_g(nswbands,pcols,pver) ! snow asymmetry parameter * tau * w
         675             :    real(r8) :: grau_tau    (nswbands,pcols,pver) ! graupel extinction optical depth
         676             :    real(r8) :: grau_tau_w  (nswbands,pcols,pver) ! graupel single scattering albedo * tau
         677             :    real(r8) :: grau_tau_w_g(nswbands,pcols,pver) ! graupel asymmetry parameter * tau * w
         678             :    real(r8) :: c_cld_tau    (nswbands,pcols,pver) ! combined cloud extinction optical depth
         679             :    real(r8) :: c_cld_tau_w  (nswbands,pcols,pver) ! combined cloud single scattering albedo * tau
         680             :    real(r8) :: c_cld_tau_w_g(nswbands,pcols,pver) ! combined cloud asymmetry parameter * w * tau
         681             :
         682             :    ! RRTMGP does not use this property in its 2-stream calculations.
         683             :    real(r8) :: sw_tau_w_f(nswbands,pcols,pver) ! Forward scattered fraction * tau * w.
         684             :
         685             :    ! Arrays for converting from CAM chunks to RRTMGP inputs.
         686      746136 :    real(r8), allocatable :: cldf(:,:)
         687      746136 :    real(r8), allocatable :: tauc(:,:,:)
         688      746136 :    real(r8), allocatable :: ssac(:,:,:)
         689      746136 :    real(r8), allocatable :: asmc(:,:,:)
         690      746136 :    real(r8), allocatable :: taucmcl(:,:,:)
         691      746136 :    real(r8), allocatable :: ssacmcl(:,:,:)
         692      746136 :    real(r8), allocatable :: asmcmcl(:,:,:)
         693      746136 :    real(r8), allocatable :: day_cld_tau(:,:,:)
         694      746136 :    real(r8), allocatable :: day_cld_tau_w(:,:,:)
         695      746136 :    real(r8), allocatable :: day_cld_tau_w_g(:,:,:)
         696             :
         697             :    character(len=128) :: errmsg
         698             :    character(len=*), parameter :: sub = 'rrtmgp_set_cloud_sw'
         699             :    !--------------------------------------------------------------------------------
         700             :
         701      746136 :    ncol = state%ncol
         702             :
         703             :    ! Combine the cloud optical properties.  These calculations are done on CAM "chunks".
         704             :
         705             :    ! gammadist liquid optics
         706      746136 :    call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, sw_tau_w_f)
         707             :    ! Mitchell ice optics
         708      746136 :    call get_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, sw_tau_w_f)
         709             :
         710 16409213784 :    cld_tau(:,:ncol,:)     =  liq_tau(:,:ncol,:)     + ice_tau(:,:ncol,:)
         711 16409213784 :    cld_tau_w(:,:ncol,:)   =  liq_tau_w(:,:ncol,:)   + ice_tau_w(:,:ncol,:)
         712 16409213784 :    cld_tau_w_g(:,:ncol,:) =  liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:)
         713             :
         714             :    ! add in snow
         715      746136 :    if (associated(cldfsnow)) then
         716      746136 :       call get_snow_optics_sw(state, pbuf, snow_tau, snow_tau_w, snow_tau_w_g, sw_tau_w_f)
         717    12458736 :       do i = 1, ncol
         718  1101730536 :          do k = 1, pver
         719  1100984400 :             if (cldfprime(i,k) > 0._r8) then
         720           0 :                c_cld_tau(:,i,k)     = ( cldfsnow(i,k)*snow_tau(:,i,k) &
         721  3504855690 :                                       + cld(i,k)*cld_tau(:,i,k) )/cldfprime(i,k)
         722           0 :                c_cld_tau_w(:,i,k)   = ( cldfsnow(i,k)*snow_tau_w(:,i,k)  &
         723  3504855690 :                                       + cld(i,k)*cld_tau_w(:,i,k) )/cldfprime(i,k)
         724           0 :                c_cld_tau_w_g(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_g(:,i,k) &
         725  3504855690 :                                       + cld(i,k)*cld_tau_w_g(:,i,k) )/cldfprime(i,k)
         726             :             else
         727 12834221310 :                c_cld_tau(:,i,k)     = 0._r8
         728 12834221310 :                c_cld_tau_w(:,i,k)   = 0._r8
         729 12834221310 :                c_cld_tau_w_g(:,i,k) = 0._r8
         730             :             end if
         731             :          end do
         732             :       end do
         733             :    else
         734           0 :       c_cld_tau(:,:ncol,:)     = cld_tau(:,:ncol,:)
         735           0 :       c_cld_tau_w(:,:ncol,:)   = cld_tau_w(:,:ncol,:)
         736           0 :       c_cld_tau_w_g(:,:ncol,:) = cld_tau_w_g(:,:ncol,:)
         737             :    end if
         738             :
         739             :    ! add in graupel
         740      746136 :    if (associated(cldfgrau) .and. graupel_in_rad) then
         741           0 :       call get_grau_optics_sw(state, pbuf, grau_tau, grau_tau_w, grau_tau_w_g, sw_tau_w_f)
         742           0 :       do i = 1, ncol
         743           0 :          do k = 1, pver
         744           0 :             if (cldfprime(i,k) > 0._r8) then
         745           0 :                c_cld_tau(:,i,k)     = ( cldfgrau(i,k)*grau_tau(:,i,k) &
         746           0 :                                       + cld(i,k)*c_cld_tau(:,i,k) )/cldfprime(i,k)
         747           0 :                c_cld_tau_w(:,i,k)   = ( cldfgrau(i,k)*grau_tau_w(:,i,k)  &
         748           0 :                                       + cld(i,k)*c_cld_tau_w(:,i,k) )/cldfprime(i,k)
         749           0 :                c_cld_tau_w_g(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_g(:,i,k) &
         750           0 :                                       + cld(i,k)*c_cld_tau_w_g(:,i,k) )/cldfprime(i,k)
         751             :             else
         752           0 :                c_cld_tau(:,i,k)     = 0._r8
         753           0 :                c_cld_tau_w(:,i,k)   = 0._r8
         754           0 :                c_cld_tau_w_g(:,i,k) = 0._r8
         755             :             end if
         756             :          end do
         757             :       end do
         758             :    end if
         759             :
         760             :    ! cloud optical properties need to be re-ordered from the RRTMG spectral bands
         761             :    ! (assumed in the optics datasets) to RRTMGP's
         762 32817681432 :    ice_tau(:,:ncol,:)       = ice_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
         763 32817681432 :    liq_tau(:,:ncol,:)       = liq_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
         764 32817681432 :    c_cld_tau(:,:ncol,:)     = c_cld_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
         765 32817681432 :    c_cld_tau_w(:,:ncol,:)   = c_cld_tau_w(rrtmg_to_rrtmgp_swbands,:ncol,:)
         766 32817681432 :    c_cld_tau_w_g(:,:ncol,:) = c_cld_tau_w_g(rrtmg_to_rrtmgp_swbands,:ncol,:)
         767      746136 :    if (associated(cldfsnow)) then
         768 32817681432 :       snow_tau(:,:ncol,:)   = snow_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
         769             :    end if
         770      746136 :    if (associated(cldfgrau) .and. graupel_in_rad) then
         771           0 :       grau_tau(:,:ncol,:)   = grau_tau(rrtmg_to_rrtmgp_swbands,:ncol,:)
         772             :    end if
         773             :
         774             :    ! Set arrays for diagnostic output.
         775             :    ! cloud optical depth fields for the visible band
         776  1159408584 :    tot_icld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)
         777  1159408584 :    liq_icld_vistau(:ncol,:) = liq_tau(idx_sw_diag,:ncol,:)
         778  1159408584 :    ice_icld_vistau(:ncol,:) = ice_tau(idx_sw_diag,:ncol,:)
         779      746136 :    if (associated(cldfsnow)) then
         780  1159408584 :       snow_icld_vistau(:ncol,:) = snow_tau(idx_sw_diag,:ncol,:)
         781             :    endif
         782      746136 :    if (associated(cldfgrau) .and. graupel_in_rad) then
         783           0 :       grau_icld_vistau(:ncol,:) = grau_tau(idx_sw_diag,:ncol,:)
         784             :    endif
         785             :
         786             :    ! multiply by total cloud fraction to get gridbox value
         787  1159408584 :    tot_cld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)*cldfprime(:ncol,:)
         788             :
         789             :    ! overwrite night columns with fillvalue
         790     6602436 :    do i = 1, Nnite
         791   550492200 :       tot_cld_vistau(IdxNite(i),:)   = fillvalue
         792   550492200 :       tot_icld_vistau(IdxNite(i),:)  = fillvalue
         793   550492200 :       liq_icld_vistau(IdxNite(i),:)  = fillvalue
         794   550492200 :       ice_icld_vistau(IdxNite(i),:)  = fillvalue
         795     5856300 :       if (associated(cldfsnow)) then
         796   550492200 :          snow_icld_vistau(IdxNite(i),:) = fillvalue
         797             :       end if
         798     6602436 :       if (associated(cldfgrau) .and. graupel_in_rad) then
         799           0 :          grau_icld_vistau(IdxNite(i),:) = fillvalue
         800             :       end if
         801             :    end do
         802             :
         803             :    ! Cloud optics for COSP
         804  1180387152 :    cld_tau_cloudsim = cld_tau(idx_sw_cloudsim,:,:)
         805  1180387152 :    snow_tau_cloudsim = snow_tau(idx_sw_cloudsim,:,:)
         806  1180387152 :    grau_tau_cloudsim = grau_tau(idx_sw_cloudsim,:,:)
         807             :
         808             :    ! if no daylight columns the cloud_sw object isn't initialized
         809      746136 :    if (nday > 0) then
         810             :
         811             :       ! number of CAM's layers in radiation calculation.  Does not include the "extra layer".
         812      389263 :       nver = pver - ktopcam + 1
         813             :
         814             :       allocate( &
         815             :          cldf(nday,nver),                      &
         816             :          day_cld_tau(nswbands,nday,nver),      &
         817             :          day_cld_tau_w(nswbands,nday,nver),    &
         818             :          day_cld_tau_w_g(nswbands,nday,nver),  &
         819             :          tauc(nswbands,nday,nver), taucmcl(nswgpts,nday,nver), &
         820             :          ssac(nswbands,nday,nver), ssacmcl(nswgpts,nday,nver), &
         821    10510101 :          asmc(nswbands,nday,nver), asmcmcl(nswgpts,nday,nver), stat=istat)
         822      389263 :       call alloc_err(istat, sub, 'cldf,..,asmcmcl', 9*nswgpts*nday*nver)
         823             :
         824             :       ! Subset "chunk" data so just the daylight columns, and the number of CAM layers in the
         825             :       ! radiation calculation are used by MCICA to produce subcolumns.
         826   581615885 :       cldf            = cldfprime(       idxday(1:nday), ktopcam:)
         827  8206518485 :       day_cld_tau     = c_cld_tau(    :, idxday(1:nday), ktopcam:)
         828  8206518485 :       day_cld_tau_w   = c_cld_tau_w(  :, idxday(1:nday), ktopcam:)
         829  8206518485 :       day_cld_tau_w_g = c_cld_tau_w_g(:, idxday(1:nday), ktopcam:)
         830             :
         831             :       ! Compute the optical properties needed for the 2-stream calculations.  These calculations
         832             :       ! are the same as the RRTMG version.
         833             :
         834             :       ! set cloud optical depth, clip @ zero
         835  8206518485 :       tauc = merge(day_cld_tau, 0.0_r8, day_cld_tau > 0.0_r8)
         836             :       ! set value of asymmetry
         837  8206518485 :       asmc = merge(day_cld_tau_w_g / max(day_cld_tau_w, tiny), 0.0_r8, day_cld_tau_w > 0.0_r8)
         838             :       ! set value of single scattering albedo
         839  8206518485 :       ssac = merge(max(day_cld_tau_w, tiny) / max(tauc, tiny), 1.0_r8 , tauc > 0.0_r8)
         840             :       ! set asymmetry to zero when tauc = 0
         841  8206518485 :       asmc = merge(asmc, 0.0_r8, tauc > 0.0_r8)
         842             :
         843             :       ! MCICA uses spectral data (on bands) to construct subcolumns (one per g-point)
         844             :       call mcica_subcol_sw( &
         845             :          kdist_sw, nswbands, nswgpts, nday, nlay, &
         846             :          nver, changeseed, pmid, cldf, tauc,     &
         847      389263 :          ssac, asmc, taucmcl, ssacmcl, asmcmcl)
         848             :
         849             :       ! Initialize object for SW cloud optical properties.
         850      389263 :       errmsg = cloud_sw%alloc_2str(nday, nlay, kdist_sw)
         851      389263 :       if (len_trim(errmsg) > 0) then
         852           0 :          call endrun(trim(sub)//': ERROR: cloud_sw%alloc_2str: '//trim(errmsg))
         853             :       end if
         854             :
         855             :       ! If there is an extra layer in the radiation then this initialization
         856             :       ! will provide the optical properties there.
         857 65797273983 :       cloud_sw%tau = 0.0_r8
         858 65797273983 :       cloud_sw%ssa = 1.0_r8
         859 65797273983 :       cloud_sw%g   = 0.0_r8
         860             :
         861             :       ! Set the properties on g-points.
         862    43986719 :       do igpt = 1,nswgpts
         863 65097381664 :          cloud_sw%g  (:,ktoprad:,igpt) = asmcmcl(igpt,:,:)
         864 65097381664 :          cloud_sw%ssa(:,ktoprad:,igpt) = ssacmcl(igpt,:,:)
         865 65097770927 :          cloud_sw%tau(:,ktoprad:,igpt) = taucmcl(igpt,:,:)
         866             :       end do
         867             :
         868             :       ! validate checks that: tau > 0, ssa is in range [0,1], and g is in range [-1,1].
         869      389263 :       errmsg = cloud_sw%validate()
         870      389263 :       if (len_trim(errmsg) > 0) then
         871           0 :          call endrun(sub//': ERROR: cloud_sw%validate: '//trim(errmsg))
         872             :       end if
         873             :
         874             :       ! delta scaling adjusts for forward scattering
         875      389263 :       errmsg = cloud_sw%delta_scale()
         876      389263 :       if (len_trim(errmsg) > 0) then
         877           0 :          call endrun(sub//': ERROR: cloud_sw%delta_scale: '//trim(errmsg))
         878             :       end if
         879             :
         880             :       ! All information is in cloud_sw, now deallocate local vars.
         881           0 :       deallocate( &
         882           0 :          cldf, tauc, ssac, asmc, &
         883           0 :          taucmcl, ssacmcl, asmcmcl,&
         884      389263 :          day_cld_tau, day_cld_tau_w, day_cld_tau_w_g )
         885             :
         886             :    end if
         887             :
         888      746136 : end subroutine rrtmgp_set_cloud_sw
         889             :
         890             : !==================================================================================================
         891             :
         892      746136 : subroutine rrtmgp_set_aer_lw(icall, state, pbuf, aer_lw)
         893             :
         894             :    ! Load LW aerosol optical properties into the RRTMGP object.
         895             :
         896             :    ! Arguments
         897             :    integer,                     intent(in) :: icall
         898             :    type(physics_state), target, intent(in) :: state
         899             :    type(physics_buffer_desc),   pointer    :: pbuf(:)
         900             :
         901             :    type(ty_optical_props_1scl), intent(inout) :: aer_lw
         902             :
         903             :    ! Local variables
         904             :    integer :: ncol
         905             :
         906             :    ! Aerosol LW absorption optical depth
         907             :    real(r8) :: aer_lw_abs (pcols,pver,nlwbands)
         908             :
         909             :    character(len=*), parameter :: sub = 'rrtmgp_set_aer_lw'
         910             :    character(len=128) :: errmsg
         911             :    !--------------------------------------------------------------------------------
         912             :
         913      746136 :    ncol = state%ncol
         914             :
         915             :    ! Get aerosol longwave optical properties.
         916      746136 :    call aer_rad_props_lw(icall, state, pbuf, aer_lw_abs)
         917             :
         918             :    ! If there is an extra layer in the radiation then this initialization
         919             :    ! will provide zero optical depths there.
         920 18750623256 :    aer_lw%tau = 0.0_r8
         921             :
         922 18551283480 :    aer_lw%tau(:ncol, ktoprad:, :) = aer_lw_abs(:ncol, ktopcam:, :)
         923             :
         924      746136 :    errmsg = aer_lw%validate()
         925      746136 :    if (len_trim(errmsg) > 0) then
         926           0 :       call endrun(sub//': ERROR: aer_lw%validate: '//trim(errmsg))
         927             :    end if
         928      746136 : end subroutine rrtmgp_set_aer_lw
         929             :
         930             : !==================================================================================================
         931             :
         932      746136 : subroutine rrtmgp_set_aer_sw( &
         933      746136 :    icall, state, pbuf, nday, idxday, nnite, idxnite, aer_sw)
         934             :
         935             :    ! Load SW aerosol optical properties into the RRTMGP object.
         936             :
         937             :    ! Arguments
         938             :    integer,                     intent(in) :: icall
         939             :    type(physics_state), target, intent(in) :: state
         940             :    type(physics_buffer_desc), pointer      :: pbuf(:)
         941             :    integer,  intent(in) :: nday
         942             :    integer,  intent(in) :: idxday(:)
         943             :    integer,  intent(in) :: nnite          ! number of night columns
         944             :    integer,  intent(in) :: idxnite(pcols) ! indices of night columns in the chunk
         945             :
         946             :    type(ty_optical_props_2str), intent(inout) :: aer_sw
         947             :
         948             :    ! local variables
         949             :    integer  :: i
         950             :
         951             :    ! The optical arrays dimensioned in the vertical as 0:pver.
         952             :    ! The index 0 is for the extra layer used in the radiation
         953             :    ! calculation.  The index ktopcam assumes the CAM vertical indices are
         954             :    ! in the range 1:pver, so using this index correctly ignores vertical
         955             :    ! index 0.  If an "extra" layer is used in the calculations, it is
         956             :    ! provided and set in the RRTMGP aerosol object aer_sw.
         957             :    real(r8) :: aer_tau    (pcols,0:pver,nswbands) ! extinction optical depth
         958             :    real(r8) :: aer_tau_w  (pcols,0:pver,nswbands) ! single scattering albedo * tau
         959             :    real(r8) :: aer_tau_w_g(pcols,0:pver,nswbands) ! asymmetry parameter * w * tau
         960             :    real(r8) :: aer_tau_w_f(pcols,0:pver,nswbands) ! forward scattered fraction * w * tau
         961             :                                                   ! aer_tau_w_f is not used by RRTMGP.
         962             :    character(len=*), parameter :: sub = 'rrtmgp_set_aer_sw'
         963             :    !--------------------------------------------------------------------------------
         964             :
         965             :    ! Get aerosol shortwave optical properties.
         966             :    ! Make outfld calls for aerosol optical property diagnostics.
         967             :    call aer_rad_props_sw( &
         968             :       icall, state, pbuf, nnite, idxnite, &
         969      746136 :       aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f)
         970             :
         971             :    ! The aer_sw object is only initialized if nday > 0.
         972      746136 :    if (nday > 0) then
         973             :
         974             :       ! aerosol optical properties need to be re-ordered from the RRTMG spectral bands
         975             :       ! (as assumed in the optics datasets) to the RRTMGP band order.
         976 17428861562 :       aer_tau(:,:,:)     = aer_tau(    :,:,rrtmg_to_rrtmgp_swbands)
         977 17428861562 :       aer_tau_w(:,:,:)   = aer_tau_w(  :,:,rrtmg_to_rrtmgp_swbands)
         978 17428861562 :       aer_tau_w_g(:,:,:) = aer_tau_w_g(:,:,rrtmg_to_rrtmgp_swbands)
         979             :
         980             :       ! If there is an extra layer in the radiation then this initialization
         981             :       ! will provide default values.
         982  8224999853 :       aer_sw%tau = 0.0_r8
         983  8224999853 :       aer_sw%ssa = 1.0_r8
         984  8224999853 :       aer_sw%g   = 0.0_r8
         985             :
         986             :       ! CAM fields are products tau, tau*ssa, tau*ssa*asy
         987             :       ! Fields expected by RRTMGP are computed by
         988             :       ! aer_sw%tau = aer_tau
         989             :       ! aer_sw%ssa = aer_tau_w / aer_tau
         990             :       ! aer_sw%g   = aer_tau_w_g / aer_taw_w
         991             :       ! aer_sw arrays have dimensions of (nday,nlay,nswbands)
         992             :
         993     6245563 :       do i = 1, nday
         994             :          ! set aerosol optical depth, clip to zero
         995  7712747100 :          aer_sw%tau(i,ktoprad:,:) = max(aer_tau(idxday(i),ktopcam:,:), 0._r8)
         996             :          ! set value of single scattering albedo
         997     5856300 :          aer_sw%ssa(i,ktoprad:,:) = merge(aer_tau_w(idxday(i),ktopcam:,:)/aer_tau(idxday(i),ktopcam:,:), &
         998  7718603400 :                                           1._r8, aer_tau(idxday(i),ktopcam:,:) > 0._r8)
         999             :          ! set value of asymmetry
        1000     5856300 :          aer_sw%g(i,ktoprad:,:) = merge(aer_tau_w_g(idxday(i),ktopcam:,:)/aer_tau_w(idxday(i),ktopcam:,:), &
        1001  7718992663 :                                         0._r8, aer_tau_w(idxday(i),ktopcam:,:) > tiny)
        1002             :       end do
        1003             :
        1004             :       ! impose limits on the components
        1005  8224999853 :       aer_sw%ssa = min(max(aer_sw%ssa, 0._r8), 1._r8)
        1006  8224999853 :       aer_sw%g   = min(max(aer_sw%g, -1._r8), 1._r8)
        1007             :
        1008             :    end if
        1009             :
        1010      746136 : end subroutine rrtmgp_set_aer_sw
        1011             :
        1012             : !==================================================================================================
        1013             :
        1014             : end module rrtmgp_inputs

Generated by: LCOV version 1.14