module SoilLittDecompMod

  !-----------------------------------------------------------------------
  ! !DESCRIPTION:
  ! Module holding routines used in litter and soil decomposition model
  ! for coupled carbon-nitrogen code.
  !
  ! !USES:
  use shr_kind_mod           , only : r8 => shr_kind_r8
  use shr_const_mod          , only : SHR_CONST_TKFRZ
  use decompMod              , only : bounds_type
  use perf_mod               , only : t_startf, t_stopf
  use perfMod_GPU
  use elm_varctl             , only : iulog, use_lch4, use_century_decomp
  use elm_varcon             , only : dzsoi_decomp
  use elm_varpar             , only : nlevdecomp, ndecomp_cascade_transitions, ndecomp_pools
  use DecompCascadeCNMod     , only : decomp_rate_constants_cn
  use DecompCascadeBGCMod    , only : decomp_rate_constants_bgc
  use NitrifDenitrifMod    , only : nitrif_denitrif
  use VerticalProfileMod     , only : decomp_vertprofiles
  use CNDecompCascadeConType , only : decomp_cascade_con
  use CNStateType            , only : cnstate_type
  !!  add phosphorus  -X. YANG
  use elm_varctl             , only : nu_com

  use PhotosynthesisType     , only : photosyns_type
  use CanopyStateType        , only : canopystate_type
  use SoilStateType          , only : soilstate_type
  use CH4Mod                 , only : ch4_type
  use cropType               , only : crop_type
  use ColumnDataType         , only : col_cs, col_cf
  use ColumnDataType         , only : col_ns, col_nf
  use ColumnDataType         , only : col_ps, col_pf
  use VegetationDataType     , only : veg_ps, veg_pf
  use ELMFatesInterfaceMod   , only : hlm_fates_interface_type
  ! clm interface & pflotran:
  use elm_varctl             , only : use_elm_interface, use_pflotran, pf_cmode
  use elm_varctl             , only : use_cn, use_fates
  use elm_instMod            , only : alm_fates
  !
  implicit none
  save
  private :: CNvariables_nan4pf  !pflotran
  !
  ! !PUBLIC MEMBER FUNCTIONS:

   public :: readSoilLittDecompParams
  public :: SoilLittDecompAlloc
  ! pflotran
  public :: SoilLittDecompAlloc2
  !
  type, public :: CNDecompParamsType
     real(r8), pointer :: dnp => null()         !denitrification proportion
  end type CNDecompParamsType

  type(CNDecompParamsType)  , public ::  CNDecompParamsInst
  !$acc declare create(CNDecompParamsInst)
  !-----------------------------------------------------------------------

contains

  !-----------------------------------------------------------------------
   subroutine readSoilLittDecompParams ( ncid )
     !
     ! !DESCRIPTION:
     ! Read parameters
     !
     ! !USES:
     use ncdio_pio    , only: file_desc_t,ncd_io
     use abortutils   , only: endrun
     use shr_log_mod  , only: errMsg => shr_log_errMsg

     ! !ARGUMENTS:
     implicit none
     type(file_desc_t),intent(inout) :: ncid   ! pio netCDF file id
     !
     ! !LOCAL VARIABLES:
     character(len=32)  :: subname = 'CNDecompParamsType'
     character(len=100) :: errCode = '-Error reading in parameters file:'
     logical            :: readv ! has variable been read in or not
     real(r8)           :: tempr ! temporary to read in constant
     character(len=100) :: tString ! temp. var for reading
     !-----------------------------------------------------------------------
     allocate(CNDecompParamsInst%dnp)
     tString='dnp'
     call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv)
     if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(__FILE__, __LINE__))
     CNDecompParamsInst%dnp=tempr

   end subroutine readSoilLittDecompParams

!-------------------------------------------------------------------------------------------------
  subroutine SoilLittDecompAlloc (bounds, num_soilc, filter_soilc,    &
                num_soilp, filter_soilp,                        &
                canopystate_vars, soilstate_vars,               &
                cnstate_vars, ch4_vars, dtime)

    !-----------------------------------------------------------------------------
    ! DESCRIPTION:
    ! Modified for clm_interface: 9/12/2015
    ! clm-bgc soil Module, can be called through clm_bgc_interface
    ! ONLY includes SOM decomposition & nitrification/denitrification
    ! CNAllocaiton is divided into 3 subroutines:
    ! (1) Allocation1_PlantNPDemand  is called in EcosystemDynNoLeaching1
    ! (2) Allocation2_ResolveNPLimit is called in SoilLittDecompAlloc (this subroutine)
    ! (3) Allocation3_PlantCNPAlloc  is called in SoilLittDecompAlloc2
    !-----------------------------------------------------------------------------

    ! !USES:
!    use AllocationMod , only: CNAllocation
      !$acc routine seq
    use AllocationMod , only: Allocation2_ResolveNPLimit ! Phase-2 of CNAllocation
    !
    ! !ARGUMENT:
    type(bounds_type)        , intent(in)    :: bounds
    integer                  , intent(in)    :: num_soilc          ! number of soil columns in filter
    integer                  , intent(in)    :: filter_soilc(:)    ! filter for soil columns
    integer                  , intent(in)    :: num_soilp          ! number of soil patches in filter
    integer                  , intent(in)    :: filter_soilp(:)    ! filter for soil patches
    type(canopystate_type)   , intent(in)    :: canopystate_vars
    type(soilstate_type)     , intent(in)    :: soilstate_vars
    type(cnstate_type)       , intent(inout) :: cnstate_vars
    type(ch4_type)           , intent(in)    :: ch4_vars

    ! add phosphorus --
