! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/triffid.F
      subroutine TRIFFID (LAND_PTS, LAND_INDEX, FORW, GAMMA, FRAC_VS
     &,                   FRAC_AGR, FRAC_MIN, FRAC_SEED, DENOM_MIN
     &,                   BF, G_LEAF, NPP, RESP_S, RESP_W, CS

     &,                   RESP_D, CS_D, DZ_GND, SCHDMIN, CS_PART
     &,                   SCHDMAX, DMIN_HOLD, PFCDENSE,TRACK_PFC
     &,                   SQU_PFC,FZN_SL

     &,                   ACTIVELT,KOVDIFF_P,KOVCON,KOVMAX,YEARKOV_SW
     &,                   ZBOT,PERMA,CS_P,RESP_P, V_SAT,KOVSAT
     &,                   AFp, AF_Ta, iAF, PREG_C

     &,                   FRAC, HT, LAI, C_VEG, CV, LIT_C, LIT_C_T)

!-----------------------------------------------------------------------
! Simulates changes in vegetation structure, areal
! coverage and the carbon contents of vegetation and soil.
! can be used to advance these variables dynamically
! (GAMMA=1/TIMESTEP) or to iterate toward  equilibrium
! (GAMMA --> 0.0, FORW=1.0).

!**********************************************************************
! this file is based on code that may have had the following copyright:
! (c) CROWN COPYRIGHT 1997, U.K. METEOROLOGICAL OFFICE.

! Permission has been granted by the authors to the public to copy
! and use this software without charge, provided that this Notice and
! any statement of authorship are reproduced on all copies. Neither the
! Crown nor the U.K. Meteorological Office makes any warranty, express
! or implied, or assumes any liability or responsibility for the use of
! this software.
!**********************************************************************
!-----------------------------------------------------------------------

      implicit none

      include "size.h"
      include "mtlm_data.h"
      include "pconst_triffid.h"

! LAND_PTS   = IN Number of points on which TRIFFID may operate.
! LAND_INDEX = IN Indices of land points on which TRIFFID may operate.

      integer LAND_PTS, LAND_INDEX(POINTS)
      integer L, N, T

! FORW      = IN Forward timestep weighting.
! FRAC_VS   = IN Total fraction of gridbox covered by veg or soil.
! GAMMA     = IN Inverse timestep (/360days).
! FRAC_AGR  = IN Fraction of agriculture.
! FRAC_MIN  = IN Minimum areal fraction for PFTs.
! FRAC_SEED = IN "Seed" fraction for PFTs.
! DENOM_MIN = IN Minimum value for the denominator of the update
!             equation. Ensures that gradient descent does not lead
!             to an unstable solution.
! BF        = Burn fraction
! G_LEAF    = IN Turnover rate for leaf and fine root biomass
!             (/360days).
! NPP       = INOUT Net primary productivity (kg C/m2/360days).
! RESP_S    = INOUT Soil respiration (kg C/m2/360days).
! RESP_W    = INOUT Wood maintenance respiration (kg C/m2/360days).
! CS        = INOUT Soil carbon (kg C/m2).
! FRAC      = INOUT Fractional cover of each Functional Type.
! HT        = INOUT Vegetation height (m).
! LAI       = INOUT Leaf area index.
! C_VEG     = OUT Total carbon content of the vegetation (kg C/m2).
! CV        = OUT Gridbox mean vegetation carbon (kg C/m2).
! LIT_C     = OUT Carbon Litter (kg C/m2/360days).
! LIT_C_T   = OUT Gridbox mean carbon litter (kg C/m2/360days).
! DCVEG     = WORK Change in vegetation carbon during the timestep
!             (kg C/m2/timestep).
! DFRAC     = WORK Change in areal fraction during the timestep
!             (/timestep).
! FRAC_FLUX = WORK PFT fraction to be used in the calculation of
!             the gridbox mean fluxes.
! LAI_BAL   = WORK Leaf area index in balanced growth state.
! LEAF      = WORK Leaf biomass (kg C/m2).
! PC_S      = WORK Net carbon flux available for spreading
!             (kg C/m2/yr).
! PHEN      = WORK Phenological state.
! ROOT      = WORK Root biomass (kg C/m2).
! WOOD      = WORK Woody biomass (kg C/m2).
! DFRAC_AGR = WORK Increment to areal fraction from agriculture.
!             (/timestep).
! DFA       = WORK Increment to areal fraction from agriculture that
!             is not burnt. (/timestep).

      real FORW, FRAC_VS(POINTS), GAMMA, FRAC_AGR(POINTS), FRAC_MIN
      real FRAC_SEED, DENOM_MIN, BF, G_LEAF(POINTS,NPFT)
      real NPP(POINTS,NPFT), RESP_S(POINTS), RESP_W(POINTS,NPFT)
      real CS(POINTS), FRAC(POINTS,NTYPE), HT(POINTS,NPFT)
      real LAI(POINTS,NPFT), C_VEG(POINTS,NPFT), CV(POINTS)
      real LIT_C(POINTS,NPFT), LIT_C_T(POINTS), DCVEG(POINTS,NPFT)
      real DFRAC(POINTS,NPFT),  FRAC_FLUX, LAI_BAL(POINTS,NPFT)
      real LEAF(POINTS,NPFT), PC_S(POINTS,NPFT), PHEN(POINTS,NPFT)
      real ROOT(POINTS,NPFT), WOOD(POINTS,NPFT)
      real DFRAC_AGR(POINTS,NPFT), DFA

