module EDCanopyStructureMod

  ! =====================================================================================
  ! Code to determine whether the canopy is closed, and which plants are either in the
  ! understorey or overstorey. This is obviosuly far too complicated for it's own good
  ! =====================================================================================

  use FatesConstantsMod     , only : r8 => fates_r8
  use FatesConstantsMod     , only : itrue, ifalse
  use FatesConstantsMod     , only : tinyr8
  use FatesConstantsMod     , only : nearzero
  use FatesConstantsMod     , only : rsnbl_math_prec
  use FatesConstantsMod     , only : nocomp_bareground
  use FatesGlobals          , only : fates_log
  use EDPftvarcon           , only : EDPftvarcon_inst
  use PRTParametersMod      , only : prt_params
  use FatesAllometryMod     , only : carea_allom
  use EDCohortDynamicsMod   , only : copy_cohort, terminate_cohorts, terminate_cohort, fuse_cohorts
  use EDCohortDynamicsMod   , only : InitPRTObject
  use EDCohortDynamicsMod   , only : InitPRTBoundaryConditions
  use FatesAllometryMod     , only : tree_lai
  use FatesAllometryMod     , only : tree_sai
  use EDtypesMod            , only : ed_site_type, ed_patch_type, ed_cohort_type
  use EDTypesMod            , only : nclmax
  use EDTypesMod            , only : nlevleaf
  use EDtypesMod            , only : AREA
  use EDLoggingMortalityMod , only : UpdateHarvestC
  use FatesGlobals          , only : endrun => fates_endrun
  use FatesInterfaceTypesMod     , only : hlm_days_per_year
  use FatesInterfaceTypesMod     , only : hlm_use_planthydro
  use FatesInterfaceTypesMod     , only : hlm_use_cohort_age_tracking
  use FatesInterfaceTypesMod     , only : hlm_use_sp
  use FatesInterfaceTypesMod     , only : numpft
  use FatesInterfaceTypesMod, only : bc_in_type
  use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage
  use PRTGenericMod,          only : leaf_organ
  use PRTGenericMod,          only : leaf_organ
  use PRTGenericMod,          only : fnrt_organ
  use PRTGenericMod,          only : sapw_organ
  use PRTGenericMod,          only : store_organ
  use PRTGenericMod,          only : repro_organ
  use PRTGenericMod,          only : struct_organ
  use PRTGenericMod,          only : SetState
  use PRTGenericMod,          only : carbon12_element

  ! CIME Globals
  use shr_log_mod           , only : errMsg => shr_log_errMsg

  implicit none
  private

  public :: canopy_structure
  public :: canopy_spread
  public :: calc_areaindex
  public :: canopy_summarization
  public :: update_hlm_dynamics
  public :: UpdateFatesAvgSnowDepth
  public :: UpdatePatchLAI
  public :: UpdateCohortLAI

  logical, parameter :: debug=.false.

  character(len=*), parameter, private :: sourcefile = &
       __FILE__

  integer :: istat           ! return status code
  character(len=255) :: smsg ! Message string for deallocation errors
  
  real(r8), parameter :: area_target_precision = 1.0E-11_r8  ! Area conservation
  ! will attempt to reduce errors
  ! below this level

  real(r8), parameter :: area_check_precision  = 1.0E-7_r8     ! Area conservation checks must
  ! be within this absolute tolerance
  real(r8), parameter :: area_check_rel_precision = 1.0E-4_r8  ! Area conservation checks must
  ! be within this relative tolerance

  real(r8), parameter :: similar_height_tol = 1.0E-3_r8    ! I think trees that differ by 1mm
  ! can be roughly considered the same right?


  ! 10/30/09: Created by Rosie Fisher
  ! 2017/2018: Modifications and updates by Ryan Knox
  ! ============================================================================