!    type(crop_type)          , intent(in)    :: crop_vars
    real(r8),   intent(in)    :: dtime
    !
    ! !LOCAL VARIABLES:
    integer :: c,j,k,l,m                                                                               !indices
    integer :: fc                                                                                      !lake filter column index
    real(r8):: p_decomp_cpool_loss(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) !potential C loss from one pool to another
    real(r8):: immob(bounds%begc:bounds%endc,1:nlevdecomp)                                             !potential N immobilization
    real(r8):: immob_p(bounds%begc:bounds%endc,1:nlevdecomp)                                             !potential P immobilization
    real(r8):: ratio                                                                                   !temporary variable
    real(r8):: dnp                                                                        !denitrification proportion
    real(r8):: cn_decomp_pools(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_pools)
    real(r8):: cp_decomp_pools(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_pools)
    real(r8):: cp_decomp_pools_new(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_pools)               !C:P ratio of new SOM
    integer, parameter :: i_atm = 0
    ! For methane code
    real(r8):: phr_vr(bounds%begc:bounds%endc,1:nlevdecomp)                                            !potential HR (gC/m3/s)
    real(r8):: hrsum(bounds%begc:bounds%endc,1:nlevdecomp)                                             !sum of HR (gC/m2/s)

    character(len=256) :: event
    !-----------------------------------------------------------------------

    associate(                                                                                           &
         cascade_donor_pool               =>    decomp_cascade_con%cascade_donor_pool                  , & ! Input:  [integer  (:)     ]  which pool is C taken from for a given decomposition step
         cascade_receiver_pool            =>    decomp_cascade_con%cascade_receiver_pool               , & ! Input:  [integer  (:)     ]  which pool is C added to for a given decomposition step
         floating_cn_ratio_decomp_pools   =>    decomp_cascade_con%floating_cn_ratio_decomp_pools      , & ! Input:  [logical  (:)     ]  TRUE => pool has fixed C:N ratio
         floating_cp_ratio_decomp_pools   =>    decomp_cascade_con%floating_cp_ratio_decomp_pools      , & ! Input:  [logical  (:)     ]  TRUE => pool has fixed C:P ratio
         initial_cn_ratio                 =>    decomp_cascade_con%initial_cn_ratio                    , & ! Input:  [real(r8) (:)     ]  c:n ratio for initialization of pools
         initial_cp_ratio                 =>    decomp_cascade_con%initial_cp_ratio                    , & ! Input:  [real(r8) (:)     ]  c:p ratio for initialization of pools

         is_cwd                           =>    decomp_cascade_con%is_cwd                              , &
!
         fpi_vr                           =>    cnstate_vars%fpi_vr_col                                , & ! Input:  [real(r8) (:,:)   ]  fraction of potential immobilization for N(no units)
         fpi_p_vr                         =>    cnstate_vars%fpi_p_vr_col                              , & ! Input:  [real(r8) (:,:)   ]  fraction of potential immobilization for P (no units)
         rf_decomp_cascade                =>    cnstate_vars%rf_decomp_cascade_col                     , & ! Input:  [real(r8) (:,:,:) ]  respired fraction in decomposition step (frac)
         pathfrac_decomp_cascade          =>    cnstate_vars%pathfrac_decomp_cascade_col               , & ! Input:  [real(r8) (:,:,:) ]  what fraction of C leaving a given pool passes through a given transition (frac)

         decomp_npools_vr                 =>    col_ns%decomp_npools_vr                , & ! Input:  [real(r8) (:,:,:) ]  (gC/m3)  vertically-resolved decomposing (litter, cwd, soil) N pools
         decomp_ppools_vr                 =>    col_ps%decomp_ppools_vr              , & ! Input:  [real(r8) (:,:,:) ]  (gC/m3)  vertically-resolved decomposing (litter, cwd, soil) P pools

         decomp_cpools_vr                 =>    col_cs%decomp_cpools_vr                  , & ! Input:  [real(r8) (:,:,:) ]  (gC/m3)  vertically-resolved decomposing (litter, cwd, soil) c pools

         w_scalar                         =>    col_cf%w_scalar                           , & ! Input:  [real(r8) (:,:)   ]  fraction by which decomposition is limited by moisture availability

         decomp_cascade_ntransfer_vr      =>    col_nf%decomp_cascade_ntransfer_vr      , & ! Output: [real(r8) (:,:,:) ]  vert-res transfer of N from donor to receiver pool along decomp. cascade (gN/m3/s)
         decomp_cascade_sminn_flux_vr     =>    col_nf%decomp_cascade_sminn_flux_vr     , & ! Output: [real(r8) (:,:,:) ]  vert-res mineral N flux for transition along decomposition cascade (gN/m3/s)
         potential_immob_vr               =>    col_nf%potential_immob_vr               , & ! Output: [real(r8) (:,:)   ]
         sminn_to_denit_decomp_cascade_vr =>    col_nf%sminn_to_denit_decomp_cascade_vr , & ! Output: [real(r8) (:,:,:) ]
         gross_nmin_vr                    =>    col_nf%gross_nmin_vr                    , & ! Output: [real(r8) (:,:)   ]
         net_nmin_vr                      =>    col_nf%net_nmin_vr                      , & ! Output: [real(r8) (:,:)   ]
         gross_nmin                       =>    col_nf%gross_nmin                       , & ! Output: [real(r8) (:)     ]  gross rate of N mineralization (gN/m2/s)
         net_nmin                         =>    col_nf%net_nmin                         , & ! Output: [real(r8) (:)     ]  net rate of N mineralization (gN/m2/s)
         ! add phosphorus
         decomp_cascade_ptransfer_vr      =>    col_pf%decomp_cascade_ptransfer_vr    , & ! Output: [real(r8) (:,:,:) ]  vert-res transfer of P from donor to receiver pool along decomp. cascade (gP/m3/s)
         decomp_cascade_sminp_flux_vr     =>    col_pf%decomp_cascade_sminp_flux_vr   , & ! Output: [real(r8) (:,:,:) ]  vert-res mineral P flux for transition along decomposition cascade (gP/m3/s)
         potential_immob_p_vr             =>    col_pf%potential_immob_p_vr           , & ! Output: [real(r8) (:,:)   ]
         gross_pmin_vr                    =>    col_pf%gross_pmin_vr                  , & ! Output: [real(r8) (:,:)   ]
         net_pmin_vr                      =>    col_pf%net_pmin_vr                    , & ! Output: [real(r8) (:,:)   ]
         gross_pmin                       =>    col_pf%gross_pmin                     , & ! Output: [real(r8) (:)     ]  gross rate of P mineralization (gP/m2/s)
         net_pmin                         =>    col_pf%net_pmin                       , & ! Output: [real(r8) (:)     ]  net rate of P mineralization (gP/m2/s)

         decomp_cascade_hr_vr             =>    col_cf%decomp_cascade_hr_vr               , & ! Output: [real(r8) (:,:,:) ]  vertically-resolved het. resp. from decomposing C pools (gC/m3/s)
         decomp_cascade_ctransfer_vr      =>    col_cf%decomp_cascade_ctransfer_vr        , & ! Output: [real(r8) (:,:,:) ]  vertically-resolved het. resp. from decomposing C pools (gC/m3/s)
         decomp_k                         =>    col_cf%decomp_k                           , & ! Output: [real(r8) (:,:,:) ]  rate constant for decomposition (1./sec)
         phr_vr                           =>    col_cf%phr_vr                             , & ! Output: [real(r8) (:,:)   ]  potential HR (gC/m3/s)
         fphr                             =>    col_cf%fphr                               , & ! Output: [real(r8) (:,:)   ]  fraction of potential SOM + LITTER heterotrophic
         pmnf_decomp_cascade              =>    col_nf%pmnf_decomp_cascade                  , &
         pmpf_decomp_cascade              =>    col_pf%pmpf_decomp_cascade                , &
         soil_n_immob_flux                =>    col_nf%soil_n_immob_flux                    , &
         soil_n_immob_flux_vr             =>    col_nf%soil_n_immob_flux_vr                 , &
         soil_n_grossmin_flux             =>    col_nf%soil_n_grossmin_flux                 , &
         soil_p_immob_flux                =>    col_pf%soil_p_immob_flux                  , &
         soil_p_immob_flux_vr             =>    col_pf%soil_p_immob_flux_vr               , &
         soil_p_grossmin_flux             =>    col_pf%soil_p_grossmin_flux               , &
         actual_immob_vr                  =>    col_nf%actual_immob_vr                  , &
         actual_immob_p_vr                =>    col_pf%actual_immob_p_vr                &
         )

      !-------------------------------------------------------------------------------------------------
      ! call decomp_rate_constants_bgc() or decomp_rate_constants_cn(): now called in EcosystemDynNoLeaching1
      !-------------------------------------------------------------------------------------------------
      ! set initial values for potential C and N fluxes
      p_decomp_cpool_loss(bounds%begc : bounds%endc, :, :) = 0._r8
      pmnf_decomp_cascade(bounds%begc : bounds%endc, :, :) = 0._r8
      pmpf_decomp_cascade(bounds%begc : bounds%endc, :, :) = 0._r8    !! initial values for potential P fluxes

      ! column loop to calculate potential decomp rates and total immobilization
      ! demand.

      !! calculate c:n ratios of applicable pools
      do l = 1, ndecomp_pools
         if ( floating_cn_ratio_decomp_pools(l) ) then
            do j = 1,nlevdecomp
               do fc = 1,num_soilc
                  c = filter_soilc(fc)
                  if ( decomp_npools_vr(c,j,l) > 0._r8 ) then
                     cn_decomp_pools(c,j,l) = decomp_cpools_vr(c,j,l) / decomp_npools_vr(c,j,l)
                  end if
               end do
            end do
         else
            do j = 1,nlevdecomp
               do fc = 1,num_soilc
                  c = filter_soilc(fc)
                  cn_decomp_pools(c,j,l) = initial_cn_ratio(l)
               end do
            end do
         end if
      end do

      !! calculate c:p ratios of applicable pools
      do l = 1, ndecomp_pools
         if ( floating_cp_ratio_decomp_pools(l) ) then
            do j = 1,nlevdecomp
               do fc = 1,num_soilc
                  c = filter_soilc(fc)
                  if ( decomp_ppools_vr(c,j,l) > 0._r8 ) then
                     cp_decomp_pools(c,j,l) = decomp_cpools_vr(c,j,l) / decomp_ppools_vr(c,j,l)
                  end if
               end do
            end do
         else
            do j = 1,nlevdecomp
               do fc = 1,num_soilc
                  c = filter_soilc(fc)
                  cp_decomp_pools(c,j,l) = initial_cp_ratio(l)
               end do
            end do
         end if
         do j = 1,nlevdecomp
            do fc = 1,num_soilc
               c = filter_soilc(fc)
               cp_decomp_pools_new(c,j,l) = initial_cp_ratio(l)
            end do
         end do
      end do

      ! calculate the non-nitrogen-limited fluxes
      ! these fluxes include the  "/ dt" term to put them on a
      ! per second basis, since the rate constants have been
      ! calculated on a per timestep basis.

      do k = 1, ndecomp_cascade_transitions
         do j = 1,nlevdecomp
            do fc = 1,num_soilc
               c = filter_soilc(fc)

               if (decomp_cpools_vr(c,j,cascade_donor_pool(k)) > 0._r8 .and. &
                    decomp_k(c,j,cascade_donor_pool(k)) > 0._r8 ) then

                 p_decomp_cpool_loss(c,j,k) = decomp_cpools_vr(c,j,cascade_donor_pool(k)) &
                       * decomp_k(c,j,cascade_donor_pool(k))  * pathfrac_decomp_cascade(c,j,k)
                  if ( .not. floating_cn_ratio_decomp_pools(cascade_receiver_pool(k)) ) then  !! not transition of cwd to litter

                     if (cascade_receiver_pool(k) /= i_atm ) then  ! not 100% respiration
                        ratio = 0._r8

                        if (decomp_npools_vr(c,j,cascade_donor_pool(k)) > 0._r8) then
                           ratio = cn_decomp_pools(c,j,cascade_receiver_pool(k))/cn_decomp_pools(c,j,cascade_donor_pool(k))
                        endif
                        pmnf_decomp_cascade(c,j,k) = (p_decomp_cpool_loss(c,j,k) * (1.0_r8 - rf_decomp_cascade(c,j,k) - ratio) &
                             / cn_decomp_pools(c,j,cascade_receiver_pool(k)) )

                     else   ! 100% respiration
                        pmnf_decomp_cascade(c,j,k) = - p_decomp_cpool_loss(c,j,k) / cn_decomp_pools(c,j,cascade_donor_pool(k))
                     endif

                  else   ! CWD -> litter
                     pmnf_decomp_cascade(c,j,k) = 0._r8
                  end if
                  !!! add phosphorus fluxes
                  if ( .not. is_cwd( (cascade_receiver_pool(k)) ) ) then  !! not transition of cwd to litter

                     if (cascade_receiver_pool(k) /= i_atm ) then  ! not 100% respiration
                        ratio = 0._r8

                        if (decomp_ppools_vr(c,j,cascade_donor_pool(k)) > 0._r8) then
