! source file: /Users/oschlies/UVIC/master/source/mtlm/triffid.F
      subroutine TRIFFID (LAND_PTS, LAND_INDEX, FORW, GAMMA, FRAC_VS
     &,                   FRAC_AGRIC, FRAC_MIN, FRAC_SEED, DENOM_MIN
     &,                   G_LEAF, NPP, RESP_S, RESP_W, CS, 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.
!**********************************************************************

! based on code by: P. Cox and R. Betts
!-----------------------------------------------------------------------

      implicit none

      include "size.h"
      include "mtlm_data.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_AGRIC = 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.
! 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).

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

!----------------------------------------------------------------------
! 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_AGRIC, FRAC_MIN, FRAC_SEED, DENOM_MIN, GAMMA
     &,           LAI_BAL, PC_S, FRAC, DFRAC)

!----------------------------------------------------------------------
! 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.0

        do N=1,NPFT
          FRAC_FLUX=FRAC(L,N)-(1.0-FORW)*DFRAC(L,N)
          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)))
          LIT_C_T(L) = LIT_C_T(L)+FRAC_FLUX*LIT_C(L,N)
        enddo
      enddo

!----------------------------------------------------------------------
! Call SOIL_C to update the soil carbon content
!----------------------------------------------------------------------
      call SOILCARB (POINTS, LAND_PTS, LAND_INDEX, FORW, GAMMA
     &,              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