contains

  ! ============================================================================
  subroutine canopy_structure( currentSite , bc_in )
    !
    ! !DESCRIPTION:
    ! create cohort instance
    !
    ! This routine allocates the 'canopy_layer' attribute to each cohort
    ! All top leaves in the same canopy layer get the same light resources.
    ! The first canopy layer is the 'canopy' or 'overstorey'. The second is the 'understorey'.
    ! More than two layers is not permitted at the moment
    ! Seeds germinating into the 3rd or higher layers are automatically removed.
    !
    ! ------Perfect Plasticity-----
    ! The idea of these canopy layers derives originally from Purves et al. 2009
    ! Their concept is that, given enoughplasticity in canopy position, size, shape and depth
    ! all of the gound area will be filled perfectly by leaves, and additional leaves will have
    ! to exist in the understorey.
    ! Purves et al. use the concept of 'Z*' to assume that the height required to attain a place in the
    ! canopy is spatially uniform. In this implementation, described in Fisher et al. (2010, New Phyt) we
    ! extent that concept to assume that position in the canopy has some random element, and that BOTH height
    ! and chance combine to determine whether trees get into the canopy.
    ! Thus, when the canopy is closed and there is excess area, some of it must be demoted
    ! If we demote -all- the trees less than a given height, there is a massive advantage in being the cohort that is
    ! the biggest when the canopy is closed.
    ! In this implementation, the amount demoted, ('weight') is a function of the height weighted by the competitive exclusion
    ! parameter (ED_val_comp_excln).

    ! Complexity in this routine results from a few things.
    ! Firstly, the complication of the demotion amount sometimes being larger than the cohort area (for a very small, short cohort)
    ! Second, occasionaly, disturbance (specifically fire) can cause the canopy layer to become less than closed,
    ! without changing the area of the patch. If this happens, then some of the plants in the lower layer need to be 'promoted' so
    ! all of the routine has to happen in both the downwards and upwards directions.
    !
    ! The order of events here is therefore:
    ! (The entire subroutine has a single outer 'patch' loop.
    ! Section 1: figure out the total area, and whether there are >1 canopy layers at all.
    !
    ! Sorts out cohorts into canopy and understorey layers...
    !
    ! !USES:

    use EDParamsMod, only : ED_val_comp_excln
    use EDTypesMod , only : min_patch_area

    !
    ! !ARGUMENTS
    type(ed_site_type) , intent(inout), target   :: currentSite
    type(bc_in_type), intent(in)                 :: bc_in

    !
    ! !LOCAL VARIABLES:
    type(ed_patch_type) , pointer :: currentPatch
    type(ed_cohort_type), pointer :: currentCohort
    integer  :: i_lyr                  ! current layer index
    integer  :: z                      ! Current number of canopy layers. (1= canopy, 2 = understorey)
    integer  :: ipft
    real(r8) :: arealayer(nclmax+2)    ! Amount of plant area currently in each canopy layer
    integer  :: patch_area_counter     ! count iterations used to solve canopy areas
    logical  :: area_not_balanced      ! logical controlling if the patch layer areas
    ! have successfully been redistributed
    integer  :: return_code            ! math checks on variables will return>0 if problems exist
    ! We only iterate because of possible imprecisions generated by the cohort
    ! termination process.  These should be super small, so at the most
    ! try to re-balance 3 times.  If that doesn't give layer areas
    ! within tolerance of canopy area, there is something wrong

    integer, parameter  :: max_patch_iterations = 10


    !----------------------------------------------------------------------
    currentPatch => currentSite%oldest_patch
    !
    ! zero site-level demotion / promotion tracking info
    currentSite%demotion_rate(:) = 0._r8
    currentSite%promotion_rate(:) = 0._r8
    currentSite%demotion_carbonflux = 0._r8
    currentSite%promotion_carbonflux = 0._r8


    !
    ! Section 1: Check  total canopy area.
    !
    do while (associated(currentPatch)) ! Patch loop

       ! ------------------------------------------------------------------------------
       ! Perform numerical checks on some cohort and patch structures
       ! ------------------------------------------------------------------------------

       ! canopy layer has a special bounds check
       currentCohort => currentPatch%tallest
       do while (associated(currentCohort))
          if( currentCohort%canopy_layer < 1 .or. currentCohort%canopy_layer > nclmax+1 ) then
             write(fates_log(),*) 'lat:',currentSite%lat
             write(fates_log(),*) 'lon:',currentSite%lon
             write(fates_log(),*) 'BOGUS CANOPY LAYER: ',currentCohort%canopy_layer
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if
          currentCohort => currentCohort%shorter
       enddo


       ! Does any layer have excess area in it? Keep going until it does not...
       patch_area_counter = 0
       area_not_balanced = .true.

       do while(area_not_balanced)

          ! ---------------------------------------------------------------------------
          ! Demotion Phase: Identify upper layers that are too full, and demote them to
          ! the layers below.
          ! ---------------------------------------------------------------------------

          ! Its possible that before we even enter this scheme
          ! some cohort numbers are very low.  Terminate them.
          call terminate_cohorts(currentSite, currentPatch, 1, 12, bc_in)

          ! Calculate how many layers we have in this canopy
          ! This also checks the understory to see if its crown
          ! area is large enough to warrant a temporary sub-understory layer
          z = NumPotentialCanopyLayers(currentPatch,currentSite%spread,include_substory=.false.)

          do i_lyr = 1,z ! Loop around the currently occupied canopy layers.
             call DemoteFromLayer(currentSite, currentPatch, i_lyr, bc_in)
          end do

          ! After demotions, we may then again have cohorts that are very very
          ! very sparse, remove them
          call terminate_cohorts(currentSite, currentPatch, 1,13,bc_in)

          call fuse_cohorts(currentSite, currentPatch, bc_in)

          ! Remove cohorts for various other reasons
          call terminate_cohorts(currentSite, currentPatch, 2,13,bc_in)


          ! ---------------------------------------------------------------------------------------
          ! Promotion Phase: Identify if any upper-layers are underful and layers below them
          ! have cohorts that can be split and promoted to the layer above.
          ! ---------------------------------------------------------------------------------------

          ! Re-calculate Number of layers without the false substory
          z = NumPotentialCanopyLayers(currentPatch,currentSite%spread,include_substory=.false.)

          ! We only promote if we have at least two layers
          if (z>1) then

             do i_lyr=1,z-1
                call PromoteIntoLayer(currentSite, currentPatch, i_lyr)
             end do

             ! Remove cohorts that are incredibly sparse
             call terminate_cohorts(currentSite, currentPatch, 1,14,bc_in)

             call fuse_cohorts(currentSite, currentPatch, bc_in)

             ! Remove cohorts for various other reasons
             call terminate_cohorts(currentSite, currentPatch, 2,14,bc_in)

          end if

          ! ---------------------------------------------------------------------------------------
          ! Check on Layer Area (if the layer differences are not small
          ! Continue trying to demote/promote. Its possible on the first pass through,
          ! that cohort fusion has nudged the areas a little bit.
          ! ---------------------------------------------------------------------------------------

          z = NumPotentialCanopyLayers(currentPatch,currentSite%spread,include_substory=.false.)
          area_not_balanced = .false.
          do i_lyr = 1,z
             call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer(i_lyr))
             if( ((arealayer(i_lyr)-currentPatch%area)/currentPatch%area > area_check_rel_precision) .or. &
                  ((arealayer(i_lyr)-currentPatch%area) > area_check_precision )  ) then
                area_not_balanced = .true.
             endif
          enddo

          ! ---------------------------------------------------------------------------------------
          ! Gracefully exit if too many iterations have gone by
          ! ---------------------------------------------------------------------------------------

          patch_area_counter = patch_area_counter + 1
          if(patch_area_counter > max_patch_iterations .and. area_not_balanced) then
             write(fates_log(),*) 'PATCH AREA CHECK NOT CLOSING'
             write(fates_log(),*) 'patch area:',currentpatch%area
             do i_lyr = 1,z
                write(fates_log(),*) 'layer: ',i_lyr,' area: ',arealayer(i_lyr)
                write(fates_log(),*) 'rel error: ',(arealayer(i_lyr)-currentPatch%area)/currentPatch%area
                write(fates_log(),*) 'abs error: ',arealayer(i_lyr)-currentPatch%area
             enddo
             write(fates_log(),*) 'lat:',currentSite%lat
             write(fates_log(),*) 'lon:',currentSite%lon
             write(fates_log(),*) 'spread:',currentSite%spread
             currentCohort => currentPatch%tallest
             do while (associated(currentCohort))
                write(fates_log(),*) 'coh ilayer:',currentCohort%canopy_layer
                write(fates_log(),*) 'coh dbh:',currentCohort%dbh
                write(fates_log(),*) 'coh pft:',currentCohort%pft
                write(fates_log(),*) 'coh n:',currentCohort%n
                write(fates_log(),*) 'coh carea:',currentCohort%c_area
                ipft=currentCohort%pft
                write(fates_log(),*) 'maxh:',prt_params%allom_dbh_maxheight(ipft)
                write(fates_log(),*) 'lmode: ',prt_params%allom_lmode(ipft)
                write(fates_log(),*) 'd2bl2: ',prt_params%allom_d2bl2(ipft)
                write(fates_log(),*) 'd2bl_ediff: ',prt_params%allom_blca_expnt_diff(ipft)
                write(fates_log(),*) 'd2ca_min: ',prt_params%allom_d2ca_coefficient_min(ipft)
                write(fates_log(),*) 'd2ca_max: ',prt_params%allom_d2ca_coefficient_max(ipft)
                currentCohort => currentCohort%shorter
             enddo
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if

       enddo ! do while(area_not_balanced)


       ! Set current canopy layer occupancy indicator.
       currentPatch%NCL_p = min(nclmax,z)

       ! -------------------------------------------------------------------------------------------
       ! if we are using "strict PPA", then calculate a z_star value as
       ! the height of the smallest tree in the canopy
       ! loop from top to bottom and locate the shortest cohort in level 1 whose shorter
       ! neighbor is in level 2 set zstar as the ehight of that shortest level 1 cohort
       ! -------------------------------------------------------------------------------------------

       if ( ED_val_comp_excln .lt. 0.0_r8) then
          currentPatch%zstar = 0._r8
          currentCohort => currentPatch%tallest
          do while (associated(currentCohort))
             if(currentCohort%canopy_layer .eq. 2)then
                if (associated(currentCohort%taller)) then
                   if (currentCohort%taller%canopy_layer .eq. 1 ) then
                      currentPatch%zstar = currentCohort%taller%hite
                   endif
                endif
             endif
             currentCohort => currentCohort%shorter
          enddo
       endif
       
       currentPatch => currentPatch%younger
    enddo !patch

    return
  end subroutine canopy_structure


  ! ==============================================================================================


  subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr,bc_in)

    use EDParamsMod, only : ED_val_comp_excln
    use SFParamsMod, only : SF_val_CWD_frac

    ! !ARGUMENTS
    type(ed_site_type), intent(inout)  :: currentSite
    type(ed_patch_type), intent(inout) :: currentPatch
    integer, intent(in)                :: i_lyr   ! Current canopy layer of interest
    type(bc_in_type), intent(in)       :: bc_in

    ! !LOCAL VARIABLES:
    type(ed_cohort_type), pointer :: currentCohort
    type(ed_cohort_type), pointer :: copyc
    type(ed_cohort_type), pointer :: nextc  ! The next cohort in line
    integer  :: i_cwd                  ! Index for CWD pool
    real(r8) :: cc_loss                ! cohort crown area loss in demotion (m2)
    real(r8) :: leaf_c             ! leaf carbon [kg]
    real(r8) :: fnrt_c             ! fineroot carbon [kg]
    real(r8) :: sapw_c             ! sapwood carbon [kg]
    real(r8) :: store_c            ! storage carbon [kg]
    real(r8) :: struct_c           ! structure carbon [kg]
    real(r8) :: scale_factor       ! for prob. exclusion - scales weight to a fraction
    real(r8) :: scale_factor_min   ! "" minimum before exeedance of 1
    real(r8) :: scale_factor_res   ! "" applied to residual areas
    real(r8) :: area_res           ! residual area to demote after weakest cohort hits max
    real(r8) :: newarea
    real(r8) :: demote_area
    real(r8) :: sumweights
    real(r8) :: sumequal              ! for rank-ordered same-size cohorts
    ! this tallies their excluded area
    real(r8) :: arealayer              ! the area of the current canopy layer
    logical  :: tied_size_with_neighbors
    real(r8) :: total_crownarea_of_tied_cohorts

    ! First, determine how much total canopy area we have in this layer
    call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer)

    demote_area = arealayer - currentPatch%area

    if ( demote_area > area_target_precision ) then

       ! Is this layer currently over-occupied?
       ! In that case, we need to work out which cohorts to demote.
       ! We go in order from shortest to tallest for ranked demotion

       sumweights  = 0.0_r8
       currentCohort => currentPatch%shortest
       do while (associated(currentCohort))
          call carea_allom(currentCohort%dbh,currentCohort%n, &
               currentSite%spread,currentCohort%pft, &
               currentCohort%crowndamage, currentCohort%c_area)

          if(debug) then
             if(currentCohort%c_area<0._r8)then
                write(fates_log(),*) 'negative c_area stage 1d: ',currentCohort%dbh,i_lyr,currentCohort%n, &
                     currentSite%spread,currentCohort%pft,currentCohort%c_area
                call endrun(msg=errMsg(sourcefile, __LINE__))
             end if
          end if

          if( currentCohort%canopy_layer == i_lyr)then

             if (ED_val_comp_excln .ge. 0.0_r8 ) then

                ! ----------------------------------------------------------
                ! Stochastic method.
                ! Weight cohort demotion by inverse size to a constant power.
                ! In this hypothesis, it is assumed that even the tallest
                ! cohorts have a chance (although smaller) of being forced
                ! to the understory.
                ! ----------------------------------------------------------

                currentCohort%excl_weight = 1._r8 / (currentCohort%hite**ED_val_comp_excln)
                sumweights = sumweights + currentCohort%excl_weight

             else

                ! -----------------------------------------------------------
                ! Rank ordered deterministic method
                ! -----------------------------------------------------------
                ! If there are cohorts that have the exact same height (which is possible, really)
                ! we don't want to unilaterally promote/demote one before the others.
                ! So we <>mote them as a unit
                ! now we need to go through and figure out how many equal-size cohorts there are.
                ! then we need to go through, add up the collective crown areas of all equal-sized
                ! and equal-canopy-layer cohorts,
                ! and then demote from each as if they were a single group

                total_crownarea_of_tied_cohorts = currentCohort%c_area

                tied_size_with_neighbors = .false.
                nextc => currentCohort%taller
                do while (associated(nextc))
                   if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then
                      if( nextc%canopy_layer .eq. currentCohort%canopy_layer ) then
                         tied_size_with_neighbors = .true.
                         total_crownarea_of_tied_cohorts = &
                              total_crownarea_of_tied_cohorts + nextc%c_area
                      end if
                   else
                      exit
                   endif
                   nextc => nextc%taller
                end do

                if ( tied_size_with_neighbors ) then

                   currentCohort%excl_weight = &
                        max(0.0_r8,min(currentCohort%c_area, &
                        (currentCohort%c_area/total_crownarea_of_tied_cohorts) * &
                        (demote_area - sumweights) ))

                   sumequal = currentCohort%excl_weight

                   nextc => currentCohort%taller
                   do while (associated(nextc))
                      if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then
                         if (nextc%canopy_layer .eq. currentCohort%canopy_layer ) then
                            ! now we know the total crown area of all equal-sized,
                            ! equal-canopy-layer cohorts
                            nextc%excl_weight = &
                                 max(0.0_r8,min(nextc%c_area, &
                                 (nextc%c_area/total_crownarea_of_tied_cohorts) * &
                                 (demote_area - sumweights) ))
                            sumequal = sumequal + nextc%excl_weight
                         end if
                      else
                         exit
                      endif
                      nextc => nextc%taller
                   end do

                   ! Update the current cohort pointer to the last similar cohort
                   ! Its ok if this is not in the right layer
                   if(associated(nextc))then
                      currentCohort => nextc%shorter
                   else
                      currentCohort => currentPatch%tallest
                   end if
                   sumweights = sumweights + sumequal

                else
                   currentCohort%excl_weight = &
                        max(min(currentCohort%c_area, demote_area - sumweights ), 0._r8)
                   sumweights = sumweights + currentCohort%excl_weight
                end if

             endif
          endif
          currentCohort => currentCohort%taller
       enddo

       ! If this is probabalistic demotion, we need to do a round of normalization.
       ! And then a few rounds where we pre-calculate the demotion areas
       ! and adjust things if the demoted area wants to be greater than
       ! what is available. The math is too hard to explain here, see
       ! the tech note section on promotion/demotion.

       if (ED_val_comp_excln .ge. 0.0_r8 ) then

          scale_factor_min  = 1.e10_r8
          scale_factor      = 0._r8
          currentCohort => currentPatch%tallest
          do while (associated(currentCohort))

             if(currentCohort%canopy_layer  ==  i_lyr) then

                currentCohort%excl_weight = currentCohort%excl_weight/sumweights
                if( 1._r8/currentCohort%excl_weight <  scale_factor_min )  &
                     scale_factor_min = 1._r8/currentCohort%excl_weight

                scale_factor = scale_factor + currentCohort%excl_weight * currentCohort%c_area

             endif
             currentCohort => currentCohort%shorter
          enddo

          ! This is the factor by which we need to multiply
          ! the demotion probabilities, so the sum result equals
          ! the total amount to demote

          scale_factor = demote_area/scale_factor

          if(scale_factor <= scale_factor_min) then

             ! Trivial case, all of the demotion fractions are less than 1.

             currentCohort => currentPatch%tallest
             do while (associated(currentCohort))
                if(currentCohort%canopy_layer  ==  i_lyr) then
                   currentCohort%excl_weight = currentCohort%c_area * currentCohort%excl_weight * scale_factor

                   if(debug) then
                      if((currentCohort%excl_weight > (currentCohort%c_area+area_target_precision)) .or. &
                           (currentCohort%excl_weight < 0._r8)  ) then
                         write(fates_log(),*) 'exclusion area too big (1)'
                         write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area
                         write(fates_log(),*) 'dbh: ',currentCohort%dbh
                         write(fates_log(),*) 'n: ',currentCohort%n
                         write(fates_log(),*) 'spread: ',currentSite%spread
                         write(fates_log(),*) 'pft: ',currentCohort%pft
                         write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight
                         write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area
                         call endrun(msg=errMsg(sourcefile, __LINE__))
                      end if
                   end if

                endif
                currentCohort => currentCohort%shorter
             enddo

          else


             ! Non-trivial case, at least 1 cohort's demotion
             ! rate would exceed its area, given the trivial scale factor

             area_res         = 0._r8
             scale_factor_res = 0._r8
             currentCohort => currentPatch%tallest
             do while (associated(currentCohort))
                if(currentCohort%canopy_layer  ==  i_lyr) then
                   area_res         = area_res + &
                        currentCohort%c_area * currentCohort%excl_weight * &
                        scale_factor_min
                   scale_factor_res = scale_factor_res + &
                        currentCohort%c_area * &
                        (1._r8 - (currentCohort%excl_weight * scale_factor_min))
                endif
                currentCohort => currentCohort%shorter
             enddo

             area_res = demote_area - area_res

             scale_factor_res = area_res / scale_factor_res

             currentCohort => currentPatch%tallest
             do while (associated(currentCohort))
                if(currentCohort%canopy_layer  ==  i_lyr) then

                   currentCohort%excl_weight = currentCohort%c_area * &
                        (currentCohort%excl_weight * scale_factor_min + &
                        (1._r8 - (currentCohort%excl_weight*scale_factor_min) ) * scale_factor_res)

                   if(debug)then
                      if((currentCohort%excl_weight > &
                           (currentCohort%c_area+area_target_precision)) .or. &
                           (currentCohort%excl_weight < 0._r8)  ) then
                         write(fates_log(),*) 'exclusion area error (2)'
                         write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area
                         write(fates_log(),*) 'currentCohort%excl_weight: ', &
                              currentCohort%excl_weight
                         write(fates_log(),*) 'excess: ', &
                              currentCohort%excl_weight - currentCohort%c_area
                         call endrun(msg=errMsg(sourcefile, __LINE__))
                      end if
                   end if

                endif
                currentCohort => currentCohort%shorter
             enddo

          end if

       end if


       ! perform a check and see if the demotions meet the demand
       sumweights = 0._r8
       currentCohort => currentPatch%tallest
       do while (associated(currentCohort))
          if(currentCohort%canopy_layer  ==  i_lyr) then
             sumweights = sumweights + currentCohort%excl_weight
          end if
          currentCohort => currentCohort%shorter
       end do

       if (abs(sumweights - demote_area) > area_check_precision ) then
          write(fates_log(),*) 'demotions dont add up'
          write(fates_log(),*) 'sum demotions: ',sumweights
          write(fates_log(),*) 'area needed to be demoted: ',demote_area
          write(fates_log(),*) 'excess: ',sumweights - demote_area
          call endrun(msg=errMsg(sourcefile, __LINE__))
       end if


       ! Weights have been calculated. Now move them to the lower layer

       currentCohort => currentPatch%tallest
       do while (associated(currentCohort))

          nextc => currentCohort%shorter

          if(currentCohort%canopy_layer == i_lyr )then

             cc_loss         = currentCohort%excl_weight
             leaf_c          = currentCohort%prt%GetState(leaf_organ,carbon12_element)
             store_c         = currentCohort%prt%GetState(store_organ,carbon12_element)
             fnrt_c          = currentCohort%prt%GetState(fnrt_organ,carbon12_element)
             sapw_c          = currentCohort%prt%GetState(sapw_organ,carbon12_element)
             struct_c        = currentCohort%prt%GetState(struct_organ,carbon12_element)

             if ( (cc_loss-currentCohort%c_area) > -nearzero .and. &
                  (cc_loss-currentCohort%c_area) < area_target_precision ) then

                ! If the whole cohort is being demoted, just change its
                ! layer index

                currentCohort%canopy_layer = i_lyr+1

                ! keep track of number and biomass of demoted cohort
                currentSite%demotion_rate(currentCohort%size_class) = &
                     currentSite%demotion_rate(currentCohort%size_class) + currentCohort%n
                currentSite%demotion_carbonflux = currentSite%demotion_carbonflux + &
                     (leaf_c + store_c + fnrt_c + sapw_c + struct_c) * currentCohort%n

             elseif( (cc_loss < currentCohort%c_area) .and. &
                  (cc_loss > area_target_precision) ) then

                ! If only part of the cohort is demoted
                ! then it must be split (little more complicated)

                ! Make a copy of the current cohort.  The copy and the original
                ! conserve total number density of the original.  The copy
                ! remains in the upper-story.  The original is the one
                ! demoted to the understory


                allocate(copyc)

                ! (keep as an example)
                ! Initialize running means
                !allocate(copyc%tveg_lpa)
                !!allocate(copyc%l2fr_ema)
                !  Note, no need to give a starter value here,
                !  that will be taken care of in copy_cohort()
                !!call copyc%l2fr_ema%InitRMean(ema_60day)
                     
                ! Initialize the PARTEH object and point to the
                ! correct boundary condition fields
                copyc%prt => null()
                call InitPRTObject(copyc%prt)

                if( hlm_use_planthydro.eq.itrue ) then
                   call InitHydrCohort(currentSite,copyc)
                endif

                call copy_cohort(currentCohort, copyc)
                call InitPRTBoundaryConditions(copyc)
                
                newarea = currentCohort%c_area - cc_loss
                copyc%n = currentCohort%n*newarea/currentCohort%c_area
                currentCohort%n = currentCohort%n - copyc%n

                copyc%canopy_layer = i_lyr !the taller cohort is the copy

                ! Demote the current cohort to the understory.
                currentCohort%canopy_layer = i_lyr + 1

                ! keep track of number and biomass of demoted cohort
                currentSite%demotion_rate(currentCohort%size_class) = &
                     currentSite%demotion_rate(currentCohort%size_class) + currentCohort%n
                currentSite%demotion_carbonflux = currentSite%demotion_carbonflux + &
                     (leaf_c + store_c + fnrt_c + sapw_c + struct_c) * currentCohort%n

                call carea_allom(copyc%dbh,copyc%n,currentSite%spread,copyc%pft, &
                     copyc%crowndamage, copyc%c_area)
                call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, &
                     currentCohort%pft,currentCohort%crowndamage, currentCohort%c_area)

                !----------- Insert copy into linked list ------------------------!
                copyc%shorter => currentCohort
                if(associated(currentCohort%taller))then
                   copyc%taller => currentCohort%taller
                   currentCohort%taller%shorter => copyc
                else
                   currentPatch%tallest => copyc
                   copyc%taller => null()
                endif
                currentCohort%taller => copyc

             elseif(cc_loss > currentCohort%c_area)then

                write(fates_log(),*) 'more area than the cohort has is being demoted'
                write(fates_log(),*) 'loss:',cc_loss
                write(fates_log(),*) 'existing area:',currentCohort%c_area
                write(fates_log(),*) 'excess: ',cc_loss - currentCohort%c_area
                call endrun(msg=errMsg(sourcefile, __LINE__))

             end if

             ! kill the ones which go into canopy layers that are not allowed
             ! USE THIS OVERRIDE IF YOU ARE FORCING A ONE COHORT SIMULATION
             ! (also make sure to turn off germination, external seed rain,
             ! (use only one PFT, and make sure disturb_frac is 0)
             ! (RGK-0822)
             !if(currentCohort%canopy_layer>1) then

             if(currentCohort%canopy_layer>nclmax )then
                ! put the litter from the terminated cohorts
                ! straight into the fragmenting pools
                call terminate_cohort(currentSite,currentPatch,currentCohort,bc_in)
                deallocate(currentCohort, stat=istat, errmsg=smsg)
                if (istat/=0) then
                   write(fates_log(),*) 'dealloc012: fail on deallocate(currentCohort):'//trim(smsg)
                   call endrun(msg=errMsg(sourcefile, __LINE__))
                endif
             else
             call carea_allom(currentCohort%dbh,currentCohort%n, &
                  currentSite%spread,currentCohort%pft,currentCohort%crowndamage, &
                  currentCohort%c_area)
             end if

          endif !canopy layer = i_ly

          ! We dont use our typical (point to smaller)
          ! here, because, we may had deallocated the existing
          ! currentCohort

          currentCohort => nextc
       enddo !currentCohort


       ! Update the area calculations of the current layer
       ! And the layer below that may or may not had recieved
       ! Demotions

       call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer)

       if ( (abs(arealayer - currentPatch%area)/arealayer > area_check_rel_precision ) .or. &
            (abs(arealayer - currentPatch%area) > area_check_precision) ) then
          write(fates_log(),*) 'demotion did not trim area within tolerance'
          write(fates_log(),*) 'arealayer:',arealayer
          write(fates_log(),*) 'patch%area:',currentPatch%area
          write(fates_log(),*) 'ilayer: ',i_lyr
          write(fates_log(),*) 'bias:',arealayer - currentPatch%area
          write(fates_log(),*) 'rel bias:',(arealayer - currentPatch%area)/arealayer
          write(fates_log(),*) 'demote_area:',demote_area
          call endrun(msg=errMsg(sourcefile, __LINE__))
       end if


    end if

    return
  end subroutine DemoteFromLayer

  ! ==============================================================================================

  subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr)

    ! -------------------------------------------------------------------------------------------
    ! Check whether the intended 'full' layers are actually filling all the space.
    ! If not, promote some fraction of cohorts upwards.
    ! THIS SECTION MIGHT BE TRIGGERED BY A FIRE OR MORTALITY EVENT, FOLLOWED BY A PATCH FUSION,
    ! SO THE TOP LAYER IS NO LONGER FULL.
    ! -------------------------------------------------------------------------------------------

    use EDParamsMod, only : ED_val_comp_excln

    ! !ARGUMENTS
    type(ed_site_type), intent(inout), target  :: currentSite
    type(ed_patch_type), intent(inout), target :: currentPatch
    integer, intent(in)                        :: i_lyr   ! Current canopy layer of interest

    ! !LOCAL VARIABLES:
    type(ed_cohort_type), pointer :: currentCohort
    type(ed_cohort_type), pointer :: copyc
    type(ed_cohort_type), pointer :: nextc   ! the next cohort, or used for looping
    ! cohorts against the current

    real(r8) :: scale_factor       ! for prob. exclusion - scales weight to a fraction
    real(r8) :: scale_factor_min   ! "" minimum before exeedance of 1
    real(r8) :: scale_factor_res   ! "" applied to residual areas
    real(r8) :: area_res           ! residual area to demote after weakest cohort hits max
    real(r8) :: promote_area
    real(r8) :: newarea
    real(r8) :: sumweights
    real(r8) :: sumequal              ! for tied cohorts, the sum of weights in
    ! their group
    real(r8) :: cc_gain                ! cohort crown area gain in promotion (m2)
    real(r8) :: arealayer_current      ! area (m2) of the current canopy layer
    real(r8) :: arealayer_below        ! area (m2) of the layer below the current layer
    real(r8) :: leaf_c             ! leaf carbon [kg]
    real(r8) :: fnrt_c             ! fineroot carbon [kg]
    real(r8) :: sapw_c             ! sapwood carbon [kg]
    real(r8) :: store_c            ! storage carbon [kg]
    real(r8) :: struct_c           ! structure carbon [kg]

    logical  :: tied_size_with_neighbors
    real(r8) :: total_crownarea_of_tied_cohorts

    call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current)
    call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr+1,arealayer_below)


    ! how much do we need to gain?
    promote_area    =  currentPatch%area - arealayer_current

    if( promote_area > area_target_precision ) then

       if(arealayer_below <= promote_area ) then

          ! ---------------------------------------------------------------------------
          ! Promote all cohorts from layer below if that whole layer has area smaller
          ! than the tolerance on the gains needed into current layer
          ! ---------------------------------------------------------------------------

          currentCohort => currentPatch%tallest
          do while (associated(currentCohort))
             !look at the cohorts in the canopy layer below...
             if(currentCohort%canopy_layer == i_lyr+1)then

                leaf_c          = currentCohort%prt%GetState(leaf_organ,carbon12_element)
                store_c         = currentCohort%prt%GetState(store_organ,carbon12_element)
                fnrt_c          = currentCohort%prt%GetState(fnrt_organ,carbon12_element)
                sapw_c          = currentCohort%prt%GetState(sapw_organ,carbon12_element)
                struct_c        = currentCohort%prt%GetState(struct_organ,carbon12_element)

                currentCohort%canopy_layer = i_lyr
                call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, &
                     currentCohort%pft,currentCohort%crowndamage, currentCohort%c_area)
                ! keep track of number and biomass of promoted cohort
                currentSite%promotion_rate(currentCohort%size_class) = &
                     currentSite%promotion_rate(currentCohort%size_class) + currentCohort%n
                currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + &
                     (leaf_c + fnrt_c + store_c + sapw_c + struct_c) * currentCohort%n

             endif
             currentCohort => currentCohort%shorter
          enddo

       else

          ! ---------------------------------------------------------------------------
          ! This is the non-trivial case where the lower layer can accomodate
          ! more than what is necessary.
          ! ---------------------------------------------------------------------------


          ! figure out with what weighting we need to promote cohorts.
          ! This is the opposite of the demotion weighting...

          sumweights = 0.0_r8
          currentCohort => currentPatch%tallest
          do while (associated(currentCohort))
             call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, &
                  currentCohort%pft,currentCohort%crowndamage,currentCohort%c_area)
             if(currentCohort%canopy_layer == i_lyr+1)then !look at the cohorts in the canopy layer below...

                if (ED_val_comp_excln .ge. 0.0_r8 ) then

                   ! ------------------------------------------------------------------
                   ! Stochastic case, as above (in demotion portion of code)
                   ! ------------------------------------------------------------------

                   currentCohort%prom_weight = currentCohort%hite**ED_val_comp_excln
                   sumweights = sumweights + currentCohort%prom_weight
                else

                   ! ------------------------------------------------------------------
                   ! Rank ordered deterministic method
                   ! If there are cohorts that have the exact same height (which is possible, really)
                   ! we don't want to unilaterally promote/demote one before the others.
                   ! So we <>mote them as a unit
                   ! now we need to go through and figure out how many equal-size cohorts there are.
                   ! then we need to go through, add up the collective crown areas of all equal-sized
                   ! and equal-canopy-layer cohorts,
                   ! and then demote from each as if they were a single group
                   ! ------------------------------------------------------------------

                   total_crownarea_of_tied_cohorts = currentCohort%c_area
                   tied_size_with_neighbors = .false.
                   nextc => currentCohort%shorter
                   do while (associated(nextc))
                      if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then
                         if( nextc%canopy_layer .eq. currentCohort%canopy_layer ) then
                            tied_size_with_neighbors = .true.
                            total_crownarea_of_tied_cohorts = &
                                 total_crownarea_of_tied_cohorts + nextc%c_area
                         end if
                      else
                         exit
                      endif
                      nextc => nextc%shorter
                   end do

                   if ( tied_size_with_neighbors ) then

                      currentCohort%prom_weight = &
                           max(0.0_r8,min(currentCohort%c_area, &
                           (currentCohort%c_area/total_crownarea_of_tied_cohorts) * &
                           (promote_area - sumweights) ))
                      sumequal = currentCohort%prom_weight

                      nextc => currentCohort%shorter
                      do while (associated(nextc))
                         if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then
                            if (nextc%canopy_layer .eq. currentCohort%canopy_layer ) then
                               ! now we know the total crown area of all equal-sized,
                               ! equal-canopy-layer cohorts
                               nextc%prom_weight = &
                                    max(0.0_r8,min(nextc%c_area, &
                                    (nextc%c_area/total_crownarea_of_tied_cohorts) * &
                                    (promote_area - sumweights) ))
                               sumequal = sumequal + nextc%prom_weight
                            end if
                         else
                            exit
                         endif
                         nextc => nextc%shorter
                      end do

                      ! Update the current cohort pointer to the last similar cohort
                      ! Its ok if this is not in the right layer
                      if(associated(nextc))then
                         currentCohort => nextc%taller
                      else
                         currentCohort => currentPatch%shortest
                      end if
                      sumweights = sumweights + sumequal

                   else
                      currentCohort%prom_weight = &
                           max(min(currentCohort%c_area, promote_area - sumweights ), 0._r8)
                      sumweights = sumweights + currentCohort%prom_weight

                   end if

                endif
             endif
             currentCohort => currentCohort%shorter
          enddo !currentCohort


          ! If this is probabalistic promotion, we need to do a round of normalization.
          ! And then a few rounds where we pre-calculate the promotion areas
          ! and adjust things if the promoted area wants to be greater than
          ! what is available.

          if (ED_val_comp_excln .ge. 0.0_r8 ) then

             scale_factor_min  = 1.e10_r8
             scale_factor      = 0._r8
             currentCohort => currentPatch%tallest
             do while (associated(currentCohort))

                if(currentCohort%canopy_layer  ==  (i_lyr+1) ) then

                   currentCohort%prom_weight = currentCohort%prom_weight/sumweights
                   if( 1._r8/currentCohort%prom_weight <  scale_factor_min )  &
                        scale_factor_min = 1._r8/currentCohort%prom_weight

                   scale_factor = scale_factor + currentCohort%prom_weight * currentCohort%c_area

                endif
                currentCohort => currentCohort%shorter
             enddo

             ! This is the factor by which we need to multiply
             ! the demotion probabilities, so the sum result equals
             ! the total amount to demote
             scale_factor = promote_area/scale_factor


             if(scale_factor <= scale_factor_min) then

                ! Trivial case, all of the demotion fractions
                ! are less than 1.

                currentCohort => currentPatch%tallest
                do while (associated(currentCohort))
                   if(currentCohort%canopy_layer  ==  (i_lyr+1) ) then
                      currentCohort%prom_weight = currentCohort%c_area * &
                           currentCohort%prom_weight * scale_factor

                      if(debug)then
                         if((currentCohort%prom_weight > &
                              (currentCohort%c_area+area_target_precision)) .or. &
                              (currentCohort%prom_weight < 0._r8)  ) then
                            write(fates_log(),*) 'promotion area too big (1)'
                            write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area
                            write(fates_log(),*) 'currentCohort%prom_weight: ', &
                                 currentCohort%prom_weight
                            write(fates_log(),*) 'excess: ', &
                                 currentCohort%prom_weight - currentCohort%c_area
                            call endrun(msg=errMsg(sourcefile, __LINE__))
                         end if
                      end if

                   endif
                   currentCohort => currentCohort%shorter
                enddo

             else

                ! Non-trivial case, at least 1 cohort's promotion
                ! rate would exceed its area, given the trivial scale factor

                area_res         = 0._r8
                scale_factor_res = 0._r8
                currentCohort => currentPatch%tallest
                do while (associated(currentCohort))
                   if(currentCohort%canopy_layer  ==  (i_lyr+1) ) then
                      area_res         = area_res + &
                           currentCohort%c_area*currentCohort%prom_weight*scale_factor_min
                      scale_factor_res = scale_factor_res + &
                           currentCohort%c_area * &
                           (1._r8 - (currentCohort%prom_weight * scale_factor_min))
                   endif
                   currentCohort => currentCohort%shorter
                enddo

                area_res = promote_area - area_res

                scale_factor_res = area_res / scale_factor_res

                currentCohort => currentPatch%tallest
                do while (associated(currentCohort))
                   if(currentCohort%canopy_layer  ==  (i_lyr+1)) then

                      currentCohort%prom_weight = currentCohort%c_area * &
                           (currentCohort%prom_weight * scale_factor_min + &
                           (1._r8 - (currentCohort%prom_weight*scale_factor_min) ) * &
                           scale_factor_res)

                      if(debug)then
                         if((currentCohort%prom_weight > &
                              (currentCohort%c_area+area_target_precision)) .or. &
                              (currentCohort%prom_weight < 0._r8)  ) then
                            write(fates_log(),*) 'promotion area error (2)'
                            write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area
                            write(fates_log(),*) 'currentCohort%prom_weight: ', &
                                 currentCohort%prom_weight
                            write(fates_log(),*) 'excess: ', &
                                 currentCohort%prom_weight - currentCohort%c_area
                            call endrun(msg=errMsg(sourcefile, __LINE__))
                         end if
                      end if

                   endif
                   currentCohort => currentCohort%shorter
                enddo

             end if

          end if


          ! lets perform a check and see if the promotions meet the demand
          sumweights = 0._r8
          currentCohort => currentPatch%tallest
          do while (associated(currentCohort))
             if(currentCohort%canopy_layer  ==  (i_lyr+1)) then
                sumweights = sumweights + currentCohort%prom_weight
             end if
             currentCohort => currentCohort%shorter
          end do

          if(debug)then
             if (abs(sumweights - promote_area) > area_check_precision ) then
                write(fates_log(),*) 'promotions dont add up'
                write(fates_log(),*) 'sum promotions: ',sumweights
                write(fates_log(),*) 'area needed to be promoted: ',promote_area
                write(fates_log(),*) 'excess: ',sumweights - promote_area
                call endrun(msg=errMsg(sourcefile, __LINE__))
             end if
          end if

          currentCohort => currentPatch%tallest
          do while (associated(currentCohort))


             !All the trees in this layer need to promote some area upwards...
             if( (currentCohort%canopy_layer == i_lyr+1) ) then

                cc_gain         = currentCohort%prom_weight
                leaf_c          = currentCohort%prt%GetState(leaf_organ,carbon12_element)
                store_c         = currentCohort%prt%GetState(store_organ,carbon12_element)
                fnrt_c          = currentCohort%prt%GetState(fnrt_organ,carbon12_element)
                sapw_c          = currentCohort%prt%GetState(sapw_organ,carbon12_element)
                struct_c        = currentCohort%prt%GetState(struct_organ,carbon12_element)

                if ( (cc_gain-currentCohort%c_area) > -nearzero .and. &
                     (cc_gain-currentCohort%c_area) < area_target_precision ) then

                   currentCohort%canopy_layer = i_lyr

                   ! keep track of number and biomass of promoted cohort
                   currentSite%promotion_rate(currentCohort%size_class) = &
                        currentSite%promotion_rate(currentCohort%size_class) + currentCohort%n

                   currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + &
                        (leaf_c + fnrt_c + store_c + sapw_c + struct_c) * currentCohort%n

                elseif ( (cc_gain < currentCohort%c_area) .and. &
                     (cc_gain > area_target_precision) ) then

                   allocate(copyc)


                   !!allocate(copyc%l2fr_ema)
                   ! Note, no need to give a starter value here,
                   ! that will be taken care of in copy_cohort()
                   !!call copyc%l2fr_ema%InitRMean(ema_60day)
                   
                   ! Initialize the PARTEH object and point to the
                   ! correct boundary condition fields
                   copyc%prt => null()
                   call InitPRTObject(copyc%prt)
                   

                   if( hlm_use_planthydro.eq.itrue ) then
                      call InitHydrCohort(CurrentSite,copyc)
                   endif

                   ! (keep as an example)
                   ! Initialize running means
                   !allocate(copyc%tveg_lpa)
                   !call copyc%tveg_lpa%InitRMean(ema_lpa,&
                   !     init_value=currentPatch%tveg_lpa%GetMean())
                   
                   call copy_cohort(currentCohort, copyc) !makes an identical copy...
                   call InitPRTBoundaryConditions(copyc)
                   
                   newarea = currentCohort%c_area - cc_gain !new area of existing cohort

                   call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, &
                        currentCohort%pft,currentCohort%crowndamage, currentCohort%c_area)

                   ! number of individuals in promoted cohort.
                   copyc%n = currentCohort%n*cc_gain/currentCohort%c_area

                   ! number of individuals in cohort remaining in understorey
                   currentCohort%n = currentCohort%n - copyc%n

                   currentCohort%canopy_layer = i_lyr + 1 ! keep current cohort in the understory.
                   copyc%canopy_layer = i_lyr             ! promote copy to the higher canopy layer.

                   ! keep track of number and biomass of promoted cohort
                   currentSite%promotion_rate(copyc%size_class) = &
                        currentSite%promotion_rate(copyc%size_class) + copyc%n

                   currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + &
                        (leaf_c + fnrt_c + store_c + sapw_c + struct_c) * copyc%n

                   call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, &
                        currentCohort%pft,currentCohort%crowndamage, currentCohort%c_area)
                   call carea_allom(copyc%dbh,copyc%n,currentSite%spread,copyc%pft,&
                        copyc%crowndamage,copyc%c_area)

                   !----------- Insert copy into linked list ------------------------!
                   copyc%shorter => currentCohort
                   if(associated(currentCohort%taller))then
                      copyc%taller => currentCohort%taller
                      currentCohort%taller%shorter => copyc
                   else
                      currentPatch%tallest => copyc
                      copyc%taller => null()
                   endif
                   currentCohort%taller => copyc

                elseif(cc_gain > currentCohort%c_area)then

                   write(fates_log(),*) 'more area than the cohort has is being promoted'
                   write(fates_log(),*) 'loss:',cc_gain
                   write(fates_log(),*) 'existing area:',currentCohort%c_area
                   call endrun(msg=errMsg(sourcefile, __LINE__))

                endif

             endif   ! if(currentCohort%canopy_layer == i_lyr+1) then
             currentCohort => currentCohort%shorter
          enddo !currentCohort

          call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current)

          if ((abs(arealayer_current - currentPatch%area)/arealayer_current > &
               area_check_rel_precision ) .or. &
               (abs(arealayer_current - currentPatch%area) > area_check_precision) ) then
             write(fates_log(),*) 'promotion did not bring area within tolerance'
             write(fates_log(),*) 'arealayer:',arealayer_current
             write(fates_log(),*) 'patch%area:',currentPatch%area
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if

       end if

    end if

    return
  end subroutine PromoteIntoLayer

  ! ============================================================================

  subroutine canopy_spread( currentSite )
    !
    ! !DESCRIPTION:
    !  Calculates the spatial spread of tree canopies based on canopy closure.
    !
    ! !USES:
    use EDTypesMod        , only : AREA
    use EDParamsMod, only : ED_val_canopy_closure_thresh
    !
    ! !ARGUMENTS
    type (ed_site_type), intent(inout), target :: currentSite
    !
    ! !LOCAL VARIABLES:
    type (ed_cohort_type), pointer :: currentCohort
    type (ed_patch_type) , pointer :: currentPatch
    real(r8) :: sitelevel_canopyarea  ! Amount of canopy in top layer at the site level
    real(r8) :: inc                   ! Arbitrary daily incremental change in canopy area
    integer  :: z
    !----------------------------------------------------------------------

    inc = 0.05_r8

    currentPatch => currentSite%oldest_patch

    sitelevel_canopyarea = 0.0_r8
    do while (associated(currentPatch))

       !calculate canopy area in each patch...
       currentCohort => currentPatch%tallest
       do while (associated(currentCohort))
          call carea_allom(currentCohort%dbh,currentCohort%n, &
               currentSite%spread,currentCohort%pft,currentCohort%crowndamage,currentCohort%c_area)
          if( ( prt_params%woody(currentCohort%pft) .eq. itrue ) .and. &
               (currentCohort%canopy_layer .eq. 1 ) ) then
             sitelevel_canopyarea = sitelevel_canopyarea + currentCohort%c_area
          endif
          currentCohort => currentCohort%shorter
       enddo

       currentPatch => currentPatch%younger

    enddo !currentPatch

    ! If the canopy area is approaching closure,
    ! squash the tree canopies and make them taller and thinner
    if( sitelevel_canopyarea/AREA .gt. ED_val_canopy_closure_thresh ) then
       currentSite%spread = currentSite%spread - inc
    else
       currentSite%spread = currentSite%spread + inc
    endif

    ! put within bounds to make sure it stays between 0 and 1
    currentSite%spread = max(min(currentSite%spread, 1._r8), 0._r8)

  end subroutine canopy_spread


  ! =====================================================================================

  subroutine canopy_summarization( nsites, sites, bc_in )

    ! ----------------------------------------------------------------------------------
    ! Much of this routine was once ed_clm_link minus all the IO and history stuff
    ! ---------------------------------------------------------------------------------

    use FatesInterfaceTypesMod    , only : hlm_use_cohort_age_tracking
    use EDPatchDynamicsMod   , only : set_patchno
    use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index
    use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index
    use EDtypesMod           , only : area
    use FatesConstantsMod    , only : itrue

    ! !ARGUMENTS
    integer                 , intent(in)            :: nsites
    type(ed_site_type)      , intent(inout), target :: sites(nsites)
    type(bc_in_type)        , intent(in)            :: bc_in(nsites)
    !
    ! !LOCAL VARIABLES:
    type (ed_patch_type)  , pointer :: currentPatch
    type (ed_cohort_type) , pointer :: currentCohort
    integer  :: s
    integer  :: ft               ! plant functional type
    integer  :: ifp              ! the number of the vegetated patch (1,2,3). In SP mode bareground patch is 0
    integer  :: patchn           ! identification number for each patch.
    real(r8) :: leaf_c           ! leaf carbon [kg]
    real(r8) :: fnrt_c           ! fineroot carbon [kg]
    real(r8) :: sapw_c           ! sapwood carbon [kg]
    real(r8) :: store_c          ! storage carbon [kg]
    real(r8) :: struct_c         ! structure carbon [kg]

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

    if ( debug ) then
       write(fates_log(),*) 'in canopy_summarization'
    endif

    do s = 1,nsites

       ! --------------------------------------------------------------------------------
       ! Set the patch indices (this is usefull mostly for communicating with a host or
       ! driving model.  Loops through all patches and sets cpatch%patchno to the integer
       ! order of oldest to youngest where the oldest is 1.
       ! --------------------------------------------------------------------------------
       call set_patchno( sites(s) )

       currentPatch => sites(s)%oldest_patch

       do while(associated(currentPatch))

          !zero cohort-summed variables.
          currentPatch%total_canopy_area = 0.0_r8
          currentPatch%total_tree_area = 0.0_r8

          !update cohort quantitie s
          currentCohort => currentPatch%shortest
          do while(associated(currentCohort))

             ft = currentCohort%pft
             leaf_c   = currentCohort%prt%GetState(leaf_organ, carbon12_element)
             sapw_c   = currentCohort%prt%GetState(sapw_organ, carbon12_element)
             struct_c = currentCohort%prt%GetState(struct_organ, carbon12_element)
             fnrt_c   = currentCohort%prt%GetState(fnrt_organ, carbon12_element)
             store_c  = currentCohort%prt%GetState(store_organ, carbon12_element)

             ! Update the cohort's index within the size bin classes
             ! Update the cohort's index within the SCPF classification system
             call sizetype_class_index(currentCohort%dbh,currentCohort%pft, &
                  currentCohort%size_class,currentCohort%size_by_pft_class)

             if (hlm_use_cohort_age_tracking .eq. itrue) then
                call coagetype_class_index(currentCohort%coage,currentCohort%pft, &
                     currentCohort%coage_class,currentCohort%coage_by_pft_class)
             end if

             if(hlm_use_sp.eq.ifalse)then
                call carea_allom(currentCohort%dbh,currentCohort%n,sites(s)%spread,&
                     currentCohort%pft,currentCohort%crowndamage, currentCohort%c_area)
             endif

             if(currentCohort%canopy_layer==1)then
                currentPatch%total_canopy_area = currentPatch%total_canopy_area + currentCohort%c_area
                if( prt_params%woody(ft) == itrue)then
                   currentPatch%total_tree_area = currentPatch%total_tree_area + currentCohort%c_area
                endif
             endif

             ! adding checks for SP and NOCOMP modes.
             if(currentPatch%nocomp_pft_label.eq.nocomp_bareground)then
                write(fates_log(),*) 'cohorts in barepatch',currentPatch%total_canopy_area,currentPatch%nocomp_pft_label
                call endrun(msg=errMsg(sourcefile, __LINE__))
             end if

             if(hlm_use_sp.eq.itrue)then

                if(associated(currentPatch%tallest%shorter))then
                   write(fates_log(),*) 'more than one cohort in SP mode',s,currentPatch%nocomp_pft_label
                   call endrun(msg=errMsg(sourcefile, __LINE__))
                end if

                if(currentPatch%total_canopy_area-currentPatch%area.gt.1.0e-16)then
                   write(fates_log(),*) 'too much canopy in summary',s, &
                        currentPatch%nocomp_pft_label, currentPatch%total_canopy_area-currentPatch%area
                   call endrun(msg=errMsg(sourcefile, __LINE__))
                end if
             end if  !sp mode

             ! Check for erroneous zero values.
             if(currentCohort%dbh <= 0._r8 .or. currentCohort%n == 0._r8)then
                write(fates_log(),*) 'FATES: dbh or n is zero in canopy_summarization', &
                     currentCohort%dbh,currentCohort%n
                call endrun(msg=errMsg(sourcefile, __LINE__))
             endif

             if(currentCohort%pft == 0.or.currentCohort%canopy_trim <= 0._r8)then
                write(fates_log(),*) 'FATES: PFT or trim is zero in canopy_summarization', &
                     currentCohort%pft,currentCohort%canopy_trim
                call endrun(msg=errMsg(sourcefile, __LINE__))
             endif
             if( (sapw_c + leaf_c + fnrt_c) <= 0._r8)then
                write(fates_log(),*) 'FATES: alive biomass is zero in canopy_summarization', &
                     sapw_c + leaf_c + fnrt_c
                call endrun(msg=errMsg(sourcefile, __LINE__))
             endif

             currentCohort => currentCohort%taller

          enddo ! ends 'do while(associated(currentCohort))

          if ( currentPatch%total_canopy_area>currentPatch%area ) then
             if ( currentPatch%total_canopy_area-currentPatch%area > 0.001_r8 ) then
                write(fates_log(),*) 'FATES: canopy area bigger than area', &
                     currentPatch%total_canopy_area ,currentPatch%area, &
                     currentPatch%total_canopy_area -currentPatch%area,&
                     currentPatch%nocomp_pft_label
                call endrun(msg=errMsg(sourcefile, __LINE__))
             end if
             currentPatch%total_canopy_area = currentPatch%area
          endif

          currentPatch => currentPatch%younger
       end do !patch loop


    
       call leaf_area_profile(sites(s))

    end do ! site loop

    return
  end subroutine canopy_summarization

  ! ====================================================================================

  subroutine UpdateFatesAvgSnowDepth(sites,bc_in)

    ! This routine updates the snow depth used in FATES to occlude vegetation
    ! Currently this average takes into account the depth of snow and the
    ! areal coverage fraction

    type(ed_site_type)      , intent(inout), target :: sites(:)
    type(bc_in_type)        , intent(in)            :: bc_in(:)

    integer  :: s

    do s = 1, size(sites,dim=1)
       sites(s)%snow_depth = bc_in(s)%snow_depth_si * bc_in(s)%frac_sno_eff_si
    end do

    return
  end subroutine UpdateFatesAvgSnowDepth


  ! =====================================================================================

  subroutine leaf_area_profile( currentSite )

    ! -----------------------------------------------------------------------------------
    ! This subroutine calculates how leaf and stem areas are distributed
    ! in vertical and horizontal space.
    !
    ! The following cohort level diagnostics are updated here:
    !
    ! currentCohort%treelai    ! LAI per unit crown area  (m2/m2)
    ! currentCohort%treesai    ! SAI per unit crown area  (m2/m2)
    ! currentCohort%NV         ! The number of discrete vegetation
    !                          ! layers needed to describe this crown
    !
    ! The following patch level diagnostics are updated here:
    !
    ! currentPatch%canopy_layer_tlai(cl)   ! total leaf area index of canopy layer
    ! currentPatch%ncan(cl,ft)             ! number of vegetation layers needed
    !                                      ! in this patch's pft/canopy-layer
    ! currentPatch%nrad(cl,ft)             ! same as ncan, but does not include
    !                                      ! layers occluded by snow
    !                                      ! CURRENTLY SAME AS NCAN
    ! currentPatch%canopy_mask(cl,ft)      ! are there canopy elements in this pft-layer?
    !                                      ! (This is redundant with nrad though...)
    ! currentPatch%tlai_profile(cl,ft,iv)  ! m2 of leaves per m2 of the PFT's footprint
    ! currentPatch%elai_profile(cl,ft,iv)  ! non-snow covered m2 of leaves per m2 of PFT footprint
    ! currentPatch%tsai_profile(cl,ft,iv)  ! m2 of stems per m2 of PFT footprint
    ! currentPatch%esai_profile(cl,ft,iv)  ! non-snow covered m2 of stems per m2 of PFT footprint
    ! currentPatch%canopy_area_profile(cl,ft,iv)  ! Fractional area of leaf layer
    !                                             ! relative to vegetated area
    ! currentPatch%layer_height_profile(cl,ft,iv) ! Elevation of layer in m
    !
    ! -----------------------------------------------------------------------------------

    ! !USES:

    use EDtypesMod           , only : area, dinc_vai, dlower_vai, hitemax, n_hite_bins

    !
    ! !ARGUMENTS
    type(ed_site_type)     , intent(inout) :: currentSite


    !
    ! !LOCAL VARIABLES:
    type (ed_patch_type)  , pointer :: currentPatch
    type (ed_cohort_type) , pointer :: currentCohort
    real(r8) :: remainder                !Thickness of layer at bottom of canopy.
    real(r8) :: fleaf                    ! fraction of cohort incepting area that is leaves.
    integer  :: ft                       ! Plant functional type index.
    integer  :: iv                       ! Vertical leaf layer index
    integer  :: cl                       ! Canopy layer index
    real(r8) :: fraction_exposed         ! how much of this layer is not covered by snow?
    real(r8) :: layer_top_hite           ! notional top height of this canopy layer (m)
    real(r8) :: layer_bottom_hite        ! notional bottom height of this canopy layer (m)
    real(r8) :: frac_canopy(N_HITE_BINS) ! amount of canopy in each height class
    real(r8) :: minh(N_HITE_BINS)        ! minimum height in height class (m)
    real(r8) :: maxh(N_HITE_BINS)        ! maximum height in height class (m)
    real(r8) :: dh                       ! vertical detph of height class (m)
    real(r8) :: min_chite                ! bottom of cohort canopy  (m)
    real(r8) :: max_chite                ! top of cohort canopy      (m)
    real(r8) :: lai                      ! leaf area per canopy area
    real(r8) :: sai                      ! stem area per canopy area

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


    ! Here we are trying to generate a profile of leaf area, indexed by 'z' and by pft
    ! We assume that each point in the canopy recieved the light attenuated by the average
    ! leaf area index above it, irrespective of PFT identity...
    ! Each leaf is defined by how deep in the canopy it is, in terms of LAI units.  (FIX(RF,032414), GB)

    currentPatch => currentSite%oldest_patch
    do while(associated(currentPatch))

       ! --------------------------------------------------------------------------------
       ! Calculate tree and canopy areas.
       ! calculate tree lai and sai.
       ! --------------------------------------------------------------------------------

       currentPatch%canopy_layer_tlai(:)        = 0._r8
       currentPatch%ncan(:,:)                   = 0
       currentPatch%nrad(:,:)                   = 0
       currentPatch%tlai_profile(:,:,:)         = 0._r8
       currentPatch%tsai_profile(:,:,:)         = 0._r8
       currentPatch%elai_profile(:,:,:)         = 0._r8
       currentPatch%esai_profile(:,:,:)         = 0._r8
       currentPatch%layer_height_profile(:,:,:) = 0._r8
       currentPatch%canopy_area_profile(:,:,:)  = 0._r8
       currentPatch%canopy_mask(:,:)            = 0
       
       ! ------------------------------------------------------------------------------
       ! It is remotely possible that in deserts we will not have any canopy
       ! area, ie not plants at all...
       ! ------------------------------------------------------------------------------

       if (currentPatch%total_canopy_area > nearzero ) then

          call UpdatePatchLAI(currentPatch)

          ! -----------------------------------------------------------------------------
          ! Standard canopy layering model.
          ! Go through all cohorts and add their leaf area
          ! and canopy area to the accumulators.
          ! -----------------------------------------------------------------------------

          currentCohort => currentPatch%shortest
          do while(associated(currentCohort))
             ft = currentCohort%pft
             cl = currentCohort%canopy_layer

             ! ----------------------------------------------------------------
             ! How much of each tree is stem area index? Assuming that there is
             ! This may indeed be zero if there is a sensecent grass
             ! ----------------------------------------------------------------
             lai = currentCohort%treelai * currentCohort%c_area/currentPatch%total_canopy_area
             sai = currentCohort%treesai * currentCohort%c_area/currentPatch%total_canopy_area
             if( (currentCohort%treelai+currentCohort%treesai) > nearzero)then
                
                ! See issue: https://github.com/NGEET/fates/issues/899
                ! fleaf = currentCohort%treelai / (currentCohort%treelai + currentCohort%treesai)
                fleaf = lai / (lai+sai)
             else
                fleaf = 0._r8
             endif

             currentPatch%nrad(cl,ft) = currentPatch%ncan(cl,ft)

             if (currentPatch%nrad(cl,ft) > nlevleaf ) then
                write(fates_log(), *) 'Number of radiative leaf layers is larger'
                write(fates_log(), *) ' than the maximum allowed.'
                write(fates_log(), *) ' cl: ',cl
                write(fates_log(), *) ' ft: ',ft
                write(fates_log(), *) ' nlevleaf: ',nlevleaf
                write(fates_log(), *) ' currentPatch%nrad(cl,ft): ', currentPatch%nrad(cl,ft)
                call endrun(msg=errMsg(sourcefile, __LINE__))
             end if


             ! --------------------------------------------------------------------------
             ! Whole layers.  Make a weighted average of the leaf area in each layer
             ! before dividing it by the total area. Fill up layer for whole layers.
             ! --------------------------------------------------------------------------

             do iv = 1,currentCohort%NV

                ! This loop builds the arrays that define the effective (not snow covered)
                ! and total (includes snow covered) area indices for leaves and stems
                ! We calculate the absolute elevation of each layer to help determine if the layer
                ! is obscured by snow.

                layer_top_hite = currentCohort%hite - &
                     ( real(iv-1,r8)/currentCohort%NV * currentCohort%hite *  &
                     prt_params%crown_depth_frac(currentCohort%pft) )

                layer_bottom_hite = currentCohort%hite - &
                     ( real(iv,r8)/currentCohort%NV * currentCohort%hite * &
                     prt_params%crown_depth_frac(currentCohort%pft) )

                fraction_exposed = 1.0_r8
                if(currentSite%snow_depth  > layer_top_hite)then
                   fraction_exposed = 0._r8
                endif
                if(currentSite%snow_depth < layer_bottom_hite)then
                   fraction_exposed = 1._r8
                endif
                if(currentSite%snow_depth >= layer_bottom_hite .and. &
                     currentSite%snow_depth <= layer_top_hite) then !only partly hidden...
                   fraction_exposed =  1._r8 - max(0._r8,(min(1.0_r8,(currentSite%snow_depth -layer_bottom_hite)/ &
                        (layer_top_hite-layer_bottom_hite ))))
                endif

                if(iv==currentCohort%NV) then
                   remainder = (currentCohort%treelai + currentCohort%treesai) - &
                        (dlower_vai(iv) - dinc_vai(iv))
                   if(remainder > dinc_vai(iv) )then
                      write(fates_log(), *)'ED: issue with remainder', &
                           currentCohort%treelai,currentCohort%treesai,dinc_vai(iv), & 
                           currentCohort%NV,remainder

                      call endrun(msg=errMsg(sourcefile, __LINE__))
                   endif
                else
                   remainder = dinc_vai(iv)
                end if

                currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) + &
                     remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area

                currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) + &
                     remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * &
                     fraction_exposed

                currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) + &
                     remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area

                currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) + &
                     remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area * &
                     fraction_exposed

                currentPatch%canopy_area_profile(cl,ft,iv) = currentPatch%canopy_area_profile(cl,ft,iv) + &
                     currentCohort%c_area/currentPatch%total_canopy_area

                currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) + &
                     (remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * &
                     (layer_top_hite+layer_bottom_hite)/2.0_r8) !average height of layer.

             end do

             currentCohort => currentCohort%taller

          enddo !cohort

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

          ! If there is an upper-story, the top canopy layer
          ! should have a value of exactly 1.0 in its top leaf layer
          ! --------------------------------------------------------------------------

          if ( (currentPatch%NCL_p > 1) .and. &
               (sum(currentPatch%canopy_area_profile(1,:,1)) < 0.9999 )) then
             write(fates_log(), *) 'FATES: canopy_area_profile was less than 1 at the canopy top'
             write(fates_log(), *) 'cl: ',1
             write(fates_log(), *) 'iv: ',1
             write(fates_log(), *) 'sum(cpatch%canopy_area_profile(1,:,1)): ', &
                  sum(currentPatch%canopy_area_profile(1,:,1))
             currentCohort => currentPatch%shortest
             do while(associated(currentCohort))
                if(currentCohort%canopy_layer==1)then
                   write(fates_log(), *) 'FATES: cohorts',currentCohort%dbh,currentCohort%c_area, &
                        currentPatch%total_canopy_area,currentPatch%area
                   write(fates_log(), *) 'ED: fracarea', currentCohort%pft, &
                        currentCohort%c_area/currentPatch%total_canopy_area
                endif
                currentCohort => currentCohort%taller
             enddo !currentCohort
             call endrun(msg=errMsg(sourcefile, __LINE__))

          end if


          ! --------------------------------------------------------------------------
          ! In the following loop we are now normalizing the effective and
          ! total area profiles to convert from units of leaf/stem area per vegetated
          ! canopy area, into leaf/stem area per area of their own radiative column
          ! which is typically the footprint of all cohorts contained in the canopy
          ! layer x pft bins.
          ! Also perform some checks on area normalization.
          ! Check the area of each leaf layer, across pfts.
          ! It should never be larger than 1 or less than 0.
          ! --------------------------------------------------------------------------

          do cl = 1,currentPatch%NCL_p
             do iv = 1,currentPatch%ncan(cl,ft)

                if( debug .and. sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then

                   write(fates_log(), *) 'FATES: A canopy_area_profile exceeded 1.0'
                   write(fates_log(), *) 'cl: ',cl
                   write(fates_log(), *) 'iv: ',iv
                   write(fates_log(), *) 'sum(cpatch%canopy_area_profile(cl,:,iv)): ', &
                        sum(currentPatch%canopy_area_profile(cl,:,iv))
                   currentCohort => currentPatch%shortest
                   do while(associated(currentCohort))
                      if(currentCohort%canopy_layer==cl)then
                         write(fates_log(), *) 'FATES: cohorts in layer cl = ',cl, &
                              currentCohort%dbh,currentCohort%c_area, &
                              currentPatch%total_canopy_area,currentPatch%area
                         write(fates_log(), *) 'ED: fracarea', currentCohort%pft, &
                              currentCohort%c_area/currentPatch%total_canopy_area
                      endif
                      currentCohort => currentCohort%taller
                   enddo !currentCohort
                   call endrun(msg=errMsg(sourcefile, __LINE__))
                end if
             end do

             do ft = 1,numpft
                do iv = 1,currentPatch%ncan(cl,ft)

                   if( currentPatch%canopy_area_profile(cl,ft,iv) > nearzero ) then

                      currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) / &
                           currentPatch%canopy_area_profile(cl,ft,iv)

                      currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) / &
                           currentPatch%canopy_area_profile(cl,ft,iv)

                      currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) / &
                           currentPatch%canopy_area_profile(cl,ft,iv)

                      currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) / &
                           currentPatch%canopy_area_profile(cl,ft,iv)
                   end if

                   if(currentPatch%tlai_profile(cl,ft,iv)>nearzero )then
                      currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) &
                           /currentPatch%tlai_profile(cl,ft,iv)
                   end if

                enddo

             enddo
          enddo

          ! --------------------------------------------------------------------------
          ! Set the mask that identifies which PFT x can-layer combinations have
          ! scattering elements in them.
          ! --------------------------------------------------------------------------

          do cl = 1,currentPatch%NCL_p
             do ft = 1,numpft
                do  iv = 1, currentPatch%nrad(cl,ft)
                   if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then
                      currentPatch%canopy_mask(cl,ft) = 1
                   endif
                end do !iv
             enddo !ft
          enddo ! loop over cl

       end if

       currentPatch => currentPatch%younger

    enddo !patch

    return
  end subroutine leaf_area_profile

  ! ======================================================================================

  subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)

    ! ----------------------------------------------------------------------------------
    ! The purpose of this routine is to package output boundary conditions related
    ! to vegetation coverage to the host land model.
    ! ----------------------------------------------------------------------------------

    use EDTypesMod        , only : ed_patch_type, ed_cohort_type, &
         ed_site_type, AREA
    use FatesInterfaceTypesMod , only : bc_out_type

    !
    ! !ARGUMENTS
    integer,            intent(in)            :: nsites
    type(ed_site_type), intent(inout), target :: sites(nsites)
    integer,            intent(in)            :: fcolumn(nsites)
    type(bc_out_type),  intent(inout)         :: bc_out(nsites)

    ! Locals
    type (ed_cohort_type) , pointer :: currentCohort
    integer :: s, ifp, c, p
    type (ed_patch_type)  , pointer :: currentPatch
    real(r8) :: bare_frac_area
    real(r8) :: total_patch_area
    real(r8) :: total_canopy_area
    real(r8) :: total_patch_leaf_stem_area
    real(r8) :: weight  ! Weighting for cohort variables in patch
    
    do s = 1,nsites

       ifp = 0
       total_patch_area = 0._r8
       total_canopy_area = 0._r8
       bc_out(s)%canopy_fraction_pa(:) = 0._r8
       bc_out(s)%dleaf_pa(:) = 0._r8
       bc_out(s)%z0m_pa(:) = 0._r8
       bc_out(s)%displa_pa(:) = 0._r8
       
       currentPatch => sites(s)%oldest_patch
       c = fcolumn(s)
       do while(associated(currentPatch))

          if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then  ! ignore the bare-ground-PFT patch entirely for these BC outs

             ifp = ifp+1

             if ( currentPatch%total_canopy_area-currentPatch%area > 0.000001_r8 ) then
                if(debug)then
                   write(fates_log(),*) 'ED: canopy area bigger than area', &
                        currentPatch%total_canopy_area ,currentPatch%area
                end if
                currentPatch%total_canopy_area = currentPatch%area
             endif

             if (associated(currentPatch%tallest)) then
                bc_out(s)%htop_pa(ifp) = currentPatch%tallest%hite
             else
                ! FIX(RF,040113) - should this be a parameter for the minimum possible vegetation height?
                bc_out(s)%htop_pa(ifp) = 0.1_r8
             endif

             bc_out(s)%hbot_pa(ifp) = max(0._r8, min(0.2_r8, bc_out(s)%htop_pa(ifp)- 1.0_r8))

             ! Use canopy-only crown area weighting for all cohorts in the patch to define the characteristic
             ! Roughness length and displacement height used by the HLM
             ! use total LAI + SAI to weight the leaft characteristic dimension
             ! Avoid this if running in satellite phenology mode
             ! ----------------------------------------------------------------------------

             if (currentPatch%total_canopy_area > nearzero) then
                currentCohort => currentPatch%shortest
                do while(associated(currentCohort))
                   if (currentCohort%canopy_layer .eq. 1) then
                      weight = min(1.0_r8,currentCohort%c_area/currentPatch%total_canopy_area)
                      bc_out(s)%z0m_pa(ifp) = bc_out(s)%z0m_pa(ifp) + &
                           EDPftvarcon_inst%z0mr(currentCohort%pft) * currentCohort%hite * weight
                      bc_out(s)%displa_pa(ifp) = bc_out(s)%displa_pa(ifp) + &
                           EDPftvarcon_inst%displar(currentCohort%pft) * currentCohort%hite * weight
                   endif
                   currentCohort => currentCohort%taller
                end do

                ! for lai, scale to total LAI + SAI in patch.  first add up all the LAI and SAI in the patch
                total_patch_leaf_stem_area = 0._r8
                currentCohort => currentPatch%shortest
                do while(associated(currentCohort))
                   total_patch_leaf_stem_area = total_patch_leaf_stem_area + &
                        (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area
                   currentCohort => currentCohort%taller
                end do

                ! make sure there is some leaf and stem area
                if (total_patch_leaf_stem_area > nearzero) then
                   currentCohort => currentPatch%shortest
                   do while(associated(currentCohort))
                      ! weight dleaf by the relative totals of leaf and stem area
                      weight = (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area / total_patch_leaf_stem_area
                      bc_out(s)%dleaf_pa(ifp) = bc_out(s)%dleaf_pa(ifp) + &
                           EDPftvarcon_inst%dleaf(currentCohort%pft) * weight
                      currentCohort => currentCohort%taller
                   end do
                else
                   ! dummy case
                   bc_out(s)%dleaf_pa(ifp)  = EDPftvarcon_inst%dleaf(1)
                endif
             else
                ! if no canopy, then use dummy values (first PFT) of aerodynamic properties
                bc_out(s)%z0m_pa(ifp)    = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp)
                bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp)
                bc_out(s)%dleaf_pa(ifp)  = EDPftvarcon_inst%dleaf(1)
             endif
             ! -----------------------------------------------------------------------------

             ! We are assuming here that grass is all located underneath tree canopies.
             ! The alternative is to assume it is all spatial distinct from tree canopies.
             ! In which case, the bare area would have to be reduced by the grass area...
             ! currentPatch%total_canopy_area/currentPatch%area is fraction of this patch cover by plants
             ! currentPatch%area/AREA is the fraction of the soil covered by this patch.

             if(currentPatch%area.gt.0.0_r8)then
                bc_out(s)%canopy_fraction_pa(ifp) = &
                     min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)*(currentPatch%area/AREA)
             else
                bc_out(s)%canopy_fraction_pa(ifp) = 0.0_r8
             endif

             bare_frac_area = (1.0_r8 - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)) * &
                  (currentPatch%area/AREA)

             total_patch_area = total_patch_area + bc_out(s)%canopy_fraction_pa(ifp) + bare_frac_area

             total_canopy_area = total_canopy_area + bc_out(s)%canopy_fraction_pa(ifp)

             bc_out(s)%nocomp_pft_label_pa(ifp) = currentPatch%nocomp_pft_label

             ! Calculate area indices for output boundary to HLM
             ! It is assumed that cpatch%canopy_area_profile and cpat%xai_profiles
             ! have been updated (ie ed_leaf_area_profile has been called since dynamics has been called)

             bc_out(s)%elai_pa(ifp) = calc_areaindex(currentPatch,'elai')
             bc_out(s)%tlai_pa(ifp) = calc_areaindex(currentPatch,'tlai')
             bc_out(s)%esai_pa(ifp) = calc_areaindex(currentPatch,'esai')
             bc_out(s)%tsai_pa(ifp) = calc_areaindex(currentPatch,'tsai')

             ! Fraction of vegetation free of snow. This is used to flag those
             ! patches which shall under-go photosynthesis
             ! INTERF-TODO: we may want to stop using frac_veg_nosno_alb and let
             ! FATES internal variables decide if photosynthesis is possible
             ! we are essentially calculating it inside FATES to tell the
             ! host to tell itself when to do things (circuitous). Just have
             ! to determine where else it is used

             if ((bc_out(s)%elai_pa(ifp) + bc_out(s)%esai_pa(ifp)) > 0._r8) then
                bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 1.0_r8
             else
                bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 0.0_r8
             end if

          else  ! nocomp or SP, and currentPatch%nocomp_pft_label .eq. 0

             total_patch_area = total_patch_area + currentPatch%area/AREA

          end if
          currentPatch => currentPatch%younger
       end do

       ! Apply patch and canopy area corrections
       ! If the difference is above reasonable math precision, apply a fix
       ! If the difference is way above reasonable math precision, gracefully exit

       if(abs(total_patch_area-1.0_r8) > rsnbl_math_prec ) then

          if(abs(total_patch_area-1.0_r8) > 1.0e-8_r8 )then
             write(fates_log(),*) 'total area is wrong in update_hlm_dynamics',total_patch_area
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if

          if(debug) then
             write(fates_log(),*) 'imprecise patch areas in update_hlm_dynamics',total_patch_area
          end if

          currentPatch => sites(s)%oldest_patch
          ifp = 0
          do while(associated(currentPatch))
             if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then ! for vegetated patches only
                ifp = ifp+1
                bc_out(s)%canopy_fraction_pa(ifp) = bc_out(s)%canopy_fraction_pa(ifp)/total_patch_area
             endif ! veg patch

             currentPatch => currentPatch%younger
          end do

       endif

       ! If running hydro, perform a final check to make sure that we
       ! have conserved water. Since this is the very end of the dynamics
       ! cycle. No water should had been added or lost to the site during dynamics.
       ! With growth and death, we may have shuffled it around.
       ! For recruitment, we initialized their water, but flagged them
       ! to not be included in the site level balance yet, for they
       ! will demand the water for their initialization on the first hydraulics time-step

       if (hlm_use_planthydro.eq.itrue) then
          call UpdateH2OVeg(sites(s),bc_out(s),bc_out(s)%plant_stored_h2o_si,1)
       end if

       ! Pass FATES Harvested C to bc_out.
       call UpdateHarvestC(sites(s),bc_out(s))

    end do

    ! This call to RecruitWaterStorage() makes an accounting of
    ! how much water is used to intialize newly recruited plants.
    ! However, it does not actually move water from the soil or create
    ! a flux, it is just accounting for diagnostics purposes.  The water
    ! will not actually be moved until the beginning of the first hydraulics
    ! call during the fast timestep sequence

    if (hlm_use_planthydro.eq.itrue) then
       call RecruitWaterStorage(nsites,sites,bc_out)
    end if

  end subroutine update_hlm_dynamics

  ! =====================================================================================

  function calc_areaindex(cpatch,ai_type) result(ai)

    ! ----------------------------------------------------------------------------------
    ! This subroutine calculates the exposed leaf area index of a patch
    ! this is the square meters of leaf per square meter of ground area
    ! It does so by integrating over the depth and functional type profile of leaf area
    ! which are per area of crown.  This value has to be scaled by crown area to convert
    ! to ground area.
    ! ----------------------------------------------------------------------------------

    ! Arguments
    type(ed_patch_type),intent(in), target :: cpatch
    character(len=*),intent(in)            :: ai_type

    integer :: cl,ft
    real(r8) :: ai
    ! TODO: THIS MIN LAI IS AN ARTIFACT FROM TESTING LONG-AGO AND SHOULD BE REMOVED
    ! THIS HAS BEEN KEPT THUS FAR TO MAINTAIN B4B IN TESTING OTHER COMMITS
    real(r8),parameter :: ai_min = 0.1_r8

    real(r8),pointer   :: ai_profile

    ai = 0._r8
    if     (trim(ai_type) == 'elai') then
       do cl = 1,cpatch%NCL_p
          do ft = 1,numpft
             ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * &
                  cpatch%elai_profile(cl,ft,1:cpatch%nrad(cl,ft)))
          enddo
       enddo
    elseif (trim(ai_type) == 'tlai') then
       do cl = 1,cpatch%NCL_p
          do ft = 1,numpft
             ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * &
                  cpatch%tlai_profile(cl,ft,1:cpatch%nrad(cl,ft)))
          enddo
       enddo

    elseif (trim(ai_type) == 'esai') then
       do cl = 1,cpatch%NCL_p
          do ft = 1,numpft
             ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * &
                  cpatch%esai_profile(cl,ft,1:cpatch%nrad(cl,ft)))
          enddo
       enddo
    elseif (trim(ai_type) == 'tsai') then
       do cl = 1,cpatch%NCL_p
          do ft = 1,numpft
             ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * &
                  cpatch%tsai_profile(cl,ft,1:cpatch%nrad(cl,ft)))
          enddo
       enddo
    else

       write(fates_log(),*) 'Unsupported area index sent to calc_areaindex'
       call endrun(msg=errMsg(sourcefile, __LINE__))
    end if

    ai = max(ai_min,ai)

    return

  end function calc_areaindex

  ! ===============================================================================================

  subroutine CanopyLayerArea(currentPatch,site_spread,layer_index,layer_area)

    ! --------------------------------------------------------------------------------------------
    ! This function calculates the total crown area footprint for a desired layer of the canopy
    ! within a patch.
    ! The return units are the same as patch%area, which is m2
    ! ---------------------------------------------------------------------------------------------

    ! Arguments
    type(ed_patch_type),intent(inout), target   :: currentPatch
    real(r8),intent(in)                         :: site_spread
    integer,intent(in)                          :: layer_index
    real(r8),intent(inout)                      :: layer_area

    type(ed_cohort_type), pointer :: currentCohort
    
    
    layer_area = 0.0_r8
    currentCohort => currentPatch%tallest
    do while (associated(currentCohort))
       call carea_allom(currentCohort%dbh,currentCohort%n,site_spread, &
            currentCohort%pft,currentCohort%crowndamage, currentCohort%c_area)
       if (currentCohort%canopy_layer .eq. layer_index) then
          layer_area = layer_area + currentCohort%c_area
       end if
       currentCohort => currentCohort%shorter
    enddo
    return
   end subroutine CanopyLayerArea
  
  ! ===============================================================================================

  subroutine UpdatePatchLAI(currentPatch)

   ! --------------------------------------------------------------------------------------------
   ! This subroutine works through the current patch cohorts and updates the canopy_layer_tlai
   ! and related variables
   ! ---------------------------------------------------------------------------------------------

   ! Uses
   use EDtypesMod, only : dlower_vai

   ! Arguments
   type(ed_patch_type),intent(inout), target   :: currentPatch

   ! Local Variables
   type(ed_cohort_type), pointer :: currentCohort
   integer  :: cl                                  ! Canopy layer index
   integer  :: ft                                  ! Plant functional type index

   ! Calculate LAI of layers above.  Because it is possible for some understory cohorts
   ! to be taller than cohorts in the top canopy layer, we must iterate through the 
   ! patch by canopy layer first.  Given that canopy_layer_tlai is a patch level variable
   ! we could iterate through each cohort in any direction as long as we go down through
   ! the canopy layers.
   
   canopyloop: do cl = 1,nclmax
      currentCohort => currentPatch%tallest
      cohortloop: do while(associated(currentCohort))

         ! Only update the current cohort tree lai if lai of the above layers have been calculated
         if (currentCohort%canopy_layer .eq. cl) then
            ft     = currentCohort%pft
            
            ! Update the cohort level lai and related variables
            call UpdateCohortLAI(currentCohort,currentPatch%canopy_layer_tlai,  &
                 currentPatch%total_canopy_area)
            
            ! Update the number of number of vegetation layers
            currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV)

            ! Update the patch canopy layer tlai (LAI per canopy area)
            currentPatch%canopy_layer_tlai(cl) = currentPatch%canopy_layer_tlai(cl) +  &
                 currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area
            
         end if
         currentCohort => currentCohort%shorter

      end do cohortloop
   end do canopyloop

  end subroutine UpdatePatchLAI
  ! ===============================================================================================
  
  subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area)
   
   ! Update LAI and related variables for a given cohort
   
   ! Uses
   use EDtypesMod, only : dlower_vai
   
   ! Arguments
   type(ed_cohort_type),intent(inout), target   :: currentCohort
   real(r8), intent(in) :: canopy_layer_tlai(nclmax)  ! total leaf area index of each canopy layer
   real(r8), intent(in) :: total_canopy_area                  ! either patch%total_canopy_area or patch%area
   
   ! Local variables
   real(r8) :: leaf_c                              ! leaf carbon [kg]
      
   ! Obtain the leaf carbon
   leaf_c = currentCohort%prt%GetState(leaf_organ,carbon12_element)

   
   ! Note that tree_lai has an internal check on the canopy locatoin
   currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, &
        currentCohort%n, currentCohort%canopy_layer,               &
        canopy_layer_tlai,currentCohort%vcmax25top )
 
   if (hlm_use_sp .eq. ifalse) then
      currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%crowndamage, &
                                       currentCohort%canopy_trim, currentCohort%efstem_coh, &
                                       currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, &
                                       canopy_layer_tlai, currentCohort%treelai , &
                                       currentCohort%vcmax25top,4)
   end if

   ! Number of actual vegetation layers in this cohort's crown
   currentCohort%nv =  count((currentCohort%treelai+currentCohort%treesai) .gt. dlower_vai(:)) + 1
   
  end subroutine UpdateCohortLAI
  
  ! ===============================================================================================

  function NumPotentialCanopyLayers(currentPatch,site_spread,include_substory) result(z)

    ! --------------------------------------------------------------------------------------------
    ! Calculate the number of canopy layers in this patch.
    ! This simple call only determines total layering by querying the cohorts
    ! which layer they are in, it doesn't do any size evaluation.
    ! It may also, optionally, account for the temporary "substory", which is the imaginary
    ! layer below the understory which will be needed to temporarily accomodate demotions from
    ! the understory in the event the understory has reached maximum allowable area.
    ! --------------------------------------------------------------------------------------------

    type(ed_patch_type),target   :: currentPatch
    real(r8),intent(in)          :: site_spread
    logical                      :: include_substory

    type(ed_cohort_type),pointer :: currentCohort

    integer :: z
    real(r8) :: c_area
    real(r8) :: arealayer

    z = 1
    currentCohort => currentPatch%tallest
    do while (associated(currentCohort))
       z = max(z,currentCohort%canopy_layer)
       currentCohort => currentCohort%shorter
    enddo

    if(include_substory)then
       arealayer = 0.0
       currentCohort => currentPatch%tallest
       do while (associated(currentCohort))
          if(currentCohort%canopy_layer == z) then
             call carea_allom(currentCohort%dbh,currentCohort%n,site_spread,currentCohort%pft, &
                  currentCohort%crowndamage, c_area)
             arealayer = arealayer + c_area
          end if
          currentCohort => currentCohort%shorter
       enddo

       ! Does the bottom layer have more than a full canopy?
       ! If so we need to make another layer.
       if(arealayer > currentPatch%area)then
          z = z + 1
          if(hlm_use_sp.eq.itrue)then
             if(debug)then
                write(fates_log(),*) 'SPmode, canopy_layer full:',arealayer,currentPatch%area
             end if
          end if

       endif
    end if

  end function NumPotentialCanopyLayers

end module EDCanopyStructureMod
