#include <misc.h>
#include <params.h>


module aer_optics 4,2

  ! Purpose:
  !   initialize aerosol spectral optical properties
  !   for use in radiation subroutines

  ! Author: D. Fillmore

  ! dust properties are from C. Zender's Mie calculations
  ! all other properties are from OPAC project
  ! Optical Properties of Aerosols and Clouds
  ! Hess M., P. Koepke and I. Schult
  ! Bull. Amer. Met. Soc. 79, 831 - 844, 1998

  use shr_kind_mod, only: r8 => shr_kind_r8
  use filenames, only: aeroptics

  implicit none

  integer, parameter, public :: idxVIS = 8     ! index to visible band
  integer, parameter, public :: nrh = 1000   ! number of relative humidity values for look-up-table
  integer, parameter, public :: nspint = 19   ! number of spectral intervals
  integer, parameter, public :: ndstsz = 4    ! number of dust size bins

  real(r8), public :: ksul(nrh, nspint)    ! sulfate specific extinction  ( m^2 g-1 )
  real(r8), public :: wsul(nrh, nspint)    ! sulfate single scattering albedo
  real(r8), public :: gsul(nrh, nspint)    ! sulfate asymmetry parameter
  real(r8), public :: kbg(nspint)          ! background specific extinction  ( m^2 g-1 )
  real(r8), public :: wbg(nspint)          ! background single scattering albedo
  real(r8), public :: gbg(nspint)          ! background asymmetry parameter
  real(r8), public :: ksslt(nrh, nspint)   ! sea-salt specific extinction  ( m^2 g-1 )
  real(r8), public :: wsslt(nrh, nspint)   ! sea-salt single scattering albedo
  real(r8), public :: gsslt(nrh, nspint)   ! sea-salt asymmetry parameter
  real(r8), public :: kcphil(nrh, nspint)  ! hydrophilic carbon specific extinction  ( m^2 g-1 )
  real(r8), public :: wcphil(nrh, nspint)  ! hydrophilic carbon single scattering albedo
  real(r8), public :: gcphil(nrh, nspint)  ! hydrophilic carbon asymmetry parameter
  real(r8), public :: kcphob(nspint)       ! hydrophobic carbon specific extinction  ( m^2 g-1 )
  real(r8), public :: wcphob(nspint)       ! hydrophobic carbon single scattering albedo
  real(r8), public :: gcphob(nspint)       ! hydrophobic carbon asymmetry parameter
  real(r8), public :: kcb(nspint)          ! black carbon specific extinction  ( m^2 g-1 )
  real(r8), public :: wcb(nspint)          ! black carbon single scattering albedo
  real(r8), public :: gcb(nspint)          ! black carbon asymmetry parameter
  real(r8), public :: kvolc(nspint)        ! volcanic specific extinction  ( m^2 g-1)
  real(r8), public :: wvolc(nspint)        ! volcanic single scattering albedo
  real(r8), public :: gvolc(nspint)        ! volcanic asymmetry parameter
  real(r8), public :: kdst(ndstsz, nspint) ! dust specific extinction  ( m^2 g-1 )
  real(r8), public :: wdst(ndstsz, nspint) ! dust single scattering albedo
  real(r8), public :: gdst(ndstsz, nspint) ! dust asymmetry parameter

  public :: aer_optics_initialize
  public :: aer_optics_log
  private :: exp_interpol

  save

contains


  subroutine aer_optics_initialize 1,96

    use shr_kind_mod, only: r8 => shr_kind_r8
    use pmgrid  ! masterproc is here
    use ioFileMod, only: getfil

#if ( defined SPMD )
    use mpishorthand
