!> \file
!! Calculates the leaf status for CTEM's pfts and leaf litter generated by
!! normal turnover of leaves and cold and drought stress; crop harvest; and
!! for grasses green leaves are converted into brown. Also updates vegetation
!! carbon pools for dcd trees and grasses accounting for nsc reallocation during leafout period.
!!
!> @author V. Arora, J. Melton, A. Asaadi
subroutine phenolgy (il1, il2, ilg, leapnow, tbar, thice, & ! In
                     thliq, THLW, THFC, ta, ipeatland, & ! In
                     anveg, iday, radl, roottemp, useTracer, & ! In
                     rmatctem, sort, fcancmx, isand, dayl_max, & ! In
                     gleafmas, gleafmas_ns, gleafmas_s, & ! In/Out
                     bleafmas, stemmass, stemmass_ns, stemmass_s, & ! In/Out
                     rootmass, rootmass_ns, rootmass_s, & ! In/Out
                     lfstatus, pandays, colddays, lfthrs, & ! In/Out
                     tracerGLeafMass, tracerBLeafMass, & ! In/Out
                     re_alloc_s2l, re_alloc_r2l, re_alloc_sr2l, & ! Out
                     gl2bl_grass_cflux, bl2ltr_grass_cflux, & ! Out
                     flhrloss, flhrloss_ns, flhrloss_s,leaflitr, leaflitr_ns, & ! Out
                     leaflitr_s, tracerLeafLitr) ! Out
  !
  !               Canadian Terrestrial Ecosystem Model (CTEM)
  !               Phenology, Leaf Turnover & Mortality Subroutine
  !
  !     28 Jun 2017   - Incorporate non-structural and structural carbohydrates
  !     A. Asaadi     
  !
  !     25 Apr 2019   - Incorporate tracer
  !     J. Melton
  !     06  Dec 2018  - Pass ilg back in as an argument + minor reorganization
  !     V. Arora        of arguments
  !
  !     7   Nov 2016  - Allow adaptions due to peatland scheme of Yuanqiao Wu.
  !
  !     14  Jan 2016  - There was a bit of hardwired code for 3 soil layers, that has been
  !                     fixed to allow >3.

  !     17  Jan 2014  - Moved parameters to global file (classicParams.f90)
  !     J. Melton
  !
  !     22  Jul 2013  - Add in module for parameters
  !     J. Melton
  !
  !     24  Sep 2012  - add in checks to prevent calculation of non-present
  !     J. Melton       pfts
  !
  !     15  Apr. 2003 - this subroutine calculates the leaf status
  !     V. Arora        for ctem's pfts and leaf litter generated by
  !                     normal turnover of leaves and cold and drought
  !                     stress. crop harvest is also modelled in this
  !                     subroutine, and for grasses green leaves are
  !                     converted into brown.
  use classicParams,    only : kn, pi, zero, kappa, eta, kappa_peat, eta_peat, &
                               lfespany, lfespany_peat, fracbofg, specsla, ignd, icc, &
                               ican, cdlsrtmx, drlsrtmx, drgta, &
                               colda, lwrthrsh, dayschk, coldlmt, &
                               coldthrs, harvthrs, flhrspan, &
                               thrprcnt, roothrsh, ctempfts, nlat, nmos, &
                               nol2pfts, min_ns2t_s, min_ns2t_r, re_alloc, &
                               classpfts, reindexPFTs
  use generalUtils,     only : findDaylength
  use ctemUtilities,    only : calcSLA

  implicit none

  integer, intent(in) :: ilg                 !< no. of grid cells in latitude circle
  integer, intent(in) :: il1                 !< il1=1
  integer, intent(in) :: il2                 !< il2=ilg
  integer, intent(in) :: iday                !< day of year
  integer, intent(in) :: isand(ilg,ignd) !<
  logical, intent(in) :: leapnow             !< true if this year is a leap year. Only used if the switch 'leap' is true.
  integer, intent(in) :: useTracer !< Switch for use of a model tracer. If useTracer is 0 then the tracer code is not used.
  !! useTracer = 1 turns on a simple tracer that tracks pools and fluxes. The simple tracer then requires that the tracer values in
  !!               the init_file and the tracerCO2file are set to meaningful values for the experiment being run.
  !! useTracer = 2 means the tracer is 14C and will then call a 14C decay scheme.
  !! useTracer = 3 means the tracer is 13C and will then call a 13C fractionation scheme.
  integer, intent(in) :: sort(icc)           !< index for correspondence between 9 pfts and the 12 values in parameters vectors
  integer, intent(in) :: ipeatland(ilg)        !< Peatland flag, non-peatlands are 0.
  real, intent(in) :: ta(ilg)                !< air temperature, k
  real, intent(in) :: tbar(ilg,ignd)         !< soil temperature, k
  real, intent(in) :: thliq(ilg,ignd)        !< liquid soil moisture content in soil layers
  real, intent(in) :: thice(ilg,ignd)        !< frozen  soil moisture content in soil layers
  real, intent(in) :: anveg(ilg,icc)         !< net photosynthesis rate of ctem's pfts, umol co2/m2.s
  real, intent(in) :: roottemp(ilg,icc)      !< root temperature, which is a function of soil temperature of course, k.
  real, intent(in) :: rmatctem(ilg,icc,ignd) !< fraction of roots in each soil layer for each pft
  real, intent(inout) :: stemmass(ilg,icc)   !< stem mass for each of the 9 ctem pfts, \f$kg c/m^2\f$
  real, intent(inout) :: stemmass_ns(ilg,icc)!< non-structural stem mass for each of the ctem pfts, \f$kg c/m^2\f$
  real, intent(in)    :: stemmass_s(ilg,icc) !< structural stem mass for each of the ctem pfts, \f$kg c/m^2\f$
  real, intent(inout) :: rootmass(ilg,icc)   !< root mass for each of the 9 ctem pfts, \f$kg c/m^2\f$
  real, intent(inout) :: rootmass_ns(ilg,icc)!< non-structural root mass for each of the 9 ctem pfts, \f$kg c/m^2\f$
  real, intent(in)    :: rootmass_s(ilg,icc) !< structural root mass for each of the 9 ctem pfts, \f$kg c/m^2\f$
  real, intent(in) :: fcancmx(ilg,icc)       !< max. fractional coverage of ctem's 9 pfts, but this can be modified
  !< by land-use change, and competition between pfts
  real, intent(in) :: THFC(ilg,ignd)      !< field capacity soil moisture content both calculated in allocate subroutine
  real, intent(in) :: THLW(ilg,ignd)       !< wilting point soil moisture content
  real, intent(in) :: radl(ilg)              !< latitude in radians
  real, intent(inout) :: gleafmas(ilg,icc)      !< green or live leaf mass in \f$kg c/m^2\f$
  real, intent(inout) :: gleafmas_ns(ilg,icc)   !< non-structural green or live leaf mass in \f$kg c/m^2\f$
  real, intent(inout) :: gleafmas_s(ilg,icc)    !< structural green or live leaf mass in \f$kg c/m^2\f$
  real, intent(inout) :: bleafmas(ilg,icc)      !< brown or dead leaf mass in \f$kg c/m^2\f$
  real, intent(inout) :: tracerGLeafMass(ilg,icc)      !< Tracer mass in the green leaf pool for each of the CTEM pfts, \f$kg c/m^2\f$
  real, intent(inout) :: tracerBLeafMass(ilg,icc)      !< Tracer mass in the brown leaf pool for each of the CTEM pfts, \f$kg c/m^2\f$
  integer, intent(inout) :: pandays(ilg,icc)    !< counter for positive net photosynthesis (an) days for initiating leaf onset
  integer, intent(inout) :: lfstatus(ilg,icc)   !< integer :: indicating leaf status or mode
  !< 1 - max. growth or onset, when all npp is allocated to leaves
  !< 2 - normal growth, when npp is allocated to leaves, stem, and root
  !< 3 - fall for dcd trees/harvest for crops, when allocation to leaves is zero.
  !< 4 - no leaves
  real, intent(in) :: dayl_max(ilg)             !< maximum daylength for at each location

  integer, intent(inout) :: colddays(ilg,2)     !< cold days counter for tracking days below a certain temperature threshold
  !< for ndl dcd and crop pfts.
  real, intent(inout) :: flhrloss(ilg,icc)      !< fall & harvest loss for bdl dcd plants and crops, respectively, \f$kg c/m^2\f$.
  real, intent(inout) :: flhrloss_ns(ilg,icc)   !< non-structural fall & harvest loss for bdl dcd plants and crops, respectively, \f$kg c/m^2\f$.
  real, intent(inout) :: flhrloss_s(ilg,icc)    !< structural fall & harvest loss for bdl dcd plants and crops, respectively, \f$kg c/m^2\f$.
  real, intent(out) :: leaflitr(ilg,icc)      !< leaf litter generated by normal turnover, cold and drought stress,
  !< and leaf fall/harvest, \f$kg c/m^2\f$
  real, intent(out) :: leaflitr_ns(ilg,icc)   !< non-structural leaf litter \f$kg c/m^2\f$
  real, intent(out) :: leaflitr_s(ilg,icc)    !< structural leaf litter \f$kg c/m^2\f$
  real, intent(out) :: tracerLeafLitr(ilg,icc)  !< leaf litter generated by normal turnover, cold and drought stress,
  !< and leaf fall/harvest, \f$tracer C units/m^2\f$
  real, intent(out) :: re_alloc_sr2l(ilg,icc) !< amount of nsc reallocated from stem and root to leaves during leaf out
  real, intent(out) :: re_alloc_s2l(ilg,icc)  !< amount of nsc reallocated from stem to leaves during leaf out
  real, intent(out) :: re_alloc_r2l(ilg,icc)  !< amount of nsc reallocated from root to leaves during leaf out
  real, intent(out) :: gl2bl_grass_cflux(ilg,icc)  !< carbon flux from gleafmas to bleafmass, to be used in nitrogen cycle
  real, intent(out) :: bl2ltr_grass_cflux(ilg,icc) !< litter carbon flux from bleafmas to litrmass, to be used in nitrogen cycle

  real :: day                    !<
  real :: daylngth(ilg)          !<
  real :: nrmlloss(ilg,icc)      !< leaf loss due to normal turnover
  real :: nrmlloss_ns(ilg,icc)   !< non-structural leaf loss due to normal turnover
  real :: nrmlloss_s(ilg,icc)    !< structural leaf loss due to normal turnover
  real :: tracerNrmlLoss(ilg,icc)!< Tracer leaf loss due to normal turnover
  real :: betadrgt(ilg,ignd)     !< (1 - drought stress)
  real :: drgtstrs(ilg,icc)      !< drought stress term
  real :: drgtlsrt(ilg,icc)      !< drought loss rate
  real :: drgtloss(ilg,icc)      !< leaf loss due to drought stress
  real :: drgtloss_ns(ilg,icc)   !< non-structural leaf loss due to drought stress
  real :: drgtloss_s(ilg,icc)    !< structural leaf loss due to drought stress
  real :: tracerDrgtLoss(ilg,icc)!< Tracer leaf loss due to drought stress
  real :: coldloss(ilg,icc)      !< leaf loss due to cold stress
  real :: coldloss_ns(ilg,icc)   !< non-structural leaf loss due to cold stress
  real :: coldloss_s(ilg,icc)    !< structural leaf loss due to cold stress
  real :: tracerColdLoss(ilg,icc)!< Tracer leaf loss due to cold stress
  real :: coldstrs(ilg,icc)      !< cold stress term
  real :: coldlsrt(ilg,icc)      !< cold loss rate
  real :: tracerFlHrLoss(ilg,icc)!< Tracer fall & harvest loss for bdl dcd plants and crops, respectively, \f$kg c/m^2\f$.
  real :: lfthrs(ilg,icc)        !< threshold lai for finding leaf status
  character(8) :: pftkind
  integer :: i, j, k, m, n
  integer :: chkmode(ilg,icc)    !< indicator for making sure that leaf status is updated
  real :: sla(icc)               !< specific leaf area
  real :: ailcg(ilg,icc)         !< green lai
  real :: ailcb(ilg,icc)         !< brown lai
  real :: frac                   !< temp var.
  real :: lastDOY                !< last day of the year 
  real :: lwrtemp                !< temp var. to set two different lwrthrsh before and after summer solstice.
  real :: daylgtloss(ilg,icc)    !< leaf loss due to decline in daylight
  real :: daylgtloss_ns(ilg,icc) !< non-structural leaf loss due to decline in daylight
  real :: daylgtloss_s(ilg,icc)  !< structural leaf loss due to decline in daylight
  real :: stemns2t, rootns2t     !< temp ratios for nonstructural to total pools
  real :: ns_s2sr,ns_r2sr        !< temp ratios for nonstructural pools: stem/(stem+root) & root/(stem+root)
  real :: tempdaylgtloss         !< temp var. used in daylgtloss 

  ! Initialize required arrays to zero
  betadrgt = 0.0
  sla = 0.0
  ailcg = 0.0
  ailcb = 0.0
  chkmode = 0.0
  leaflitr = 0.0
  leaflitr_ns = 0.0
  leaflitr_s = 0.0
  tracerLeafLitr = 0.
  nrmlloss = 0.0
  nrmlloss_ns = 0.0
  nrmlloss_s = 0.0
  drgtstrs = 0.0
  drgtlsrt = 0.0
  drgtloss = 0.0
  drgtloss_ns = 0.0
  drgtloss_s = 0.0
  coldstrs = 0.0
  coldlsrt = 0.0
  coldloss = 0.0
  coldloss_ns = 0.0
  coldloss_s = 0.0
  lfthrs = 0.0
  ! flhrloss = 0.0
  ! flhrloss_ns = 0.0
  ! flhrloss_s = 0.0
  daylgtloss = 0.0
  daylgtloss_ns = 0.0
  daylgtloss_s = 0.0
  re_alloc_sr2l = 0.0
  re_alloc_s2l = 0.0
  re_alloc_r2l = 0.0
  gl2bl_grass_cflux = 0.0
  bl2ltr_grass_cflux = 0.0
  tempdaylgtloss = 0.0
  tracerFlHrLoss = 0.
  tracerNrmlLoss = 0.
  tracerDrgtLoss = 0.
  tracerColdLoss = 0.
  ! Initialization ends

  !------------------------------------------------------------------
  !> Set the number of days in this year.
  if (.not. leapnow) then
    lastDOY = 365.
  else 
    lastDOY = 366.
  end if 
  
  !> Convert green leaf mass into leaf area index using specific leaf
  !! area \f$(sla, m^2 /kg c)\f$ estimated using leaf life span. see bio2str
  !! subroutine for more details.
  do j = 1,icc
    ! if (ipeatland(i) /= 1 .and. ipeatland(i) /= 2) then ! Uplands
    !   sla(j) = 25.0 * (lfespany(sort(j)) ** ( -0.50))
    ! else ! peatlands
    !   sla(j) = 25.0 * (lfespany_peat(sort(j)) ** ( -0.50))
    ! end if 
    ! if (specsla(sort(j)) > zero) sla(j) = specsla(sort(j))
    n = sort(j)
    do i = il1,il2
      sla(j) = calcSLA(ipeatland(i),sort(j))
      if (fcancmx(i,j) > 0.0) then
        ailcg(i,j) = sla(j) * gleafmas(i,j)
        ailcb(i,j) = sla(j) * bleafmas(i,j) * fracbofg

        !> also find threshold lai as a function of stem+root biomass
        !! which is used to determine leaf status
        if (ipeatland(i) /= 1 .and. ipeatland(i) /= 2) then !uplands
          lfthrs(i,j) = ((stemmass_s(i,j) + rootmass_s(i,j)) &
                          / eta(n)) ** (1.0 / kappa(n))
        else ! peatlands 
          lfthrs(i,j) = ((stemmass_s(i,j) + rootmass_s(i,j)) &
                          / eta_peat(n)) ** (1.0 / kappa_peat(n))
        end if 
        lfthrs(i,j) = (thrprcnt(n) / 100.0) * sla(j) * lfthrs(i,j)

        !> using green leaf area index (ailcg) determine the leaf status for
        !! each pft. loops 190 and 200 thus initialize lfstatus, if this
        !! this information is not passed specifically as an initialization quantity.
        if (lfstatus(i,j) == 0) then
          if (ailcg(i,j) <= zero) then
            lfstatus(i,j) = 4                      ! no leaves
          else if (ailcg(i,j) > lfthrs(i,j)) then
            lfstatus(i,j) = 2                      ! normal growth
          else
            lfstatus(i,j) = 4                   ! treat this as no leaves
          end if                                  ! so that we start growing
        end if                                   ! if possible
      end if ! fcancmx
    end do ! loop 180
  end do ! loop 170

  !> Knowing lfstatus (after initialization above or using value from
  !! from previous time step) we decide if we stay in a given leaf
  !! mode or we move to some other mode.
  !!
  !! we start with the "no leaves" mode
  !! ----------------------------------
  !!
  !! add one to pandays(i,j) if daily an is positive, otherwise set it to zero.
  do j = 1,icc
    n = sort(j)
    do i = il1,il2
      if (fcancmx(i,j) > 0.0) then
        if (anveg(i,j) > zero) then
          pandays(i,j) = pandays(i,j) + 1
          if (pandays(i,j) > dayschk(n)) then
            pandays(i,j) = dayschk(n)
          end if
        else
          pandays(i,j) = 0
        end if
      end if
    end do ! loop 230
  end do ! loop 220

  !> if in "no leaves" mode check if an has been positive over last
  !! dayschk(j) days to move into "max. growth" mode. if not we stay
  !! in "no leaves" mode. also set the chkmode(i,j) switch to 1.
  do j = 1,icc
    n = sort(j)
    do i = il1,il2
      if (fcancmx(i,j) > 0.0) then
        if (chkmode(i,j) == 0 .and. lfstatus(i,j) == 4) then
          if (pandays(i,j) >= dayschk(n)) then
            lfstatus(i,j) = 1        ! switch to "max. growth" mode
            chkmode(i,j) = 1         ! mode checked, no more checks further down
          else
            lfstatus(i,j) = 4        ! stay in "no leaves" mode
            chkmode(i,j) = 1         ! mode checked, no more checks further down
          end if
        end if
      end if
    end do ! loop 250
  end do ! loop 240

  !> find day length using day of year and latitude. this is to be used for
  !! initiating leaf offset for broad leaf dcd trees.
  day = real(iday)
  do i = il1,il2
    daylngth(i) = findDaylength(day,radl(i))
  end do

  !! similar to the way we count no. of days when an is positive, we
  !! find no. of days when temperature is below -5 c. we need this to
  !! determine if we go into "leaf fall" mode for needle leaf dcd trees.
  !!
  !! also estimate no. of days below 8 C. we use these days to decide
  !! if its cold enough to harvest crops.
  do k = 1,2
    do i = il1,il2
      if (ta(i) < (coldthrs(k) + 273.16)) then
        colddays(i,k) = colddays(i,k) + 1
        if (colddays(i,k) > coldlmt(k)) then
          colddays(i,k) = coldlmt(k)
        end if
      else
        colddays(i,k) = 0
      end if
    end do ! loop 290
  end do ! loop 280

  !> Even if pandays criteria has been satisfied do not go into max. growth
  !! mode if environmental conditions are such that they will force leaf
  !! fall or harvest mode.
  do i = il1,il2
    do j = 1,icc
      if (fcancmx(i,j) > 0.0) then
        if (lfstatus(i,j) == 1 .and. chkmode(i,j) == 1) then
          pftkind = ctempfts(j)
          select case (pftkind)

          case ('NdlDcdTr') ! Needle leaf deciduous tree

            if (ta(i) < (coldthrs(1) + 273.16)) lfstatus(i,j) = 4

          case ('BdlDCoTr') ! Broadleaf deciduous cold tree

            if (roottemp(i,j) < (roothrsh + 273.16) .or. &
                (daylngth(i) < 11.0 .and. &
                roottemp(i,j) < (11.15 + 273.16))) lfstatus(i,j) = 4

          case ('BdlDDrTr') ! Broadleaf deciduous drought dry tree

            if (roottemp(i,j) < (roothrsh + 273.16) .or. &
                (daylngth(i) < 11.0 .and. &
                roottemp(i,j) < (11.15 + 273.16))) lfstatus(i,j) = 4

          case ('CropC3  ')

            if (ta(i) < (coldthrs(2) + 273.16)) lfstatus(i,j) = 4

          case ('CropC4  ')

            if (ta(i) < (coldthrs(2) + 273.16)) lfstatus(i,j) = 4

          case ('BdlDCoSh') ! Broadleaf deciduous shrub

            !           if (roottemp(i,j)<(roothrsh+273.16).or.
            !     &    (daylngth(i)<11.0.and.roottemp(i,j)<(11.15+273.16))) then
            if (ta(i) < (coldthrs(1) + 273.16)) lfstatus(i,j) = 4

          case ('GrassC3 ','Sedge   ')

            if (ta(i) < (coldthrs(2) + 273.16)) lfstatus(i,j) = 4

          case ('GrassC4 ')

            if (ta(i) < (coldthrs(2) + 273.16)) lfstatus(i,j) = 4

          case ('NdlEvgTr','BdlEvgTr','BdlEvgSh')

            ! Evergreen, leave unchanged.

          case default

            print * ,'Unknown CTEM PFT in phenolgy ',pftkind
            call errorHandler('phenolgy', - 1)

          end select
        end if
      end if
    end do
  end do ! loop 270

  !> if in "max growth" mode
  !! ------------------------
  !!
  !> Re-allocating non-structural carbon from stem and root ns pools (for dcd trees and grasses)
  !! during early spring (max growth mode) in order to push out leaf faster following
  !! observation. For grasses, stemmass is equal to zero thus a different loop is used.
  do j = 1,icc
       pftkind = ctempfts(j)
       do i = il1,il2
        if (lfstatus(i,j) == 1 .and. ailcg(i,j) < lfthrs(i,j) .and. &
                                             fcancmx(i,j) > 0.0) then
 
         ! For now, we use the following condition to allow reallocation only before
         ! summer solstice and aviod it on occasions in the fall season when lfstatus=1.
         ! This model bug (lfstatus=1 in fall) might need to be fixed in future.
         if ((radl(i) > 0. .and. iday < 172 ) .or. &
             (radl(i) < 0. .and. (iday >= 172 .and. iday < 356))) then
 
          select case (pftkind)

          case ('NdlDcdTr','BdlDCoTr','BdlDDrTr','BdlDCoSh')
          
            if (stemmass(i,j) > 0.0 .and. rootmass(i,j) > 0.0) then

               stemns2t = stemmass_ns(i,j) / stemmass(i,j)
               rootns2t = rootmass_ns(i,j) / rootmass(i,j)

               if (stemns2t < 1.0 .and. rootns2t < 1.0) then
                 ns_s2sr = stemmass_ns(i,j) / (stemmass_ns(i,j) + rootmass_ns(i,j))
                 ns_r2sr = rootmass_ns(i,j) / (stemmass_ns(i,j) + rootmass_ns(i,j))
                 re_alloc_sr2l(i,j) = re_alloc(sort(j)) * lfthrs(i,j) &
                                      * 1000. / (sla(j)*15.)

                 if (stemns2t > min_ns2t_s(sort(j)) .and. rootns2t > min_ns2t_r(sort(j))) then
                  re_alloc_s2l(i,j) = re_alloc_sr2l(i,j) * ns_s2sr
                  re_alloc_r2l(i,j) = re_alloc_sr2l(i,j) * ns_r2sr
                 else if (stemns2t > min_ns2t_s(sort(j)) .and. &
                          rootns2t <= min_ns2t_r(sort(j))) then
                  re_alloc_s2l(i,j) = re_alloc_sr2l(i,j)
                  re_alloc_r2l(i,j) = 0.0
                 else if (stemns2t <= min_ns2t_s(sort(j)) .and. &
                          rootns2t > min_ns2t_r(sort(j))) then
                  re_alloc_s2l(i,j) = 0.0
                  re_alloc_r2l(i,j) = re_alloc_sr2l(i,j)
                 else if (stemns2t <= min_ns2t_s(sort(j)) .and. &
                          rootns2t <= min_ns2t_r(sort(j))) then
                  re_alloc_sr2l(i,j)= 0.0
                  re_alloc_s2l(i,j) = 0.0
                  re_alloc_r2l(i,j) = 0.0
                 end if

               end if

            end if
          
          case ('GrassC3 ','GrassC4 ','Sedge   ')
    
            if (rootmass(i,j) > 0.0) then

               rootns2t = rootmass_ns(i,j) / rootmass(i,j)
               
               if (rootns2t /= 1.0) then
 
                 re_alloc_sr2l(i,j) = re_alloc(sort(j))* lfthrs(i,j)*1000./(sla(j)*15.)
 
                 if (rootns2t .gt. min_ns2t_r(sort(j))) then
                   re_alloc_s2l(i,j)  = 0.0
                   re_alloc_r2l(i,j)  = re_alloc_sr2l(i,j)
                 else
                   re_alloc_sr2l(i,j) = 0.0
                   re_alloc_s2l(i,j)  = 0.0
                   re_alloc_r2l(i,j)  = 0.0
                 end if

               end if

           end if

         case ('NdlEvgTr','BdlEvgTr','CropC3  ','CropC4  ','BdlEvgSh')

            ! Evergreen and crops leave unchanged.
          case default

            print * ,'Unknown CTEM PFT in phenolgy ',pftkind
            call errorHandler('phenolgy', - 2)

          end select
         end if !before summer solstice
        end if ! lfstatus == 1

        !! Update the non-structural pools and refind the total pool sizes
        gleafmas_ns(i,j) = gleafmas_ns(i,j) + re_alloc_sr2l(i,j) / 1000.
        !gleafmas(i,j)    = gleafmas(i,j)    + re_alloc_sr2l(i,j) / 1000.
        gleafmas(i,j)    = gleafmas_ns(i,j) + gleafmas_s(i,j)

        stemmass_ns(i,j) = stemmass_ns(i,j) - re_alloc_s2l(i,j) / 1000.
        !stemmass(i,j)    = stemmass(i,j)    - re_alloc_s2l(i,j) / 1000.
        stemmass(i,j)    = stemmass_ns(i,j) + stemmass_s(i,j)

        rootmass_ns(i,j) = rootmass_ns(i,j) - re_alloc_r2l(i,j) / 1000.
        !rootmass(i,j)    = rootmass(i,j)    - re_alloc_r2l(i,j) / 1000.
        rootmass(i,j)    = rootmass_ns(i,j) + rootmass_s(i,j)
 
        !! Undo the reallocation if it leads to negative pools
        if (gleafmas_ns(i,j) < 0.0 .or. gleafmas(i,j) < 0.0 .or.  &
            stemmass_ns(i,j) < 0.0 .or. stemmass(i,j) < 0.0 .or.  &
            rootmass_ns(i,j) < 0.0 .or. rootmass(i,j) < 0.0) then
            gleafmas_ns(i,j) = gleafmas_ns(i,j) - re_alloc_sr2l(i,j) / 1000.
            !gleafmas(i,j)    = gleafmas(i,j)    - re_alloc_sr2l(i,j)/ 1000.
             gleafmas(i,j)    = gleafmas_ns(i,j) + gleafmas_s(i,j)
            stemmass_ns(i,j) = stemmass_ns(i,j) + re_alloc_s2l(i,j) / 1000.
            !stemmass(i,j)    = stemmass(i,j)    + re_alloc_s2l(i,j) / 1000.
            stemmass(i,j)    = stemmass_ns(i,j) + stemmass_s(i,j)
            rootmass_ns(i,j) = rootmass_ns(i,j) + re_alloc_r2l(i,j) / 1000.
            !rootmass(i,j)    = rootmass(i,j)    + re_alloc_r2l(i,j) / 1000.
            rootmass(i,j)    = rootmass_ns(i,j) + rootmass_s(i,j)
            re_alloc_sr2l(i,j) = 0.0
            re_alloc_s2l(i,j)  = 0.0
            re_alloc_r2l(i,j)  = 0.0
         end if
 
       end do ! loop 290
  end do ! loop 280
  !!
  !! if mode hasn't been checked and we are in "max. growth" mode, then
  !! check if we are above pft-dependent lai threshold. if lai is more
  !! then this threshold we move into "normal growth" mode, otherwise
  !! we stay in "max growth" mode so that leaves can grow at their
  !! max. climate-dependent rate
  do j = 1,icc
    do i = il1,il2
      if (fcancmx(i,j) > 0.0) then
        if (chkmode(i,j) == 0 .and. lfstatus(i,j) == 1) then
          if (ailcg(i,j) >= lfthrs(i,j)) then
            lfstatus(i,j) = 2        ! switch to "normal growth" mode
            chkmode(i,j) = 1
          else if (ailcg(i,j) <= zero) then
            lfstatus(i,j) = 4        ! switch to "no leaves" mode
            chkmode(i,j) = 1
            pandays(i,j) = 0
          else
            lfstatus(i,j) = 1        ! stay in "max. growth" mode
            chkmode(i,j) = 1
          end if

          !> for dcd trees we also need to go into "leaf fall" mode
          !! directly from "max. growth" mode.
          pftkind = ctempfts(j)
          select case (pftkind)
          case ('NdlDcdTr')

            if (ailcg(i,j) < lfthrs(i,j) .and. colddays(i,1) >= coldlmt(1) .and. &
                ailcg(i,j) > zero) then
              lfstatus(i,j) = 3
              chkmode(i,j) = 1
            end if

          case ('BdlDCoTr')

            if (ailcg(i,j) > zero .and. ((daylngth(i) < 11.0 .and. roottemp(i,j) < (11.15 + 273.16)) & ! FLAG put 11.15 in param !
                .or. roottemp(i,j) < (roothrsh + 273.16))) then
              lfstatus(i,j) = 3
              chkmode(i,j) = 1
              flhrloss(i,j) = gleafmas(i,j) * (1.0 / flhrspan(2))
              flhrloss_ns(i,j) = gleafmas_ns(i,j) * (1.0 / flhrspan(2))
              flhrloss_s(i,j)  = gleafmas_s(i,j)  * (1.0 / flhrspan(2))
              if (useTracer > 0) tracerFlHrLoss(i,j) = tracerGLeafMass(i,j) * (1.0 / flhrspan(2))

            end if

          case ('BdlDDrTr')

            if (ailcg(i,j) > zero .and. ((daylngth(i) < 11.0 .and. roottemp(i,j) < (11.15 + 273.16)) &
                .or. roottemp(i,j) < (roothrsh + 273.16))) then
              lfstatus(i,j) = 3        ! go into "leaf fall" mode
              chkmode(i,j) = 1
            end if

          case ('BdlDCoSh') ! Shrub

            if (ailcg(i,j) > zero .and. ((daylngth(i) < 11.0 .and. roottemp(i,j) < (11.15 + 273.16)) & ! FLAG param to move out.
                .or. roottemp(i,j) < (roothrsh + 273.16))) then
              lfstatus(i,j) = 3        ! go into "leaf fall" mode
              chkmode(i,j) = 1
              flhrloss(i,j) = gleafmas(i,j) * (1.0 / flhrspan(2))
              flhrloss_ns(i,j) = gleafmas_ns(i,j) * (1.0 / flhrspan(2))
              flhrloss_s(i,j)  = gleafmas_s(i,j)  * (1.0 / flhrspan(2))
              if (useTracer > 0) tracerFlHrLoss(i,j) = tracerGLeafMass(i,j) * (1.0 / flhrspan(2))

            end if

          case ('NdlEvgTr' , 'BdlEvgTr','CropC3  ','CropC4  ','GrassC3 ','GrassC4 ','BdlEvgSh','Sedge   ')

            ! Evergreen,crop or grass so do nothing.

          case default

            print * ,'Unknown CTEM PFT in phenolgy ',pftkind
            call errorHandler('phenolgy', - 3)

          end select
        end if
      end if
    end do ! loop 310
  end do ! loop 300
  !>
  !! if in "normal growth" mode then go through every pft individually
  !! and follow set of rules to determine if we go into "fall/harvest" mode
  !!
  do i =  il1,il2
    do j = 1,icc
      if (fcancmx(i,j) > 0.0) then
        if (chkmode(i,j) == 0 .and. lfstatus(i,j) == 2) then
          pftkind = ctempfts(j)
          select case (pftkind)

          case ('NdlEvgTr') ! needle leaf evg tree

            if (ailcg(i,j) < lfthrs(i,j) .and. ailcg(i,j) > zero) then
              lfstatus(i,j) = 1         ! go back to "max. growth" mode
              chkmode(i,j) = 1
            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4         ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
            else
              lfstatus(i,j) = 2         ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case ('NdlDcdTr') ! Needle leaf deciduous tree

            if (ailcg(i,j) < lfthrs(i,j) .and. colddays(i,1) > coldlmt(1) &
                .and. ailcg(i,j) > zero) then
              lfstatus(i,j) = 3         ! go into "leaf fall" mode
              chkmode(i,j) = 1
            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4         ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
            else
              lfstatus(i,j) = 2         ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case ('BdlEvgTr') ! broad leaf evg tree

            if (ailcg(i,j) < lfthrs(i,j) .and. ailcg(i,j) > zero) then
              lfstatus(i,j) = 1         ! go back to "max. growth" mode
              chkmode(i,j) = 1
            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4         ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
            else
              lfstatus(i,j) = 2         ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case ('BdlDCoTr') ! Broadleaf deciduous cold tree

            ! we use daylength and roottemp to initiate leaf offset
            if (ailcg(i,j) > zero .and. ((daylngth(i) < 11.0 .and. &
                roottemp(i,j) < (11.15 + 273.16)) .or. &
                roottemp(i,j) < (roothrsh + 273.16))) then
              lfstatus(i,j) = 3         ! go into "leaf fall" mode
              chkmode(i,j) = 1
              flhrloss(i,j) = gleafmas(i,j) * (1.0 / flhrspan(2))
              flhrloss_ns(i,j) = gleafmas_ns(i,j) * (1.0 / flhrspan(2))
              flhrloss_s(i,j)  = gleafmas_s(i,j)  * (1.0 / flhrspan(2))
              if (useTracer > 0) tracerFlHrLoss(i,j) = tracerGLeafMass(i,j) * (1.0 / flhrspan(2))

            else if (ailcg(i,j) > zero .and. ailcg(i,j) < lfthrs(i,j)) then
              lfstatus(i,j) = 1         ! switch to "max. growth" mode
              chkmode(i,j) = 1
            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4         ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
              flhrloss(i,j) = 0.0
              flhrloss_ns(i,j) = 0.0
              flhrloss_s(i,j)  = 0.0
              if (useTracer > 0) tracerFlHrLoss(i,j) = 0.
            else
              lfstatus(i,j) = 2         ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case ('BdlDDrTr') ! Broadleaf deciduous drought dry tree

            !! We still use daylength and roottemp to initiate leaf offset,
            !! for the pathological cases of dry dcd trees being further
            !! away from the equator then we can imagine. other wise leaf
            !! loss will occur due to drought anyway.
            if (ailcg(i,j) > zero .and. &
                ((daylngth(i) < 11.0 .and. roottemp(i,j) < (11.15 + 273.16)) .or. &
                roottemp(i,j) < (roothrsh + 273.16))) then
              lfstatus(i,j) = 3         ! go into "leaf fall" mode
              chkmode(i,j) =  1
            else if (ailcg(i,j) > zero .and. ailcg(i,j) < lfthrs(i,j)) then
              lfstatus(i,j) = 1         ! switch to "max. growth" mode
              chkmode(i,j) = 1
            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4         ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
            else
              lfstatus(i,j) = 2         ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case ('CropC3  ','CropC4  ') ! Crops

            !! "normal growth" to "fall/harvest" transition for crops is based on
            !! specified lai. we harvest if lai of crops reaches a threshold.
            !! if lai doesn't reach this threshold (say due to a bad year)
            !! we harvest anyway if it starts getting cold, otherwise we don't harvest.
            n = sort(j)
            if (ailcg(i,j) >= harvthrs(n)) then
              lfstatus(i,j) = 3        ! go into "harvest" mode
              chkmode(i,j) = 1
              flhrloss(i,j) = gleafmas(i,j) * (1.0 / flhrspan(1))
              flhrloss_ns(i,j) = gleafmas_ns(i,j) * (1.0 / flhrspan(1))
              flhrloss_s(i,j)  = gleafmas_s(i,j)  * (1.0 / flhrspan(1))
              if (useTracer > 0) tracerFlHrLoss(i,j) = tracerGLeafMass(i,j) * (1.0 / flhrspan(1))

            else if (ailcg(i,j) > zero .and. colddays(i,2) >= coldlmt(2)) then
              lfstatus(i,j) = 3        ! go into "harvest" mode
              chkmode(i,j) = 1          ! regardless of lai
              flhrloss(i,j) = gleafmas(i,j) * (1.0 / flhrspan(1))
              flhrloss_ns(i,j) = gleafmas_ns(i,j) * (1.0 / flhrspan(1))
              flhrloss_s(i,j)  = gleafmas_s(i,j)  * (1.0 / flhrspan(1))
              if (useTracer > 0) tracerFlHrLoss(i,j) = tracerGLeafMass(i,j) * (1.0 / flhrspan(1))

            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4        ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
              flhrloss(i,j) = 0.0
              flhrloss_ns(i,j) = 0.0
              flhrloss_s(i,j)  = 0.0
              if (useTracer > 0) tracerFlHrLoss(i,j) = 0.

            else
              lfstatus(i,j) = 2        ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case ('GrassC3 ','GrassC4 ') ! Grasses

            !! "normal growth" to "max. growth" transition for grasses
            if (ailcg(i,j) < lfthrs(i,j) .and. ailcg(i,j) > zero) then
              lfstatus(i,j) = 1        ! switch back to "max. growth" mode
              chkmode(i,j) = 1
            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4        ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
            else
              lfstatus(i,j) = 2        ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case ('BdlEvgSh') ! EVG-SHRUB, treated the same as needel leaf EVG

            if (ailcg(i,j) > lfthrs(i,j) .and. ailcg(i,j) > zero) then
              lfstatus(i,j) = 1         ! go back to "max. growth" mode
              chkmode(i,j) = 1
            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4         ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
            else
              lfstatus(i,j) = 2         ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case ('BdlDCoSh') ! DCD-SHRUB     treated the same as needle leaf dcd

            if (ailcg(i,j) < lfthrs(i,j) .and. colddays(i,1) >= coldlmt(1) .and. &
                ailcg(i,j) > zero) then
              lfstatus(i,j) = 3         ! go into "leaf fall" mode
              chkmode(i,j) = 1
            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4         ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
            else
              lfstatus(i,j) = 2         ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case ('Sedge   ') ! Sedges

            if (ailcg(i,j) < lfthrs(i,j) .and. ailcg(i,j) > zero) then
              lfstatus(i,j) = 1        ! switch back to "max. growth" mode
              chkmode(i,j) = 1
            else if (ailcg(i,j) <= zero) then
              lfstatus(i,j) = 4        ! switch to "no leaves" mode
              chkmode(i,j) = 1
              pandays(i,j) = 0
            else
              lfstatus(i,j) = 2        ! stay in "normal growth" mode
              chkmode(i,j) = 1
            end if

          case default

            print * ,'Unknown CTEM PFT in phenolgy ',pftkind
            call errorHandler('phenolgy', - 4)

          end select
        end if
      end if
    end do
  end do ! loop 320


  !> If in "fall/harvest" mode
  !! --------------------------
  !!
  !! grasses and evg trees do not come into this mode, because they want
  !! to stay green if possible. this mode is activated for dcd plants and
  !! crops. once in this mode dcd trees loose their leaves and crops are
  !! harvested. ndl dcd trees keep loosing their leaves at rate determined
  !! by cold stress, bdl dcd trees loose their leaves at a specified
  !! rate, and crops are harvested over a period of  15 days. dcd trees
  !! and crops stay in "leaf fall/harvest" model until all green leaves
  !! are gone at which time they switch into "no leaves" mode, and then
  !! wait for the climate to become favourable to go into "max. growth" mode
  do i = il1,il2
    do j = 1,icc
      pftkind = ctempfts(j)
      select case (pftkind)

      case ('NdlDcdTr','BdlDCoTr','BdlDDrTr','BdlDCoSh','CropC3  ','CropC4  ') ! All deciduous PFTs

        if (fcancmx(i,j) > 0.0 .and. chkmode(i,j) == 0 .and. lfstatus(i,j) == 3) then
          if (ailcg(i,j) <= 0.01) then
            lfstatus(i,j) = 4            ! go into "no leaves" mode
            chkmode(i,j) = 1
            pandays(i,j) = 0
            flhrloss(i,j) = 0.0
            flhrloss_ns(i,j) = 0.0
            flhrloss_s(i,j)  = 0.0
            if (useTracer > 0) tracerFlHrLoss(i,j) = 0.

          else
            if (pftkind == 'NdlDcdTr' .or. pftkind == 'BdlDCoSh') then !
              if (pandays(i,j) >= dayschk(j) .and. ta(i) > (coldthrs(1) + 273.16)) then
                if (ailcg(i,j) < lfthrs(i,j)) then
                  lfstatus(i,j) = 1      ! go into "max. growth" mode
                  chkmode(i,j) = 1
                else
                  lfstatus(i,j) = 2      ! go into "normal growth" mode
                  chkmode(i,j) = 1
                end if
              else
                lfstatus(i,j) = 3        ! stay in "fall/harvest" mode
                chkmode(i,j) = 1
              end if
            else if (pftkind == 'BdlDCoTr' .or. pftkind == 'BdlDDrTr') then
              if ((pandays(i,j) >= dayschk(j)) .and. ((roottemp(i,j) > (roothrsh + 273.16)) .and. &
                  (daylngth(i) > 11.0))) then
                if (ailcg(i,j) < lfthrs(i,j)) then
                  lfstatus(i,j) = 1      ! go into "max. growth" mode
                  chkmode(i,j) = 1
                else
                  lfstatus(i,j) = 2      ! go into "normal growth" mode
                  chkmode(i,j) = 1
                end if
              else
                lfstatus(i,j) = 3        ! stay in "fall/harvest" mode
                chkmode(i,j) = 1
              end if
            else if (pftkind == 'CropC3  ' .or. pftkind == 'CropC4  ') then  ! Crops
              lfstatus(i,j) = 3          ! stay in "fall/harvest" mode
              chkmode(i,j) = 1
            end if ! pftkind
          end if ! ailcg
        end if ! fcancmax etc.

      case ('NdlEvgTr' ,'BdlEvgTr','GrassC3 ','GrassC4 ','Sedge   ','BdlEvgSh')

        ! Evergreens and grasses do not enter this mode.

      case default

        print * ,'Unknown CTEM PFT in phenolgy ',pftkind
        call errorHandler('phenolgy', - 5)

      end select
    end do ! loop 400
  end do ! loop 410

  !> Check that leaf status of all vegetation types in all grid cells has been updated
  do j = 1,icc
    do i = il1,il2
      if (fcancmx(i,j) > 0.0) then
        if (chkmode(i,j) == 0) then
          write(6,2000) i,j