! SCHDMAX    = IN maximum depth of active soil carbon

      integer SCHDMAX

! DZ_GND    = IN thickness of each ground layer (m)
! CS_PART   = IN partitioning of Litter into subsurface layers
! SCHDMIN   = IN active carbon depth tracer (layer #)
! PFCDENSE  = IN permafrost carbon density (kg C/m3)
! FZN_SL    = IN frozen soil layers (layer #)
! DMIN_HOLD = INOUT holds SCHDMIN from last triffid time-step
! SQU_PFC   = INOUT sequestered permafrost carbon (kg C/m2)
! TRACK_PFC = INOUT tracks permafrost carbon (kg C/m2)
! RESP_D    = INOUT Soil respiration at depth (kg C/m3/360days).
! CS_D      = INOUT Soil carbon at depth (kg C/m3).

      real DZ_GND(NGND), RESP_D(POINTS,NGND), CS_D(POINTS,NGND)
      real SCHDMIN(POINTS), CS_PART(SCHDMAX+1,NGND)
      real DMIN_HOLD(POINTS), PFCDENSE, TRACK_PFC(POINTS)
      real SQU_PFC(POINTS), FZN_SL(POINTS)

! KOVDIFF_P  = IN Parameters for carbon diffusion (m-2)
! KOVCON     = IN Constant of diffution (m2a-1)
! KOVMAX   = IN Maximum depth of mixing KOVMAX*active layer thickness
! YEARKOV_SW = INOUT Swithch for year count
! CS_P = INOUT Depth varying permafrost carbon density (kg m-3)
! RESP_P    = INOUT permafrost carbon respiration at depth (kg C/m3/360days).
! V_SAT     =  IN Saturated soil porosity
! AFp = INOUT Available fraction permafrost carbon
! AF_Ta = INOUT Available fraction transmutation accumulator
! iAF = Initial available fraction
! PREG_C = Soil carbon in permafrost region both in permafrost and active layer

      integer YEARKOV_SW
      real KOVDIFF_P(SCHDMAX,3),KOVCON,KOVMAX,ZBOT(NGND)
      real ACTIVELT(POINTS),PERMA(POINTS),CS_P(POINTS,NGND)
      real RESP_P(POINTS,NGND), V_SAT(POINTS,NGND),KOVSAT
      real AFp(POINTS,NGND),AF_Ta(POINTS,NGND),iAF,PREG_C(POINTS)

!----------------------------------------------------------------------
! Loop through Functional Types
!----------------------------------------------------------------------
      do N=1,NPFT

!----------------------------------------------------------------------
! Loop through TRIFFID points
!----------------------------------------------------------------------
        do T=1,LAND_PTS
          L=LAND_INDEX(T)

!----------------------------------------------------------------------
! Diagnose the balanced-growth leaf area index and the associated leaf,
! wood, root and total vegetation carbon
!----------------------------------------------------------------------
          LAI_BAL(L,N) = (A_WS(N)*ETA_SL(N)*HT(L,N)
     &                   /A_WL(N))**(1.0/(B_WL(N)-1))
          LEAF(L,N) = SIGL(N)*LAI_BAL(L,N)
          ROOT(L,N) = LEAF(L,N)
          WOOD(L,N) = A_WL(N)*(LAI_BAL(L,N)**B_WL(N))
          C_VEG(L,N) = LEAF(L,N) + ROOT(L,N) + WOOD(L,N)

!----------------------------------------------------------------------
! Diagnose the phenological state
!----------------------------------------------------------------------
          PHEN(L,N) = LAI(L,N)/LAI_BAL(L,N)
        enddo

!----------------------------------------------------------------------
! Update vegetation carbon contents
!----------------------------------------------------------------------
        call VEGCARB (LAND_PTS, LAND_INDEX, N, FORW, GAMMA, DENOM_MIN
     &,               G_LEAF(1,N), NPP(1,N), RESP_W(1,N), LEAF(1,N)
     &,               ROOT(1,N), WOOD(1,N), DCVEG(1,N), PC_S(1,N))

      enddo

!-----------------------------------------------------------------------
! Diagnose the new value of Canopy Height, Leaf Area Index and Total
! Vegetation Carbon
!-----------------------------------------------------------------------
      do N=1,NPFT

        do T=1,LAND_PTS
          L=LAND_INDEX(T)

          HT(L,N) = WOOD(L,N)/(A_WS(N)*ETA_SL(N))
     &              *(A_WL(N)/WOOD(L,N))**(1.0/B_WL(N))
          LAI_BAL(L,N) = LEAF(L,N)/SIGL(N)
          LAI(L,N) = PHEN(L,N)*LAI_BAL(L,N)
          C_VEG(L,N) = LEAF(L,N) + ROOT(L,N) + WOOD(L,N)

        enddo

      enddo

!----------------------------------------------------------------------
! Update the areal coverage of each functional type
!----------------------------------------------------------------------
      call LOTKA (LAND_PTS, LAND_INDEX, C_VEG, FORW, FRAC_VS
     &,           FRAC_AGR, FRAC_MIN, FRAC_SEED, DENOM_MIN, GAMMA
     &,           LAI_BAL, PC_S, FRAC, DFRAC, DFRAC_AGR)

!----------------------------------------------------------------------
! Diagnose the litter fall from the carbon balance of each vegetation
! type
!----------------------------------------------------------------------
      do T=1,LAND_PTS
        L=LAND_INDEX(T)
        LIT_C_T(L) = 0.
        do N=1,NPFT

          DFA = DFRAC_AGR(L,N)*(1. - BF)
          FRAC_FLUX = FRAC(L,N) - (1. - FORW)*(DFRAC(L,N) + DFA)
          LIT_C(L,N) = NPP(L,N) - GAMMA/FRAC_FLUX*(C_VEG(L,N)*FRAC(L,N)
     &               - (C_VEG(L,N) - DCVEG(L,N))*(FRAC(L,N)
     &               - DFRAC(L,N) - DFA))

! NM 27.05.19 disable negative Litter fluxes
          LIT_C(L,N) = max(c0,LIT_C(L,N))

          LIT_C_T(L) = LIT_C_T(L) + FRAC_FLUX*LIT_C(L,N)
        enddo
      enddo

!      print*, "LIT_C_T = ", LIT_C_T

!----------------------------------------------------------------------
! Call SOILCARB to update the soil carbon content
!----------------------------------------------------------------------
      call SOILCARB (POINTS, LAND_PTS, LAND_INDEX, FORW, GAMMA

     &,              RESP_D, CS_D, DZ_GND, NGND, SCHDMIN
     &,              CS_PART, SCHDMAX, DMIN_HOLD, PFCDENSE
     &,              TRACK_PFC, SQU_PFC, FZN_SL

     &,              ACTIVELT,KOVDIFF_P,KOVCON,KOVMAX,YEARKOV_SW
     &,              ZBOT,PERMA,CS_P,RESP_P,V_SAT,KOVSAT
     &,              AFp, AF_Ta, iAF, PREG_C

     &,              DENOM_MIN, LIT_C_T, RESP_S, CS)

!----------------------------------------------------------------------
! Diagnose the gridbox mean vegetation carbon
!----------------------------------------------------------------------
      do T=1,LAND_PTS
        L=LAND_INDEX(T)
        CV(L) = 0.0
        do N=1,NPFT
          CV(L) = CV(L) + FRAC(L,N)*C_VEG(L,N)
        enddo
      enddo

      return
      end