#endif
    implicit none

    include 'netcdf.inc'


    integer :: nrh_opac  ! number of relative humidity values for OPAC data
    integer :: nbnd      ! number of spectral bands, should be identical to nspint
    real(r8), parameter :: wgt_sscm = 6.0 / 7.0

    real(r8), allocatable :: rh_opac(:)        ! OPAC relative humidity grid
    real(r8), allocatable :: ksul_opac(:,:)    ! sulfate  extinction
    real(r8), allocatable :: wsul_opac(:,:)    !          single scattering albedo
    real(r8), allocatable :: gsul_opac(:,:)    !          asymmetry parameter
    real(r8), allocatable :: ksslt_opac(:,:)   ! sea-salt
    real(r8), allocatable :: wsslt_opac(:,:)
    real(r8), allocatable :: gsslt_opac(:,:)
    real(r8), allocatable :: kssam_opac(:,:)   ! sea-salt accumulation mode
    real(r8), allocatable :: wssam_opac(:,:)
    real(r8), allocatable :: gssam_opac(:,:)
    real(r8), allocatable :: ksscm_opac(:,:)   ! sea-salt coarse mode
    real(r8), allocatable :: wsscm_opac(:,:)
    real(r8), allocatable :: gsscm_opac(:,:)
    real(r8), allocatable :: kcphil_opac(:,:)  ! hydrophilic organic carbon
    real(r8), allocatable :: wcphil_opac(:,:)
    real(r8), allocatable :: gcphil_opac(:,:)

    character(len=256) :: locfn            ! local file

    integer :: krh_opac  ! rh index for OPAC rh grid
    integer :: krh       ! another rh index
    integer :: ksz       ! dust size bin index
    integer :: kbnd      ! band index

    real(r8) :: rh   ! local relative humidity variable

    integer :: nc_id
    integer :: dst_id, rh_id, bnd_id
    integer :: rh_opac_id
    integer :: ksul_opac_id, wsul_opac_id, gsul_opac_id
    integer :: kssam_opac_id, wssam_opac_id, gssam_opac_id
    integer :: ksscm_opac_id, wsscm_opac_id, gsscm_opac_id
    integer :: kcphil_opac_id, wcphil_opac_id, gcphil_opac_id
    integer :: kcb_id, wcb_id, gcb_id
    integer :: kvolc_id, wvolc_id, gvolc_id
    integer :: kdst_id, wdst_id, gdst_id
    integer :: kbg_id, wbg_id, gbg_id

    if ( masterproc ) then

    write (6, '(2x, a)') '_______________________________________________________'
    write (6, '(2x, a)') '_______ initializing aerosol optical properties _______'
    write (6, '(2x, a)') '_______________________________________________________'

    call getfil(aeroptics, locfn, 0)

    call wrap_open(locfn, 0, nc_id)
    call wrap_inq_dimid(nc_id, 'dst', dst_id)
    call wrap_inq_dimid(nc_id, 'rh', rh_id)
    call wrap_inq_dimid(nc_id, 'wvl_cam', bnd_id)

    call wrap_inq_dimlen(nc_id, rh_id, nrh_opac)
    call wrap_inq_dimlen(nc_id, bnd_id, nbnd)

    write (6, '(2x, a, x, i2)') 'aer_optics_initialize : nrh_opac =', nrh_opac
    write (6, '(2x, a, x, i2)') 'aer_optics_initialize : nbnd =', nbnd

    ! allocate local arrays
    allocate(rh_opac(nrh_opac))
    allocate(ksul_opac(nrh_opac, nbnd))
    allocate(wsul_opac(nrh_opac, nbnd))
    allocate(gsul_opac(nrh_opac, nbnd))
    allocate(ksslt_opac(nrh_opac, nbnd))
    allocate(wsslt_opac(nrh_opac, nbnd))
    allocate(gsslt_opac(nrh_opac, nbnd))
    allocate(kssam_opac(nrh_opac, nbnd))
    allocate(wssam_opac(nrh_opac, nbnd))
    allocate(gssam_opac(nrh_opac, nbnd))
    allocate(ksscm_opac(nrh_opac, nbnd))
    allocate(wsscm_opac(nrh_opac, nbnd))
    allocate(gsscm_opac(nrh_opac, nbnd))
    allocate(kcphil_opac(nrh_opac, nbnd))
    allocate(wcphil_opac(nrh_opac, nbnd))
    allocate(gcphil_opac(nrh_opac, nbnd))

    call wrap_inq_varid(nc_id, 'rh', rh_opac_id)
    call wrap_inq_varid(nc_id, 'ext_cam_sul', ksul_opac_id)
    call wrap_inq_varid(nc_id, 'ssa_cam_sul', wsul_opac_id)
    call wrap_inq_varid(nc_id, 'asm_cam_sul', gsul_opac_id)
    call wrap_inq_varid(nc_id, 'ext_cam_ssam', kssam_opac_id)
    call wrap_inq_varid(nc_id, 'ssa_cam_ssam', wssam_opac_id)
    call wrap_inq_varid(nc_id, 'asm_cam_ssam', gssam_opac_id)
    call wrap_inq_varid(nc_id, 'ext_cam_sscm', ksscm_opac_id)
    call wrap_inq_varid(nc_id, 'ssa_cam_sscm', wsscm_opac_id)
    call wrap_inq_varid(nc_id, 'asm_cam_sscm', gsscm_opac_id)
    call wrap_inq_varid(nc_id, 'ext_cam_waso', kcphil_opac_id)
    call wrap_inq_varid(nc_id, 'ssa_cam_waso', wcphil_opac_id)
    call wrap_inq_varid(nc_id, 'asm_cam_waso', gcphil_opac_id)
    call wrap_inq_varid(nc_id, 'ext_cam_soot', kcb_id)
    call wrap_inq_varid(nc_id, 'ssa_cam_soot', wcb_id)
    call wrap_inq_varid(nc_id, 'asm_cam_soot', gcb_id)
    call wrap_inq_varid(nc_id, 'ext_cam_volc', kvolc_id)
    call wrap_inq_varid(nc_id, 'ssa_cam_volc', wvolc_id)
    call wrap_inq_varid(nc_id, 'asm_cam_volc', gvolc_id)
    call wrap_inq_varid(nc_id, 'ext_cam_dust', kdst_id)
    call wrap_inq_varid(nc_id, 'ssa_cam_dust', wdst_id)
    call wrap_inq_varid(nc_id, 'asm_cam_dust', gdst_id)
    call wrap_inq_varid(nc_id, 'ext_cam_bg', kbg_id)
    call wrap_inq_varid(nc_id, 'ssa_cam_bg', wbg_id)
    call wrap_inq_varid(nc_id, 'asm_cam_bg', gbg_id)

    call wrap_get_var_realx(nc_id, rh_opac_id, rh_opac)
    call wrap_get_var_realx(nc_id, ksul_opac_id, ksul_opac)
    call wrap_get_var_realx(nc_id, wsul_opac_id, wsul_opac)
    call wrap_get_var_realx(nc_id, gsul_opac_id, gsul_opac)
    call wrap_get_var_realx(nc_id, kssam_opac_id, kssam_opac)
    call wrap_get_var_realx(nc_id, wssam_opac_id, wssam_opac)
    call wrap_get_var_realx(nc_id, gssam_opac_id, gssam_opac)
    call wrap_get_var_realx(nc_id, ksscm_opac_id, ksscm_opac)
    call wrap_get_var_realx(nc_id, wsscm_opac_id, wsscm_opac)
    call wrap_get_var_realx(nc_id, gsscm_opac_id, gsscm_opac)
    call wrap_get_var_realx(nc_id, kcphil_opac_id, kcphil_opac)
    call wrap_get_var_realx(nc_id, wcphil_opac_id, wcphil_opac)
    call wrap_get_var_realx(nc_id, gcphil_opac_id, gcphil_opac)
    call wrap_get_var_realx(nc_id, kcb_id, kcb)
    call wrap_get_var_realx(nc_id, wcb_id, wcb)
    call wrap_get_var_realx(nc_id, gcb_id, gcb)
    call wrap_get_var_realx(nc_id, kvolc_id, kvolc)
    call wrap_get_var_realx(nc_id, wvolc_id, wvolc)
    call wrap_get_var_realx(nc_id, gvolc_id, gvolc)
    call wrap_get_var_realx(nc_id, kdst_id, kdst)
    call wrap_get_var_realx(nc_id, wdst_id, wdst)
    call wrap_get_var_realx(nc_id, gdst_id, gdst)
    call wrap_get_var_realx(nc_id, kbg_id, kbg)
    call wrap_get_var_realx(nc_id, wbg_id, wbg)
    call wrap_get_var_realx(nc_id, gbg_id, gbg)

    ! map OPAC aerosol species onto CAM aerosol species
    ! CAM name             OPAC name
    ! sul   or SO4         = suso                  sulfate soluble
    ! sslt  or SSLT        = 1/7 ssam + 6/7 sscm   sea-salt accumulation/coagulation mode
    ! cphil or CPHI        = waso                  water soluble (carbon)
    ! cphob or CPHO        = waso @ rh = 0
    ! cb    or BCPHI/BCPHO = soot

    ksslt_opac(:,:) = (1.0 - wgt_sscm) * kssam_opac(:,:) + wgt_sscm * ksscm_opac(:,:)

    wsslt_opac(:,:) = ( (1.0 - wgt_sscm) * kssam_opac(:,:) * wssam_opac(:,:) &
                  + wgt_sscm * ksscm_opac(:,:) * wsscm_opac(:,:) ) &
                  / ksslt_opac(:,:)

    gsslt_opac(:,:) = ( (1.0 - wgt_sscm) * kssam_opac(:,:) * wssam_opac(:,:) * gssam_opac(:,:) &
                  + wgt_sscm * ksscm_opac(:,:) * wsscm_opac(:,:) * gsscm_opac(:,:) ) &
                   / ( ksslt_opac(:,:) * wsslt_opac(:,:) )

    kcphob(:) = kcphil_opac(1,:)
    wcphob(:) = wcphil_opac(1,:)
    gcphob(:) = gcphil_opac(1,:)

    ! interpolate optical properties of hygrospopic aerosol species
    !   onto a uniform relative humidity grid

    do krh = 1, nrh
      rh = 1.0_r8 / nrh * (krh - 1)
      do kbnd = 1, nbnd
        ksul(krh, kbnd) = exp_interpol( rh_opac, &
          ksul_opac(:, kbnd) / ksul_opac(1, kbnd), rh ) * ksul_opac(1, kbnd)
        wsul(krh, kbnd) = lin_interpol( rh_opac, &
          wsul_opac(:, kbnd) / wsul_opac(1, kbnd), rh ) * wsul_opac(1, kbnd)
        gsul(krh, kbnd) = lin_interpol( rh_opac, &
          gsul_opac(:, kbnd) / gsul_opac(1, kbnd), rh ) * gsul_opac(1, kbnd)
        ksslt(krh, kbnd) = exp_interpol( rh_opac, &
          ksslt_opac(:, kbnd) / ksslt_opac(1, kbnd), rh ) * ksslt_opac(1, kbnd)
        wsslt(krh, kbnd) = lin_interpol( rh_opac, &
          wsslt_opac(:, kbnd) / wsslt_opac(1, kbnd), rh ) * wsslt_opac(1, kbnd)
        gsslt(krh, kbnd) = lin_interpol( rh_opac, &
          gsslt_opac(:, kbnd) / gsslt_opac(1, kbnd), rh ) * gsslt_opac(1, kbnd)
        kcphil(krh, kbnd) = exp_interpol( rh_opac, &
          kcphil_opac(:, kbnd) / kcphil_opac(1, kbnd), rh ) * kcphil_opac(1, kbnd)
        wcphil(krh, kbnd) = lin_interpol( rh_opac, &
          wcphil_opac(:, kbnd) / wcphil_opac(1, kbnd), rh ) * wcphil_opac(1, kbnd)
        gcphil(krh, kbnd) = lin_interpol( rh_opac, &
          gcphil_opac(:, kbnd) / gcphil_opac(1, kbnd), rh )  * gcphil_opac(1, kbnd)
      end do
    end do

    deallocate(rh_opac)
    deallocate(ksul_opac)
    deallocate(wsul_opac)
    deallocate(gsul_opac)
    deallocate(ksslt_opac)
    deallocate(wsslt_opac)
    deallocate(gsslt_opac)
    deallocate(kssam_opac)
    deallocate(wssam_opac)
    deallocate(gssam_opac)
    deallocate(ksscm_opac)
    deallocate(wsscm_opac)
    deallocate(gsscm_opac)
    deallocate(kcphil_opac)
    deallocate(wcphil_opac)
    deallocate(gcphil_opac)

    ! write optical constants to log file for debugging

    write (6, '(2x, a)') '_______ hygroscopic growth in visible band _______'
    call aer_optics_log_rh( 'SO4', ksul(:, 8), wsul(:, 8), gsul(:, 8) )
    call aer_optics_log_rh( 'SSLT', ksslt(:, 8), wsslt(:, 8), gsslt(:, 8) )
    call aer_optics_log_rh( 'OCPHI', kcphil(:, 8), wcphil(:, 8), gcphil(:, 8) )

    end if ! end if ( masterproc )

    ! broadcast aerosol optical constants to all nodes
