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