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

              Line data    Source code
           1             : module ebert_curry_ice_optics
           2             :
           3             :
           4             : use shr_kind_mod,     only: r8 => shr_kind_r8
           5             : use physconst,        only: gravit
           6             : use ppgrid,           only: pcols, pver
           7             : use physics_types,    only: physics_state
           8             : use physics_buffer,   only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx
           9             : use constituents,     only: cnst_get_ind
          10             : use radconstants,     only: nswbands, nlwbands, get_sw_spectral_boundaries
          11             : use cam_abortutils,   only: endrun
          12             :
          13             : implicit none
          14             : private
          15             : save
          16             :
          17             : public :: &
          18             :    ec_rad_props_init,        &
          19             :    ec_ice_optics_sw,       &
          20             :    ec_ice_get_rad_props_lw
          21             :
          22             :
          23             : real(r8), public, parameter:: scalefactor = 1._r8 !500._r8/917._r8
          24             :
          25             : ! indices into pbuf
          26             : integer :: iciwp_idx   = 0
          27             : integer :: iclwp_idx   = 0
          28             : integer :: cld_idx     = 0
          29             : integer :: rei_idx     = 0
          30             :
          31             : ! indices into constituents for old optics
          32             : integer :: ixcldice  ! cloud ice water index
          33             : integer :: ixcldliq  ! cloud liquid water index
          34             :
          35             :
          36             : !==============================================================================
          37             : contains
          38             : !==============================================================================
          39             :
          40        1536 : subroutine ec_rad_props_init()
          41             :
          42             :    integer :: err
          43             :
          44        1536 :    iciwp_idx  = pbuf_get_index('ICIWP',errcode=err)
          45        1536 :    iclwp_idx  = pbuf_get_index('ICLWP',errcode=err)
          46        1536 :    cld_idx    = pbuf_get_index('CLD')
          47        1536 :    rei_idx    = pbuf_get_index('REI')
          48             :
          49             :    ! old optics
          50        1536 :    call cnst_get_ind('CLDICE', ixcldice)
          51        1536 :    call cnst_get_ind('CLDLIQ', ixcldliq)
          52             :
          53        1536 : end subroutine ec_rad_props_init
          54             :
          55             : !==============================================================================
          56             :
          57           0 : subroutine ec_ice_optics_sw   (state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp)
          58             :
          59             :    type(physics_state),    intent(in) :: state
          60             :    type(physics_buffer_desc), pointer :: pbuf(:)
          61             :
          62             :    real(r8),intent(out) :: ice_tau    (nswbands,pcols,pver) ! extinction optical depth
          63             :    real(r8),intent(out) :: ice_tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
          64             :    real(r8),intent(out) :: ice_tau_w_g(nswbands,pcols,pver) ! asymmetry parameter * tau * w
          65             :    real(r8),intent(out) :: ice_tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w
          66             :    logical, intent(in) :: oldicewp
          67             :
          68           0 :    real(r8), pointer, dimension(:,:) :: rei
          69           0 :    real(r8), pointer, dimension(:,:) :: cldn
          70           0 :    real(r8), pointer, dimension(:,:) :: tmpptr
          71             :    real(r8), dimension(pcols,pver) :: cicewp
          72             :    real(r8), dimension(nswbands) :: wavmin
          73             :    real(r8), dimension(nswbands) :: wavmax
          74             :    !
          75             :    ! ice water coefficients (Ebert and Curry,1992, JGR, 97, 3831-3836)
          76             :    real(r8) :: abari(4) = &     ! a coefficient for extinction optical depth
          77             :       (/ 3.448e-03_r8, 3.448e-03_r8,3.448e-03_r8,3.448e-03_r8/)
          78             :    real(r8) :: bbari(4) = &     ! b coefficient for extinction optical depth
          79             :       (/ 2.431_r8    , 2.431_r8    ,2.431_r8    ,2.431_r8    /)
          80             :    real(r8) :: cbari(4) = &     ! c coefficient for single scat albedo
          81             :       (/ 1.00e-05_r8 , 1.10e-04_r8 ,1.861e-02_r8,.46658_r8   /)
          82             :    real(r8) :: dbari(4) = &     ! d coefficient for single scat albedo
          83             :       (/ 0.0_r8      , 1.405e-05_r8,8.328e-04_r8,2.05e-05_r8 /)
          84             :    real(r8) :: ebari(4) = &     ! e coefficient for asymmetry parameter
          85             :       (/ 0.7661_r8   , 0.7730_r8   ,0.794_r8    ,0.9595_r8   /)
          86             :    real(r8) :: fbari(4) = &     ! f coefficient for asymmetry parameter
          87             :       (/ 5.851e-04_r8, 5.665e-04_r8,7.267e-04_r8,1.076e-04_r8/)
          88             :
          89             :    real(r8) :: abarii           ! A coefficient for current spectral band
          90             :    real(r8) :: bbarii           ! B coefficient for current spectral band
          91             :    real(r8) :: cbarii           ! C coefficient for current spectral band
          92             :    real(r8) :: dbarii           ! D coefficient for current spectral band
          93             :    real(r8) :: ebarii           ! E coefficient for current spectral band
          94             :    real(r8) :: fbarii           ! F coefficient for current spectral band
          95             :
          96             :    ! Minimum cloud amount (as a fraction of the grid-box area) to
          97             :    ! distinguish from clear sky
          98             :    real(r8), parameter :: cldmin = 1.0e-80_r8
          99             :
         100             :    ! Decimal precision of cloud amount (0 -> preserve full resolution;
         101             :    ! 10^-n -> preserve n digits of cloud amount)
         102             :    real(r8), parameter :: cldeps = 0.0_r8
         103             :
         104             :    integer :: ns, i, k, indxsl, lchnk, Nday
         105             :    integer :: itim_old
         106             :    real(r8) :: tmp1i, tmp2i, tmp3i, g
         107             :
         108           0 :    Nday = state%ncol
         109           0 :    lchnk = state%lchnk
         110             :
         111           0 :    itim_old = pbuf_old_tim_idx()
         112           0 :    call pbuf_get_field(pbuf, cld_idx,cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
         113           0 :    call pbuf_get_field(pbuf, rei_idx,rei)
         114             :
         115           0 :    if(oldicewp) then
         116           0 :      do k=1,pver
         117           0 :         do i = 1,Nday
         118           0 :            cicewp(i,k) = 1000.0_r8*state%q(i,k,ixcldice)*state%pdel(i,k) /(gravit* max(0.01_r8,cldn(i,k)))
         119             :         end do
         120             :      end do
         121             :    else
         122           0 :      if (iciwp_idx<=0) then
         123           0 :         call endrun('ec_ice_optics_sw: oldicewp must be set to true since ICIWP was not found in pbuf')
         124             :      endif
         125           0 :      call pbuf_get_field(pbuf, iciwp_idx, tmpptr)
         126           0 :      cicewp(1:pcols,1:pver) =  1000.0_r8*tmpptr(1:pcols,1:pver)
         127             :    endif
         128             :
         129           0 :    call get_sw_spectral_boundaries(wavmin,wavmax,'microns')
         130             :
         131           0 :    do ns = 1, nswbands
         132             :
         133           0 :       if(wavmax(ns) <= 0.7_r8) then
         134             :          indxsl = 1
         135           0 :       else if(wavmax(ns) <= 1.25_r8) then
         136             :          indxsl = 2
         137           0 :       else if(wavmax(ns) <= 2.38_r8) then
         138             :          indxsl = 3
         139           0 :       else if(wavmax(ns) > 2.38_r8) then
         140             :          indxsl = 4
         141             :       end if
         142             :
         143           0 :       abarii = abari(indxsl)
         144           0 :       bbarii = bbari(indxsl)
         145           0 :       cbarii = cbari(indxsl)
         146           0 :       dbarii = dbari(indxsl)
         147           0 :       ebarii = ebari(indxsl)
         148           0 :       fbarii = fbari(indxsl)
         149             :
         150           0 :       do k=1,pver
         151           0 :          do i=1,Nday
         152             :
         153             :             ! note that optical properties for ice valid only
         154             :             ! in range of 13 > rei > 130 micron (Ebert and Curry 92)
         155           0 :             if (cldn(i,k) >= cldmin .and. cldn(i,k) >= cldeps) then
         156           0 :                tmp1i = abarii + bbarii/max(13._r8,min(scalefactor*rei(i,k),130._r8))
         157           0 :                ice_tau(ns,i,k) = cicewp(i,k)*tmp1i
         158             :             else
         159           0 :                ice_tau(ns,i,k) = 0.0_r8
         160             :             endif
         161             :
         162           0 :             tmp2i = 1._r8 - cbarii - dbarii*min(max(13._r8,scalefactor*rei(i,k)),130._r8)
         163           0 :             tmp3i = fbarii*min(max(13._r8,scalefactor*rei(i,k)),130._r8)
         164             :             ! Do not let single scatter albedo be 1.  Delta-eddington solution
         165             :             ! for non-conservative case has different analytic form from solution
         166             :             ! for conservative case, and raddedmx is written for non-conservative case.
         167           0 :             ice_tau_w(ns,i,k) = ice_tau(ns,i,k) * min(tmp2i,.999999_r8)
         168           0 :             g = ebarii + tmp3i
         169           0 :             ice_tau_w_g(ns,i,k) = ice_tau_w(ns,i,k) * g
         170           0 :             ice_tau_w_f(ns,i,k) = ice_tau_w(ns,i,k) * g * g
         171             :
         172             :          end do ! End do i=1,Nday
         173             :       end do    ! End do k=1,pver
         174             :    end do ! nswbands
         175             :
         176           0 : end subroutine ec_ice_optics_sw
         177             :
         178             : !==============================================================================
         179             :
         180           0 : subroutine ec_ice_get_rad_props_lw(state, pbuf, abs_od, oldicewp)
         181             :
         182             :    type(physics_state), intent(in)   :: state
         183             :    type(physics_buffer_desc),pointer :: pbuf(:)
         184             :     real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)
         185             :    logical, intent(in) :: oldicewp
         186             :
         187             :    real(r8) :: gicewp(pcols,pver)
         188             :    real(r8) :: gliqwp(pcols,pver)
         189             :    real(r8) :: cicewp(pcols,pver)
         190             :    real(r8) :: cliqwp(pcols,pver)
         191             :    real(r8) :: ficemr(pcols,pver)
         192             :    real(r8) :: cwp(pcols,pver)
         193             :    real(r8) :: cldtau(pcols,pver)
         194             :
         195           0 :    real(r8), pointer, dimension(:,:) :: cldn
         196           0 :    real(r8), pointer, dimension(:,:) :: rei
         197             :    integer :: ncol, itim_old, lwband, i, k, lchnk
         198             :
         199             :     real(r8) :: kabs, kabsi
         200             :
         201             :     real(r8) kabsl                  ! longwave liquid absorption coeff (m**2/g)
         202             :     parameter (kabsl = 0.090361_r8)
         203             :
         204           0 :    real(r8), pointer, dimension(:,:) :: iclwpth, iciwpth
         205             :
         206             :
         207           0 :    ncol = state%ncol
         208           0 :    lchnk = state%lchnk
         209             :
         210           0 :    itim_old  =  pbuf_old_tim_idx()
         211           0 :    call pbuf_get_field(pbuf, rei_idx,   rei)
         212           0 :    call pbuf_get_field(pbuf, cld_idx,   cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
         213             :
         214             :
         215           0 :    if(oldicewp) then
         216           0 :      do k=1,pver
         217           0 :          do i = 1,ncol
         218           0 :             gicewp(i,k) = state%q(i,k,ixcldice)*state%pdel(i,k)/gravit*1000.0_r8  ! Grid box ice water path.
         219           0 :             gliqwp(i,k) = state%q(i,k,ixcldliq)*state%pdel(i,k)/gravit*1000.0_r8  ! Grid box liquid water path.
         220           0 :             cicewp(i,k) = gicewp(i,k) / max(0.01_r8,cldn(i,k))                 ! In-cloud ice water path.
         221           0 :             cliqwp(i,k) = gliqwp(i,k) / max(0.01_r8,cldn(i,k))                 ! In-cloud liquid water path.
         222             :             ficemr(i,k) = state%q(i,k,ixcldice) /                 &
         223           0 :                  max(1.e-10_r8,(state%q(i,k,ixcldice)+state%q(i,k,ixcldliq)))
         224             :          end do
         225             :      end do
         226           0 :      cwp(:ncol,:pver) = cicewp(:ncol,:pver) + cliqwp(:ncol,:pver)
         227             :    else
         228           0 :       if (iclwp_idx<=0 .or. iciwp_idx<=0) then
         229           0 :          call endrun('ec_ice_get_rad_props_lw: oldicewp must be set to true since ICIWP and/or ICLWP were not found in pbuf')
         230             :       endif
         231           0 :       call pbuf_get_field(pbuf, iclwp_idx, iclwpth)
         232           0 :       call pbuf_get_field(pbuf, iciwp_idx, iciwpth)
         233           0 :       do k=1,pver
         234           0 :          do i = 1,ncol
         235           0 :             cwp(i,k) = 1000.0_r8 *iciwpth(i,k) + 1000.0_r8 *iclwpth(i,k)
         236           0 :             ficemr(i,k) = 1000.0_r8*iciwpth(i,k)/(max(1.e-18_r8,cwp(i,k)))
         237             :          end do
         238             :       end do
         239             :    endif
         240             :
         241           0 :    do k=1,pver
         242           0 :        do i=1,ncol
         243             :
         244             :           ! Note from Andrew Conley:
         245             :           !  Optics for RK no longer supported, This is constructed to get
         246             :           !  close to bit for bit.  Otherwise we could simply use ice water path
         247             :           !note that optical properties for ice valid only
         248             :           !in range of 13 > rei > 130 micron (Ebert and Curry 92)
         249           0 :           kabsi = 0.005_r8 + 1._r8/min(max(13._r8,scalefactor*rei(i,k)),130._r8)
         250           0 :           kabs =  kabsi*ficemr(i,k) ! kabsl*(1._r8-ficemr(i,k)) + kabsi*ficemr(i,k)
         251             :           !emis(i,k) = 1._r8 - exp(-1.66_r8*kabs*clwp(i,k))
         252           0 :           cldtau(i,k) = kabs*cwp(i,k)
         253             :        end do
         254             :    end do
         255             :
         256           0 :    do lwband = 1,nlwbands
         257           0 :       abs_od(lwband,1:ncol,1:pver)=cldtau(1:ncol,1:pver)
         258             :    enddo
         259             :
         260           0 : end subroutine ec_ice_get_rad_props_lw
         261             :
         262             : !==============================================================================
         263             :
         264             : end module ebert_curry_ice_optics

Generated by: LCOV version 1.14