#if ( defined SPMD )
    call mpibcast(ksul, nrh * nspint, mpir8, 0, mpicom)
    call mpibcast(wsul, nrh * nspint, mpir8, 0, mpicom)
    call mpibcast(gsul, nrh * nspint, mpir8, 0, mpicom)
    call mpibcast(ksslt, nrh * nspint, mpir8, 0, mpicom)
    call mpibcast(wsslt, nrh * nspint, mpir8, 0, mpicom)
    call mpibcast(gsslt, nrh * nspint, mpir8, 0, mpicom)
    call mpibcast(kcphil, nrh * nspint, mpir8, 0, mpicom)
    call mpibcast(wcphil, nrh * nspint, mpir8, 0, mpicom)
    call mpibcast(gcphil, nrh * nspint, mpir8, 0, mpicom)
    call mpibcast(kcphob, nspint, mpir8, 0, mpicom)
    call mpibcast(wcphob, nspint, mpir8, 0, mpicom)
    call mpibcast(gcphob, nspint, mpir8, 0, mpicom)
    call mpibcast(kcb, nspint, mpir8, 0, mpicom)
    call mpibcast(wcb, nspint, mpir8, 0, mpicom)
    call mpibcast(gcb, nspint, mpir8, 0, mpicom)
    call mpibcast(kvolc, nspint, mpir8, 0, mpicom)
    call mpibcast(wvolc, nspint, mpir8, 0, mpicom)
    call mpibcast(gvolc, nspint, mpir8, 0, mpicom)
    call mpibcast(kdst, ndstsz * nspint, mpir8, 0, mpicom)
    call mpibcast(wdst, ndstsz * nspint, mpir8, 0, mpicom)
    call mpibcast(gdst, ndstsz * nspint, mpir8, 0, mpicom)
    call mpibcast(kbg, nspint, mpir8, 0, mpicom)
    call mpibcast(wbg, nspint, mpir8, 0, mpicom)
    call mpibcast(gbg, nspint, mpir8, 0, mpicom)
