LCOV - code coverage report
Current view: | top level - chemistry/aerosol - refractive_aerosol_optics_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 139 | 156 | 89.1 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 6 | 7 | 85.7 %

              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     5969088 :   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    11938176 :     real(r8) :: vol(ncol)             ! volume concentration of aerosol species (m3/kg)
          91    11938176 :     real(r8) :: dryvol(ncol)          ! volume concentration of aerosol mode (m3/kg)
          92             :     real(r8) :: specdens              ! species density (kg/m3)
          93     5969088 :     real(r8), pointer :: specmmr(:,:) ! species mass mixing ratio
          94             :     real(r8) :: logsigma              ! geometric standard deviation of number distribution
          95             :
          96    11938176 :     real(r8) :: dgnumwet(ncol,nlev)   ! aerosol wet number mode diameter (m)
          97    11938176 :     real(r8) :: qaerwat(ncol,nlev)    ! aerosol water (g/g)
          98             :
          99             :     real(r8), parameter :: rh2odens = 1._r8/rhoh2o
         100             :
         101     5969088 :     allocate(newobj, stat=ierr)
         102     5969088 :     if (ierr/=0) then
         103     5969088 :        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     5969088 :          absplw=newobj%absplw, ncoef=newobj%ncoef, prefr=newobj%prefr, prefi=newobj%prefi)
         113             :
         114    23876352 :     allocate(newobj%watervol(ncol,nlev),stat=ierr)
         115     5969088 :     if (ierr/=0) then
         116           0 :        nullify(newobj)
         117           0 :        return
         118             :     end if
         119    17907264 :     allocate(newobj%wetvol(ncol,nlev),stat=ierr)
         120     5969088 :     if (ierr/=0) then
         121           0 :        nullify(newobj)
         122           0 :        return
         123             :     end if
         124    29845440 :     allocate(newobj%cheb(newobj%ncoef,ncol,nlev),stat=ierr)
         125     5969088 :     if (ierr/=0) then
         126           0 :        nullify(newobj)
         127           0 :        return
         128             :     end if
         129    17907264 :     allocate(newobj%radsurf(ncol,nlev),stat=ierr)
         130     5969088 :     if (ierr/=0) then
         131           0 :        nullify(newobj)
         132           0 :        return
         133             :     end if
         134    17907264 :     allocate(newobj%logradsurf(ncol,nlev),stat=ierr)
         135     5969088 :     if (ierr/=0) then
         136           0 :        nullify(newobj)
         137           0 :        return
         138             :     end if
         139             :
         140    17907264 :     allocate(newobj%crefwlw(nlw),stat=ierr)
         141     5969088 :     if (ierr/=0) then
         142           0 :        nullify(newobj)
         143           0 :        return
         144             :     end if
         145   101474496 :     newobj%crefwlw(:) = crefwlw(:)
         146             :
         147    17907264 :     allocate(newobj%crefwsw(nsw),stat=ierr)
         148     5969088 :     if (ierr/=0) then
         149           0 :        nullify(newobj)
         150           0 :        return
         151             :     end if
         152    89536320 :     newobj%crefwsw(:) = crefwsw(:)
         153             :
         154     5969088 :     call aero_state%water_uptake(aero_props, ilist, ibin,  ncol, nlev, dgnumwet, qaerwat)
         155             :
         156     5969088 :     nspec = aero_props%nspecies(ilist,ibin)
         157             :
         158     5969088 :     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     5969088 :                                newobj%radsurf, newobj%logradsurf, newobj%cheb)
         163             :
         164   561094272 :     do ilev = 1, nlev
         165  9269299584 :        dryvol(:ncol) = 0._r8
         166  2642813712 :        do ispec = 1, nspec
         167  2081719440 :           call aero_state%get_ambient_mmr(ilist,ispec,ibin,specmmr)
         168  2081719440 :           call aero_props%get(ibin, ispec, list_ndx=ilist, density=specdens)
         169             :
         170 35314998624 :           do icol = 1, ncol
         171 32678154000 :              vol(icol) = specmmr(icol,ilev)/specdens
         172 32678154000 :              dryvol(icol) = dryvol(icol) + vol(icol)
         173             :
         174 32678154000 :              newobj%watervol(icol,ilev) = qaerwat(icol,ilev)*rh2odens
         175 32678154000 :              newobj%wetvol(icol,ilev) = newobj%watervol(icol,ilev) + dryvol(icol)
         176 34759873440 :              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     5969088 :     newobj%aero_state => aero_state
         185     5969088 :     newobj%aero_props => aero_props
         186     5969088 :     newobj%ilist = ilist
         187     5969088 :     newobj%ibin = ibin
         188             :
         189     5969088 :   end function constructor
         190             :
         191             :   !------------------------------------------------------------------------------
         192             :   ! returns short wave aerosol optics properties
         193             :   !------------------------------------------------------------------------------
         194  3885876288 :   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  7771752576 :     real(r8) :: refr(ncol)      ! real part of refractive index
         206  7771752576 :     real(r8) :: refi(ncol)      ! imaginary part of refractive index
         207  7771752576 :     real(r8) :: cext(self%ncoef,ncol), cabs(self%ncoef,ncol), casm(self%ncoef,ncol)
         208             :
         209  7771752576 :     complex(r8) :: crefin(ncol) ! complex refractive index
         210             :     integer :: icol,icoef
         211             :
         212  7771752576 :     type(table_interp_wghts) :: wghtsr(ncol)
         213  7771752576 :     type(table_interp_wghts) :: wghtsi(ncol)
         214             :
         215 64885097088 :     crefin(:ncol) = self%aero_state%refractive_index_sw(ncol, ilev, self%ilist, self%ibin, iwav, self%aero_props)
         216             :
         217 64885097088 :     do icol = 1, ncol
         218 60999220800 :        crefin(icol) = crefin(icol) + self%watervol(icol,ilev)*self%crefwsw(iwav)
         219 60999220800 :        crefin(icol) = crefin(icol)/max(self%wetvol(icol,ilev),1.e-60_r8)
         220 60999220800 :        refr(icol) = real(crefin(icol))
         221 64885097088 :        refi(icol) = abs(aimag(crefin(icol)))
         222             :     end do
         223             :
         224             :     ! interpolate coefficients linear in refractive index
         225             :
         226  3885876288 :     wghtsr = table_interp_calcwghts( self%prefr, self%refrtabsw(:,iwav), ncol, refr(:ncol) )
         227  3885876288 :     wghtsi = table_interp_calcwghts( self%prefi, self%refitabsw(:,iwav), ncol, refi(:ncol) )
         228             :
         229  3885876288 :     cext(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%extpsw(:,:,:,iwav))
         230  3885876288 :     cabs(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%abspsw(:,:,:,iwav))
         231  3885876288 :     casm(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%asmpsw(:,:,:,iwav))
         232             :
         233 64885097088 :     do icol = 1,ncol
         234             :
         235 60999220800 :        if (self%logradsurf(icol,ilev) <= xrmax) then
         236 60999220800 :           pext(icol) = 0.5_r8*cext(1,icol)
         237 >30499*10^7 :           do icoef = 2, self%ncoef
         238 >30499*10^7 :              pext(icol) = pext(icol) + self%cheb(icoef,icol,ilev)*cext(icoef,icol)
         239             :           enddo
         240 60999220800 :           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 60999220800 :        pext(icol) = pext(icol)*self%wetvol(icol,ilev)*rhoh2o
         247 60999220800 :        pabs(icol) = 0.5_r8*cabs(1,icol)
         248 60999220800 :        pasm(icol) = 0.5_r8*casm(1,icol)
         249 >30499*10^7 :        do icoef = 2, self%ncoef
         250 >24399*10^7 :           pabs(icol) = pabs(icol) + self%cheb(icoef,icol,ilev)*cabs(icoef,icol)
         251 >30499*10^7 :           pasm(icol) = pasm(icol) + self%cheb(icoef,icol,ilev)*casm(icoef,icol)
         252             :        enddo
         253 60999220800 :        pabs(icol) = pabs(icol)*self%wetvol(icol,ilev)*rhoh2o
         254 60999220800 :        pabs(icol) = max(0._r8,pabs(icol))
         255 60999220800 :        pabs(icol) = min(pext(icol),pabs(icol))
         256             :
         257 64885097088 :        palb(icol) = 1._r8-pabs(icol)/max(pext(icol),1.e-40_r8)
         258             :
         259             :     end do
         260             :
         261  3885876288 :   end subroutine sw_props
         262             :
         263             :   !------------------------------------------------------------------------------
         264             :   ! returns long wave aerosol optics properties
         265             :   !------------------------------------------------------------------------------
         266  4441001472 :   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  8882002944 :     real(r8) :: refr(ncol)      ! real part of refractive index
         275  8882002944 :     real(r8) :: refi(ncol)      ! imaginary part of refractive index
         276  8882002944 :     real(r8) :: cabs(self%ncoef,ncol)
         277             :
         278  8882002944 :     complex(r8) :: crefin(ncol) ! complex refractive index
         279             :     integer :: icol, icoef
         280             :
         281  8882002944 :     type(table_interp_wghts) :: wghtsr(ncol)
         282  8882002944 :     type(table_interp_wghts) :: wghtsi(ncol)
         283             :
         284 74154396672 :     crefin(:ncol) = self%aero_state%refractive_index_lw(ncol, ilev, self%ilist, self%ibin, iwav, self%aero_props)
         285             :
         286 74154396672 :     do icol = 1, ncol
         287 69713395200 :        crefin(icol) = crefin(icol) + self%watervol(icol,ilev)*self%crefwlw(iwav)
         288 69713395200 :        crefin(icol) = crefin(icol)/max(self%wetvol(icol,ilev), 1.e-40_r8)
         289             :
         290 69713395200 :        refr(icol) = real(crefin(icol))
         291 74154396672 :        refi(icol) = aimag(crefin(icol))
         292             :
         293             :     end do
         294             :
         295             :     ! interpolate coefficients linear in refractive index
         296             :
         297  4441001472 :     wghtsr = table_interp_calcwghts( self%prefr, self%refrtablw(:,iwav), ncol, refr(:ncol) )
         298  4441001472 :     wghtsi = table_interp_calcwghts( self%prefi, self%refitablw(:,iwav), ncol, refi(:ncol) )
         299             :
         300  4441001472 :     cabs(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%absplw(:,:,:,iwav))
         301             :
         302 74154396672 :     do icol = 1,ncol
         303 69713395200 :        pabs(icol) = 0.5_r8*cabs(1,icol)
         304 >34856*10^7 :        do icoef = 2, self%ncoef
         305 >34856*10^7 :           pabs(icol) = pabs(icol) + self%cheb(icoef,icol,ilev)*cabs(icoef,icol)
         306             :        end do
         307 69713395200 :        pabs(icol) = pabs(icol)*self%wetvol(icol,ilev)*rhoh2o
         308 74154396672 :        pabs(icol) = max(0._r8,pabs(icol))
         309             :     end do
         310             :
         311  4441001472 :   end subroutine lw_props
         312             :
         313             :
         314             :   !------------------------------------------------------------------------------
         315             :   !------------------------------------------------------------------------------
         316     5969088 :   subroutine destructor(self)
         317             :
         318             :     type(refractive_aerosol_optics), intent(inout) :: self
         319             :
         320     5969088 :     deallocate(self%watervol)
         321     5969088 :     deallocate(self%wetvol)
         322     5969088 :     deallocate(self%cheb)
         323     5969088 :     deallocate(self%radsurf)
         324     5969088 :     deallocate(self%logradsurf)
         325     5969088 :     deallocate(self%crefwsw)
         326     5969088 :     deallocate(self%crefwlw)
         327             :
         328     5969088 :     nullify(self%aero_state)
         329     5969088 :     nullify(self%aero_props)
         330     5969088 :     nullify(self%extpsw)
         331     5969088 :     nullify(self%abspsw)
         332     5969088 :     nullify(self%asmpsw)
         333     5969088 :     nullify(self%absplw)
         334     5969088 :     nullify(self%refrtabsw)
         335     5969088 :     nullify(self%refitabsw)
         336     5969088 :     nullify(self%refrtablw)
         337     5969088 :     nullify(self%refitablw)
         338             :
         339     5969088 :   end subroutine destructor
         340             :
         341             :
         342             :   ! Private routines
         343             :   !===============================================================================
         344             :
         345             :   !===============================================================================
         346             :
         347     5969088 :   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    11938176 :     real(r8) :: xrad(ncol) ! normalized aerosol radius
         359             :
         360             :     !-------------------------------------------------------------------------------
         361             :
         362     5969088 :     explnsigma = exp(2.0_r8*alnsg_amode*alnsg_amode)
         363             :
         364   561094272 :     do k = 1, nlev
         365  9275268672 :        do i = 1, ncol
         366             :           ! convert from number mode diameter to surface area
         367  8714174400 :           radsurf(i,k) = max(0.5_r8*dgnumwet(i,k)*explnsigma,radmin)
         368  8714174400 :           logradsurf(i,k) = log(radsurf(i,k))
         369             :           ! normalize size parameter
         370  8714174400 :           xrad(i) = max(logradsurf(i,k),xrmin)
         371  8714174400 :           xrad(i) = min(xrad(i),xrmax)
         372  8714174400 :           xrad(i) = (2._r8*xrad(i)-xrmax-xrmin)/(xrmax-xrmin)
         373             :           ! chebyshev polynomials
         374  8714174400 :           cheb(1,i,k) = 1._r8
         375  8714174400 :           cheb(2,i,k) = xrad(i)
         376 35411822784 :           do nc = 3, ncoef
         377 34856697600 :              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     5969088 :   end subroutine modal_size_parameters
         383             :
         384    23876352 : end module refractive_aerosol_optics_mod

Generated by: LCOV version 1.14