LCOV - code coverage report
Current view: | top level - chemistry/aerosol - refractive_aerosol_optics_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 156 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 7 | 0.0 %

              Line data    Source code
           1             : module refractive_aerosol_optics_mod
           2             :   use shr_kind_mod, only: r8 => shr_kind_r8
           3             :   use aerosol_optics_mod, only: aerosol_optics
           4             :   use physconst,  only: rhoh2o
           5             :   use aerosol_state_mod, only: aerosol_state
           6             :   use aerosol_properties_mod, only: aerosol_properties
           7             :
           8             :   use table_interp_mod, only: table_interp, table_interp_wghts, table_interp_calcwghts
           9             :
          10             :   implicit none
          11             :
          12             :   private
          13             :   public :: refractive_aerosol_optics
          14             :
          15             :   !> refractive_aerosol_optics
          16             :   !! Table look up implementation of aerosol_optics to parameterize aerosol radiative properties in terms of
          17             :   !! surface mode wet radius and wet refractive index using chebychev polynomials
          18             :   type, extends(aerosol_optics) :: refractive_aerosol_optics
          19             :
          20             :      integer :: ibin, ilist
          21             :      class(aerosol_state), pointer :: aero_state      ! aerosol_state object
          22             :      class(aerosol_properties), pointer :: aero_props ! aerosol_properties object
          23             :
          24             :      real(r8), allocatable :: watervol(:,:)   ! volume concentration of water in each mode (m3/kg)
          25             :      real(r8), allocatable :: wetvol(:,:)     ! volume concentration of wet mode (m3/kg)
          26             :      real(r8), allocatable :: cheb(:,:,:)     ! chebychev polynomials
          27             :      real(r8), allocatable :: radsurf(:,:)    ! aerosol surface mode radius
          28             :      real(r8), allocatable :: logradsurf(:,:) ! log(aerosol surface mode radius)
          29             :
          30             :      ! refractive index for water read in read_water_refindex
          31             :      complex(r8), allocatable :: crefwsw(:) ! complex refractive index for water visible
          32             :      complex(r8), allocatable :: crefwlw(:) ! complex refractive index for water infrared
          33             :
          34             :      real(r8), pointer :: extpsw(:,:,:,:) => null() ! specific extinction
          35             :      real(r8), pointer :: abspsw(:,:,:,:) => null() ! specific absorption
          36             :      real(r8), pointer :: asmpsw(:,:,:,:) => null() ! asymmetry factor
          37             :      real(r8), pointer :: absplw(:,:,:,:) => null() ! specific absorption
          38             :
          39             :      real(r8), pointer :: refrtabsw(:,:) => null()  ! table of real refractive indices for aerosols
          40             :      real(r8), pointer :: refitabsw(:,:) => null()  ! table of imag refractive indices for aerosols
          41             :      real(r8), pointer :: refrtablw(:,:) => null()  ! table of real refractive indices for aerosols
          42             :      real(r8), pointer :: refitablw(:,:) => null()  ! table of imag refractive indices for aerosols
          43             :
          44             :      ! Dimension sizes in coefficient arrays used to parameterize aerosol radiative properties
          45             :      ! in terms of refractive index and wet radius
          46             :      integer :: ncoef = -1  ! number of chebychev coeficients
          47             :      integer :: prefr = -1  ! number of real refractive indices
          48             :      integer :: prefi = -1  ! number of imaginary refractive indices
          49             :
          50             :    contains
          51             :
          52             :      procedure :: sw_props
          53             :      procedure :: lw_props
          54             :
          55             :      final :: destructor
          56             :
          57             :   end type refractive_aerosol_optics
          58             :
          59             :   interface refractive_aerosol_optics
          60             :      procedure :: constructor
          61             :   end interface refractive_aerosol_optics
          62             :
          63             :   ! radius limits (m)
          64             :   real(r8), parameter :: radmin = 0.01e-6_r8 ! min aerosol surface mode radius (m)
          65             :   real(r8), parameter :: radmax = 25.e-6_r8  ! max aerosol surface mode radius (m)
          66             :   real(r8), parameter :: xrmin=log(radmin)   ! min log(aerosol surface mode radius)
          67             :   real(r8), parameter :: xrmax=log(radmax)   ! max log(aerosol surface mode radius)
          68             :
          69             : contains
          70             :
          71             :   !------------------------------------------------------------------------------
          72             :   !------------------------------------------------------------------------------
          73           0 :   function constructor(aero_props, aero_state, ilist, ibin, ncol, nlev, nsw, nlw, crefwsw, crefwlw) &
          74             :        result(newobj)
          75             :
          76             :     class(aerosol_properties),intent(in), target :: aero_props   ! aerosol_properties object
          77             :     class(aerosol_state),intent(in), target :: aero_state        ! aerosol_state object
          78             :     integer, intent(in) :: ilist  ! climate or a diagnostic list number
          79             :     integer, intent(in) :: ibin   ! bin number
          80             :     integer, intent(in) :: ncol   ! number of columns
          81             :     integer, intent(in) :: nlev   ! number of levels
          82             :     integer, intent(in) :: nsw    ! number of short wave lengths
          83             :     integer, intent(in) :: nlw    ! number of long wave lengths
          84             :     complex(r8), intent(in) :: crefwsw(nsw) ! complex refractive index for water visible
          85             :     complex(r8), intent(in) :: crefwlw(nlw) ! complex refractive index for water infrared
          86             :
          87             :     type(refractive_aerosol_optics), pointer :: newobj
          88             :
          89             :     integer :: ierr, icol, ilev, ispec, nspec
          90           0 :     real(r8) :: vol(ncol)             ! volume concentration of aerosol species (m3/kg)
          91           0 :     real(r8) :: dryvol(ncol)          ! volume concentration of aerosol mode (m3/kg)
          92             :     real(r8) :: specdens              ! species density (kg/m3)
          93           0 :     real(r8), pointer :: specmmr(:,:) ! species mass mixing ratio
          94             :     real(r8) :: logsigma              ! geometric standard deviation of number distribution
          95             :
          96           0 :     real(r8) :: dgnumwet(ncol,nlev)   ! aerosol wet number mode diameter (m)
          97           0 :     real(r8) :: qaerwat(ncol,nlev)    ! aerosol water (g/g)
          98             :
          99             :     real(r8), parameter :: rh2odens = 1._r8/rhoh2o
         100             :
         101           0 :     allocate(newobj, stat=ierr)
         102           0 :     if (ierr/=0) then
         103           0 :        nullify(newobj)
         104             :        return
         105             :     end if
         106             :
         107             :     ! get mode properties
         108             :     call aero_props%optics_params(ilist, ibin, &
         109             :          refrtabsw=newobj%refrtabsw, refitabsw=newobj%refitabsw, &
         110             :          refrtablw=newobj%refrtablw, refitablw=newobj%refitablw,&
         111             :          extpsw=newobj%extpsw, abspsw=newobj%abspsw, asmpsw=newobj%asmpsw, &
         112           0 :          absplw=newobj%absplw, ncoef=newobj%ncoef, prefr=newobj%prefr, prefi=newobj%prefi)
         113             :
         114           0 :     allocate(newobj%watervol(ncol,nlev),stat=ierr)
         115           0 :     if (ierr/=0) then
         116           0 :        nullify(newobj)
         117           0 :        return
         118             :     end if
         119           0 :     allocate(newobj%wetvol(ncol,nlev),stat=ierr)
         120           0 :     if (ierr/=0) then
         121           0 :        nullify(newobj)
         122           0 :        return
         123             :     end if
         124           0 :     allocate(newobj%cheb(newobj%ncoef,ncol,nlev),stat=ierr)
         125           0 :     if (ierr/=0) then
         126           0 :        nullify(newobj)
         127           0 :        return
         128             :     end if
         129           0 :     allocate(newobj%radsurf(ncol,nlev),stat=ierr)
         130           0 :     if (ierr/=0) then
         131           0 :        nullify(newobj)
         132           0 :        return
         133             :     end if
         134           0 :     allocate(newobj%logradsurf(ncol,nlev),stat=ierr)
         135           0 :     if (ierr/=0) then
         136           0 :        nullify(newobj)
         137           0 :        return
         138             :     end if
         139             :
         140           0 :     allocate(newobj%crefwlw(nlw),stat=ierr)
         141           0 :     if (ierr/=0) then
         142           0 :        nullify(newobj)
         143           0 :        return
         144             :     end if
         145           0 :     newobj%crefwlw(:) = crefwlw(:)
         146             :
         147           0 :     allocate(newobj%crefwsw(nsw),stat=ierr)
         148           0 :     if (ierr/=0) then
         149           0 :        nullify(newobj)
         150           0 :        return
         151             :     end if
         152           0 :     newobj%crefwsw(:) = crefwsw(:)
         153             :
         154           0 :     call aero_state%water_uptake(aero_props, ilist, ibin,  ncol, nlev, dgnumwet, qaerwat)
         155             :
         156           0 :     nspec = aero_props%nspecies(ilist,ibin)
         157             :
         158           0 :     logsigma=aero_props%alogsig(ilist,ibin)
         159             :
         160             :     ! calc size parameter for all columns
         161             :     call modal_size_parameters(newobj%ncoef, ncol, nlev, logsigma, dgnumwet, &
         162           0 :                                newobj%radsurf, newobj%logradsurf, newobj%cheb)
         163             :
         164           0 :     do ilev = 1, nlev
         165           0 :        dryvol(:ncol) = 0._r8
         166           0 :        do ispec = 1, nspec
         167           0 :           call aero_state%get_ambient_mmr(ilist,ispec,ibin,specmmr)
         168           0 :           call aero_props%get(ibin, ispec, list_ndx=ilist, density=specdens)
         169             :
         170           0 :           do icol = 1, ncol
         171           0 :              vol(icol) = specmmr(icol,ilev)/specdens
         172           0 :              dryvol(icol) = dryvol(icol) + vol(icol)
         173             :
         174           0 :              newobj%watervol(icol,ilev) = qaerwat(icol,ilev)*rh2odens
         175           0 :              newobj%wetvol(icol,ilev) = newobj%watervol(icol,ilev) + dryvol(icol)
         176           0 :              if (newobj%watervol(icol,ilev) < 0._r8) then
         177           0 :                 newobj%watervol(icol,ilev) = 0._r8
         178           0 :                 newobj%wetvol(icol,ilev) = dryvol(icol)
         179             :              end if
         180             :           end do
         181             :        end do
         182             :     end do
         183             :
         184           0 :     newobj%aero_state => aero_state
         185           0 :     newobj%aero_props => aero_props
         186           0 :     newobj%ilist = ilist
         187           0 :     newobj%ibin = ibin
         188             :
         189           0 :   end function constructor
         190             :
         191             :   !------------------------------------------------------------------------------
         192             :   ! returns short wave aerosol optics properties
         193             :   !------------------------------------------------------------------------------
         194           0 :   subroutine sw_props(self, ncol, ilev, iwav, pext, pabs, palb, pasm)
         195             :
         196             :     class(refractive_aerosol_optics), intent(in) :: self
         197             :     integer, intent(in) :: ncol        ! number of columns
         198             :     integer, intent(in) :: ilev        ! vertical level index
         199             :     integer, intent(in) :: iwav        ! wave length index
         200             :     real(r8),intent(out) :: pext(ncol) ! parameterized specific extinction (m2/kg)
         201             :     real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg)
         202             :     real(r8),intent(out) :: palb(ncol) ! parameterized asymmetry factor
         203             :     real(r8),intent(out) :: pasm(ncol) ! parameterized single scattering albedo
         204             :
         205           0 :     real(r8) :: refr(ncol)      ! real part of refractive index
         206           0 :     real(r8) :: refi(ncol)      ! imaginary part of refractive index
         207           0 :     real(r8) :: cext(self%ncoef,ncol), cabs(self%ncoef,ncol), casm(self%ncoef,ncol)
         208             :
         209           0 :     complex(r8) :: crefin(ncol) ! complex refractive index
         210             :     integer :: icol,icoef
         211             :
         212           0 :     type(table_interp_wghts) :: wghtsr(ncol)
         213           0 :     type(table_interp_wghts) :: wghtsi(ncol)
         214             :
         215           0 :     crefin(:ncol) = self%aero_state%refractive_index_sw(ncol, ilev, self%ilist, self%ibin, iwav, self%aero_props)
         216             :
         217           0 :     do icol = 1, ncol
         218           0 :        crefin(icol) = crefin(icol) + self%watervol(icol,ilev)*self%crefwsw(iwav)
         219           0 :        crefin(icol) = crefin(icol)/max(self%wetvol(icol,ilev),1.e-60_r8)
         220           0 :        refr(icol) = real(crefin(icol))
         221           0 :        refi(icol) = abs(aimag(crefin(icol)))
         222             :     end do
         223             :
         224             :     ! interpolate coefficients linear in refractive index
         225             :
         226           0 :     wghtsr = table_interp_calcwghts( self%prefr, self%refrtabsw(:,iwav), ncol, refr(:ncol) )
         227           0 :     wghtsi = table_interp_calcwghts( self%prefi, self%refitabsw(:,iwav), ncol, refi(:ncol) )
         228             :
         229           0 :     cext(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%extpsw(:,:,:,iwav))
         230           0 :     cabs(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%abspsw(:,:,:,iwav))
         231           0 :     casm(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%asmpsw(:,:,:,iwav))
         232             :
         233           0 :     do icol = 1,ncol
         234             :
         235           0 :        if (self%logradsurf(icol,ilev) <= xrmax) then
         236           0 :           pext(icol) = 0.5_r8*cext(1,icol)
         237           0 :           do icoef = 2, self%ncoef
         238           0 :              pext(icol) = pext(icol) + self%cheb(icoef,icol,ilev)*cext(icoef,icol)
         239             :           enddo
         240           0 :           pext(icol) = exp(pext(icol))
         241             :        else
         242           0 :           pext(icol) = 1.5_r8/(self%radsurf(icol,ilev)*rhoh2o) ! geometric optics
         243             :        endif
         244             :
         245             :        ! convert from m2/kg water to m2/kg aerosol
         246           0 :        pext(icol) = pext(icol)*self%wetvol(icol,ilev)*rhoh2o
         247           0 :        pabs(icol) = 0.5_r8*cabs(1,icol)
         248           0 :        pasm(icol) = 0.5_r8*casm(1,icol)
         249           0 :        do icoef = 2, self%ncoef
         250           0 :           pabs(icol) = pabs(icol) + self%cheb(icoef,icol,ilev)*cabs(icoef,icol)
         251           0 :           pasm(icol) = pasm(icol) + self%cheb(icoef,icol,ilev)*casm(icoef,icol)
         252             :        enddo
         253           0 :        pabs(icol) = pabs(icol)*self%wetvol(icol,ilev)*rhoh2o
         254           0 :        pabs(icol) = max(0._r8,pabs(icol))
         255           0 :        pabs(icol) = min(pext(icol),pabs(icol))
         256             :
         257           0 :        palb(icol) = 1._r8-pabs(icol)/max(pext(icol),1.e-40_r8)
         258             :
         259             :     end do
         260             :
         261           0 :   end subroutine sw_props
         262             :
         263             :   !------------------------------------------------------------------------------
         264             :   ! returns long wave aerosol optics properties
         265             :   !------------------------------------------------------------------------------
         266           0 :   subroutine lw_props(self, ncol, ilev, iwav, pabs)
         267             :
         268             :     class(refractive_aerosol_optics), intent(in) :: self
         269             :     integer, intent(in) :: ncol        ! number of columns
         270             :     integer, intent(in) :: ilev        ! vertical level index
         271             :     integer, intent(in) :: iwav        ! wave length index
         272             :     real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg)
         273             :
         274           0 :     real(r8) :: refr(ncol)      ! real part of refractive index
         275           0 :     real(r8) :: refi(ncol)      ! imaginary part of refractive index
         276           0 :     real(r8) :: cabs(self%ncoef,ncol)
         277             :
         278           0 :     complex(r8) :: crefin(ncol) ! complex refractive index
         279             :     integer :: icol, icoef
         280             :
         281           0 :     type(table_interp_wghts) :: wghtsr(ncol)
         282           0 :     type(table_interp_wghts) :: wghtsi(ncol)
         283             :
         284           0 :     crefin(:ncol) = self%aero_state%refractive_index_lw(ncol, ilev, self%ilist, self%ibin, iwav, self%aero_props)
         285             :
         286           0 :     do icol = 1, ncol
         287           0 :        crefin(icol) = crefin(icol) + self%watervol(icol,ilev)*self%crefwlw(iwav)
         288           0 :        crefin(icol) = crefin(icol)/max(self%wetvol(icol,ilev), 1.e-40_r8)
         289             :
         290           0 :        refr(icol) = real(crefin(icol))
         291           0 :        refi(icol) = aimag(crefin(icol))
         292             :
         293             :     end do
         294             :
         295             :     ! interpolate coefficients linear in refractive index
         296             :
         297           0 :     wghtsr = table_interp_calcwghts( self%prefr, self%refrtablw(:,iwav), ncol, refr(:ncol) )
         298           0 :     wghtsi = table_interp_calcwghts( self%prefi, self%refitablw(:,iwav), ncol, refi(:ncol) )
         299             :
         300           0 :     cabs(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%absplw(:,:,:,iwav))
         301             :
         302           0 :     do icol = 1,ncol
         303           0 :        pabs(icol) = 0.5_r8*cabs(1,icol)
         304           0 :        do icoef = 2, self%ncoef
         305           0 :           pabs(icol) = pabs(icol) + self%cheb(icoef,icol,ilev)*cabs(icoef,icol)
         306             :        end do
         307           0 :        pabs(icol) = pabs(icol)*self%wetvol(icol,ilev)*rhoh2o
         308           0 :        pabs(icol) = max(0._r8,pabs(icol))
         309             :     end do
         310             :
         311           0 :   end subroutine lw_props
         312             :
         313             :
         314             :   !------------------------------------------------------------------------------
         315             :   !------------------------------------------------------------------------------
         316           0 :   subroutine destructor(self)
         317             :
         318             :     type(refractive_aerosol_optics), intent(inout) :: self
         319             :
         320           0 :     deallocate(self%watervol)
         321           0 :     deallocate(self%wetvol)
         322           0 :     deallocate(self%cheb)
         323           0 :     deallocate(self%radsurf)
         324           0 :     deallocate(self%logradsurf)
         325           0 :     deallocate(self%crefwsw)
         326           0 :     deallocate(self%crefwlw)
         327             :
         328           0 :     nullify(self%aero_state)
         329           0 :     nullify(self%aero_props)
         330           0 :     nullify(self%extpsw)
         331           0 :     nullify(self%abspsw)
         332           0 :     nullify(self%asmpsw)
         333           0 :     nullify(self%absplw)
         334           0 :     nullify(self%refrtabsw)
         335           0 :     nullify(self%refitabsw)
         336           0 :     nullify(self%refrtablw)
         337           0 :     nullify(self%refitablw)
         338             :
         339           0 :   end subroutine destructor
         340             :
         341             :
         342             :   ! Private routines
         343             :   !===============================================================================
         344             :
         345             :   !===============================================================================
         346             :
         347           0 :   subroutine modal_size_parameters(ncoef,ncol,nlev, alnsg_amode, dgnumwet, radsurf, logradsurf, cheb)
         348             :
         349             :     integer,  intent(in)  :: ncoef,ncol,nlev
         350             :     real(r8), intent(in)  :: alnsg_amode     ! geometric standard deviation of number distribution
         351             :     real(r8), intent(in)  :: dgnumwet(:,:)   ! aerosol wet number mode diameter (m)
         352             :     real(r8), intent(out) :: radsurf(:,:)    ! aerosol surface mode radius
         353             :     real(r8), intent(out) :: logradsurf(:,:) ! log(aerosol surface mode radius)
         354             :     real(r8), intent(out) :: cheb(:,:,:)
         355             :
         356             :     integer  :: i, k, nc
         357             :     real(r8) :: explnsigma
         358           0 :     real(r8) :: xrad(ncol) ! normalized aerosol radius
         359             :
         360             :     !-------------------------------------------------------------------------------
         361             :
         362           0 :     explnsigma = exp(2.0_r8*alnsg_amode*alnsg_amode)
         363             :
         364           0 :     do k = 1, nlev
         365           0 :        do i = 1, ncol
         366             :           ! convert from number mode diameter to surface area
         367           0 :           radsurf(i,k) = max(0.5_r8*dgnumwet(i,k)*explnsigma,radmin)
         368           0 :           logradsurf(i,k) = log(radsurf(i,k))
         369             :           ! normalize size parameter
         370           0 :           xrad(i) = max(logradsurf(i,k),xrmin)
         371           0 :           xrad(i) = min(xrad(i),xrmax)
         372           0 :           xrad(i) = (2._r8*xrad(i)-xrmax-xrmin)/(xrmax-xrmin)
         373             :           ! chebyshev polynomials
         374           0 :           cheb(1,i,k) = 1._r8
         375           0 :           cheb(2,i,k) = xrad(i)
         376           0 :           do nc = 3, ncoef
         377           0 :              cheb(nc,i,k) = 2._r8*xrad(i)*cheb(nc-1,i,k)-cheb(nc-2,i,k)
         378             :           end do
         379             :        end do
         380             :     end do
         381             :
         382           0 :   end subroutine modal_size_parameters
         383             :
         384           0 : end module refractive_aerosol_optics_mod

Generated by: LCOV version 1.14