!                           ratio = cp_decomp_pools(c,j,cascade_receiver_pool(k))/cp_decomp_pools(c,j,cascade_donor_pool(k))
                           ratio = cp_decomp_pools_new(c,j,cascade_receiver_pool(k))/cp_decomp_pools(c,j,cascade_donor_pool(k))
                        endif

                        pmpf_decomp_cascade(c,j,k) = (p_decomp_cpool_loss(c,j,k) * (1.0_r8 - rf_decomp_cascade(c,j,k) - ratio) &
                             / cp_decomp_pools_new(c,j,cascade_receiver_pool(k)) )

                     else   ! 100% respiration
                        pmpf_decomp_cascade(c,j,k) = - p_decomp_cpool_loss(c,j,k) / cp_decomp_pools(c,j,cascade_donor_pool(k))
                     endif

                  else   ! CWD -> litter
                     pmpf_decomp_cascade(c,j,k) = 0._r8
                  end if

               end if
            end do
         end do
      end do

      ! Sum up all the potential immobilization fluxes (positive pmnf flux)
      ! and all the mineralization fluxes (negative pmnf flux)
      do j = 1,nlevdecomp
         do fc = 1,num_soilc
            c = filter_soilc(fc)
            immob(c,j) = 0._r8
            immob_p(c,j) = 0._r8
            gross_nmin_vr(c,j) = 0._r8
            gross_pmin_vr(c,j) = 0._r8
         end do
      end do
      do k = 1, ndecomp_cascade_transitions
         do j = 1,nlevdecomp
            do fc = 1,num_soilc
               c = filter_soilc(fc)
               if (pmnf_decomp_cascade(c,j,k) > 0._r8) then
                  immob(c,j) = immob(c,j) + pmnf_decomp_cascade(c,j,k)
               else
                  gross_nmin_vr(c,j) = gross_nmin_vr(c,j) - pmnf_decomp_cascade(c,j,k)
               end if
               if (pmpf_decomp_cascade(c,j,k) > 0._r8) then
                  immob_p(c,j) = immob_p(c,j) + pmpf_decomp_cascade(c,j,k)
               else
                  gross_pmin_vr(c,j) = gross_pmin_vr(c,j) - pmpf_decomp_cascade(c,j,k)
               end if
            end do
         end do
      end do

      do j = 1,nlevdecomp
         do fc = 1,num_soilc
            c = filter_soilc(fc)
            potential_immob_vr(c,j) = immob(c,j)
            potential_immob_p_vr(c,j) = immob_p(c,j)
         end do
      end do

      ! Add up potential hr for methane calculations
      do j = 1,nlevdecomp
         do fc = 1,num_soilc
            c = filter_soilc(fc)
            phr_vr(c,j) = 0._r8
         end do
      end do
      do k = 1, ndecomp_cascade_transitions
         do j = 1,nlevdecomp
            do fc = 1,num_soilc
               c = filter_soilc(fc)
               phr_vr(c,j) = phr_vr(c,j) + rf_decomp_cascade(c,j,k) * p_decomp_cpool_loss(c,j,k)
            end do
         end do
      end do