2000      format(' at (i) = (',i3,'),pft = ',i2,' lfstatus not updated')
          call errorHandler('phenolgy', - 6)
        end if
      end if
    end do ! loop 412
  end do ! loop 411

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

  !> Having decided leaf status for every pft, we now calculate normal
  !! leaf turnover, cold and drought stress mortality, and for bdl dcd
  !! plants we also calculate specified loss rate if they are in "leaf fall"
  !! mode, and for crops we calculate harvest loss, if they are in "harvest" mode.
  !!
  !! all these loss calculations will yield leaf litter in \f$kg c/m^2\f$ for the given day for all pfts

  ! Normal leaf turn over
  do j = 1,icc
    n = sort(j)
    do i = il1,il2
      if (fcancmx(i,j) > 0.0) then
        if (ipeatland(i) /= 1 .and. ipeatland(i) /= 2) then ! Uplands
          nrmlloss(i,j) = gleafmas(i,j) * (1.0 - exp( -1.0 / (lastDOY * lfespany(n))))
          nrmlloss_ns(i,j) = gleafmas_ns(i,j) * (1.0 - exp( -1.0 / (lastDOY * lfespany(n))))
          nrmlloss_s(i,j)  = gleafmas_s(i,j)  * (1.0 - exp( -1.0 / (lastDOY * lfespany(n))))
          if (useTracer > 0) tracerNrmlLoss(i,j) = tracerGLeafMass(i,j) * (1.0 - exp( -1.0 / (lastDOY * lfespany(n))))
        else ! peatlands 
          nrmlloss(i,j) = gleafmas(i,j) * (1.0 - exp( -1.0 / (lastDOY * lfespany_peat(n))))
          nrmlloss_ns(i,j) = gleafmas_ns(i,j) * (1.0 - exp( -1.0 / (lastDOY * lfespany_peat(n))))
          nrmlloss_s(i,j)  = gleafmas_s(i,j)  * (1.0 - exp( -1.0 / (lastDOY * lfespany_peat(n))))
          if (useTracer > 0) tracerNrmlLoss(i,j) = tracerGLeafMass(i,j) * (1.0 - exp( -1.0 / (lastDOY * lfespany_peat(n))))
        end if 
      end if   ! fcancmx
    end do ! loop 430
  end do ! loop 420

  !> For drought stress related mortality we need field capacity and wilting point soil
  !! moisture contents, which we calculated in allocate subroutine
  do j = 1,ignd
    do i = il1,il2
      if (isand(i,j) /= - 3) then
        ! for non-bedrock estimate (1-drought stress)
        ! (bedrock keeps the initialization value of 0)
        if (thliq(i,j) <= THLW(i,j)) then
          betadrgt(i,j) = 0.0
        else if (thliq(i,j) > (THLW(i,j) - thice(i,j)) .and. (thliq(i,j) < (THFC(i,j) - thice(i,j)))) then
          betadrgt(i,j) = thliq(i,j) - THLW(i,j) - thice(i,j)
          betadrgt(i,j) = betadrgt(i,j) / (THFC(i,j) - thice(i,j) - THLW(i,j))
        else
          betadrgt(i,j) = 1.0
        end if
        betadrgt(i,j) = max(0.0,min(1.0,betadrgt(i,j)))
      end if
    end do ! loop 460
  end do ! loop 450

  !> Estimate drought stress term averaged over the rooting depth for each pft
  do j = 1,icc
    n = sort(j)
    pftkind = ctempfts(j)    
    do i = il1,il2
      if (fcancmx(i,j) > 0.0) then
        drgtstrs(i,j) =  sum((1.0 - betadrgt(i,:)) * rmatctem(i,j,:))
        drgtstrs(i,j) = drgtstrs(i,j) / sum(rmatctem(i,j,:))
        !    ---------YW May 04,2015 ---------------------------------------/
        !          drgtstrs(i,j)=max(0.0,min(1.0,drgtstrs(i,j))) ! FLAG why is this commented out? JM Nov 2016.

        !> Using this drought stress term and our two vegetation-dependent
        !! parameters we find leaf loss rate associated with drought

        ! drought related leaf loss rate
        drgtlsrt(i,j) = drlsrtmx(n) * drgtstrs(i,j) ** drgta(n)

        !> Estimate leaf loss in \f$kg c/m^2\f$ due to drought stress
        drgtloss(i,j) = gleafmas(i,j) * (1.0 - exp( - drgtlsrt(i,j)))
        drgtloss_ns(i,j) = gleafmas_ns(i,j) * (1.0 - exp( - drgtlsrt(i,j)))
        drgtloss_s(i,j)  = gleafmas_s(i,j)  * (1.0 - exp( - drgtlsrt(i,j)))

        if (useTracer > 0) tracerDrgtLoss(i,j) = tracerGLeafMass(i,j) * (1.0 - exp( - drgtlsrt(i,j)))

        !! We use a lower value for lwrthrsh before summer solstice which includes leafout period in spring.
        !! This helps to get a more prominent shifting of the seasonal LAI at the leafout period.

        lwrtemp = lwrthrsh(n)

        select case (pftkind)

        case ('BdlDCoTr','BdlDCoSh')

           if ((radl(i) > 0. .and. iday <= 172) .or. (radl(i) < 0. .and. &
               (iday >= 356 .or. iday <= 172))) then ! before summer solstice
               lwrtemp = lwrthrsh(n) - 4.0
           end if

        case ('NdlEvgTr','NdlDcdTr','BdlEvgTr','BdlDDrTr','CropC3  ','CropC4  ','GrassC3 ','GrassC4 ',&
              'BdlEvgSh','Sedge   ')

           ! other pfts, do nothing.

        case default
           print * ,'Unknown CTEM PFT in phenolgy ',pftkind
           call errorHandler('phenolgy', - 7)

        end select
        
        !> Similar to drgtstrs we find coldstrs for each pft. we assume that
        !! max. cold stress related leaf loss occurs when temperature is 5 c
        !! or more below pft's threshold
        if (ta(i) <= (lwrthrsh(n) - 5.0 + 273.16)) then
          coldstrs(i,j) = 1.0
        else if (ta(i) > (lwrthrsh(n) - 5.0 + 273.16) .and. &
                     ta(i) < (lwrthrsh(n) + 273.16)) then
          coldstrs(i,j) = 1.0 - ((ta(i) - (lwrthrsh(n) - 5.0 + 273.16)) / 5.0)
        else
          coldstrs(i,j) = 0.0
        end if
        coldstrs(i,j) = max(0.0,min(1.0,coldstrs(i,j)))

        !> Using this cold stress term and our two vegetation-dependent
        !> parameters we find leaf loss rate associated with cold
        !> cold related leaf loss rate
        coldlsrt(i,j) = cdlsrtmx(n) * coldstrs(i,j) ** colda(n)

        !> estimate leaf loss in \f$kg c/m^2\f$ due to cold stress
        coldloss(i,j) = gleafmas(i,j) * (1.0 - exp( - coldlsrt(i,j)))
        coldloss_ns(i,j) = gleafmas_ns(i,j) * (1.0 - exp( - coldlsrt(i,j)))
        coldloss_s(i,j)  = gleafmas_s(i,j)  * (1.0 - exp( - coldlsrt(i,j)))

        if (useTracer > 0) tracerColdLoss(i,j) = tracerGLeafMass(i,j) * (1.0 - exp( - coldlsrt(i,j)))

        !> For NdlDcdTr & BdlDCoTr trees trigger leaf loss as a function of day length as well.
        !! The idea is to get rid of leaves at a maximum rate by the time day length hits 9.5 hours
        !! after summer solstice. The rationale is that even in the Cretaceous period
        !! when atmos. CO2 was ~1000 ppm, there were deciduous trees that shed their leaves in polar
        !! regions. Since it was so warm back then it wasn't the cold stress, it was day length that
        !! gave deciduous trees their deciduous nature.

        select case (pftkind)

        case ('NdlDcdTr','BdlDCoTr','BdlDCoSh')

          if ((radl(i) > 0. .and. iday >= 172) .or. &  ! NH, after summer solstice
              (radl(i) < 0. .and. (iday >= 356 .or. iday <= 172))) then ! SH, after summer solstice
   
             tempdaylgtloss = 0.15 * max(0., (1. - max(((daylngth(i) - 8.0) &
                                           / (dayl_max(i) - 8.0)), 0.0))**1.0)
             daylgtloss(i,j)    = tempdaylgtloss * gleafmas(i,j)
             daylgtloss_ns(i,j) = tempdaylgtloss * gleafmas_ns(i,j)
             daylgtloss_s(i,j)  = tempdaylgtloss * gleafmas_s(i,j)

          end if

        case ('NdlEvgTr','BdlEvgTr','BdlDDrTr','CropC3  ','CropC4  ','GrassC3 ','GrassC4 ',&
              'BdlEvgSh','Sedge   ')

          ! other pfts, do nothing

        case default
          print * ,'Unknown CTEM PFT in phenolgy ',pftkind
          call errorHandler('phenolgy', - 8)

        end select
      end if
    end do ! loop 490
  end do ! loop 480

  !> Now that we have all types of leaf losses (due to normal turnover,
  !! cold and drought stress, and fall/harvest) we take the losses
  !! for grasses and use those to turn live green grass into dead
  !! brown grass. we then find the leaf litter from the brown grass
  !! which will then go into the litter pool.
  do j = 1,icc
    pftkind = ctempfts(j)
    select case (pftkind)

    case ('GrassC3 ','GrassC4 ','Sedge   ')

      n = sort(j)
      do i = il1,il2
        if (fcancmx(i,j) > 0.0) then
          !gleafmas(i,j) = gleafmas(i,j) - nrmlloss(i,j) - drgtloss(i,j) - coldloss(i,j) - daylgtloss(i,j)
          gleafmas_ns(i,j) = gleafmas_ns(i,j) - nrmlloss_ns(i,j) - drgtloss_ns(i,j) &
                            - coldloss_ns(i,j)- daylgtloss_ns(i,j)
          gleafmas_s(i,j)  = gleafmas_s(i,j)  - nrmlloss_s(i,j)  - drgtloss_s(i,j) &
                            - coldloss_s(i,j) - daylgtloss_s(i,j)
          gleafmas(i,j) = gleafmas_ns(i,j) + gleafmas_s(i,j)
          
          if (useTracer > 0) tracerGLeafMass(i,j) = tracerGLeafMass(i,j) - tracerNrmlLoss(i,j) &
              - tracerDrgtLoss(i,j) - tracerColdLoss(i,j)

          if (gleafmas(i,j) < 0.0) then

            gl2bl_grass_cflux(i,j) = nrmlloss(i,j) + drgtloss(i,j) + coldloss(i,j) + gleafmas(i,j) &
                                     + daylgtloss(i,j)
            bleafmas(i,j) = bleafmas(i,j) + gl2bl_grass_cflux(i,j)
            gleafmas(i,j) = 0.0
            gleafmas_ns(i,j) = 0.0
            gleafmas_s(i,j)  = 0.0

            if (useTracer > 0) then
              tracerBLeafMass(i,j) = tracerBLeafMass(i,j) + tracerNrmlLoss(i,j) &
                                     + tracerDrgtLoss(i,j) &
                                     + tracerColdLoss(i,j) + tracerGLeafMass(i,j)
              tracerGLeafMass(i,j) = 0.
            end if

          else
            gl2bl_grass_cflux(i,j) = nrmlloss(i,j) + drgtloss(i,j) + coldloss(i,j) &
                                     + daylgtloss(i,j)
            bleafmas(i,j) = bleafmas(i,j) + gl2bl_grass_cflux(i,j)
            if (useTracer > 0) tracerBLeafMass(i,j) = tracerBLeafMass(i,j) + tracerNrmlLoss(i,j) &
                + tracerDrgtLoss(i,j) + tracerColdLoss(i,j)
          end if
          
          nrmlloss(i,j) = 0.0
          nrmlloss_ns(i,j) = 0.0
          nrmlloss_s(i,j)  = 0.0
          drgtloss(i,j) = 0.0
          drgtloss_ns(i,j) = 0.0
          drgtloss_s(i,j)  = 0.0
          coldloss(i,j) = 0.0
          coldloss_ns(i,j) = 0.0
          coldloss_s(i,j)  = 0.0
          daylgtloss(i,j)    = 0.0
          daylgtloss_ns(i,j) = 0.0
          daylgtloss_s(i,j)  = 0.0

          if (useTracer > 0) then
            tracerNrmlLoss(i,j) = 0.
            tracerDrgtLoss(i,j) = 0.
            tracerColdLoss(i,j) = 0.
          end if

          !> we assume life span of brown grass is 10% that of green grass
          !! but this is an adjustable parameter.
          if (ipeatland(i) /= 1 .and. ipeatland(i) /= 2) then ! Uplands
            nrmlloss(i,j) = bleafmas(i,j) * (1.0 - exp( - 1.0 / (0.10 * lastDOY * lfespany(n))))
            if (useTracer > 0) tracerNrmlLoss(i,j) = tracerBLeafMass(i,j) &
                       * (1.0 - exp( - 1.0 / (0.10 * lastDOY * lfespany(n))))
          else ! peatlands 
            nrmlloss(i,j) = bleafmas(i,j) * (1.0 - exp( - 1.0 / (0.10 * lastDOY * lfespany_peat(n))))
            if (useTracer > 0) tracerNrmlLoss(i,j) = tracerBLeafMass(i,j) &
                       * (1.0 - exp( - 1.0 / (0.10 * lastDOY * lfespany_peat(n))))
          end if 
          bl2ltr_grass_cflux(i,j) = nrmlloss(i,j)
        end if
      end do

    case ('NdlEvgTr' , 'NdlDcdTr', 'BdlEvgTr','BdlDCoTr', 'BdlDDrTr', &
          'CropC3  ', 'CropC4  ','BdlEvgSh','BdlDCoSh')

      ! Do nothing for the non-grasses

    case default

      print * ,'Unknown CTEM PFT in phenolgy ',pftkind
      call errorHandler('phenolgy', - 6)

    end select
  end do ! loop 620

  !> combine nrmlloss, drgtloss, and coldloss together to get total
  !! leaf litter generated, which is what we use to update gleafmass,
  !! except for grasses, for which we have already updated gleafmass.
  do j = 1,icc
    do i = il1,il2
      if (fcancmx(i,j) > 0.0) then
        leaflitr(i,j) = nrmlloss(i,j) + drgtloss(i,j) &
                        + coldloss(i,j) + flhrloss(i,j) + daylgtloss(i,j)
        leaflitr_ns(i,j) = nrmlloss_ns(i,j) + drgtloss_ns(i,j) &
                           + coldloss_ns(i,j) + flhrloss_ns(i,j) + daylgtloss_ns(i,j)
        leaflitr_s(i,j) = nrmlloss_s(i,j) + drgtloss_s(i,j) &
                          + coldloss_s(i,j) + flhrloss_s(i,j) + daylgtloss_s(i,j)
        if (useTracer > 0) tracerLeafLitr(i,j) = tracerNrmlLoss(i,j) &
            + tracerDrgtLoss(i,j) + tracerColdLoss(i,j) + tracerFlHrLoss(i,j)

      end if
    end do ! loop 660
  end do ! loop 650

  return
