LCOV - code coverage report
Current view: | top level - physics/cam - cloud_rad_props.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 218 | 260 | 83.8 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 12 | 15 | 80.0 %

              Line data    Source code
           1             : module cloud_rad_props
           2             :
           3             : !------------------------------------------------------------------------------------------------
           4             : !------------------------------------------------------------------------------------------------
           5             :
           6             : use shr_kind_mod,     only: r8 => shr_kind_r8
           7             : use ppgrid,           only: pcols, pver, pverp
           8             : use physics_types,    only: physics_state
           9             : use physics_buffer,   only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx
          10             : use constituents,     only: cnst_get_ind
          11             : use radconstants,     only: nswbands, nlwbands, idx_sw_diag
          12             : use rad_constituents, only: iceopticsfile, liqopticsfile
          13             : use oldcloud_optics,  only: oldcloud_init, oldcloud_lw, &
          14             :                             old_liq_get_rad_props_lw, old_ice_get_rad_props_lw
          15             :
          16             : use slingo_liq_optics,      only: slingo_rad_props_init
          17             : use ebert_curry_ice_optics, only: ec_rad_props_init, scalefactor
          18             :
          19             : use interpolate_data, only: interp_type, lininterp_init, lininterp, &
          20             :                             extrap_method_bndry, lininterp_finish
          21             :
          22             : use cam_logfile,      only: iulog
          23             : use cam_abortutils,   only: endrun
          24             :
          25             :
          26             : implicit none
          27             : private
          28             : save
          29             :
          30             : public :: &
          31             :    cloud_rad_props_init,          &
          32             :    cloud_rad_props_get_lw,        & ! return LW optical props for old cloud optics
          33             :    get_ice_optics_sw,             & ! return Mitchell SW ice radiative properties
          34             :    ice_cloud_get_rad_props_lw,    & ! return Mitchell LW ice radiative properties
          35             :    get_liquid_optics_sw,          & ! return Conley SW radiative properties
          36             :    liquid_cloud_get_rad_props_lw, & ! return Conley LW radiative properties
          37             :    get_snow_optics_sw,            &
          38             :    snow_cloud_get_rad_props_lw,   &
          39             :    get_grau_optics_sw,            &
          40             :    grau_cloud_get_rad_props_lw
          41             :
          42             :
          43             : integer :: nmu, nlambda
          44             : real(r8), allocatable :: g_mu(:)           ! mu samples on grid
          45             : real(r8), allocatable :: g_lambda(:,:)     ! lambda scale samples on grid
          46             : real(r8), allocatable :: ext_sw_liq(:,:,:)
          47             : real(r8), allocatable :: ssa_sw_liq(:,:,:)
          48             : real(r8), allocatable :: asm_sw_liq(:,:,:)
          49             : real(r8), allocatable :: abs_lw_liq(:,:,:)
          50             :
          51             : integer :: n_g_d
          52             : real(r8), allocatable :: g_d_eff(:)        ! radiative effective diameter samples on grid
          53             : real(r8), allocatable :: ext_sw_ice(:,:)
          54             : real(r8), allocatable :: ssa_sw_ice(:,:)
          55             : real(r8), allocatable :: asm_sw_ice(:,:)
          56             : real(r8), allocatable :: abs_lw_ice(:,:)
          57             :
          58             : ! indexes into pbuf for optical parameters of MG clouds
          59             : integer :: i_dei=0
          60             : integer :: i_mu=0
          61             : integer :: i_lambda=0
          62             : integer :: i_iciwp=0
          63             : integer :: i_iclwp=0
          64             : integer :: i_des=0
          65             : integer :: i_icswp=0
          66             : integer :: i_degrau=0
          67             : integer :: i_icgrauwp=0
          68             :
          69             : ! indexes into constituents for old optics
          70             : integer :: &
          71             :    ixcldice,           & ! cloud ice water index
          72             :    ixcldliq              ! cloud liquid water index
          73             :
          74             : real(r8), parameter :: tiny = 1.e-80_r8
          75             :
          76             : !==============================================================================
          77             : contains
          78             : !==============================================================================
          79             :
          80        1536 : subroutine cloud_rad_props_init()
          81             :
          82             :    use netcdf
          83             :    use spmd_utils,     only: masterproc
          84             :    use ioFileMod,      only: getfil
          85             :    use error_messages, only: handle_ncerr
          86             : #if ( defined SPMD )
          87             :    use mpishorthand
          88             : #endif
          89             :
          90             :    character(len=256) :: liquidfile
          91             :    character(len=256) :: icefile
          92             :    character(len=256) :: locfn
          93             :
          94             :    integer :: ncid, dimid, f_nlwbands, f_nswbands, ierr
          95             :    integer :: vdimids(NF90_MAX_VAR_DIMS), ndims, templen
          96             :    ! liquid clouds
          97             :    integer :: mudimid, lambdadimid
          98             :    integer :: mu_id, lambda_id, ext_sw_liq_id, ssa_sw_liq_id, asm_sw_liq_id, abs_lw_liq_id
          99             :
         100             :    ! ice clouds
         101             :    integer :: d_dimid ! diameters
         102             :    integer :: d_id, ext_sw_ice_id, ssa_sw_ice_id, asm_sw_ice_id, abs_lw_ice_id
         103             :
         104             :    integer :: err
         105             :    character(len=*), parameter :: sub = 'cloud_rad_props_init'
         106             :
         107        1536 :    liquidfile = liqopticsfile
         108        1536 :    icefile = iceopticsfile
         109             :
         110        1536 :    call slingo_rad_props_init
         111        1536 :    call ec_rad_props_init
         112        1536 :    call oldcloud_init
         113             :
         114        1536 :    i_dei    = pbuf_get_index('DEI',errcode=err)
         115        1536 :    i_mu     = pbuf_get_index('MU',errcode=err)
         116        1536 :    i_lambda = pbuf_get_index('LAMBDAC',errcode=err)
         117        1536 :    i_iciwp  = pbuf_get_index('ICIWP',errcode=err)
         118        1536 :    i_iclwp  = pbuf_get_index('ICLWP',errcode=err)
         119        1536 :    i_des    = pbuf_get_index('DES',errcode=err)
         120        1536 :    i_icswp  = pbuf_get_index('ICSWP',errcode=err)
         121        1536 :    i_icgrauwp  = pbuf_get_index('ICGRAUWP',errcode=err) ! Available when using MG3
         122        1536 :    i_degrau    = pbuf_get_index('DEGRAU',errcode=err)   ! Available when using MG3
         123             :
         124             :    ! old optics
         125        1536 :    call cnst_get_ind('CLDICE', ixcldice)
         126        1536 :    call cnst_get_ind('CLDLIQ', ixcldliq)
         127             :
         128             :    ! read liquid cloud optics
         129        1536 :    if (masterproc) then
         130           2 :       call getfil( trim(liquidfile), locfn, 0)
         131           2 :       call handle_ncerr( nf90_open(locfn, NF90_NOWRITE, ncid), 'liquid optics file missing')
         132           2 :       write(iulog,*)' reading liquid cloud optics from file ',locfn
         133             :
         134           2 :       call handle_ncerr(nf90_inq_dimid( ncid, 'lw_band', dimid), 'getting lw_band dim')
         135           2 :       call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nlwbands), 'getting n lw bands')
         136           2 :       if (f_nlwbands /= nlwbands) call endrun(sub//': number of lw bands does not match')
         137             :
         138           2 :       call handle_ncerr(nf90_inq_dimid( ncid, 'sw_band', dimid), 'getting sw_band_dim')
         139           2 :       call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nswbands), 'getting n sw bands')
         140           2 :       if (f_nswbands /= nswbands) call endrun(sub//': number of sw bands does not match')
         141             :
         142           2 :       call handle_ncerr(nf90_inq_dimid( ncid, 'mu', mudimid), 'getting mu dim')
         143           2 :       call handle_ncerr(nf90_inquire_dimension( ncid, mudimid, len=nmu), 'getting n mu samples')
         144             :
         145           2 :       call handle_ncerr(nf90_inq_dimid( ncid, 'lambda_scale', lambdadimid), 'getting lambda dim')
         146           2 :       call handle_ncerr(nf90_inquire_dimension( ncid, lambdadimid, len=nlambda), 'getting n lambda samples')
         147             :    end if ! if (masterproc)
         148             :
         149             : #if ( defined SPMD )
         150        1536 :    call mpibcast(nmu, 1, mpiint, 0, mpicom, ierr)
         151        1536 :    call mpibcast(nlambda, 1, mpiint, 0, mpicom, ierr)
         152             : #endif
         153             :
         154        4608 :    allocate(g_mu(nmu))
         155        6144 :    allocate(g_lambda(nmu,nlambda))
         156        7680 :    allocate(ext_sw_liq(nmu,nlambda,nswbands) )
         157        4608 :    allocate(ssa_sw_liq(nmu,nlambda,nswbands))
         158        4608 :    allocate(asm_sw_liq(nmu,nlambda,nswbands))
         159        7680 :    allocate(abs_lw_liq(nmu,nlambda,nlwbands))
         160             :
         161        1536 :    if (masterproc) then
         162             :       call handle_ncerr( nf90_inq_varid(ncid, 'mu', mu_id),&
         163           2 :          'cloud optics mu get')
         164             :       call handle_ncerr( nf90_get_var(ncid, mu_id, g_mu),&
         165           2 :          'read cloud optics mu values')
         166             :
         167             :       call handle_ncerr( nf90_inq_varid(ncid, 'lambda', lambda_id),&
         168           2 :          'cloud optics lambda get')
         169             :       call handle_ncerr( nf90_get_var(ncid, lambda_id, g_lambda),&
         170           2 :          'read cloud optics lambda values')
         171             :
         172             :       call handle_ncerr( nf90_inq_varid(ncid, 'k_ext_sw', ext_sw_liq_id),&
         173           2 :          'cloud optics ext_sw_liq get')
         174             :       call handle_ncerr( nf90_get_var(ncid, ext_sw_liq_id, ext_sw_liq),&
         175           2 :          'read cloud optics ext_sw_liq values')
         176             :
         177             :       call handle_ncerr( nf90_inq_varid(ncid, 'ssa_sw', ssa_sw_liq_id),&
         178           2 :          'cloud optics ssa_sw_liq get')
         179             :       call handle_ncerr( nf90_get_var(ncid, ssa_sw_liq_id, ssa_sw_liq),&
         180           2 :          'read cloud optics ssa_sw_liq values')
         181             :
         182             :       call handle_ncerr( nf90_inq_varid(ncid, 'asm_sw', asm_sw_liq_id),&
         183           2 :          'cloud optics asm_sw_liq get')
         184             :       call handle_ncerr( nf90_get_var(ncid, asm_sw_liq_id, asm_sw_liq),&
         185           2 :          'read cloud optics asm_sw_liq values')
         186             :
         187             :       call handle_ncerr( nf90_inq_varid(ncid, 'k_abs_lw', abs_lw_liq_id),&
         188           2 :          'cloud optics abs_lw_liq get')
         189             :       call handle_ncerr( nf90_get_var(ncid, abs_lw_liq_id, abs_lw_liq),&
         190           2 :          'read cloud optics abs_lw_liq values')
         191             :
         192           2 :       call handle_ncerr( nf90_close(ncid), 'liquid optics file missing')
         193             :    end if ! if masterproc
         194             :
         195             : #if ( defined SPMD )
         196        1536 :     call mpibcast(g_mu, nmu, mpir8, 0, mpicom, ierr)
         197        1536 :     call mpibcast(g_lambda, nmu*nlambda, mpir8, 0, mpicom, ierr)
         198        1536 :     call mpibcast(ext_sw_liq, nmu*nlambda*nswbands, mpir8, 0, mpicom, ierr)
         199        1536 :     call mpibcast(ssa_sw_liq, nmu*nlambda*nswbands, mpir8, 0, mpicom, ierr)
         200        1536 :     call mpibcast(asm_sw_liq, nmu*nlambda*nswbands, mpir8, 0, mpicom, ierr)
         201        1536 :     call mpibcast(abs_lw_liq, nmu*nlambda*nlwbands, mpir8, 0, mpicom, ierr)
         202             : #endif
         203             :    ! Convert kext from m^2/Volume to m^2/Kg
         204    22602240 :    ext_sw_liq = ext_sw_liq / 0.9970449e3_r8
         205    25830912 :    abs_lw_liq = abs_lw_liq / 0.9970449e3_r8
         206             :
         207             :    ! read ice cloud optics
         208        1536 :    if (masterproc) then
         209           2 :       call getfil( trim(icefile), locfn, 0)
         210           2 :       call handle_ncerr( nf90_open(locfn, NF90_NOWRITE, ncid), 'ice optics file missing')
         211           2 :       write(iulog,*)' reading ice cloud optics from file ',locfn
         212           2 :       call handle_ncerr(nf90_inq_dimid( ncid, 'lw_band', dimid), 'getting lw_band dim')
         213           2 :       call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nlwbands), 'getting n lw bands')
         214           2 :       if (f_nlwbands /= nlwbands) then
         215           0 :          call endrun(sub//': number of lw bands does not match')
         216             :       end if
         217           2 :       call handle_ncerr(nf90_inq_dimid( ncid, 'sw_band', dimid), 'getting sw_band_dim')
         218           2 :       call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nswbands), 'getting n sw bands')
         219           2 :       if (f_nswbands /= nswbands) then
         220           0 :          call endrun(sub//': number of sw bands does not match')
         221             :       end if
         222           2 :       call handle_ncerr(nf90_inq_dimid( ncid, 'd_eff', d_dimid), 'getting deff dim')
         223           2 :       call handle_ncerr(nf90_inquire_dimension( ncid, d_dimid, len=n_g_d), 'getting n deff samples')
         224             :    end if ! if (masterproc)
         225             :
         226             : #if ( defined SPMD )
         227        1536 :    call mpibcast(n_g_d, 1, mpiint, 0, mpicom, ierr)
         228             : !   call mpibcast(nswbands, 1, mpiint, 0, mpicom, ierr)
         229             : !   call mpibcast(nlwbands, 1, mpiint, 0, mpicom, ierr)
         230             : #endif
         231             :
         232        4608 :    allocate(g_d_eff(n_g_d))
         233        4608 :    allocate(ext_sw_ice(n_g_d,nswbands))
         234        3072 :    allocate(ssa_sw_ice(n_g_d,nswbands))
         235        3072 :    allocate(asm_sw_ice(n_g_d,nswbands))
         236        4608 :    allocate(abs_lw_ice(n_g_d,nlwbands))
         237             :
         238        1536 :    if (masterproc) then
         239             :       call handle_ncerr( nf90_inq_varid(ncid, 'd_eff', d_id),&
         240           2 :          'cloud optics deff get')
         241             :       call handle_ncerr( nf90_get_var(ncid, d_id, g_d_eff),&
         242           2 :          'read cloud optics deff values')
         243             :
         244             :       call handle_ncerr( nf90_inq_varid(ncid, 'sw_ext', ext_sw_ice_id),&
         245           2 :          'cloud optics ext_sw_ice get')
         246             :       call handle_ncerr(nf90_inquire_variable ( ncid, ext_sw_ice_id, ndims=ndims, dimids=vdimids),&
         247           2 :          'checking dimensions of ext_sw_ice')
         248             :       call handle_ncerr(nf90_inquire_dimension( ncid, vdimids(1), len=templen),&
         249           2 :          'getting first dimension sw_ext')
         250             :       call handle_ncerr(nf90_inquire_dimension( ncid, vdimids(2), len=templen),&
         251           2 :          'getting first dimension sw_ext')
         252             :       call handle_ncerr( nf90_get_var(ncid, ext_sw_ice_id, ext_sw_ice),&
         253           2 :          'read cloud optics ext_sw_ice values')
         254             :
         255             :       call handle_ncerr( nf90_inq_varid(ncid, 'sw_ssa', ssa_sw_ice_id),&
         256           2 :          'cloud optics ssa_sw_ice get')
         257             :       call handle_ncerr( nf90_get_var(ncid, ssa_sw_ice_id, ssa_sw_ice),&
         258           2 :          'read cloud optics ssa_sw_ice values')
         259             :
         260             :       call handle_ncerr( nf90_inq_varid(ncid, 'sw_asm', asm_sw_ice_id),&
         261           2 :          'cloud optics asm_sw_ice get')
         262             :       call handle_ncerr( nf90_get_var(ncid, asm_sw_ice_id, asm_sw_ice),&
         263           2 :          'read cloud optics asm_sw_ice values')
         264             :
         265             :       call handle_ncerr( nf90_inq_varid(ncid, 'lw_abs', abs_lw_ice_id),&
         266           2 :          'cloud optics abs_lw_ice get')
         267             :       call handle_ncerr( nf90_get_var(ncid, abs_lw_ice_id, abs_lw_ice),&
         268           2 :          'read cloud optics abs_lw_ice values')
         269             :
         270           2 :       call handle_ncerr( nf90_close(ncid), 'ice optics file missing')
         271             :    end if ! if masterproc
         272             :
         273             : #if ( defined SPMD )
         274        1536 :    call mpibcast(g_d_eff, n_g_d, mpir8, 0, mpicom, ierr)
         275        1536 :    call mpibcast(ext_sw_ice, n_g_d*nswbands, mpir8, 0, mpicom, ierr)
         276        1536 :    call mpibcast(ssa_sw_ice, n_g_d*nswbands, mpir8, 0, mpicom, ierr)
         277        1536 :    call mpibcast(asm_sw_ice, n_g_d*nswbands, mpir8, 0, mpicom, ierr)
         278        1536 :    call mpibcast(abs_lw_ice, n_g_d*nlwbands, mpir8, 0, mpicom, ierr)
         279             : #endif
         280             :
         281        1536 :    return
         282             :
         283             : end subroutine cloud_rad_props_init
         284             :
         285             : !==============================================================================
         286             :
         287           0 : subroutine cloud_rad_props_get_lw(state, pbuf, cld_abs_od, oldliq, oldice, oldcloud)
         288             :
         289             :    ! Purpose: Compute cloud longwave absorption optical depth
         290             :
         291             :    ! Arguments
         292             :    type(physics_state), intent(in)  :: state
         293             :    type(physics_buffer_desc),pointer:: pbuf(:)
         294             :    real(r8),            intent(out) :: cld_abs_od(nlwbands,pcols,pver) ! [fraction] absorption optical depth, per layer
         295             :    logical, optional,   intent(in)  :: oldliq  ! use old liquid optics
         296             :    logical, optional,   intent(in)  :: oldice  ! use old ice optics
         297             :    logical, optional,   intent(in)  :: oldcloud  ! use old optics for both (b4b)
         298             :
         299             :    ! Local variables
         300             :    integer :: ncol        ! number of columns
         301             :
         302             :    ! rad properties for liquid clouds
         303             :    real(r8) :: liq_tau_abs_od(nlwbands,pcols,pver) ! liquid cloud absorption optical depth
         304             :
         305             :    ! rad properties for ice clouds
         306             :    real(r8) :: ice_tau_abs_od(nlwbands,pcols,pver) ! ice cloud absorption optical depth
         307             :    !-----------------------------------------------------------------------------
         308             :
         309           0 :    ncol = state%ncol
         310             :
         311           0 :    cld_abs_od = 0._r8
         312             :
         313           0 :    if(present(oldcloud))then
         314           0 :       if(oldcloud) then
         315           0 :          call oldcloud_lw(state,pbuf,cld_abs_od,oldwp=.false.)
         316           0 :          return
         317             :       endif
         318             :    endif
         319             :
         320           0 :    if(present(oldliq))then
         321           0 :       if(oldliq) then
         322           0 :          call old_liq_get_rad_props_lw(state, pbuf, liq_tau_abs_od, oldliqwp=.false.)
         323             :       else
         324           0 :          call liquid_cloud_get_rad_props_lw(state, pbuf, liq_tau_abs_od)
         325             :       endif
         326             :    else
         327           0 :       call liquid_cloud_get_rad_props_lw(state, pbuf, liq_tau_abs_od)
         328             :    endif
         329             :
         330           0 :    if(present(oldice))then
         331           0 :       if(oldice) then
         332           0 :          call old_ice_get_rad_props_lw(state, pbuf, ice_tau_abs_od, oldicewp=.false.)
         333             :       else
         334           0 :          call ice_cloud_get_rad_props_lw(state, pbuf, ice_tau_abs_od)
         335             :       endif
         336             :    else
         337           0 :       call ice_cloud_get_rad_props_lw(state, pbuf, ice_tau_abs_od)
         338             :    endif
         339             :
         340           0 :    cld_abs_od(:,1:ncol,:) = liq_tau_abs_od(:,1:ncol,:) + ice_tau_abs_od(:,1:ncol,:)
         341             :
         342             : end subroutine cloud_rad_props_get_lw
         343             :
         344             : !==============================================================================
         345             :
         346      746136 : subroutine get_ice_optics_sw(state, pbuf, tau, tau_w, tau_w_g, tau_w_f)
         347             :    type(physics_state), intent(in)   :: state
         348             :    type(physics_buffer_desc),pointer :: pbuf(:)
         349             :
         350             :    real(r8),intent(out) :: tau    (nswbands,pcols,pver) ! extinction optical depth
         351             :    real(r8),intent(out) :: tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
         352             :    real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! asymmetry parameter * tau * w
         353             :    real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w
         354             :
         355      746136 :    real(r8), pointer :: iciwpth(:,:), dei(:,:)
         356             :
         357             :    ! Get relevant pbuf fields, and interpolate optical properties from
         358             :    ! the lookup tables.
         359      746136 :    call pbuf_get_field(pbuf, i_iciwp, iciwpth)
         360      746136 :    call pbuf_get_field(pbuf, i_dei,   dei)
         361             :
         362             :    call interpolate_ice_optics_sw(state%ncol, iciwpth, dei, tau, tau_w, &
         363      746136 :         tau_w_g, tau_w_f)
         364             :
         365      746136 : end subroutine get_ice_optics_sw
         366             :
         367             : !==============================================================================
         368             :
         369      746136 : subroutine get_snow_optics_sw(state, pbuf, tau, tau_w, tau_w_g, tau_w_f)
         370             :    type(physics_state), intent(in)   :: state
         371             :    type(physics_buffer_desc),pointer :: pbuf(:)
         372             :
         373             :    real(r8),intent(out) :: tau    (nswbands,pcols,pver) ! extinction optical depth
         374             :    real(r8),intent(out) :: tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
         375             :    real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! asymmetry parameter * tau * w
         376             :    real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w
         377             :
         378      746136 :    real(r8), pointer :: icswpth(:,:), des(:,:)
         379             :
         380             :    ! This does the same thing as get_ice_optics_sw, except with a different
         381             :    ! water path and effective diameter.
         382      746136 :    call pbuf_get_field(pbuf, i_icswp, icswpth)
         383      746136 :    call pbuf_get_field(pbuf, i_des,   des)
         384             :
         385             :    call interpolate_ice_optics_sw(state%ncol, icswpth, des, tau, tau_w, &
         386      746136 :         tau_w_g, tau_w_f)
         387             :
         388      746136 : end subroutine get_snow_optics_sw
         389             :
         390             : !==============================================================================
         391             :
         392           0 : subroutine get_grau_optics_sw(state, pbuf, tau, tau_w, tau_w_g, tau_w_f)
         393             :    type(physics_state), intent(in)   :: state
         394             :    type(physics_buffer_desc),pointer :: pbuf(:)
         395             :
         396             :    real(r8),intent(out) :: tau    (nswbands,pcols,pver) ! extinction optical depth
         397             :    real(r8),intent(out) :: tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
         398             :    real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! asymmetry parameter * tau * w
         399             :    real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w
         400             :
         401           0 :    real(r8), pointer :: icgrauwpth(:,:), degrau(:,:)
         402             :
         403             :    integer :: i,k
         404             :    character(len=*), parameter :: sub = 'get_grau_optics_sw'
         405             :
         406             :    ! This does the same thing as get_ice_optics_sw, except with a different
         407             :    ! water path and effective diameter.
         408           0 :    if((i_icgrauwp > 0) .and. (i_degrau > 0)) then
         409             :
         410           0 :       call pbuf_get_field(pbuf, i_icgrauwp, icgrauwpth)
         411           0 :       call pbuf_get_field(pbuf, i_degrau,   degrau)
         412             :
         413             :       call interpolate_ice_optics_sw(state%ncol, icgrauwpth, degrau, tau, tau_w, &
         414           0 :            tau_w_g, tau_w_f)
         415           0 :       do i = 1, pcols
         416           0 :          do k = 1, pver
         417           0 :             if (tau(idx_sw_diag,i,k).gt.100._r8) then
         418           0 :                write(iulog,*) 'WARNING: SW Graupel Tau > 100  (i,k,icgrauwpth,degrau,tau):'
         419           0 :                write(iulog,*) i,k,icgrauwpth(i,k), degrau(i,k), tau(idx_sw_diag,i,k)
         420             :             end if
         421             :          enddo
         422             :       enddo
         423             :
         424             :    else
         425           0 :       call endrun(sub//': ERROR: Get_grau_optics_sw called when graupel properties not supported')
         426             :    end if
         427             :
         428           0 : end subroutine get_grau_optics_sw
         429             :
         430             : !==============================================================================
         431             :
         432      746136 : subroutine get_liquid_optics_sw(state, pbuf, tau, tau_w, tau_w_g, tau_w_f)
         433             :    type(physics_state), intent(in)   :: state
         434             :    type(physics_buffer_desc),pointer :: pbuf(:)
         435             :
         436             :    real(r8),intent(out) :: tau    (nswbands,pcols,pver) ! extinction optical depth
         437             :    real(r8),intent(out) :: tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
         438             :    real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! asymmetry parameter * tau * w
         439             :    real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w
         440             :
         441      746136 :    real(r8), pointer, dimension(:,:) :: lamc, pgam, iclwpth
         442             :    real(r8), dimension(pcols,pver) :: kext
         443             :    integer i,k,swband, ncol
         444             :
         445      746136 :    ncol = state%ncol
         446             :
         447             :
         448      746136 :    call pbuf_get_field(pbuf, i_lambda,  lamc)
         449      746136 :    call pbuf_get_field(pbuf, i_mu,      pgam)
         450      746136 :    call pbuf_get_field(pbuf, i_iclwp,   iclwpth)
         451             :
         452    70136784 :    do k = 1,pver
         453  1159408584 :       do i = 1,ncol
         454  1158662448 :          if(lamc(i,k) > 0._r8) then ! This seems to be clue from microphysics of no cloud
         455           0 :             call gam_liquid_sw(iclwpth(i,k), lamc(i,k), pgam(i,k), &
         456   858219788 :                 tau(1:nswbands,i,k), tau_w(1:nswbands,i,k), tau_w_g(1:nswbands,i,k), tau_w_f(1:nswbands,i,k))
         457             :          else
         458  3465780180 :             tau(1:nswbands,i,k) = 0._r8
         459  3465780180 :             tau_w(1:nswbands,i,k) = 0._r8
         460  3465780180 :             tau_w_g(1:nswbands,i,k) = 0._r8
         461  3465780180 :             tau_w_f(1:nswbands,i,k) = 0._r8
         462             :          endif
         463             :       enddo
         464             :    enddo
         465             :
         466      746136 : end subroutine get_liquid_optics_sw
         467             :
         468             : !==============================================================================
         469             :
         470      746136 : subroutine liquid_cloud_get_rad_props_lw(state, pbuf, abs_od)
         471             :    type(physics_state), intent(in)    :: state
         472             :    type(physics_buffer_desc),pointer  :: pbuf(:)
         473             :    real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)
         474             :
         475             :    integer :: ncol
         476      746136 :    real(r8), pointer, dimension(:,:) :: lamc, pgam, iclwpth
         477             :
         478             :    integer lwband, i, k
         479             :
         480      746136 :    abs_od = 0._r8
         481             :
         482      746136 :    ncol = state%ncol
         483             :
         484      746136 :    call pbuf_get_field(pbuf, i_lambda,  lamc)
         485      746136 :    call pbuf_get_field(pbuf, i_mu,      pgam)
         486      746136 :    call pbuf_get_field(pbuf, i_iclwp,   iclwpth)
         487             :
         488    70136784 :    do k = 1,pver
         489  1159408584 :       do i = 1,ncol
         490  1158662448 :          if(lamc(i,k) > 0._r8) then ! This seems to be the clue for no cloud from microphysics formulation
         491   858219788 :             call gam_liquid_lw(iclwpth(i,k), lamc(i,k), pgam(i,k), abs_od(1:nlwbands,i,k))
         492             :          else
         493  3927884204 :             abs_od(1:nlwbands,i,k) = 0._r8
         494             :          endif
         495             :       enddo
         496             :    enddo
         497             :
         498      746136 : end subroutine liquid_cloud_get_rad_props_lw
         499             : !==============================================================================
         500             :
         501      746136 : subroutine snow_cloud_get_rad_props_lw(state, pbuf, abs_od)
         502             :    type(physics_state), intent(in)    :: state
         503             :    type(physics_buffer_desc), pointer :: pbuf(:)
         504             :    real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)
         505             :
         506      746136 :    real(r8), pointer :: icswpth(:,:), des(:,:)
         507             :
         508             :    ! This does the same thing as ice_cloud_get_rad_props_lw, except with a
         509             :    ! different water path and effective diameter.
         510      746136 :    call pbuf_get_field(pbuf, i_icswp, icswpth)
         511      746136 :    call pbuf_get_field(pbuf, i_des,   des)
         512             :
         513      746136 :    call interpolate_ice_optics_lw(state%ncol,icswpth, des, abs_od)
         514             :
         515      746136 : end subroutine snow_cloud_get_rad_props_lw
         516             :
         517             :
         518             : !==============================================================================
         519             :
         520           0 : subroutine grau_cloud_get_rad_props_lw(state, pbuf, abs_od)
         521             :    type(physics_state), intent(in)    :: state
         522             :    type(physics_buffer_desc), pointer :: pbuf(:)
         523             :    real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)
         524             :
         525           0 :    real(r8), pointer :: icgrauwpth(:,:), degrau(:,:)
         526             :    character(len=*), parameter :: sub = 'grau_cloud_get_rad_props_lw'
         527             :
         528             :    ! This does the same thing as ice_cloud_get_rad_props_lw, except with a
         529             :    ! different water path and effective diameter.
         530           0 :    if((i_icgrauwp > 0) .and. (i_degrau > 0)) then
         531           0 :       call pbuf_get_field(pbuf, i_icgrauwp, icgrauwpth)
         532           0 :       call pbuf_get_field(pbuf, i_degrau,   degrau)
         533             :
         534           0 :       call interpolate_ice_optics_lw(state%ncol,icgrauwpth, degrau, abs_od)
         535             :    else
         536             :       call endrun(sub//': ERROR: Grau_cloud_get_rad_props_lw called when graupel &
         537           0 :            &properties not supported')
         538             :    end if
         539             :
         540           0 : end subroutine grau_cloud_get_rad_props_lw
         541             :
         542             : !==============================================================================
         543             :
         544      746136 : subroutine ice_cloud_get_rad_props_lw(state, pbuf, abs_od)
         545             :    type(physics_state), intent(in)     :: state
         546             :    type(physics_buffer_desc), pointer  :: pbuf(:)
         547             :    real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)
         548             :
         549      746136 :    real(r8), pointer :: iciwpth(:,:), dei(:,:)
         550             :
         551             :    ! Get relevant pbuf fields, and interpolate optical properties from
         552             :    ! the lookup tables.
         553      746136 :    call pbuf_get_field(pbuf, i_iciwp, iciwpth)
         554      746136 :    call pbuf_get_field(pbuf, i_dei,   dei)
         555             :
         556      746136 :    call interpolate_ice_optics_lw(state%ncol,iciwpth, dei, abs_od)
         557             :
         558      746136 : end subroutine ice_cloud_get_rad_props_lw
         559             :
         560             : !==============================================================================
         561             : ! Private methods
         562             : !==============================================================================
         563             :
         564     1492272 : subroutine interpolate_ice_optics_sw(ncol, iciwpth, dei, tau, tau_w, &
         565             :      tau_w_g, tau_w_f)
         566             :
         567             :   integer, intent(in) :: ncol
         568             :   real(r8), intent(in) :: iciwpth(pcols,pver)
         569             :   real(r8), intent(in) :: dei(pcols,pver)
         570             :
         571             :   real(r8),intent(out) :: tau    (nswbands,pcols,pver) ! extinction optical depth
         572             :   real(r8),intent(out) :: tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
         573             :   real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! asymmetry parameter * tau * w
         574             :   real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w
         575             :
         576             :   type(interp_type) :: dei_wgts
         577             :
         578             :   integer :: i, k, swband
         579             :   real(r8) :: ext(nswbands), ssa(nswbands), asm(nswbands)
         580             :
         581   140273568 :   do k = 1,pver
         582  2318817168 :      do i = 1,ncol
         583  2317324896 :         if( iciwpth(i,k) < tiny .or. dei(i,k) == 0._r8) then
         584             :            ! if ice water path is too small, OD := 0
         585 29847898740 :            tau    (:,i,k) = 0._r8
         586 29847898740 :            tau_w  (:,i,k) = 0._r8
         587 29847898740 :            tau_w_g(:,i,k) = 0._r8
         588 29847898740 :            tau_w_f(:,i,k) = 0._r8
         589             :         else
         590             :            ! for each cell interpolate to find weights in g_d_eff grid.
         591             :            call lininterp_init(g_d_eff, n_g_d, dei(i:i,k), 1, &
         592   188683684 :                 extrap_method_bndry, dei_wgts)
         593             :            ! interpolate into grid and extract radiative properties
         594  2830255260 :            do swband = 1, nswbands
         595             :               call lininterp(ext_sw_ice(:,swband), n_g_d, &
         596  2641571576 :                    ext(swband:swband), 1, dei_wgts)
         597             :               call lininterp(ssa_sw_ice(:,swband), n_g_d, &
         598  2641571576 :                    ssa(swband:swband), 1, dei_wgts)
         599             :               call lininterp(asm_sw_ice(:,swband), n_g_d, &
         600  2830255260 :                    asm(swband:swband), 1, dei_wgts)
         601             :            end do
         602  2830255260 :            tau    (:,i,k) = iciwpth(i,k) * ext
         603  2830255260 :            tau_w  (:,i,k) = tau(:,i,k) * ssa
         604  2830255260 :            tau_w_g(:,i,k) = tau_w(:,i,k) * asm
         605  2830255260 :            tau_w_f(:,i,k) = tau_w_g(:,i,k) * asm
         606   188683684 :            call lininterp_finish(dei_wgts)
         607             :         endif
         608             :      enddo
         609             :   enddo
         610             :
         611     1492272 : end subroutine interpolate_ice_optics_sw
         612             :
         613             : !==============================================================================
         614             :
         615     1492272 : subroutine interpolate_ice_optics_lw(ncol, iciwpth, dei, abs_od)
         616             :
         617             :   integer, intent(in) :: ncol
         618             :   real(r8), intent(in) :: iciwpth(pcols,pver)
         619             :   real(r8), intent(in) :: dei(pcols,pver)
         620             :
         621             :   real(r8),intent(out) :: abs_od(nlwbands,pcols,pver)
         622             :
         623             :   type(interp_type) :: dei_wgts
         624             :
         625             :   integer :: i, k, lwband
         626             :   real(r8) :: absor(nlwbands)
         627             :
         628   140273568 :   do k = 1,pver
         629  2318817168 :      do i = 1,ncol
         630             :         ! if ice water path is too small, OD := 0
         631  2317324896 :         if( iciwpth(i,k) < tiny .or. dei(i,k) == 0._r8) then
         632 33827618572 :            abs_od (:,i,k) = 0._r8
         633             :         else
         634             :            ! for each cell interpolate to find weights in g_d_eff grid.
         635             :            call lininterp_init(g_d_eff, n_g_d, dei(i:i,k), 1, &
         636   188683684 :                 extrap_method_bndry, dei_wgts)
         637             :            ! interpolate into grid and extract radiative properties
         638  3207622628 :            do lwband = 1, nlwbands
         639             :               call lininterp(abs_lw_ice(:,lwband), n_g_d, &
         640  3207622628 :                    absor(lwband:lwband), 1, dei_wgts)
         641             :            enddo
         642  3207622628 :            abs_od(:,i,k) = iciwpth(i,k) * absor
         643  3207622628 :            where(abs_od(:,i,k) > 50.0_r8) abs_od(:,i,k) = 50.0_r8
         644   188683684 :            call lininterp_finish(dei_wgts)
         645             :         endif
         646             :      enddo
         647             :   enddo
         648             :
         649     1492272 : end subroutine interpolate_ice_optics_lw
         650             :
         651             : !==============================================================================
         652             :
         653   858219788 : subroutine gam_liquid_lw(clwptn, lamc, pgam, abs_od)
         654             :   real(r8), intent(in) :: clwptn ! cloud water liquid path new (in cloud) (in g/m^2)?
         655             :   real(r8), intent(in) :: lamc   ! prognosed value of lambda for cloud
         656             :   real(r8), intent(in) :: pgam   ! prognosed value of mu for cloud
         657             :   real(r8), intent(out) :: abs_od(1:nlwbands)
         658             :
         659             :   integer :: lwband ! sw band index
         660             :
         661             :   type(interp_type) :: mu_wgts
         662             :   type(interp_type) :: lambda_wgts
         663             :
         664   858219788 :   if (clwptn < tiny) then
         665   791514480 :     abs_od = 0._r8
         666   791514480 :     return
         667             :   endif
         668             :
         669    66705308 :   call get_mu_lambda_weights(lamc, pgam, mu_wgts, lambda_wgts)
         670             :
         671  1133990236 :   do lwband = 1, nlwbands
         672             :      call lininterp(abs_lw_liq(:,:,lwband), nmu, nlambda, &
         673  1133990236 :           abs_od(lwband:lwband), 1, mu_wgts, lambda_wgts)
         674             :   enddo
         675             :
         676  1133990236 :   abs_od = clwptn * abs_od
         677             :
         678    66705308 :   call lininterp_finish(mu_wgts)
         679    66705308 :   call lininterp_finish(lambda_wgts)
         680             :
         681             : end subroutine gam_liquid_lw
         682             :
         683             : !==============================================================================
         684             :
         685   858219788 : subroutine gam_liquid_sw(clwptn, lamc, pgam, tau, tau_w, tau_w_g, tau_w_f)
         686             :   real(r8), intent(in) :: clwptn ! cloud water liquid path new (in cloud) (in g/m^2)?
         687             :   real(r8), intent(in) :: lamc   ! prognosed value of lambda for cloud
         688             :   real(r8), intent(in) :: pgam   ! prognosed value of mu for cloud
         689             :   real(r8), intent(out) :: tau(1:nswbands), tau_w(1:nswbands), tau_w_f(1:nswbands), tau_w_g(1:nswbands)
         690             :
         691             :   integer :: swband ! sw band index
         692             :
         693             :   real(r8) :: ext(nswbands), ssa(nswbands), asm(nswbands)
         694             :
         695             :   type(interp_type) :: mu_wgts
         696             :   type(interp_type) :: lambda_wgts
         697             :
         698   858219788 :   if (clwptn < tiny) then
         699   791514480 :     tau = 0._r8
         700   791514480 :     tau_w = 0._r8
         701   791514480 :     tau_w_g = 0._r8
         702   791514480 :     tau_w_f = 0._r8
         703   791514480 :     return
         704             :   endif
         705             :
         706    66705308 :   call get_mu_lambda_weights(lamc, pgam, mu_wgts, lambda_wgts)
         707             :
         708  1000579620 :   do swband = 1, nswbands
         709             :      call lininterp(ext_sw_liq(:,:,swband), nmu, nlambda, &
         710   933874312 :           ext(swband:swband), 1, mu_wgts, lambda_wgts)
         711             :      call lininterp(ssa_sw_liq(:,:,swband), nmu, nlambda, &
         712   933874312 :           ssa(swband:swband), 1, mu_wgts, lambda_wgts)
         713             :      call lininterp(asm_sw_liq(:,:,swband), nmu, nlambda, &
         714  1000579620 :           asm(swband:swband), 1, mu_wgts, lambda_wgts)
         715             :   enddo
         716             :
         717             :   ! compute radiative properties
         718  1000579620 :   tau = clwptn * ext
         719  1000579620 :   tau_w = tau * ssa
         720  1000579620 :   tau_w_g = tau_w * asm
         721  1000579620 :   tau_w_f = tau_w_g * asm
         722             :
         723    66705308 :   call lininterp_finish(mu_wgts)
         724    66705308 :   call lininterp_finish(lambda_wgts)
         725             :
         726             : end subroutine gam_liquid_sw
         727             :
         728             : !==============================================================================
         729             :
         730   133410616 : subroutine get_mu_lambda_weights(lamc, pgam, mu_wgts, lambda_wgts)
         731             :   real(r8), intent(in) :: lamc   ! prognosed value of lambda for cloud
         732             :   real(r8), intent(in) :: pgam   ! prognosed value of mu for cloud
         733             :   ! Output interpolation weights. Caller is responsible for freeing these.
         734             :   type(interp_type), intent(out) :: mu_wgts
         735             :   type(interp_type), intent(out) :: lambda_wgts
         736             :
         737             :   integer :: ilambda
         738   266821232 :   real(r8) :: g_lambda_interp(nlambda)
         739             :
         740             :   ! Make interpolation weights for mu.
         741             :   ! (Put pgam in a temporary array for this purpose.)
         742   266821232 :   call lininterp_init(g_mu, nmu, [pgam], 1, extrap_method_bndry, mu_wgts)
         743             :
         744             :   ! Use mu weights to interpolate to a row in the lambda table.
         745  6803941416 :   do ilambda = 1, nlambda
         746             :      call lininterp(g_lambda(:,ilambda), nmu, &
         747  6803941416 :           g_lambda_interp(ilambda:ilambda), 1, mu_wgts)
         748             :   end do
         749             :
         750             :   ! Make interpolation weights for lambda.
         751             :   call lininterp_init(g_lambda_interp, nlambda, [lamc], 1, &
         752   266821232 :        extrap_method_bndry, lambda_wgts)
         753             :
         754   133410616 : end subroutine get_mu_lambda_weights
         755             :
         756             : !==============================================================================
         757             :
         758             : end module cloud_rad_props

Generated by: LCOV version 1.14