!-------------------------------------------------------------------------------------------------
! 'call decomp_vertprofiles()' (calc nfixation_prof) is moved to EcosystemDynNoLeaching1
! 'nfixation_prof' is used in 'calc_nuptake_prof' & 'calc_puptake_prof', which are called in Allocation1,2,3
!-------------------------------------------------------------------------------------------------
      call nitrif_denitrif(bounds, &
           num_soilc, filter_soilc, &
           soilstate_vars,  ch4_vars )

      ! now that potential N immobilization is known, call allocation
      ! to resolve the competition between plants and soil heterotrophs
      ! for available soil mineral N resource.
      ! in addition, calculate fpi_vr, fpi_p_vr, & fgp
      event = 'CNAllocation - phase-2'
      call t_start_lnd(event)
      call Allocation2_ResolveNPLimit(bounds,                       &
               num_soilc, filter_soilc, num_soilp, filter_soilp,    &
               cnstate_vars,                                        &
               soilstate_vars, dtime,                               &
               alm_fates)
      call t_stop_lnd(event)


      ! column loop to calculate actual immobilization and decomp rates, following
      ! resolution of plant/heterotroph  competition for mineral N


          !-------------------------------------------------------------------------------------------------
          ! delete c:n,c:p ratios calculation, they have been calculated at the beginning of this subroutine
          !-------------------------------------------------------------------------------------------------

          ! upon return from CNAllocation, the fraction of potential immobilization
      ! has been set (cnstate_vars%fpi_vr_col). now finish the decomp calculations.
      ! Only the immobilization steps are limited by fpi_vr (pmnf > 0)
      ! Also calculate denitrification losses as a simple proportion
      ! of mineralization flux.

      do k = 1, ndecomp_cascade_transitions
         do j = 1,nlevdecomp
            do fc = 1,num_soilc
               c = filter_soilc(fc)

               if (decomp_cpools_vr(c,j,cascade_donor_pool(k)) > 0._r8) then
                  if ( pmnf_decomp_cascade(c,j,k) > 0._r8 .and. pmpf_decomp_cascade(c,j,k) > 0._r8 ) then    ! N and P co-limitation
                     p_decomp_cpool_loss(c,j,k) = p_decomp_cpool_loss(c,j,k) * min( fpi_vr(c,j),fpi_p_vr(c,j) )
                     pmnf_decomp_cascade(c,j,k) = pmnf_decomp_cascade(c,j,k) * min( fpi_vr(c,j),fpi_p_vr(c,j) )
                     pmpf_decomp_cascade(c,j,k) = pmpf_decomp_cascade(c,j,k) * min( fpi_vr(c,j),fpi_p_vr(c,j) )   !!! immobilization step
                  elseif ( pmnf_decomp_cascade(c,j,k) > 0._r8 .and. pmpf_decomp_cascade(c,j,k) <= 0._r8 ) then  ! N limitation 
                     p_decomp_cpool_loss(c,j,k) = p_decomp_cpool_loss(c,j,k) * fpi_vr(c,j)
                     pmnf_decomp_cascade(c,j,k) = pmnf_decomp_cascade(c,j,k) * fpi_vr(c,j)
                     pmpf_decomp_cascade(c,j,k) = pmpf_decomp_cascade(c,j,k) * fpi_vr(c,j) !!! immobilization step
                  elseif ( pmnf_decomp_cascade(c,j,k) <= 0._r8 .and. pmpf_decomp_cascade(c,j,k) >  0._r8 ) then  ! P limitation 
                     p_decomp_cpool_loss(c,j,k) = p_decomp_cpool_loss(c,j,k) * fpi_p_vr(c,j)
                     pmnf_decomp_cascade(c,j,k) = pmnf_decomp_cascade(c,j,k) * fpi_p_vr(c,j)
                     pmpf_decomp_cascade(c,j,k) = pmpf_decomp_cascade(c,j,k) * fpi_p_vr(c,j) !!! immobilization step
                  end if
                  decomp_cascade_hr_vr(c,j,k) = rf_decomp_cascade(c,j,k) * p_decomp_cpool_loss(c,j,k)
                  decomp_cascade_ctransfer_vr(c,j,k) = (1._r8 - rf_decomp_cascade(c,j,k)) * p_decomp_cpool_loss(c,j,k)
                  if (decomp_npools_vr(c,j,cascade_donor_pool(k)) > 0._r8 .and. cascade_receiver_pool(k) /= i_atm) then
                     decomp_cascade_ntransfer_vr(c,j,k) = p_decomp_cpool_loss(c,j,k) / cn_decomp_pools(c,j,cascade_donor_pool(k))
                  else
                     decomp_cascade_ntransfer_vr(c,j,k) = 0._r8
                  endif
                  !!! phosphorus fluxes
                  if (decomp_ppools_vr(c,j,cascade_donor_pool(k)) > 0._r8 .and. cascade_receiver_pool(k) /= i_atm) then
                     decomp_cascade_ptransfer_vr(c,j,k) = p_decomp_cpool_loss(c,j,k) / cp_decomp_pools(c,j,cascade_donor_pool(k))
                  else
                     decomp_cascade_ptransfer_vr(c,j,k) = 0._r8
                  endif
                  if ( cascade_receiver_pool(k) /= 0 ) then
                     decomp_cascade_sminn_flux_vr(c,j,k) = pmnf_decomp_cascade(c,j,k)
                     decomp_cascade_sminp_flux_vr(c,j,k) = pmpf_decomp_cascade(c,j,k)
                  else  ! keep sign convention negative for terminal pools
                     decomp_cascade_sminn_flux_vr(c,j,k) = - pmnf_decomp_cascade(c,j,k)
                     decomp_cascade_sminp_flux_vr(c,j,k) = - pmpf_decomp_cascade(c,j,k)
                  endif
                  net_nmin_vr(c,j) = net_nmin_vr(c,j) - pmnf_decomp_cascade(c,j,k)
                  net_pmin_vr(c,j) = net_pmin_vr(c,j) - pmpf_decomp_cascade(c,j,k)
               else
                  decomp_cascade_ntransfer_vr(c,j,k) = 0._r8
                  decomp_cascade_ptransfer_vr(c,j,k) = 0._r8
                  decomp_cascade_sminn_flux_vr(c,j,k) = 0._r8
                  decomp_cascade_sminp_flux_vr(c,j,k) = 0._r8
               end if

            end do
         end do
      end do

      if (nu_com .eq. 'RD') then
         do fc = 1,num_soilc
            c = filter_soilc(fc)
            do j = 1,nlevdecomp
                gross_nmin_vr(c,j) = 0.0_r8
                gross_pmin_vr(c,j) = 0.0_r8
            end do
         end do
         do k = 1, ndecomp_cascade_transitions
            do j = 1,nlevdecomp
               do fc = 1,num_soilc
                  c = filter_soilc(fc)
                  if (pmnf_decomp_cascade(c,j,k) <= 0._r8) then
                      gross_nmin_vr(c,j) = gross_nmin_vr(c,j) - 1.0_r8*pmnf_decomp_cascade(c,j,k)
                  end if
                  if (pmpf_decomp_cascade(c,j,k) <= 0._r8) then
                      gross_pmin_vr(c,j) = gross_pmin_vr(c,j) - 1.0_r8*pmpf_decomp_cascade(c,j,k)
                  end if
                end do
             end do
          end do
      end if

      if (nu_com .ne. 'RD') then
      do fc = 1,num_soilc
          c = filter_soilc(fc)
          soil_n_immob_flux(c) =0.0_r8
          soil_p_immob_flux(c) = 0.0_r8
          soil_n_grossmin_flux(c) = 0.0_r8
          soil_p_grossmin_flux(c) = 0.0_r8
          do j = 1,nlevdecomp
              soil_n_immob_flux_vr(c,j) = 0.0_r8
              soil_p_immob_flux_vr(c,j) = 0.0_r8
              gross_nmin_vr(c,j) = 0.0_r8
              gross_pmin_vr(c,j) = 0.0_r8
          end do
      end do
      do k = 1, ndecomp_cascade_transitions
         do j = 1,nlevdecomp
            do fc = 1,num_soilc
               c = filter_soilc(fc)
                   if (pmnf_decomp_cascade(c,j,k) > 0._r8) then
                   soil_n_immob_flux(c) = soil_n_immob_flux(c) + pmnf_decomp_cascade(c,j,k)*dzsoi_decomp(j)
                   soil_n_immob_flux_vr(c,j) = soil_n_immob_flux_vr(c,j) + pmnf_decomp_cascade(c,j,k)
               else
                   soil_n_grossmin_flux(c) = soil_n_grossmin_flux(c) -1.0_r8*pmnf_decomp_cascade(c,j,k)*dzsoi_decomp(j)
                   gross_nmin_vr(c,j) = gross_nmin_vr(c,j) - 1.0_r8*pmnf_decomp_cascade(c,j,k)
               end if
               if (pmpf_decomp_cascade(c,j,k) > 0._r8) then
                   soil_p_immob_flux(c) = soil_p_immob_flux(c) + pmpf_decomp_cascade(c,j,k)*dzsoi_decomp(j)
                   soil_p_immob_flux_vr(c,j) = soil_p_immob_flux_vr(c,j) + pmpf_decomp_cascade(c,j,k)
               else
                   soil_p_grossmin_flux(c) = soil_p_grossmin_flux(c) -1.0_r8*pmpf_decomp_cascade(c,j,k)*dzsoi_decomp(j)
                   gross_pmin_vr(c,j) = gross_pmin_vr(c,j) - 1.0_r8*pmpf_decomp_cascade(c,j,k)
               end if
             end do
          end do
       end do
       end if

      if (use_lch4) then
         ! Calculate total fraction of potential HR, for methane code
         do j = 1,nlevdecomp
            do fc = 1,num_soilc
               c = filter_soilc(fc)
               hrsum(c,j) = 0._r8
            end do
         end do
         do k = 1, ndecomp_cascade_transitions
            do j = 1,nlevdecomp
               do fc = 1,num_soilc
                  c = filter_soilc(fc)
                  hrsum(c,j) = hrsum(c,j) + rf_decomp_cascade(c,j,k) * p_decomp_cpool_loss(c,j,k)
               end do
            end do
         end do

         ! Nitrogen limitation / (low)-moisture limitation
         do j = 1,nlevdecomp
            do fc = 1,num_soilc
               c = filter_soilc(fc)
               if (phr_vr(c,j) > 0._r8) then
                  fphr(c,j) = hrsum(c,j) / phr_vr(c,j) * w_scalar(c,j)
                  fphr(c,j) = max(fphr(c,j), 0.01_r8) ! Prevent overflow errors for 0 respiration
               else
                  fphr(c,j) = 1._r8
               end if
            end do
         end do
      end if

      ! vertically integrate net and gross mineralization fluxes for diagnostic output
      ! moved to SoilLittDecompAlloc2
