#include <misc.h>
#include <preproc.h>


module VOCEmissionMod 1,2

!-----------------------------------------------------------------------
!BOP
!
! !MODULE: VOCEmissionMod
!
! !DESCRIPTION:
! Volatile organic compound emission
!
! !USES:
  use clm_varctl, only: iulog
  use abortutils, only: endrun
!
! !PUBLIC TYPES:
  implicit none
  save
!
! !PUBLIC MEMBER FUNCTIONS:
  public :: VOCEmission
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!EOP
!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: VOCEmission
!
! !INTERFACE:

  subroutine VOCEmission (lbp, ubp, num_soilp, filter_soilp ) 1,21
!
! ! NEW DESCRIPTION
! Volatile organic compound emission
! This code simulates volatile organic compound emissions following:
! 1. Isoprene: Guenther et al., 2006 description of MEGAN emissions
!     following equations 2-9, 16-17, 20
! 2. Monoterpenes/OVOCs/ORVOCs/CO: algorithm presented in Guenther, A., 
!    1999: Modeling Biogenic Volatile Organic Compound Emissions to the 
!    Atmosphere. In Reactive Hydrocarbons in the Atmosphere, Ch. 3
!    With updates from MEGAN online user's guide 
!    ( http://acd.ucar.edu/~guenther/MEGAN/MEGANusersguide.pdf)
! This model relies on the assumption that 90% of isoprene and monoterpene
! emissions originate from canopy foliage:
!    E= epsilon * gamma * rho
! VOC flux (E) [ugC m-2 h-1] is calculated from baseline emission
! factors (epsilon) [ugC m-2 h-1] which are mapped for each PFT (isoprene)
! or constant for each PFT (others).  Note that for constant EFs the units
! of [ugC g-1 h-1] must be multiplied by the source density factor.
! The emission activity factor (gamma) [unitless] for isoprene includes 
! dependence on PPFT, temperature, LAI, leaf age and soil moisture.  
! The canopy environment constant was calculated offline for CLM+CAM at 
! standard conditions.
! The emission activity factor for the other emissions depends on temperature.
! We assume that the escape efficiency (rho) here is unity following
! Guenther et al., 2006.
! Subroutine written to operate at the patch level.
! IN FINAL IMPLEMENTATION, REMEMBER:
! 1. may wish to call this routine only as freq. as rad. calculations
! 2. may wish to place epsilon values directly in pft-physiology file
! Output: vocflx(nvoc) !VOC flux [ug C m-2 h-1]
!
!
! !USES:
    use shr_kind_mod , only : r8 => shr_kind_r8
    use clm_atmlnd   , only : clm_a2l
    use clmtype
    use clm_varpar   , only : nvoc, numpft
    use clm_varpar   , only : nvoc, numpft
    use clm_atmlnd   , only : clm_a2l
    use shr_const_mod, only : SHR_CONST_RGAS
    use clm_varcon   , only : denice
    use clm_varpar   , only : nlevsoi
    use pftvarcon    , only : ndllf_evr_tmp_tree,  ndllf_evr_brl_tree,    &
                              ndllf_dcd_brl_tree,  nbrdlf_evr_trp_tree,   &
                              nbrdlf_evr_tmp_tree, nbrdlf_dcd_brl_shrub,  &
                              nbrdlf_dcd_trp_tree, nbrdlf_dcd_tmp_tree,   &
                              nbrdlf_dcd_brl_tree, nbrdlf_evr_shrub,      &
                              nc3_arctic_grass,    nwheat,                &
                              ncorn,               nc4_grass,             &
                              noveg
!
! !ARGUMENTS:
    implicit none
    integer, intent(in) :: lbp, ubp                    ! pft bounds
    integer, intent(in) :: num_soilp                   ! number of columns in soil pft filter
    integer, intent(in) :: filter_soilp(num_soilp)     ! pft filter for soil
!
! !CALLED FROM:
!
! !REVISION HISTORY:
! Author: Sam Levis
! 2/1/02, Peter Thornton: migration to new data structure
! 4/15/06, Colette L. Heald: modify for updated MEGAN model (Guenther et al., 2006)
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in arguments
!
    integer , pointer :: pgridcell(:)     ! gridcell index of corresponding pft
    integer , pointer :: pcolumn(:)       ! column index of corresponding pft
    integer , pointer :: ivt(:)           ! pft vegetation type for current
    real(r8), pointer :: t_veg(:)         ! pft vegetation temperature (Kelvin)
    real(r8), pointer :: fsun(:)          ! sunlit fraction of canopy
    real(r8), pointer :: elai(:)          ! one-sided leaf area index with burying by snow
    real(r8), pointer :: clayfrac(:)      ! fraction of soil that is clay
    real(r8), pointer :: sandfrac(:)      ! fraction of soil that is sand
    real(r8), pointer :: forc_solad(:,:)  ! direct beam radiation (visible only)
    real(r8), pointer :: forc_solai(:,:)  ! diffuse radiation     (visible only)
    real(r8), pointer :: sla(:)           ! specific leaf area [m2 leaf g-1 C]
    real(r8), pointer :: h2osoi_vol(:,:)  ! volumetric soil water (m3/m3)
    real(r8), pointer :: h2osoi_ice(:,:)  ! ice soil content (kg/m3)
    real(r8), pointer :: dz(:,:)          ! depth of layer (m)
    real(r8), pointer :: coszen(:)        ! cosine of solar zenith angle
    real(r8), pointer :: efisop(:,:)      ! emission factors for isoprene for each pft [ug C m-2 h-1]
    real(r8), pointer :: elai_p(:)        ! one-sided leaf area index from previous timestep
    real(r8), pointer :: t_veg24(:)       ! avg pft vegetation temperature for last 24 hrs
    real(r8), pointer :: t_veg240(:)      ! avg pft vegetation temperature for last 240 hrs
    real(r8), pointer :: fsun24(:)        ! sunlit fraction of canopy last 24 hrs
    real(r8), pointer :: fsun240(:)       ! sunlit fraction of canopy last 240 hrs
    real(r8), pointer :: forc_solad24(:)  ! direct beam radiation last 24hrs  (visible only)
    real(r8), pointer :: forc_solai24(:)  ! diffuse radiation  last 24hrs     (visible only)
    real(r8), pointer :: forc_solad240(:) ! direct beam radiation last 240hrs (visible only)
    real(r8), pointer :: forc_solai240(:) ! diffuse radiation  last 240hrs    (visible only)
    real(r8), pointer :: bsw(:,:)         ! Clapp and Hornberger "b" (nlevgrnd)
    real(r8), pointer :: watsat(:,:)      ! volumetric soil water at saturation (porosity) (nlevgrnd)
    real(r8), pointer :: sucsat(:,:)      ! minimum soil suction (mm) (nlevgrnd)

    real(r8), parameter :: smpmax = 2.57e5_r8 ! maximum soil matrix potential
!
! local pointers to original implicit out arrays
!
    real(r8), pointer :: vocflx(:,:)      ! VOC flux [ug C m-2 h-1]
    real(r8), pointer :: vocflx_tot(:)    ! VOC flux [ug C m-2 h-1]
    real(r8), pointer :: vocflx_1(:)      ! VOC flux(1) [ug C m-2 h-1]
    real(r8), pointer :: vocflx_2(:)      ! VOC flux(2) [ug C m-2 h-1]
    real(r8), pointer :: vocflx_3(:)      ! VOC flux(3) [ug C m-2 h-1]
    real(r8), pointer :: vocflx_4(:)      ! VOC flux(4) [ug C m-2 h-1]
    real(r8), pointer :: vocflx_5(:)      ! VOC flux(5) [ug C m-2 h-1]
    real(r8), pointer :: Eopt_out(:)     
    real(r8), pointer :: topt_out(:)
    real(r8), pointer :: alpha_out(:)
    real(r8), pointer :: cp_out(:)
    real(r8), pointer :: paru_out(:)
    real(r8), pointer :: par24u_out(:)
    real(r8), pointer :: par240u_out(:)
    real(r8), pointer :: para_out(:)
    real(r8), pointer :: par24a_out(:)
    real(r8), pointer :: par240a_out(:)
    real(r8), pointer :: gamma_out(:)
    real(r8), pointer :: gammaT_out(:)
    real(r8), pointer :: gammaP_out(:)
    real(r8), pointer :: gammaL_out(:)
    real(r8), pointer :: gammaA_out(:)
    real(r8), pointer :: gammaS_out(:)
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
    integer  :: fp,p,g,c,n,j            ! indices
    integer  :: ct_bad
    real(r8) :: epsilon(lbp:ubp)        ! emission factor [ugC m-2 h-1]
    real(r8) :: par                     ! temporary
    real(r8) :: par24                   ! temporary
    real(r8) :: par240                  ! temporary
    real(r8) :: density                 ! source density factor [g dry wgt foliar mass/m2 ground]
    real(r8) :: gamma(lbp:ubp)          ! activity factor (accounting for light, T, age, LAI conditions)
    real(r8) :: gamma_p                 ! activity factor for PPFD
    real(r8) :: gamma_l                 ! activity factor for PPFD & LAI
    real(r8) :: gamma_t                 ! activity factor for temperature
    real(r8) :: gamma_a                 ! activity factor for leaf age
    real(r8) :: gamma_sm                ! activity factor for soil moisture
    real(r8) :: x                       ! temporary 
    real(r8) :: Eopt                    ! temporary 
    real(r8) :: topt                    ! temporary 
    real(r8) :: cp                      ! temporary
    real(r8) :: alpha                   ! temporary
    real(r8) :: elai_prev               ! lai for previous timestep
    real(r8) :: fnew, fgro, fmat, fsen  ! fractions of leaves at different phenological stages
    real(r8) :: nl                      ! temporary number of soil levels
    real(r8) :: theta_ice               ! water content in ice in m3/m3
    real(r8) :: wilt                    ! wilting point in m3/m3
    real(r8) :: theta1                  ! temporary
!
! Constants
!
    real(r8), parameter :: R   = SHR_CONST_RGAS*0.001_r8   ! univ. gas constant [J K-1 mol-1]
    real(r8), parameter :: scale_mw =0.882_r8              ! conversion factor for isoprene -> carbon
    real(r8), parameter :: alpha_fix = 0.001_r8            ! empirical coefficient
    real(r8), parameter :: cp_fix = 1.21_r8                ! empirical coefficient
    real(r8), parameter :: ct1 = 95.0_r8                   ! empirical coefficient (70 in User's Guide)
    real(r8), parameter :: ct2 = 230.0_r8                  ! empirical coefficient  (200 in User's Guide)
    real(r8), parameter :: ct3 = 0.00831_r8                ! empirical coefficient (0.0083 in User's Guide)
    real(r8), parameter :: topt_fix = 317._r8              ! std temperature [K]
    real(r8), parameter :: Eopt_fix = 2.26_r8              ! empirical coefficient
    real(r8), parameter :: tstd = 303.15_r8                ! std temperature [K]
    real(r8), parameter :: bet = 0.09_r8                   ! beta empirical coefficient [K-1]
    real(r8), parameter :: clai1 = 0.49_r8                 ! empirical coefficient
    real(r8), parameter :: clai2 = 0.2_r8                  ! empirical coefficient
    real(r8), parameter :: clai3 = 5.0_r8                  ! empirical coefficient
    real(r8), parameter :: Anew = 0.01_r8                  ! relative emission factor for new plants
    real(r8), parameter :: Agro = 0.5_r8                   ! relative emission factor for new plants
    real(r8), parameter :: Amat = 1.0_r8                   ! relative emission factor for new plants
    real(r8), parameter :: Asen = 0.33_r8                  ! relative emission factor for new plants
    real(r8), parameter :: cce = 0.40_r8                   ! factor to set emissions to unity @ std
    real(r8), parameter :: cce1 = 0.47_r8                  ! same as Cce but for non-accumulated vars
    real(r8), parameter :: ca1 = 0.004_r8                  ! empirical coefficent for alpha
    real(r8), parameter :: ca2 = 0.0005_r8                 ! empirical coefficent for alpha
    real(r8), parameter :: ca3 = 0.0468_r8                 ! empirical coefficent for cp
    real(r8), parameter :: par0_sun = 200._r8              ! std conditions for past 24 hrs [umol/m2/s]
    real(r8), parameter :: par0_shade = 50._r8             ! std conditions for past 24 hrs [umol/m2/s]
    real(r8), parameter :: co1 = 313._r8                   ! empirical coefficient
    real(r8), parameter :: co2 = 0.6_r8                    ! empirical coefficient
    real(r8), parameter :: co3 = 2.034_r8                  ! empirical coefficient
    real(r8), parameter :: co4 = 0.05_r8                   ! empirical coefficient
    real(r8), parameter :: tstd0 = 297_r8                  ! std temperature [K]
    real(r8), parameter :: deltheta1=0.06_r8               ! empirical coefficient
    real(r8), parameter :: scaling_to_500_Tg = 5._r8/7._r8 ! J-F Larmaque's empirical scaling factor

!
! These are the values from version of genesis-ibis / 1000.
! CN calculates its own sla [m2 leaf g-1 C]
! Divide by 2 in the equation to get dry weight foliar mass from grams carbon
!
    real(r8) :: hardwire_sla(0:numpft)
    real(r8) :: slarea(lbp:ubp)           ! Specific leaf areas [m2 leaf g-1 C]
    real(r8) :: hardwire_droot(0:numpft)  ! Root depth [m]
!-----------------------------------------------------------------------

    ! Assign local pointers to derived type members (gridcell-level)
    forc_solad => clm_a2l%forc_solad
    forc_solai => clm_a2l%forc_solai
    efisop     => clm3%g%gve%efisop

    ! Assign local pointers to derived subtypes components (column-level)
    h2osoi_vol       => clm3%g%l%c%cws%h2osoi_vol
    h2osoi_ice       => clm3%g%l%c%cws%h2osoi_ice
    dz               => clm3%g%l%c%cps%dz
    bsw              => clm3%g%l%c%cps%bsw
    watsat           => clm3%g%l%c%cps%watsat
    sucsat           => clm3%g%l%c%cps%sucsat

    ! Assign local pointers to derived subtypes components (pft-level)
    pgridcell        => clm3%g%l%c%p%gridcell
    pcolumn          => clm3%g%l%c%p%column
    ivt              => clm3%g%l%c%p%itype
    t_veg            => clm3%g%l%c%p%pes%t_veg
    fsun             => clm3%g%l%c%p%pps%fsun
    elai             => clm3%g%l%c%p%pps%elai
    clayfrac         => clm3%g%l%c%p%pps%clayfrac
    sandfrac         => clm3%g%l%c%p%pps%sandfrac
    vocflx           => clm3%g%l%c%p%pvf%vocflx
    vocflx_tot       => clm3%g%l%c%p%pvf%vocflx_tot
    vocflx_1         => clm3%g%l%c%p%pvf%vocflx_1
    vocflx_2         => clm3%g%l%c%p%pvf%vocflx_2
    vocflx_3         => clm3%g%l%c%p%pvf%vocflx_3
    vocflx_4         => clm3%g%l%c%p%pvf%vocflx_4
    vocflx_5         => clm3%g%l%c%p%pvf%vocflx_5
    Eopt_out         => clm3%g%l%c%p%pvf%Eopt_out
    topt_out         => clm3%g%l%c%p%pvf%topt_out
    alpha_out        => clm3%g%l%c%p%pvf%alpha_out
    cp_out           => clm3%g%l%c%p%pvf%cp_out
    paru_out         => clm3%g%l%c%p%pvf%paru_out
    par24u_out       => clm3%g%l%c%p%pvf%par24u_out
    par240u_out      => clm3%g%l%c%p%pvf%par240u_out
    para_out         => clm3%g%l%c%p%pvf%para_out
    par24a_out       => clm3%g%l%c%p%pvf%par24a_out
    par240a_out      => clm3%g%l%c%p%pvf%par240a_out
    gammaL_out       => clm3%g%l%c%p%pvf%gammaL_out
    gammaT_out       => clm3%g%l%c%p%pvf%gammaT_out
    gammaP_out       => clm3%g%l%c%p%pvf%gammaP_out
    gammaA_out       => clm3%g%l%c%p%pvf%gammaA_out
    gammaS_out       => clm3%g%l%c%p%pvf%gammaS_out
    gamma_out        => clm3%g%l%c%p%pvf%gamma_out
    sla              => clm3%g%l%c%p%pps%slasha

    t_veg24          => clm3%g%l%c%p%pvs%t_veg24
    t_veg240         => clm3%g%l%c%p%pvs%t_veg240
    forc_solad24     => clm3%g%l%c%p%pvs%fsd24
    forc_solad240    => clm3%g%l%c%p%pvs%fsd240
    forc_solai24     => clm3%g%l%c%p%pvs%fsi24
    forc_solai240    => clm3%g%l%c%p%pvs%fsi240
    fsun24           => clm3%g%l%c%p%pvs%fsun24
    fsun240          => clm3%g%l%c%p%pvs%fsun240
    elai_p           => clm3%g%l%c%p%pvs%elai_p

    hardwire_sla(noveg)                                    = 0._r8     ! bare-soil

    hardwire_sla(ndllf_evr_tmp_tree)                       = 0.0125_r8 !needleleaf
    hardwire_sla(ndllf_evr_brl_tree)                       = 0.0125_r8 !Gordon Bonan suggests NET = 0.0076
    hardwire_sla(ndllf_dcd_brl_tree)                       = 0.0125_r8 !Gordon Bonan suggests NDT = 0.0200

    hardwire_sla(nbrdlf_evr_trp_tree)                      = 0.0250_r8 !broadleaf
    hardwire_sla(nbrdlf_evr_tmp_tree)                      = 0.0250_r8 !Gordon Bonan suggests BET = 0.0178
    hardwire_sla(nbrdlf_dcd_trp_tree)                      = 0.0250_r8 !Gordon Bonan suggests BDT = 0.0274
    hardwire_sla(nbrdlf_dcd_tmp_tree:nbrdlf_dcd_brl_shrub) = 0.0250_r8 

    hardwire_sla(nc3_arctic_grass:numpft)                  = 0.0200_r8 !grass/crop

! root depth (m) (defined based on Zeng et al., 2001, cf Guenther 2006)

    hardwire_droot(noveg)                                     = 0._r8   ! bare-soil
    hardwire_droot(ndllf_evr_tmp_tree:ndllf_evr_brl_tree)     = 1.8_r8  ! evergreen tree
    hardwire_droot(ndllf_dcd_brl_tree)                        = 2.0_r8  ! needleleaf deciduous boreal tree
    hardwire_droot(nbrdlf_evr_trp_tree:nbrdlf_evr_tmp_tree)   = 3.0_r8  ! broadleaf evergreen tree
    hardwire_droot(nbrdlf_dcd_trp_tree:nbrdlf_dcd_brl_tree)   = 2.0_r8  ! broadleaf deciduous tree
    hardwire_droot(nbrdlf_evr_shrub:nbrdlf_dcd_brl_shrub)     = 2.5_r8  ! shrub
    hardwire_droot(nc3_arctic_grass:numpft)                   = 1.5_r8  ! grass/crop

! initialize variables which get passed to the atmosphere
    vocflx(lbp:ubp, :)=0._r8

    ! Determine specific leaf array
    do fp = 1,num_soilp
       p = filter_soilp(fp)

       slarea(p) = hardwire_sla(ivt(p))

    end do


    ! Begin loop through voc species
    !_______________________________________________________________________________

    do n = 1, nvoc
       select case (n)

       case(1)	

          do fp = 1,num_soilp
             p = filter_soilp(fp)
             g = pgridcell(p)


             ! epsilon: use gridded values for 6 PFTs specified by MEGAN following
             ! -------  Guenther et al. (2006).  Map the numpft CLM PFTs to these 6.
             !          Units: [ug C m-2 h-1] (convert input files from units of 
             !                 [ug isop m-2 h-1])
    	     epsilon(p) = 0._r8

             ! isoprenes:
             if (     ivt(p) == ndllf_evr_tmp_tree  &
             .or.     ivt(p) == ndllf_evr_brl_tree) then     !fineleaf evergreen
                	epsilon(p) = efisop(2,g)*scale_mw
             else if (ivt(p) == ndllf_dcd_brl_tree) then     !fineleaf deciduous
                	epsilon(p) = efisop(3,g)*scale_mw
             else if (ivt(p) >= nbrdlf_evr_trp_tree &
             .and.    ivt(p) <= nbrdlf_dcd_brl_tree) then    !broadleaf trees
                	epsilon(p) = efisop(1,g)*scale_mw
             else if (ivt(p) >= nbrdlf_evr_shrub &
             .and.    ivt(p) <= nbrdlf_dcd_brl_shrub) then   !shrubs
                	epsilon(p) = efisop(4,g)*scale_mw
             else if (ivt(p) >= nc3_arctic_grass &
             .and.    ivt(p) <= nc4_grass) then              !grass
                	epsilon(p) = efisop(5,g)*scale_mw
             else if (ivt(p) >= ncorn) then                  !crops
                	epsilon(p) =efisop(6,g)*scale_mw
             end if

          end do

       case(2)

          do fp = 1,num_soilp
             p = filter_soilp(fp)
             g = pgridcell(p)

             ! epsilon: use values from table 3 in Guenther (1997) which originate in
             ! -------  Guenther et al. (1995). In the comments below, I mention the pft
             !          category as described in table 3. Some values were taken directly
             !          from Guenther et al. (1995). Units: [ugC g-1 h-1]
             !          Values were updated on 1/2002 (Guenther, personal communication)

             ! monoterpenes:
             epsilon(p) = 0._r8
             ! monoterpenes:
             if (     ivt(p) >= ndllf_evr_tmp_tree &
             .and.    ivt(p) <= ndllf_evr_brl_tree) then     !needleleaf evergreen
                epsilon(p) = 2.0_r8
             else if (ivt(p) == ndllf_dcd_brl_tree) then     !needleleaf deciduous
                epsilon(p) = 1.6_r8
             else if (ivt(p) >= nbrdlf_evr_trp_tree  &
             .and.    ivt(p) <= nbrdlf_dcd_brl_tree) then    !broadleaf everg trop
                epsilon(p) = 0.4_r8
             else if (ivt(p) >= nbrdlf_evr_shrub &
             .and.    ivt(p) <= nbrdlf_dcd_brl_shrub) then   !other woody veg
                epsilon(p) = 0.8_r8
             else if (ivt(p) >= nc3_arctic_grass &
             .and.    ivt(p) <= numpft) then                 !grass & crop
                epsilon(p) = 0.1_r8
             end if
          end do

       case (3)
          do fp = 1,num_soilp
             p = filter_soilp(fp)
             g = pgridcell(p)

             ! other VOCs (OVOCs)
             epsilon(p) = 1.0_r8                 !Guenther (personal communication)
          end do

       case (4)
          do fp = 1,num_soilp
             p = filter_soilp(fp)
             g = pgridcell(p)

             ! other reactive VOCs (ORVOCs)
             epsilon(p) = 1.0_r8                 !Guenther (personal communication)
          end do

       case (5)
          do fp = 1,num_soilp
             p = filter_soilp(fp)
             g = pgridcell(p)

             ! CO
             epsilon(p) = 0.3_r8                 !Guenther (personal communication)
          end do


       case default

          write(iulog,*)'only nvocs up to index 5 are currently supported'
          call endrun()

       end select
       
       
       ct_bad=0

       select case (n)

       case (1)

          do fp = 1,num_soilp
             p = filter_soilp(fp)
             g = pgridcell(p)
             c = pcolumn(p)


             ! gamma: Activity factor. Units [dimensionless]
             ! =====  For isoprene include activity factors for LAI,PPFD, T, leaf age, and soil moisture

             ! Activity factor for LAI (Guenther et al., 2006)
             !------------------------
             ! Guenther et al., 2006 eq 3
             if ( (fsun240(p) > 0.0_r8) .and. (fsun240(p) < 1.e30_r8) ) then 
                 gamma_l = cce * elai(p)
             else
                 gamma_l = cce1 * elai(p)
             end if
	     gammaL_out(p)=gamma_l

             ! Activity factor for PPFD (Guenther et al., 2006)
             !-------------------------
	     ! With distinction between sunlit and shaded leafs, weight scalings by
             ! fsun and fshade 
             ! Scale total incident par by fraction of sunlit leaves (added on 1/2002)
             ! multiply w/m2 by 4.6 to get umol/m2/s for par (added 8/14/02)

             ! fvitt -- forc_solad240, forc_solai240 can be zero when CLM finidat is specified
             !          which will cause par240 to be zero and produce NaNs via log(par240)
             ! dml   -- fsun240 can be equal to or greater than one before 10 day averages are
             !           set on startup or if a new pft comes online during land cover change.
             !           Avoid this problem by only doing calculations with fsun240 when fsun240 is
             !           between 0 and 1
             if ( (fsun240(p) > 0._r8) .and. (fsun240(p) < 1._r8) .and.  (forc_solad240(p) > 0._r8) &
             .and. (forc_solai240(p) > 0._r8)) then
                ! With alpha and cp calculated based on eq 6 and 7:
                ! Note indexing for accumulated variables is all at pft level
                ! SUN:
                par = (forc_solad(g,1) + fsun(p) * forc_solai(g,1)) * 4.6_r8
                par24 = (forc_solad24(p) + fsun24(p) * forc_solai24(p)) * 4.6_r8
                par240 = (forc_solad240(p) + fsun240(p) * forc_solai240(p)) * 4.6_r8
                alpha = ca1 - ca2 * log(par240)
                cp = ca3 * exp(ca2 * (par24-par0_sun))*par240**(0.6_r8)
                gamma_p = fsun(p) * ( cp * alpha*par * (1._r8 + alpha*alpha*par*par)**(-0.5_r8) )
	        paru_out(p)=par
		par24u_out(p)=par24
                par240u_out(p)=par240
                ! SHADE:
                par = ((1._r8 - fsun(p)) * forc_solai(g,1)) * 4.6_r8
                par24 = ((1._r8 - fsun24(p)) * forc_solai24(p)) * 4.6_r8
                par240 = ((1._r8 - fsun240(p)) * forc_solai240(p)) * 4.6_r8
                alpha = ca1 - ca2 * log(par240)
                cp = ca3 * exp(ca2 * (par24-par0_shade))*par240**(0.6_r8)
                par = ((1._r8 - fsun(p)) * forc_solai(g,1)) * 4.6_r8
                gamma_p = gamma_p + (1-fsun(p)) * (cp*alpha*par*(1._r8 + alpha*alpha*par*par)**(-0.5_r8))
                para_out(p)=par
		par24a_out(p)=par24
 		par240a_out(p)=par240
             else
                ! With fixed alpha and cp (from MEGAN User's Guide):
                ! SUN: direct + diffuse  
                par = (forc_solad(g,1) + fsun(p) * forc_solai(g,1)) * 4.6_r8
                alpha = alpha_fix
                cp = cp_fix
                gamma_p = fsun(p) * ( cp * alpha*par * (1._r8 + alpha*alpha*par*par)**(-0.5_r8) )
		paru_out(p)=par
	        par24u_out(p)=-999
	        par240u_out(p)=-999
                ! SHADE: diffuse 
                par = ((1._r8 - fsun(p)) * forc_solai(g,1)) * 4.6_r8
                gamma_p = gamma_p + (1-fsun(p)) * (cp*alpha*par*(1._r8 + alpha*alpha*par*par)**(-0.5_r8))
		para_out(p)=par
                par24a_out(p)=-999
                par240a_out(p)=-999
             end if 
             alpha_out(p)=alpha
             cp_out(p)=cp
             gammaP_out(p)=gamma_p


             ! Activity factor for temperature (Guenther et al., 2006)
             !--------------------------------
             if ( (t_veg240(p) > 0.0_r8) .and. (t_veg240(p) < 1.e30_r8) ) then 
                ! topt and Eopt from eq 8 and 9:
                topt = co1 + (co2 * (t_veg240(p)-tstd0))
                Eopt = co3 * exp (co4 * (t_veg24(p)-tstd0)) * exp(co4 * (t_veg240(p) -tstd0))
	     else
                topt = topt_fix
                Eopt = Eopt_fix
             endif 
             x = ( (1._r8/topt) - (1._r8/(t_veg(p))) ) / ct3
             gamma_t = Eopt * ( ct2 * exp(ct1 * x)/(ct2 - ct1 * (1._r8 - exp(ct2 * x))) )
             topt_out(p)=topt
             Eopt_out(p)=Eopt
             gammaT_out(p)=gamma_t


             ! Activity factor for leaf age (Guenther et al., 2006)
             !-----------------------------
             ! If not CNDV elai is constant therefore gamma_a=1.0
             ! gamma_a set to unity for evergreens (PFTs 1, 2, 4, 5)
             ! Note that we assume here that the time step is shorter than the number of 
             !days after budbreak required to induce isoprene emissions (ti=12 days) and 
             ! the number of days after budbreak to reach peak emission (tm=28 days)
	     if ( (ivt(p) == ndllf_dcd_brl_tree) .or. (ivt(p) >= nbrdlf_dcd_trp_tree) ) then  ! non-evergreen

                if ( (elai_p(p) > 0.0_r8) .and. (elai_p(p) < 1.e30_r8) )then 
                   elai_prev = 2._r8*elai_p(p)-elai(p)  ! have accumulated average lai over last timestep
                   if (elai_prev == elai(p)) then
                      fnew = 0.0_r8
                      fgro = 0.0_r8
                      fmat = 1.0_r8
                      fsen = 0.0_r8
                   else if (elai_prev > elai(p)) then
                      fnew = 0.0_r8
                      fgro = 0.0_r8
                      fmat = 1.0_r8 - (elai_prev - elai(p))/elai_prev
                      fsen = (elai_prev - elai(p))/elai_prev
                   else if (elai_prev < elai(p)) then
                      fnew = 1 - (elai_prev / elai(p))
                      fgro = 0.0_r8
                      fmat = (elai_prev / elai(p))
                      fsen = 0.0_r8
                   end if             
                
                   gamma_a = fnew * Anew + fgro * Agro + fmat * Amat + fsen * Asen
	        else
                   gamma_a = 1.0_r8
                end if

             else
                gamma_a = 1.0_r8
             end if
             gammaA_out(p)=gamma_a


             ! Activity factor for soil moisture (Guenther et al., 2006) 
             !----------------------------------
             ! Calculate the mean scaling factor throughout the root depth.
             ! wilting point potential is in units of matric potential (mm) 
             ! (1 J/Kg = 0.001 MPa, approx = 0.1 m)
             ! convert to volumetric soil water using equation 7.118 of the CLM4 Technical Note
             if ((clayfrac(p) > 0) .and. (sandfrac(p) > 0)) then 
               gamma_sm = 0._r8
	       nl=0._r8

               do j = 1,nlevsoi
	         if  (sum(dz(c,1:j)) < hardwire_droot(ivt(p)))  then
                   theta_ice = h2osoi_ice(c,j)/(dz(c,j)*denice)
                   wilt = ((smpmax/sucsat(c,j))**(-1._r8/bsw(c,j))) * (watsat(c,j) - theta_ice)
                   theta1 = wilt + deltheta1
                   if (h2osoi_vol(c,j) >= theta1) then 
             	      gamma_sm = gamma_sm + 1._r8
                   else if ( (h2osoi_vol(c,j) > wilt) .and. (h2osoi_vol(c,j) < theta1) ) then
		      gamma_sm = gamma_sm + ( h2osoi_vol(c,j) - wilt ) / deltheta1
                   else
		      gamma_sm = gamma_sm + 0._r8
                   end if
		   nl=nl+1._r8
                 end if
 	       end do 

	       if (nl > 0) then
	         gamma_sm = gamma_sm/nl
	       endif
             else
	       gamma_sm = 1.0_r8
             end if
             gammaS_out(p)=gamma_sm


             ! Calculate total scaling factor
             !--------------------------------
	     gamma(p) = gamma_l * gamma_p * gamma_t * gamma_a * gamma_sm
             if ( (gamma(p) >=0.0_r8) .and. (gamma(p)< 100._r8) ) then
                gamma_out(p)=gamma(p)
	     else
                gamma_out(p)=gamma(p)
                write(6,*) 'clh GAMMA: ',gamma(p),gamma_l,gamma_p,gamma_t,gamma_a,gamma_sm
             end if

          end do

       case (2,3,4,5)

          do fp = 1,num_soilp
             p = filter_soilp(fp)
             g = pgridcell(p)

             ! gamma: Activity factor. Units [dimensionless]
             ! -----  For monoterpenes, OVOCs, ORVOCs, CO include simple activity factors 
             !        for LAI and T only (Guenther et al., 1995)
             gamma_t = exp(bet * (t_veg(p) - tstd))
	     gamma(p)=gamma_t

          end do

       end select

       do fp = 1,num_soilp
          p = filter_soilp(fp)
          g = pgridcell(p)

          ! density: Source density factor [g dry weight foliar mass m-2 ground]
          ! -------  Other than isoprene, need to convert EF units from 
          ! [ug g-1 h-1] to [ug m-2 h-1]
          if (ivt(p) > noveg) then
             density = elai(p) / (slarea(p) * 0.5_r8)
          else
             density = 0._r8
          end if

          ! calculate the voc flux
          ! ----------------------
	  select case (n)

          case(1)

              vocflx(p,n) = epsilon(p) * gamma(p) * scaling_to_500_Tg

          case(2,3,4,5)

              vocflx(p,n) = epsilon(p) * gamma(p) * density

          end select


       end do   ! end pft loop

    end do   ! end voc species loop
    !_______________________________________________________________________________

    ! Calculate total voc flux and individual components for history output

    do fp = 1,num_soilp
       p = filter_soilp(fp)
       vocflx_tot(p) = 0._r8
    end do
    do n = 1, nvoc
       do fp = 1,num_soilp
          p = filter_soilp(fp)
          vocflx_tot(p) = vocflx_tot(p) + vocflx(p,n)
       end do
    end do
    do fp = 1,num_soilp
       p = filter_soilp(fp)
       g = pgridcell(p)
       vocflx_1(p) = vocflx(p,1)
       vocflx_2(p) = vocflx(p,2)
       vocflx_3(p) = vocflx(p,3)
       vocflx_4(p) = vocflx(p,4)
       vocflx_5(p) = vocflx(p,5)
    end do

  end subroutine VOCEmission

end module VOCEmissionMod