LCOV - code coverage report
Current view: | top level - physics/cam - slingo_liq_optics.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 9 | 84 | 10.7 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 3 | 33.3 %

              Line data    Source code
           1             : module slingo_liq_optics
           2             :
           3             : !------------------------------------------------------------------------------------------------
           4             : !  Implements Slingo Optics for MG/RRTMG for liquid clouds and
           5             : !  a copy of the old cloud routine for reference
           6             : !------------------------------------------------------------------------------------------------
           7             :
           8             : use shr_kind_mod,     only: r8 => shr_kind_r8
           9             : use physconst,        only: gravit
          10             : use ppgrid,           only: pcols, pver, pverp
          11             : use physics_types,    only: physics_state
          12             : use physics_buffer,   only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx
          13             : use radconstants,     only: nswbands, nlwbands, get_sw_spectral_boundaries
          14             : use cam_abortutils,   only: endrun
          15             :
          16             : implicit none
          17             : private
          18             : save
          19             :
          20             : public :: &
          21             :    slingo_rad_props_init,        &
          22             :    slingo_liq_get_rad_props_lw, &
          23             :    slingo_liq_optics_sw
          24             :
          25             : ! Minimum cloud amount (as a fraction of the grid-box area) to
          26             : ! distinguish from clear sky
          27             : real(r8), parameter :: cldmin = 1.0e-80_r8
          28             :
          29             : ! Decimal precision of cloud amount (0 -> preserve full resolution;
          30             : ! 10^-n -> preserve n digits of cloud amount)
          31             : real(r8), parameter :: cldeps = 0.0_r8
          32             :
          33             : ! indexes into pbuf for optical parameters of MG clouds
          34             : integer :: iclwp_idx  = 0
          35             : integer :: iciwp_idx  = 0
          36             : integer :: cld_idx    = 0
          37             : integer :: rel_idx  = 0
          38             : integer :: rei_idx  = 0
          39             :
          40             : ! indexes into constituents for old optics
          41             : integer :: &
          42             :    ixcldliq,   &         ! cloud liquid water index
          43             :    ixcldice              ! cloud liquid water index
          44             :
          45             : !==============================================================================
          46             : contains
          47             : !==============================================================================
          48             :
          49        1536 : subroutine slingo_rad_props_init()
          50             :
          51             : !   use cam_history, only: addfld
          52             :    use netcdf
          53             :    use spmd_utils,     only: masterproc
          54             :    use ioFileMod,      only: getfil
          55             :    use cam_logfile,    only: iulog
          56             :    use error_messages, only: handle_ncerr
          57             : #if ( defined SPMD )
          58             :    use mpishorthand
          59             : #endif
          60             :    use constituents,   only: cnst_get_ind
          61             :
          62             :    integer :: err
          63             :
          64        1536 :    iciwp_idx  = pbuf_get_index('ICIWP',errcode=err)
          65        1536 :    iclwp_idx  = pbuf_get_index('ICLWP',errcode=err)
          66        1536 :    cld_idx    = pbuf_get_index('CLD')
          67        1536 :    rel_idx    = pbuf_get_index('REL')
          68        1536 :    rei_idx    = pbuf_get_index('REI')
          69             :
          70             :    ! old optics
          71        1536 :    call cnst_get_ind('CLDLIQ', ixcldliq)
          72        1536 :    call cnst_get_ind('CLDICE', ixcldice)
          73             :
          74        1536 : end subroutine slingo_rad_props_init
          75             :
          76             : !==============================================================================
          77             :
          78           0 : subroutine slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp)
          79             :
          80             :    type(physics_state), intent(in) :: state
          81             :    type(physics_buffer_desc), pointer :: pbuf(:)
          82             :
          83             :    real(r8),intent(out) :: liq_tau    (nswbands,pcols,pver) ! extinction optical depth
          84             :    real(r8),intent(out) :: liq_tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
          85             :    real(r8),intent(out) :: liq_tau_w_g(nswbands,pcols,pver) ! asymmetry parameter * tau * w
          86             :    real(r8),intent(out) :: liq_tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w
          87             :    logical, intent(in) :: oldliqwp
          88             :
          89           0 :    real(r8), pointer, dimension(:,:) :: rel
          90           0 :    real(r8), pointer, dimension(:,:) :: cldn
          91           0 :    real(r8), pointer, dimension(:,:) :: tmpptr
          92             :    real(r8), dimension(pcols,pver) :: cliqwp
          93             :    real(r8), dimension(nswbands) :: wavmin
          94             :    real(r8), dimension(nswbands) :: wavmax
          95             :
          96             :    ! A. Slingo's data for cloud particle radiative properties (from 'A GCM
          97             :    ! Parameterization for the Shortwave Properties of Water Clouds' JAS
          98             :    ! vol. 46 may 1989 pp 1419-1427)
          99             :    real(r8) :: abarl(4) = &  ! A coefficient for extinction optical depth
         100             :       (/ 2.817e-02_r8, 2.682e-02_r8,2.264e-02_r8,1.281e-02_r8/)
         101             :    real(r8) :: bbarl(4) = &  ! B coefficient for extinction optical depth
         102             :       (/ 1.305_r8    , 1.346_r8    ,1.454_r8    ,1.641_r8    /)
         103             :    real(r8) :: cbarl(4) = &  ! C coefficient for single scat albedo
         104             :       (/-5.62e-08_r8 ,-6.94e-06_r8 ,4.64e-04_r8 ,0.201_r8    /)
         105             :    real(r8) :: dbarl(4) = &  ! D coefficient for single  scat albedo
         106             :       (/ 1.63e-07_r8 , 2.35e-05_r8 ,1.24e-03_r8 ,7.56e-03_r8 /)
         107             :    real(r8) :: ebarl(4) = &  ! E coefficient for asymmetry parameter
         108             :       (/ 0.829_r8    , 0.794_r8    ,0.754_r8    ,0.826_r8    /)
         109             :    real(r8) :: fbarl(4) = &  ! F coefficient for asymmetry parameter
         110             :       (/ 2.482e-03_r8, 4.226e-03_r8,6.560e-03_r8,4.353e-03_r8/)
         111             :
         112             :    real(r8) :: abarli        ! A coefficient for current spectral band
         113             :    real(r8) :: bbarli        ! B coefficient for current spectral band
         114             :    real(r8) :: cbarli        ! C coefficient for current spectral band
         115             :    real(r8) :: dbarli        ! D coefficient for current spectral band
         116             :    real(r8) :: ebarli        ! E coefficient for current spectral band
         117             :    real(r8) :: fbarli        ! F coefficient for current spectral band
         118             :
         119             :    ! Caution... A. Slingo recommends no less than 4.0 micro-meters nor
         120             :    ! greater than 20 micro-meters
         121             :
         122             :    integer :: ns, i, k, indxsl, Nday
         123             :    integer :: i_rel, lchnk, icld, itim_old
         124             :    real(r8) :: tmp1l, tmp2l, tmp3l, g
         125             :    real(r8) :: kext(pcols,pver)
         126             :    real(r8), pointer, dimension(:,:) :: iclwpth
         127             :
         128           0 :    Nday = state%ncol
         129           0 :    lchnk = state%lchnk
         130             :
         131           0 :    itim_old = pbuf_old_tim_idx()
         132           0 :    call pbuf_get_field(pbuf, cld_idx, cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
         133           0 :    call pbuf_get_field(pbuf, rel_idx, rel)
         134             :
         135           0 :    if (oldliqwp) then
         136           0 :      do k=1,pver
         137           0 :         do i = 1,Nday
         138           0 :            cliqwp(i,k) = state%q(i,k,ixcldliq)*state%pdel(i,k)/(gravit*max(0.01_r8,cldn(i,k)))
         139             :         end do
         140             :      end do
         141             :    else
         142           0 :      if (iclwp_idx<=0) then
         143           0 :         call endrun('slingo_liq_optics_sw: oldliqwp must be set to true since ICLWP was not found in pbuf')
         144             :      endif
         145             :      ! The following is the eventual target specification for in cloud liquid water path.
         146           0 :      call pbuf_get_field(pbuf, iclwp_idx, tmpptr)
         147           0 :      cliqwp = tmpptr
         148             :    endif
         149             :
         150           0 :    call get_sw_spectral_boundaries(wavmin,wavmax,'microns')
         151             :
         152           0 :    do ns = 1, nswbands
         153             :       ! Set index for cloud particle properties based on the wavelength,
         154             :       ! according to A. Slingo (1989) equations 1-3:
         155             :       ! Use index 1 (0.25 to 0.69 micrometers) for visible
         156             :       ! Use index 2 (0.69 - 1.19 micrometers) for near-infrared
         157             :       ! Use index 3 (1.19 to 2.38 micrometers) for near-infrared
         158             :       ! Use index 4 (2.38 to 4.00 micrometers) for near-infrared
         159           0 :       if(wavmax(ns) <= 0.7_r8) then
         160             :          indxsl = 1
         161           0 :       else if(wavmax(ns) <= 1.25_r8) then
         162             :          indxsl = 2
         163           0 :       else if(wavmax(ns) <= 2.38_r8) then
         164             :          indxsl = 3
         165           0 :       else if(wavmax(ns) > 2.38_r8) then
         166             :          indxsl = 4
         167             :       end if
         168             :
         169             :       ! Set cloud extinction optical depth, single scatter albedo,
         170             :       ! asymmetry parameter, and forward scattered fraction:
         171           0 :       abarli = abarl(indxsl)
         172           0 :       bbarli = bbarl(indxsl)
         173           0 :       cbarli = cbarl(indxsl)
         174           0 :       dbarli = dbarl(indxsl)
         175           0 :       ebarli = ebarl(indxsl)
         176           0 :       fbarli = fbarl(indxsl)
         177             :
         178           0 :       do k=1,pver
         179           0 :          do i=1,Nday
         180             :
         181             :             ! note that optical properties for liquid valid only
         182             :             ! in range of 4.2 > rel > 16 micron (Slingo 89)
         183           0 :             if (cldn(i,k) >= cldmin .and. cldn(i,k) >= cldeps) then
         184           0 :                tmp1l = abarli + bbarli/min(max(4.2_r8,rel(i,k)),16._r8)
         185           0 :                liq_tau(ns,i,k) = 1000._r8*cliqwp(i,k)*tmp1l
         186             :             else
         187           0 :                liq_tau(ns,i,k) = 0.0_r8
         188             :             endif
         189             :
         190           0 :             tmp2l = 1._r8 - cbarli - dbarli*min(max(4.2_r8,rel(i,k)),16._r8)
         191           0 :             tmp3l = fbarli*min(max(4.2_r8,rel(i,k)),16._r8)
         192             :             ! Do not let single scatter albedo be 1.  Delta-eddington solution
         193             :             ! for non-conservative case has different analytic form from solution
         194             :             ! for conservative case, and raddedmx is written for non-conservative case.
         195           0 :             liq_tau_w(ns,i,k) = liq_tau(ns,i,k) * min(tmp2l,.999999_r8)
         196           0 :             g = ebarli + tmp3l
         197           0 :             liq_tau_w_g(ns,i,k) = liq_tau_w(ns,i,k) * g
         198           0 :             liq_tau_w_f(ns,i,k) = liq_tau_w(ns,i,k) * g * g
         199             :
         200             :          end do ! End do i=1,Nday
         201             :       end do    ! End do k=1,pver
         202             :    end do ! nswbands
         203             :
         204           0 : end subroutine slingo_liq_optics_sw
         205             :
         206           0 : subroutine slingo_liq_get_rad_props_lw(state, pbuf, abs_od, oldliqwp)
         207             :
         208             :    type(physics_state), intent(in)  :: state
         209             :    type(physics_buffer_desc),pointer  :: pbuf(:)
         210             :    real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)
         211             :    logical, intent(in) :: oldliqwp
         212             :
         213             :    real(r8) :: gicewp(pcols,pver)
         214             :    real(r8) :: gliqwp(pcols,pver)
         215             :    real(r8) :: cicewp(pcols,pver)
         216             :    real(r8) :: cliqwp(pcols,pver)
         217             :    real(r8) :: ficemr(pcols,pver)
         218             :    real(r8) :: cwp(pcols,pver)
         219             :    real(r8) :: cldtau(pcols,pver)
         220             :
         221           0 :    real(r8), pointer, dimension(:,:) :: cldn
         222           0 :    real(r8), pointer, dimension(:,:) :: rei
         223             :    integer :: ncol, icld, itim_old, i_rei, lwband, i, k, lchnk
         224             :
         225             :     real(r8) :: kabs, kabsi
         226             :     real(r8) kabsl                  ! longwave liquid absorption coeff (m**2/g)
         227             :     parameter (kabsl = 0.090361_r8)
         228             :
         229           0 :    real(r8), pointer, dimension(:,:) :: iclwpth, iciwpth
         230             :
         231           0 :     ncol=state%ncol
         232           0 :    lchnk = state%lchnk
         233             :
         234           0 :    itim_old  =  pbuf_old_tim_idx()
         235           0 :    call pbuf_get_field(pbuf, rei_idx,   rei)
         236           0 :    call pbuf_get_field(pbuf, cld_idx,   cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
         237             :
         238           0 :    if (oldliqwp) then
         239           0 :      do k=1,pver
         240           0 :          do i = 1,ncol
         241           0 :             gicewp(i,k) = state%q(i,k,ixcldice)*state%pdel(i,k)/gravit*1000.0_r8  ! Grid box ice water path.
         242           0 :             gliqwp(i,k) = state%q(i,k,ixcldliq)*state%pdel(i,k)/gravit*1000.0_r8  ! Grid box liquid water path.
         243           0 :             cicewp(i,k) = gicewp(i,k) / max(0.01_r8,cldn(i,k))                 ! In-cloud ice water path.
         244           0 :             cliqwp(i,k) = gliqwp(i,k) / max(0.01_r8,cldn(i,k))                 ! In-cloud liquid water path.
         245             :             ficemr(i,k) = state%q(i,k,ixcldice) /                 &
         246           0 :                  max(1.e-10_r8,(state%q(i,k,ixcldice)+state%q(i,k,ixcldliq)))
         247             :          end do
         248             :      end do
         249           0 :      cwp(:ncol,:pver) = cicewp(:ncol,:pver) + cliqwp(:ncol,:pver)
         250             :    else
         251           0 :      if (iclwp_idx<=0 .or. iciwp_idx<=0) then
         252           0 :         call endrun('slingo_liq_get_rad_props_lw: oldliqwp must be set to true since ICIWP and/or ICLWP were not found in pbuf')
         253             :      endif
         254           0 :      call pbuf_get_field(pbuf, iclwp_idx, iclwpth)
         255           0 :      call pbuf_get_field(pbuf, iciwp_idx, iciwpth)
         256           0 :      do k=1,pver
         257           0 :          do i = 1,ncol
         258           0 :               cwp   (i,k) = 1000.0_r8 * iclwpth(i,k) + 1000.0_r8 * iciwpth(i, k)
         259           0 :               ficemr(i,k) = 1000.0_r8 * iciwpth(i,k)/(max(1.e-18_r8, cwp(i,k)))
         260             :          end do
         261             :      end do
         262             :    endif
         263             :
         264             :
         265           0 :    do k=1,pver
         266           0 :        do i=1,ncol
         267             :
         268             :           ! Note from Andrew Conley:
         269             :           !  Optics for RK no longer supported, This is constructed to get
         270             :           !  close to bit for bit.  Otherwise we could simply use liquid water path
         271             :           !note that optical properties for ice valid only
         272             :           !in range of 13 > rei > 130 micron (Ebert and Curry 92)
         273           0 :           kabs = kabsl*(1._r8-ficemr(i,k))
         274           0 :           cldtau(i,k) = kabs*cwp(i,k)
         275             :        end do
         276             :    end do
         277             : !
         278           0 :    do lwband = 1,nlwbands
         279           0 :       abs_od(lwband,1:ncol,1:pver)=cldtau(1:ncol,1:pver)
         280             :    enddo
         281             :
         282           0 : end subroutine slingo_liq_get_rad_props_lw
         283             :
         284             : end module slingo_liq_optics

Generated by: LCOV version 1.14