!      do j = 1,nlevdecomp
!         do fc = 1,num_soilc
!            c = filter_soilc(fc)
!            net_nmin(c) = net_nmin(c) + net_nmin_vr(c,j) * dzsoi_decomp(j)
!            gross_nmin(c) = gross_nmin(c) + gross_nmin_vr(c,j) * dzsoi_decomp(j)
!            net_pmin(c) = net_pmin(c) + net_pmin_vr(c,j) * dzsoi_decomp(j)
!            gross_pmin(c) = gross_pmin(c) + gross_pmin_vr(c,j) * dzsoi_decomp(j)
!         end do
!      end do

    end associate

  end subroutine SoilLittDecompAlloc

!-------------------------------------------------------------------------------------------------

  subroutine SoilLittDecompAlloc2 (bounds, num_soilc, filter_soilc, num_soilp, filter_soilp,   &
       photosyns_vars, canopystate_vars, soilstate_vars,          &
       cnstate_vars, ch4_vars, crop_vars, atm2lnd_vars, dt)
    !-----------------------------------------------------------------------------
    ! DESCRIPTION:
    ! bgc interface & pflotran:
    ! (1) Simplified codes of SoilLittDecompAlloc subroutine for coupling with pflotran
    ! (2) call Allocation3_PlantCNPAlloc
    ! (3) calculate net_nmin(c), gross_nmin(c), net_pmin(c), gross_pmin(c)
    !-----------------------------------------------------------------------------

    ! !USES:
      !$acc routine seq
    use AllocationMod , only: Allocation3_PlantCNPAlloc ! Phase-3 of CNAllocation
    use atm2lndType     , only: atm2lnd_type

    !
    ! !ARGUMENT:
    type(bounds_type)        , intent(in)    :: bounds
    integer                  , intent(in)    :: num_soilc          ! number of soil columns in filter
    integer                  , intent(in)    :: filter_soilc(:)    ! filter for soil columns
    integer                  , intent(in)    :: num_soilp          ! number of soil patches in filter
    integer                  , intent(in)    :: filter_soilp(:)    ! filter for soil patches
    type(photosyns_type)     , intent(in)    :: photosyns_vars
    type(canopystate_type)   , intent(in)    :: canopystate_vars
    type(soilstate_type)     , intent(in)    :: soilstate_vars
    type(cnstate_type)       , intent(inout) :: cnstate_vars
    type(ch4_type)           , intent(in)    :: ch4_vars
    type(crop_type)          , intent(inout) :: crop_vars
    type(atm2lnd_type)       , intent(in)    :: atm2lnd_vars
    real(r8), intent(in) :: dt                                           ! time step (seconds)

    !
    ! !LOCAL VARIABLES:
    integer :: fc, c, j                                     ! indices