end subroutine phenolgy
!> \file
!> The leaf phenology parametrization used in CTEM v. 1.0 is described in detail by \cite Arora2005-6b1 . Changes between
!! version 1.0 and 2.0 are limited to parameter values and the parametrization is briefly described here. There are four
!! different leaf phenological states in which vegetation can be at a given instant: (i) no leaves or dormant,
!! (ii) maximum growth, (iii) normal growth and (iv) leaf fall or harvest. PFTs may go through only some, or all,
!! of these phenological states depending on their deciduousness. A broadleaf cold deciduous tree, for example,
!! transitions through all these four states in a year. In winter, the broadleaf cold deciduous trees are in the
!! no leaves/dormant state; favourable climatic conditions in spring trigger leaf growth and the tree enters the
!! maximum leaf growth state when all the NPP is allocated to leaves to accelerate leaf out; when the LAI reaches a
!! threshold (described below) the tree enters the normal leaf growth state and NPP is also allocated to stem and root
!! components; finally the arrival of autumn triggers leaf fall and the trees go into the leaf fall mode where no carbon
!! is allocated to leaves (but it continues for roots and stems). When all the leaves have been shed, the trees go into
!! the no leaves or dormant state again and the cycle is repeated the next year. The evergreen tree PFTs and the grass
!! PFTs do not enter the leaf fall state and maintain a leaf canopy as long as environmental conditions are favourable.
!! Although drought and cold stress cause accelerated leaf loss compared to the normal leaf turnover from these PFTs,
!! they do not explicitly go into the leaf fall mode where the intent is to lose all leaves in a specified amount of time.
!!
!! The leaf phenological state transitions are dependent upon environmental conditions. In particular, the transition
!! from no leaves/dormant state to the maximum growth state is based on the carbon-gain approach. CTEM uses
!! \f$\textit{virtual}\f$ leaves to assess favourable meteorological conditions for leaf out. The virtual leaves
!! photosynthesize and respire in a manner similar to normal leaves except the carbon gain or loss is not taken into
!! account in vegetation's carbon balance. A positive net leaf photosynthesis rate (\f$ G_{canopy, net} \f$, Eq. \f$\ref{Gnet}\f$)
!! for the virtual leaves over seven consecutive days indicates the arrival of favourable growth conditions and triggers
!! leaf onset and the associated transition from the no leaves/dormant state to the maximum leaf growth state, when the
!! entire positive NPP is allocated to leaves (\f$a_{fL} = 1\f$, \f$a_{fS} = a_{fR} = 0\f$). When LAI reaches \f${LAI}_{thrs}\f$
!! then the vegetation switches to the normal growth mode and positive NPP is allocated to all three vegetation
!! components -- leaves, stem and roots (\f$a_{fL}\f$, \f$a_{fS}\f$, \f$a_{fR}\f$ > 0). \f${LAI}_{thrs}\f$ is calculated as
!!
!! \f[
!! {LAI}_{thrs} = L_f [ {SLA} (\frac{C_S + C_R}{\eta} ) ^ {1/\kappa} ].
!! \f]
!!
!! The PFT-specific \f$L_f\f$ term (see also classicParams.f90) calculates \f${LAI}_{thrs}\f$ to be typically between
!! 40 and 50 % of the maximum LAI that a given amount of stem and root biomass can support (based on the terms in the square
!! brackets and Eq. (\f$\ref{propwoody}\f$). \f$SLA\f$ is the specific leaf area (Eq.\f$ \ref{sla}\f$).  This rule for transition
!! from a maximum to a normal growth state is also used for evergreen tree PFTs and grass PFTs. Similar to
!! \f${LAI}_{thrs}\f$, the LAI of virtual leaves is \f$7.5\, {\%}\f$ of the maximum LAI a given amount of root
!! and stem biomass can support for tree and crop PFTsand \f$2.5\, {\%}\f$ for grass PFTs. In
!! addition, the LAI of virtual leaves is constrained to be, at least, \f$0.3\, m^2\, m^{-2}\f$
!! for tree PFTs and \f$0.2\, m^2\, m^{-2}\f$ for crop and grass PFTs.
!!
!! The transition from the normal growth state to the leaf fall state is triggered by unfavourable environmental conditions
!! and shorter day length. Broadleaf deciduous trees transition to the leaf fall state when either: (i) day length is less
!! than \f$11\, h\f$ and the rooting zone temperature drops below \f$11.15\, C\f$ or (ii) when the rooting zone temperature
!! drops below \f$8\, C\f$ regardless of the day length. Needleleaf deciduous tress begin leaf fall after seven consecutive
!! days with daily mean air temperature below \f$-5\, C\f$. Leaf fall occurs over a period of 15 days. In the leaf fall
!! state, the vegetation continues carbon allocation to its root and stem components, but not to leaves (\f$a_{fL} = 0\f$,
!! \f$a_{fS} + a_{fR} = 1\f$). Evergreen trees and grasses do not enter the leaf fall state and neither do the
!! broadleaf drought deciduous trees. The implication for the latter PFT is that if the climate changes and the
!! dry season becomes shorter, then the trees will keep their leaves on for a longer period of time since broadleaf
!! drought deciduous trees lose leaves due to soil moisture stress (described below).
!!
!! The model vegetation is able to transition between the different leaf phenological states in response to changing
!! conditions. For example, a leaf out in spring for broadleaf cold deciduous trees can be interrupted by a cold event
!! when the vegetation goes into a leaf fall state until the return of more favourable conditions.
!!
!! Leaf litter generation is caused by normal turnover of leaves (\f$\Omega_N\f$, \f$day^{-1}\f$) and also
!! due to cold (\f$ \Omega_C \f$, \f$day^{-1}\f$) and drought (\f$\Omega_D\f$, \f$day^{-1}\f$) stress,
!! both of which contribute to seasonality of LAI. For example, the leaf loss associated with drought and reduced
!! photosynthesis during the dry season are the principal causes of the seasonality of LAI for the broadleaf drought deciduous tree PFT.
!!
!! The conversion of leaf carbon to leaf litter (\f$D_L\f$, \f$kg\, C\, m^{-2}\, day^{-1}\f$) is expressed as
!!
!! \f[ D_L = C_L[1 - \exp(-\Omega_N - \Omega_C - \Omega_{D})], \f]
!!
!! where (\f$\Omega_{N, C, D}\f$, \f$day^{-1}\f$) are the leaf loss rates associated with normal turnover of leaves and
!! the cold and drought stress. The rate of normal turnover of leaves is governed by PFT-specific leaf lifespan
!! (\f$\tau_L\f$, \f$yr\f$) as \f$\Omega_N= 1/(number of days in year) \tau_L\f$ (see also classicParams.f90}
!! for PFT specific values of \f$\tau_L\f$).  The leaf loss rate associated with cold stress
!! (\f$\Omega_C\f$) is calculated as
!!
!! \f[ \label{gamma_cold} \Omega_C = \Omega_{C, max}L_{cold}^3, \f]
!!
!! where \f$\Omega_{C, max}\f$ (\f$day^{-1}\f$, see also classicParams.f90) is the maximum cold stress loss rate.
!! \f$L_{cold}\f$ is a scalar that varies between 0 and 1 as
!!
!! \f[ \label{cldls} L_{cold} = \begin{cases} 1, \quad T_a < \left(T_{cold}^{leaf} -
!! 5\right) \\ 1 - \frac{T_a - \left(T_{cold}^{leaf} - 5\right)}{5}, \\ \quad T_{cold}^{leaf}
!! > T_a > (T_{cold}^{leaf} - 5) \\ 0, \quad T_a > T_{cold}^{leaf} , \\ \end{cases} \f]
!!
!! where \f$T_{cold}^{leaf}\f$ is a PFT-specific temperature threshold below which a PFT experiences damage
!! to its leaves promoting leaf loss (see also classicParams.f90) and \f$T_a\f$ is the daily mean air
!! temperature (\f$C\f$).  The leaf loss rate due to drought stress is calculated in a similar manner
!!
!! \f[ \label{gamma_dry} \Omega_{D} = \Omega_{D, max}\, (1-\phi_{root})^3, \f]
!!
!! where \f$\Omega_{D, max}\f$ (\f$day^{-1}\f$, see also classicParams.f90) is the maximum drought
!! stress loss rate and \f$\phi_{root}\f$ (Eq. \ref{degsoilsat}) is the degree of soil saturation in the rooting zone.