#endif

  end subroutine aer_optics_initialize


  subroutine aer_optics_log(name, ext, ssa, asm),1

    ! Purpose:
    !   write aerosol optical constants to log file

    ! Author: D. Fillmore

    use shr_kind_mod, only: r8 => shr_kind_r8

    implicit none

    character(len=*), intent(in) :: name
    real(r8), intent(in) :: ext(nspint)
    real(r8), intent(in) :: ssa(nspint)
    real(r8), intent(in) :: asm(nspint)

    integer :: kbnd

    write (6, '(2x, a)') name
    write (6, '(2x, a, 4x, a, 4x, a, 4x, a)') 'SW band', 'ext (m^2 g-1)', ' ssa', ' asm'
    do kbnd = 1, nspint
      write (6, '(2x, i7, 4x, f13.2, 4x, f4.2, 4x, f4.2)') kbnd, ext(kbnd), ssa(kbnd), asm(kbnd)
    end do

  end subroutine aer_optics_log



  subroutine aer_optics_log_rh(name, ext, ssa, asm) 3,1

    ! Purpose:
    !   write out aerosol optical properties
    !   for a set of test rh values
    !   to test hygroscopic growth interpolation

    ! Author: D. Fillmore

    use shr_kind_mod, only: r8 => shr_kind_r8

    implicit none

    character(len=*), intent(in) :: name
    real(r8), intent(in) :: ext(nrh)
    real(r8), intent(in) :: ssa(nrh)
    real(r8), intent(in) :: asm(nrh)

    integer :: krh_test
    integer, parameter :: nrh_test = 36
    integer :: krh
    real(r8) :: rh
    real(r8) :: rh_test(nrh_test)
    real(r8) :: exti
    real(r8) :: ssai
    real(r8) :: asmi
    real(r8) :: wrh

    do krh_test = 1, nrh_test
      rh_test(krh_test) = sqrt(sqrt(sqrt(sqrt(((krh_test - 1.0) / (nrh_test - 1))))))
    enddo
    write (6, '(2x, a)') name
    write (6, '(2x, a, 4x, a, 4x, a, 4x, a)') '   rh', 'ext (m^2 g-1)', '  ssa', '  asm'

    ! loop through test rh values
    do krh_test = 1, nrh_test
      ! find corresponding rh index
      rh = rh_test(krh_test)
      krh = min(floor( (rh) * nrh ) + 1, nrh - 1)
      wrh = (rh) *nrh - krh
      exti = ext(krh + 1) * (wrh + 1) - ext(krh) * wrh
      ssai = ssa(krh + 1) * (wrh + 1) - ssa(krh) * wrh
      asmi = asm(krh + 1) * (wrh + 1) - asm(krh) * wrh
      write (6, '(2x, f5.3, 4x, f13.3, 4x, f5.3, 4x, f5.3)') rh_test(krh_test), exti, ssai, asmi
      ! write (6, '(2x, f5.3, 4x, f13.3, 4x, f5.3, 4x, f5.3)') rh_test(krh_test), ext(krh), ssa(krh), asm(krh)
    end do

  end subroutine aer_optics_log_rh



  function exp_interpol(x, f, y) result(g) 3,1

    ! Purpose:
    !   interpolates f(x) to point y
    !   assuming f(x) = f(x0) exp a(x - x0)
    !   where a = ( ln f(x1) - ln f(x0) ) / (x1 - x0)
    !   x0 <= x <= x1
    !   assumes x is monotonically increasing

    ! Author: D. Fillmore

    use shr_kind_mod, only: r8 => shr_kind_r8

    implicit none

    real(r8), intent(in), dimension(:) :: x  ! grid points
    real(r8), intent(in), dimension(:) :: f  ! grid function values
    real(r8), intent(in) :: y                ! interpolation point
    real(r8) :: g                            ! interpolated function value

    integer :: k  ! interpolation point index
    integer :: n  ! length of x
    real(r8) :: a

    n = size(x)

    ! find k such that x(k) < y =< x(k+1)
    ! set k = 1 if y <= x(1)  and  k = n-1 if y > x(n)

    if (y <= x(1)) then
      k = 1
    else if (y >= x(n)) then
      k = n - 1
    else
      k = 1
      do while (y > x(k+1) .and. k < n)
        k = k + 1
      end do
    end if

    ! interpolate
    a = (  log( f(k+1) / f(k) )  ) / ( x(k+1) - x(k) )
    g = f(k) * exp( a * (y - x(k)) )

  end function exp_interpol


  function lin_interpol(x, f, y) result(g) 6,1

    ! Purpose:
    !   interpolates f(x) to point y
    !   assuming f(x) = f(x0) + a * (x - x0)
    !   where a = ( f(x1) - f(x0) ) / (x1 - x0)
    !   x0 <= x <= x1
    !   assumes x is monotonically increasing

    ! Author: D. Fillmore

    use shr_kind_mod, only: r8 => shr_kind_r8

    implicit none

    real(r8), intent(in), dimension(:) :: x  ! grid points
    real(r8), intent(in), dimension(:) :: f  ! grid function values
    real(r8), intent(in) :: y                ! interpolation point
    real(r8) :: g                            ! interpolated function value

    integer :: k  ! interpolation point index
    integer :: n  ! length of x
    real(r8) :: a

    n = size(x)

    ! find k such that x(k) < y =< x(k+1)
    ! set k = 1 if y <= x(1)  and  k = n-1 if y > x(n)

    if (y <= x(1)) then
      k = 1
    else if (y >= x(n)) then
      k = n - 1
    else
      k = 1
      do while (y > x(k+1) .and. k < n)
        k = k + 1
      end do
    end if

    ! interpolate
    a = (  f(k+1) - f(k) ) / ( x(k+1) - x(k) )
    g = f(k) + a * (y - x(k))

  end function lin_interpol

end module aer_optics