!    real(r8):: col_plant_ndemand(bounds%begc:bounds%endc)   ! column-level vertically-integrated plant N demand (gN/m2/s)
    real(r8) :: smin_nh4_to_plant_vr_loc(bounds%begc:bounds%endc,1:nlevdecomp)
    real(r8) :: smin_no3_to_plant_vr_loc(bounds%begc:bounds%endc,1:nlevdecomp)

    ! For methane code
    real(r8):: hrsum(bounds%begc:bounds%endc,1:nlevdecomp)                                             !sum of HR (gC/m2/s)

    character(len=256) :: event
    !-----------------------------------------------------------------------

    associate(                                                                                      &
         gross_nmin_vr                    =>    col_nf%gross_nmin_vr                    , & ! Output: [real(r8) (:,:)   ]
         net_nmin_vr                      =>    col_nf%net_nmin_vr                      , & ! Output: [real(r8) (:,:)   ]
         gross_nmin                       =>    col_nf%gross_nmin                       , & ! Output: [real(r8) (:)     ]  gross rate of N mineralization (gN/m2/s)
         net_nmin                         =>    col_nf%net_nmin                         , & ! Output: [real(r8) (:)     ]  net rate of N mineralization (gN/m2/s)
         !phosphorus
         gross_pmin_vr                    =>    col_pf%gross_pmin_vr                  , & ! Output: [real(r8) (:,:)   ]
         net_pmin_vr                      =>    col_pf%net_pmin_vr                    , & ! Output: [real(r8) (:,:)   ]
         gross_pmin                       =>    col_pf%gross_pmin                     , & ! Output: [real(r8) (:)     ]  gross rate of P mineralization (gP/m2/s)
         net_pmin                         =>    col_pf%net_pmin                       , & ! Output: [real(r8) (:)     ]  net rate of P mineralization (gP/m2/s)


         fpi_vr                           =>    cnstate_vars%fpi_vr_col                                , & ! Output:  [real(r8) (:,:)   ]  fraction of potential immobilization (no units)
         fpi                              =>    cnstate_vars%fpi_col                                   , & ! Output: [real(r8) (:)   ]  fraction of potential immobilization (no units)
         potential_immob_vr               =>    col_nf%potential_immob_vr               , & ! Input:
         actual_immob_vr                  =>    col_nf%actual_immob_vr                  , & ! Input:
         potential_immob                  =>    col_nf%potential_immob                  , & ! Output: [real(r8) (:)   ]
         actual_immob                     =>    col_nf%actual_immob                     , & ! Output: [real(r8) (:)   ]

         fpg                              =>    cnstate_vars%fpg_col                                   , & ! Output: [real(r8) (:)   ]  fraction of potential gpp (no units)
         sminn_to_plant                   =>    col_nf%sminn_to_plant                   , & ! Output: [real(r8) (:)     ]  col N uptake (gN/m2/s)
         sminn_to_plant_vr                =>    col_nf%sminn_to_plant_vr                , & ! Input:  [real(r8) (:,:)    ]  vertically-resolved N uptake (gN/m3/s)

         smin_no3_to_plant_vr             =>    col_nf%smin_no3_to_plant_vr             , & ! Output: [real(r8) (:,:) ]
         smin_nh4_to_plant_vr             =>    col_nf%smin_nh4_to_plant_vr             , & ! Output: [real(r8) (:,:) ]

         col_plant_ndemand_vr             =>    col_nf%plant_ndemand_vr                 , & ! Input:  [real(r8) (:)     ]  col N uptake (gN/m2/s)

         plant_ndemand_col                =>    col_nf%plant_ndemand                    , & ! Output:  [real(r8) (:,:) ]

         w_scalar                         =>    col_cf%w_scalar                           , & ! Input:  [real(r8) (:,:)   ]  fraction by which decomposition is limited by moisture availability
         decomp_cascade_hr_vr             =>    col_cf%decomp_cascade_hr_vr               , & ! Output: [real(r8) (:,:,:) ]  vertically-resolved het. resp. from decomposing C pools (gC/m3/s)
         decomp_cascade_ctransfer_vr      =>    col_cf%decomp_cascade_ctransfer_vr        , & ! Output: [real(r8) (:,:,:) ]  vertically-resolved het. resp. from decomposing C pools (gC/m3/s)
         decomp_k                         =>    col_cf%decomp_k                           , & ! Output: [real(r8) (:,:,:) ]  rate constant for decomposition (1./sec)
         hr_vr                            =>    col_cf%hr_vr                              , & ! Output: [real(r8) (:,:)   ]  potential HR (gC/m3/s)
         phr_vr                           =>    col_cf%phr_vr                             , & ! Output: [real(r8) (:,:)   ]  potential HR (gC/m3/s)
         fphr                             =>    col_cf%fphr                               , & ! Output: [real(r8) (:,:)   ]  fraction of potential SOM + LITTER heterotrophic

         smin_no3_vr                      =>    col_ns%smin_no3_vr                     , &
         smin_nh4_vr                      =>    col_ns%smin_nh4_vr                       &
         )

      ! set time steps
            !------------------------------------------------------------------
            ! 'call decomp_vertprofiles()' moved to EcosystemDynNoLeaching1
            !------------------------------------------------------------------
            smin_nh4_to_plant_vr_loc(:,:) = 0._r8
            smin_no3_to_plant_vr_loc(:,:) = 0._r8


      ! MUST have already updated needed bgc variables from PFLOTRAN by this point
      if(use_elm_interface.and.use_pflotran.and.pf_cmode) then
         ! fpg calculation
         do fc=1,num_soilc
            c = filter_soilc(fc)
            sminn_to_plant(c)       = 0._r8
            do j = 1, nlevdecomp           ! sum up actual and potential column-level N fluxes to plant
               sminn_to_plant(c)    = sminn_to_plant(c) + sminn_to_plant_vr(c,j) * dzsoi_decomp(j)
            end do
         end do
         do fc=1,num_soilc
            c = filter_soilc(fc)
            ! calculate the fraction of potential growth that can be
            ! acheived with the N available to plants
            if (plant_ndemand_col(c) > 0.0_r8) then
               fpg(c) = max(0._r8,sminn_to_plant(c)) / plant_ndemand_col(c)
               fpg(c) = min(1._r8, fpg(c))
            else
               fpg(c) = 1.0_r8
            end if
         end do

         ! fpi calculation
         do fc=1,num_soilc
            c = filter_soilc(fc)
            potential_immob(c) = 0._r8
            actual_immob(c)    = 0._r8
            do j = 1, nlevdecomp
               if (potential_immob_vr(c,j) > 0.0_r8) then
                  fpi_vr(c,j) = actual_immob_vr(c,j) / potential_immob_vr(c,j)
                  potential_immob(c) = potential_immob(c) + potential_immob_vr(c,j)*dzsoi_decomp(j)
                  actual_immob(c) = actual_immob(c) + actual_immob_vr(c,j)*dzsoi_decomp(j)
               else
                  fpi_vr(c,j) = 0.0_r8
               end if
            end do
         end do
         do fc=1,num_soilc
            c = filter_soilc(fc)
            if (potential_immob(c) > 0.0_r8) then
               fpi(c) = max(0._r8,actual_immob(c)) / potential_immob(c)
               fpi(c) = min(1._r8, fpi(c))
            else
               fpi(c) = 1.0_r8
            end if
         end do


         if (use_lch4) then
            ! Add up potential hr for methane calculations
            ! potential hr is not available from PFLOTRAN, so here temporarily as actual hr
            ! in the end, this methane module will be moving into PFLOTRAN as well

            do j = 1,nlevdecomp
               do fc = 1, num_soilc
                  c = filter_soilc(fc)
                  phr_vr(c,j) = hr_vr(c,j)
               end do
            end do

            ! Calculate total fraction of potential HR, for methane code
            do j = 1,nlevdecomp
               do fc = 1,num_soilc
                  c = filter_soilc(fc)
                  hrsum(c,j) = hr_vr(c,j)
               end do
            end do

            ! Nitrogen limitation / (low)-moisture limitation
            do j = 1,nlevdecomp
               do fc = 1,num_soilc
                  c = filter_soilc(fc)
                  if (phr_vr(c,j) > 0._r8) then
                     fphr(c,j) = hrsum(c,j) / phr_vr(c,j) * w_scalar(c,j)
                     fphr(c,j) = max(fphr(c,j), 0.01_r8) ! Prevent overflow errors for 0 respiration
                  else
                     fphr(c,j) = 1._r8
                  end if
               end do
            end do
         end if

         ! needs to zero CLM-CNP variables NOT available from pflotran bgc coupling
         call CNvariables_nan4pf(bounds, num_soilc, filter_soilc, &
                        num_soilp, filter_soilp)

         ! save variables before updating
         do fc = 1,num_soilc
            c = filter_soilc(fc)
            do j = 1,nlevdecomp
                smin_no3_to_plant_vr_loc(c,j) = smin_no3_to_plant_vr(c,j)
                smin_nh4_to_plant_vr_loc(c,j) = smin_nh4_to_plant_vr(c,j)
            end do
         end do

      end if !if(use_elm_interface.and.use_pflotran.and.pf_cmode)

      !------------------------------------------------------------------
      ! phase-3 Allocation for plants
      if(.not.use_fates)then
        event = 'CNAllocation - phase-3'
        call t_start_lnd(event)
        call Allocation3_PlantCNPAlloc (bounds                      , &
                  num_soilc, filter_soilc, num_soilp, filter_soilp    , &
                  canopystate_vars                                    , &
                  cnstate_vars, crop_vars, dt)
        call t_stop_lnd(event)
      end if
      !------------------------------------------------------------------

    if(use_pflotran.and.pf_cmode) then
    ! in Allocation3_PlantCNPAlloc():
    ! smin_nh4_to_plant_vr(c,j), smin_no3_to_plant_vr(c,j), sminn_to_plant_vr(c,j) may be adjusted
    ! therefore, we need to update smin_no3_vr(c,j) & smin_nh4_vr(c,j)
      do fc = 1,num_soilc
           c = filter_soilc(fc)
           do j = 1,nlevdecomp
               smin_no3_vr(c,j) = smin_no3_vr(c,j) - (smin_no3_to_plant_vr(c,j) - smin_no3_to_plant_vr_loc(c,j))*dt
               smin_nh4_vr(c,j) = smin_nh4_vr(c,j) - (smin_nh4_to_plant_vr(c,j) - smin_nh4_to_plant_vr_loc(c,j))*dt
               smin_no3_vr(c,j) = max(0._r8, smin_no3_vr(c,j))
               smin_nh4_vr(c,j) = max(0._r8, smin_nh4_vr(c,j))
            end do
      end do
    end if !(use_pflotran.and.pf_cmode)
    !------------------------------------------------------------------
      ! vertically integrate net and gross mineralization fluxes for diagnostic output
      do fc=1,num_soilc
         c = filter_soilc(fc)
         net_nmin(c)    = 0._r8
         gross_nmin(c)  = 0._r8
         net_pmin(c)    = 0._r8
         gross_pmin(c)  = 0._r8
      end do
      do j = 1,nlevdecomp
         do fc = 1,num_soilc
            c = filter_soilc(fc)
            net_nmin(c) = net_nmin(c) + net_nmin_vr(c,j) * dzsoi_decomp(j)
            gross_nmin(c) = gross_nmin(c) + gross_nmin_vr(c,j) * dzsoi_decomp(j)
            !phosphorus
            net_pmin(c) = net_pmin(c) + net_pmin_vr(c,j) * dzsoi_decomp(j)
            gross_pmin(c) = gross_pmin(c) + gross_pmin_vr(c,j) * dzsoi_decomp(j)
         end do
      end do
    end associate
  end subroutine SoilLittDecompAlloc2

  !-------------------------------------------------------------------------------------------------
  !
  subroutine CNvariables_nan4pf (bounds, num_soilc, filter_soilc, num_soilp, filter_soilp)
  !
  !DESCRIPTION:
  !  CN variables not available from PFLOTRAN, some of which may be output and may cause issues,
  !  if not properly set.
  !
  !USES:
    !$acc routine seq
    use elm_varctl   , only: carbon_only, carbonnitrogen_only
    use elm_varpar   , only: nlevdecomp, ndecomp_cascade_transitions
   !
   !ARGUMENTS:
    type(bounds_type)        , intent(in)    :: bounds
    integer                  , intent(in)    :: num_soilc          ! number of soil columns in filter
    integer                  , intent(in)    :: filter_soilc(:)    ! filter for soil columns
    integer                  , intent(in)    :: num_soilp          ! number of soil patches in filter
    integer                  , intent(in)    :: filter_soilp(:)    ! filter for soil patches
   !
   !CALLED FROM:
   !
   !LOCAL VARIABLES:
   integer :: c,j,k,fc          !indices
   !
   !-----------------------------------------------------------------------
   associate (&
         decomp_cascade_hr_vr             =>    col_cf%decomp_cascade_hr_vr               , & ! Output: [real(r8) (:,:,:) ]  vertically-resolved het. resp. from decomposing C pools (gC/m3/s)
         decomp_cascade_ctransfer_vr      =>    col_cf%decomp_cascade_ctransfer_vr        , & ! Output: [real(r8) (:,:,:) ]  vertically-resolved het. resp. from decomposing C pools (gC/m3/s)
         decomp_cascade_ntransfer_vr      =>    col_nf%decomp_cascade_ntransfer_vr        & ! Output: [real(r8) (:,:,:) ]  vert-res transfer of N from donor to receiver pool along decomp. cascade (gN/m3/s)
   )
   ! set zeros for those variables NOT available from PFLOTRAN
   ! (TODO) will check 'zero' or '-9999' is better.

   do fc = 1,num_soilc
      c = filter_soilc(fc)
      do j = 1,nlevdecomp
         do k = 1, ndecomp_cascade_transitions
            decomp_cascade_ctransfer_vr(c,j,k) = 0._r8
            decomp_cascade_ntransfer_vr(c,j,k) = 0._r8
            decomp_cascade_hr_vr(c,j,k)        = 0._r8
         end do
      end do

   end do

   ! pflotran not yet support phosphous cycle
   if ( carbon_only .or.  carbonnitrogen_only  ) then
      call veg_ps%SetValues(num_patch=num_soilp,  filter_patch=filter_soilp,  value_patch=0._r8)
      call col_ps%SetValues(num_column=num_soilc, filter_column=filter_soilc, value_column=0._r8)

      call veg_pf%setvalues( num_patch=num_soilp,  filter_patch=filter_soilp,  value_patch=0._r8)
      call col_pf%setvalues( num_column=num_soilc, filter_column=filter_soilc, value_column=0._r8)
   end if

  end associate
  end subroutine CNvariables_nan4pf

  !-------------------------------------------------------------------------------------------------

end module SoilLittDecompMod
