! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/mtlm.F
      subroutine MTLM (is, ie, js, je)

!-----------------------------------------------------------------------
! land-surface and vegetation model
! The land model is a modified form of the MOSES and TRIFFID

!**********************************************************************
! 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 "csbc.h"
      include "mtlm.h"
      include "mtlm_data.h"
      include "switch.h"

      integer ie, is, je, js, NOTCROP(NPFT)
      integer I, J, K, L, N, KITER             ! Loop counters

!-----------------------------------------------------------------------
! Driving variables
!-----------------------------------------------------------------------
! LW_OUT_P = WORK Surface longwave (W/m**2)
! LW_OUT_S = WORK Surface longwave (W/m**2)
! PAR      = WORK Photosynthetically active
! SWN_P    = WORK Net shortwave for each PFT
! SWN_S    = WORK Net shortwave for bare soil
! RTSB     = WORK Inverse of the timestep for burning

      real LW_OUT_P(POINTS,NPFT), LW_OUT_S(POINTS), PAR(POINTS)
      real SWN_P(POINTS,NPFT), SWN_S(POINTS)
      real RTSB

!-----------------------------------------------------------------------
! Other variables which need to be saved between calls
!-----------------------------------------------------------------------
! ALBLAND_SOIL = WORK surface albedo
! ALBLAND_VEG  = WORK surface albedo

      real ALBLAND_SOIL(POINTS), ALBLAND_VEG(POINTS,NPFT)

!-----------------------------------------------------------------------
! Work Carbon variables.
!-----------------------------------------------------------------------
! FTIME_PHEN = Weighting factor for accumulations.
! GC         = Canopy conductance (m/s). phenology (/yr).
! LAI_BAL    = Balanced growth LAI.
! LIT_C      = Carbon litter (kg C/m2/yr).
! RESP_P     = Plant respiration rate (kg C/m2/s).
! RESP_W     = Wood maintenance respiration rate (kg C/m2/s).

      real FTIME_PHEN, GC(POINTS,NPFT), LAI_BAL(POINTS,NPFT)

      real RESP_P(POINTS,NPFT), RESP_W(POINTS,NPFT)

!-----------------------------------------------------------------------
! Work Hydrology variables
!-----------------------------------------------------------------------
! ESOIL      = Soil evaporation (kg/m2/s).
! ETRAN      = Transpiration (kg/m2/s).
! G_T        = ground heat flux on veg tiles (W/m2).
! G          = ground heat flux (W/m2).
! GSOIL      = Bare soil ground heat flux (W/m2).
! SH_T       = Sensible heat flux on veg tiles (W/m2
! HSOIL      = Bare soil sensible heat flux (W/m2).
! LESOIL     = Latent heat of soil evaporation (W/m2).
! LETRAN     = Latent heat of transpiration (W/m2).
! RADNET_T   = Net radiation on veg tiles (W/m2).
! RADNET     = Net Radiation (W/m2).
! RADNETSOIL = Bare soil net radiation (W/m2).
! STHU       = Unfrozen soil moisture as a fraction of saturation.
! V_ROOT     = Liquid volumetric soil moisture concentration in the
!              rootzone (m3 H2O/m3 soil).

      real ESOIL(POINTS), ETRAN(POINTS,NPFT), G_T(POINTS,NPFT)
      real G(POINTS), GSOIL(POINTS), SH_T(POINTS,NPFT), HSOIL(POINTS)
      real LESOIL(POINTS), LETRAN(POINTS,NPFT), RADNET_T(POINTS,NPFT)
      real RADNET(POINTS), RADNETSOIL(POINTS), STHU(POINTS)
      real V_ROOT(POINTS)

!-----------------------------------------------------------------------
! Other surface quantities
!-----------------------------------------------------------------------
! RA     = Aerodynamic resistance (s/m).
! RS     = Surface resistance (s/m).
! RSSOIL = Surface resistance for bare soil (s/m).
! V_SAT  = Volumetric soil moisture concentration at saturation
!          (m3 H2O/m3 soil).
! V_WILT = Volumetric soil moisture concentration below which stomata
!          close (m3 H2O/m3 soil).

      real RA(POINTS), RS(POINTS,NPFT), RSSOIL(POINTS)

!-----------------------------------------------------------------------
! Local parameters
!-----------------------------------------------------------------------
! ITER_EQ = Number of TRIFFID iterations for gradient descent to
!           equilibrium.

      integer ITER_EQ
      parameter (ITER_EQ=10)

! DENOM_MIN  = Minimum value for the denominator of the update equation.
!              Ensures that gradient descent does not lead to an
!              unstable solution.
! GAMMA_EQ   = Inverse timestep for gradient descent to equilibrium
!              (/yr).
! CP         = Specific heat of dry air at constant pressure (J/kg/K).
! EPSILON    = Ratio of molecular weights of water and dry air.
! KARMAN     = Von Karman's constant.
! TM         = Melting point of fresh water (K).
! ZERODEGC   = Zero Celsius (K).
! MSAT       = Saturated soil moisture (kg/m2).
! RSS        = Surface resistance for bare soil (s/m
! SATCON     = Saturated hydraulic conductivity (kg/m2/s).
! V_CRIT     = Volumetric soil moisture concentration above which
!              stomata are not sensitive to soil water (m3 H2O/m3 soil).
! VSAT       = Volumetric soil moisture concentration at saturation
!              (m3 H2O/m3 soil).
! VWILT      = Volumetric soil moisture concentration below which
!              stomata close (m3 H2O/m3 soil).
! Z0_SOIL    = Roughness length for bare soil (m).
! Z1_REF     = Reference level height (m).
! DTIME_PHEN = Timestep for phenology (/yr).
! DT         = Timestep for updating state
! FORW       = Forward timestep weighting.
! GAMMA      = Inverse timestep (/yr).
! RHO_WATER  = density of pure water (Kg/m3)
! PI         = PI

      REAL HCAIR, HCICE, HCWAT, THWAT, THICE, STH, HCSAT
      PARAMETER (HCAIR=0.025,HCICE=2.24,HCWAT=0.56)

      real DENOM_MIN ,GAMMA_EQ, CP, EPSILON, KARMAN, TM, RSS, SATCON
      real Z1_REF, DTIME_PHEN, DT, FORW, GAMMA, RHO_WATER, PI

      parameter (DENOM_MIN=1.0E-6, GAMMA_EQ=1.0E-1, CP=1005.0)
      parameter (EPSILON=0.62198, KARMAN=0.4, TM=273.15)
      parameter (RSS=100.0, SATCON=0.0005)
      parameter (Z1_REF=10.0)
      parameter (RHO_WATER = 1000.0, PI=3.14159265358979323846)

!-----------------------------------------------------------------------
! Surface physical variables
!-----------------------------------------------------------------------
! ALB_P    = Surface albedo (PFT)
! ALB_S    = Surface albedo (soil)
! RSUSOIL  = Surface resistance for bare soil evaporation (s/m).
! RSFSOIL  = Surface resistance for bare soil sublimation (s/m).

      real ALB_S(POINTS), ALB_P(POINTS,NPFT)
      real RSUSOIL(POINTS), RSFSOIL(POINTS)

!~ Vegetation (carbon) variables ~
! FRAC_OLD   = Previous vegetation fractions
! VEG_FRAC_OLD = Previous vegetated fraction

      real FRAC_OLD(POINTS, NPFT), VEG_FRAC_OLD(POINTS)

!-----------------------------------------------------------------------
! Soil variables
!-----------------------------------------------------------------------
! ET_P       = Evaporation from plants (transpiration) (kg/m2/s).
! ET_S       = Evaporation from bare soil (kg/m2/s)
! ET_L       = Transpiration from each soil layer
! G_S        = Bare soil ground heat flux (W/m2).
! G_P        = ground heat flux on veg tiles (W/m2).
! L_RATIO    = Layer ratio for evapotranspiration withdrawal
! LH_P       = Latent heat of transpiration (W/m2).
! LH_S       = Latent heat of soil evaporation (W/m2).
! MH_FIX     = Adjustment to ground heat fluxes in the case of negative melt heat
! MH_P       = Energy available to melt snow (PFT) (W/m2).
! MH_S       = Energy available to melt snow (SOIL)(W/m2)
! MH         = Energy available to melt snow (TOTAL GRIDBOX) (W/m2)
! SH_P       = Sensible heat flux on veg tiles (W/m2)
! SH_S       = Bare soil sensible heat flux (W/m2).
! SUBSNOW_P  = Sublimation from lying snow on PFT TILES (kg/m2/s)
! SUBSNOW_S  = Sublimation from lying snow on soil tiles (kg/m2/s)

      real ET_P(POINTS, NPFT), ET_S(POINTS), ET_L(POINTS,NGND)
      real G_P(POINTS,NPFT), G_S(POINTS), SH_P(POINTS,NPFT)
      real SH_S(POINTS), LH_S(POINTS), LH_P(POINTS,NPFT)
      real MH(POINTS), MH_P(POINTS, NPFT), MH_S(POINTS), MH_FIX(POINTS)
      real L_RATIO(POINTS,NPFT,NGND), SUBSNOW_P(POINTS,NPFT)
      real SUBSNOW_S(POINTS), MSAT(POINTS,NGND)
      real TMAX, MU, MELTABLE

!-----------------------------------------------------------------------
! Factors for accumulation variables
!-----------------------------------------------------------------------
      FTIME = TIMESTEP/REAL(SEC_DAY*DAY_TRIF)
      FTIME_PHEN = TIMESTEP/REAL(SEC_DAY*DAY_PHEN)

!-----------------------------------------------------------------------
! Calculate timestep values of the driving data.
!-----------------------------------------------------------------------
      LAND_COUNTER = LAND_COUNTER + 1
      ISTEP = mod(LAND_COUNTER-1,STEP_DAY) + 1

      L_TRIF = .false.
      L_PHEN = .false.
      if (ISTEP .eq. STEP_DAY) then
        L = LAND_COUNTER/STEP_DAY
        if (INT_VEG .and. MOD(L,DAY_TRIF) .eq. 0.) L_TRIF=.true.
        if (MOD(L,DAY_PHEN) .eq. 0.) L_PHEN=.true.
      endif

      TIMEDAY = (REAL(ISTEP) - 0.5)*SEC_DAY/STEP_DAY

      do I=1,LAND_PTS
        L = LAND_INDEX(I)
        SAT_D(L) = T_C(L) + 0.5*DTEMP_DAY(L)
     &         *COS(2*PI*(TIMEDAY-TIME_MAX(L))/SEC_DAY)

        LW(L) = SIGMA*T_C(L)*T_C(L)*T_C(L)*T_C(L)
     &        *(4.0*SAT_D(L)/T_C(L)-3.0)

        SW(L) = SUN(L,ISTEP)
      enddo

      call QSAT (POINTS, LAND_PTS, LAND_INDEX, EPSILON, ZERODEGC
     &,          QS, SAT_D, PSTAR)

!CDIR NODEP
      do I=1,LAND_PTS

        L = LAND_INDEX(I)
        if (RH_C(L) .gt. 0.) then
          Q(L) = QS(L)*RH_C(L)
        else
          Q(L) = 0.
        endif

!-----------------------------------------------------------------------
! Approximate PAR, assume atmosphere is neutral and setup vector fields.
!-----------------------------------------------------------------------
        PAR(L) = 0.45 * SW(L)
        Z0S(L) = Z0_SOIL

        RA(L) = 0.0
        HT_S(L) = 0.0
        SNOW_HT(L) = LYING_SNOW(L)/RHO_S

!----------------------------------------------------------------------
! Calculate the soil moisture stress factor - bare soil
!----------------------------------------------------------------------
! Bare soil points will experience evaporation

        RSUSOIL(L) = 0.
        RSFSOIL(L) = 0.
        if (V_SAT(L,1) .gt. 0) then
          if (SU(L,1) .gt. 0) then
            RSUSOIL(L) = RSS*(V_CRIT(L,1)/(SU(L,1)*V_SAT(L,1)))**2
          endif
          if (SF(L,1) .gt. 0) then
            RSFSOIL(L) = RSS*(V_CRIT(L,1)
     &                   /(SF(L,1)*RHO_W/RHO_I*V_SAT(L,1)))**2
          endif
        endif

        if (RSUSOIL(L) .gt. 10**6) RSUSOIL(L) = 10**6
        if (RSFSOIL(L) .gt. 10**6) RSFSOIL(L) = 10**6

        do N = 1, NSOIL(L)
          MSAT(L,N) = RHO_W*DZ_GND(N)*V_SAT(L,N)
        enddo

      enddo

!-----------------------------------------------------------------------
! Calculate the soil moisture stress factor.
!-----------------------------------------------------------------------
      call STRESSFACTOR (POINTS, LAND_PTS,LAND_INDEX, NPFT, NGND
     &,                  NSOIL, F_ROOT,  FSMC_PFT, V_CRIT, V_WILT
     &,                  DZ_GND, V_SAT, SU, ROOTDEP, PERMA, L_RATIO
     &,                  ZTOP, ZBOT)

!-----------------------------------------------------------------------
! Calculate the soil respiration
!-----------------------------------------------------------------------
      call MICROBE (POINTS, LAND_PTS, LAND_INDEX, NGND, DZ_GND
     &,             NSOIL,  CS_D, SU, V_SAT, V_WILT
     &,             C_WEIGHT, TGND, RESP_S, ZTOP, ZBOT, RESP_SD

     &,             RESP_SP, CS_P, KAPSp, AFp, PPtm, AF_Ta, iAF

     &             )

!-----------------------------------------------------------------------
! Calculate the bare soil evaporation.
!-----------------------------------------------------------------------
      call SWRAD (POINTS, LAND_PTS, LAND_INDEX, ALBSOIL, ALBSNOW
     &,           LYING_SNOW, SW, TSTAR_S, ZERODEGC, SWN_S, ALB_S)

      do N=1,NPFT
!-----------------------------------------------------------------------
! Calculate the albedo of the plant functional types
!-----------------------------------------------------------------------
        call SWRAD (POINTS, LAND_PTS, LAND_INDEX, ALBSNF(1,N)
     &,             ALBSNC(1,N), LYING_SNOW, SW, TSTAR_P(1,N), ZERODEGC
     &,             SWN_P(1,N), ALB_P(1,N))

!-----------------------------------------------------------------------
! Calculate the canopy resistance and the resistance factor.
!-----------------------------------------------------------------------
        call SF_STOM (LAND_PTS, LAND_INDEX, N, CO2, FSMC_PFT(1,N)
     &,               HT(1,N), PAR, LAI(1,N), PSTAR, Q, RA, SAT_D
     &,               ZERODEGC, EPCO2, EPSILON, GPP(1,N), NPP(1,N)
     &,               RESP_P(1,N), RESP_W(1,N), GC(1,N))

        do I=1,LAND_PTS
          L = LAND_INDEX(I)
          RS(L,N) = 1.0 / GC(L,N)
        enddo

!-----------------------------------------------------------------------
! Calculate the leaf turnover rate.
!-----------------------------------------------------------------------
        call LEAF_LIT (LAND_PTS, LAND_INDEX, N
     &,                FSMC_PFT(1,N), TSTAR_P(1,N), G_LEAF(1,N))

!-----------------------------------------------------------------------
! Update leaf phenological state
!-----------------------------------------------------------------------
        do I=1,LAND_PTS
          L = LAND_INDEX(I)
          G_LEAF_DAY(L,N) = G_LEAF_DAY(L,N) + G_LEAF(L,N)*FTIME_PHEN
        enddo
      enddo ! PFT LOOP

      if (L_PHEN) then

        if (DAY_PHEN .gt. 1) then
          DTIME_PHEN = real(DAY_PHEN)*segtim/DAY_YEAR
        else
          DTIME_PHEN = real(DAY_PHEN)/DAY_YEAR
        endif

        do N=1,NPFT
          call PHENOL (LAND_PTS, LAND_INDEX, N, G_LEAF_DAY(1,N)
     &,                HT(1,N), DTIME_PHEN, G_LEAF_PHEN(1,N), LAI(1,N))
          do I=1,LAND_PTS
            L = LAND_INDEX(I)
            G_LEAF_DAY(L,N) = 0.0
          enddo
        enddo
      endif ! End of PHENOL call

      DT = TIMESTEP*segtim

!-----------------------------------------------------------------------
! Solve the surface energy balance.
!-----------------------------------------------------------------------

      do I=1,LAND_PTS
        L = LAND_INDEX(I)
        IF (SU(L,1) .GT. 0.) THEN
          THWAT = V_SAT(L,1)*SU(L,1)/(SU(L,1) + SF(L,1))
        ELSE
          THWAT = 0.
        ENDIF
        IF (SF(L,1) .GT. 0.0) THEN
          THICE = V_SAT(L,1)*SF(L,1)/(SU(L,1) + SF(L,1))
        ELSE
          THICE = 0.
        ENDIF
        STH = SU(L,1) + SF(L,1)
        HCSAT = HCON_D(L,1)*(HCWAT**THWAT)*(HCICE**THICE)
     &          /(HCAIR**V_SAT(L,1))
        HCON_B(L,1) = ((HCSAT - HCON_D(L,1))*STH + HCON_D(L,1))
      enddo

      do N=1,NPFT
        call PENMON (POINTS, LAND_PTS, LAND_INDEX, NGND, DT
     &,              RHO_S, DZ_GND, HCON_B(1,1), RS(1,N),RSFSOIL
     &,              Z0(1,N), LW, SWN_P(1,N), PSTAR, Q, SAT_D, TGND
     &,              WIND, LHC, LHF, EPSILON, SIGMA, ZERODEGC
     &,              ET_P(1,N), SUBSNOW_P(1,N), SUBSOIL, LH_P(1,N)
     &,              SH_P(1,N), G_P(1,N), TSTAR_P(1,N), LW_OUT_P(1,N)
     &,              SU, SF, MH_P(1,N), SNOW_HT, HT(1,N),HCON_SNOW)
      enddo ! FT Loop

!-----------------------------------------------------------------------
! Solve the soil surface energy balance.
!-----------------------------------------------------------------------
      call PENMON (POINTS ,LAND_PTS, LAND_INDEX, NGND, DT, RHO_S
     &,            DZ_GND, HCON_B(1,1), RSUSOIL, RSFSOIL, Z0S, LW
     &,            SWN_S, PSTAR, Q, SAT_D, TGND, WIND, LHC, LHF
     &,            EPSILON, SIGMA, ZERODEGC, ET_S, SUBSNOW_S, SUBSOIL
     &,            LH_S, SH_S, G_S, TSTAR_S, LW_OUT_S, SU, SF,MH_S
     &,            SNOW_HT, HT_S,HCON_SNOW)

!----------------------------------------------------------------------
! Work out the evapotranspiration to come out of each layer
! Note that this accounts for the gridbox mean evapotranspiration and
! latent heat fluxes
!----------------------------------------------------------------------
      ET(:) = 0.
      LE(:) = 0.
      ET_L(:,:) = 0.
      SUBSNOW(:) = 0.

      do I=1,LAND_PTS
        L = LAND_INDEX(I)
        do N=1,NPFT
          ET(L) = ET(L) + FRAC(L,N) * ET_P(L,N)
          LE(L) = LE(L) + FRAC(L,N) * LH_P(L,N)
          SUBSNOW(L) = SUBSNOW(L) + SUBSNOW_P(L,N)*FRAC(L,N)
          do K=1,NSOIL(L)
            ET_L(L,K) = ET_L(L,K)
     &                + FRAC(L,N)*L_RATIO(L,N,K)*ET_P(L,N)
          enddo
        enddo
        ET(L) = ET(L) + (1. - VEG_FRAC(L))*(ET_S(L))
        ET_L(L,1) = ET_L(L,1) + (1. -VEG_FRAC(L))*ET_S(L)
        LE(L) = LE(L) + (1. - VEG_FRAC(L))* LH_S(L)
        SUBSOIL(L) = (1. - VEG_FRAC(L))*SUBSOIL(L)

!---------------------------------------------------------------------------------
! Check to ensure we are not evaporating water that's not there
!---------------------------------------------------------------------------------
        do K=1,NSOIL(L)
          MU = SU(L,K)*MSAT(L,K)
          if (ET_L(L,K) .gt. 0.25*MU/DT) then
            ET_L(L,K+1) = ET_L(L,K+1) + ET_L(L,K)
     &                     - (0.25*MU/DT)
            ET_L(L,K) = 0.25*MU/DT
          endif
        enddo
        SUBSNOW(L) = SUBSNOW(L) + (1. - VEG_FRAC(L))*SUBSNOW_S(L)
        MELTABLE = LYING_SNOW(L)/DT + SNOW(L)
        MH_FIX(L) = 0.
        if (SUBSNOW(L) .gt. MELTABLE) then
          MH_FIX(L) = (SUBSNOW(L) - MELTABLE)*(LHF+LHC)
          LE(L) = LE(L) - MH_FIX(L)
          SUBSNOW(L) = MELTABLE
        endif

      enddo

!-----------------------------------------------------------------------
! Calculate gridbox mean fluxes and surface temperature.
!-----------------------------------------------------------------------
!!  !CDIR NODEP
      do I=1,LAND_PTS
        L = LAND_INDEX(I)
        SH(L) = 0.0
        G(L) = 0.0
        TSTAR_GB(L) = 0.0
        ALBLAND(L) = 0.0
        LW_OUT(L) = 0.0
        SWN(L) = 0.0
        MH(L) = 0.0
!!  !CDIR unroll=5
        do N=1,NPFT
          FRAC_OLD(L,N) = FRAC(L,N)
          TSTAR_GB(L) = TSTAR_GB(L) + FRAC(L,N)*TSTAR_P(L,N)
          SH(L) = SH(L) + FRAC(L,N)*SH_P(L,N)
          G(L) = G(L) + FRAC(L,N)*G_P(L,N)
          ALBLAND(L) = ALBLAND(L) + FRAC(L,N)*ALB_P(L,N)
          LW_OUT(L) = LW_OUT(L) + FRAC(L,N)*LW_OUT_P(L,N)
          SWN(L) = SWN(L) + FRAC(L,N)*SWN_P(L,N)
          MH(L) = MH(L) + FRAC(L,N)*MH_P(L,N)
        enddo
        VEG_FRAC_OLD(L) = VEG_FRAC(L)
        ALBLAND(L) = ALBLAND(L) + (1 - VEG_FRAC(L))*ALB_S(L)
        LW_OUT(L) = LW_OUT(L) + (1 - VEG_FRAC(L))*LW_OUT_S(L)
        SWN(L) = SWN(L) + (1 - VEG_FRAC(L))*SWN_S(L)
        TSTAR_GB(L) = TSTAR_GB(L) + (1 - VEG_FRAC(L)) * TSTAR_S(L)
        SH(L) = SH(L) + (1 - VEG_FRAC(L)) * SH_S(L)
        G(L) = G(L) + (1 - VEG_FRAC(L)) * G_S(L) + MH_FIX(L)
        MH(L) = MH(L) + (1-VEG_FRAC(L))*MH_S(L)
      enddo

!-----------------------------------------------------------------------
! Update the land surface state
!-----------------------------------------------------------------------

      call MTLM_STATE (DT, G, ET_L, MH)

!     set DT back to the original timestep. This is only done for
!     consistency with previous code. The value of DT is irrelevant
!     after this point (as long as it is not zero) since it is only
!     used to multiply fluxes but then divided out when again when
!     calculating flux averages.
      DT = TIMESTEP

!---------------------------------------------------------------------
! Form the mean for variables given to the atmospheric and ocean model
!---------------------------------------------------------------------
      do j=2,jmt-1
        do i=2,imt-1
          L = land_map(i,j)
          if (L .ne. 0) then
            sbc(i,j,iro) = sbc(i,j,iro) + SURF_ROFF(L)*DT*0.1
            sbc(i,j,isca) = sbc(i,j,isca) + (1. - ALBLAND(L))*DT
            sbc(i,j,ievap) = sbc(i,j,ievap)

     &                     + 0.1*(ET(L) + SUBSNOW(L) + SUBSOIL(L))*DT

            sbc(i,j,ilwr) = sbc(i,j,ilwr) - 1000.*(LW_OUT(L) - SW_C(L)
     &                    + SWN(L))*DT
            sbc(i,j,isens) = sbc(i,j,isens) + 1000.*SH(L)*DT

            do N=1,NPFT
              sbc(i,j,inpp) = sbc(i,j,inpp)
     &                      + NPP(L,N)*FRAC(L,N)*DT
            enddo
            sbc(i,j,isr) = sbc(i,j,isr) + RESP_S(L)*DT

          endif
        enddo
      enddo
      atlnd = atlnd + DT

!-----------------------------------------------------------------------
! Accumulate TRIFFID driving variables
!-----------------------------------------------------------------------
      do I=1,LAND_PTS
        L = LAND_INDEX(I)
	RESP_S_DR(L) = RESP_S_DR(L)+RESP_S(L)*FTIME*SEC_YEAR

	do N=1,NGND
           RESP_S_DD(L,N) = RESP_S_DD(L,N)+RESP_SD(L,N)*FTIME*SEC_YEAR

           RESP_S_DP(L,N) = RESP_S_DP(L,N)+RESP_SP(L,N)*FTIME*SEC_YEAR
	   AF_DTa(L,N) = AF_DTa(L,N)+AF_Ta(L,N)*FTIME*SEC_YEAR

	enddo

        do N=1,NPFT
          NPP_DR(L,N) = NPP_DR(L,N) + NPP(L,N)*FTIME*SEC_YEAR
          G_LEAF_DR(L,N) = G_LEAF_DR(L,N) + G_LEAF_PHEN(L,N)*FTIME
          RESP_W_DR(L,N) = RESP_W_DR(L,N) + RESP_W(L,N)*FTIME*SEC_YEAR
        enddo
      enddo

!----------------------------------------------------------------------
! Update the vegetation areal coverages, structural parameters,
! and soil carbon.
!----------------------------------------------------------------------
      if (L_TRIF) then

        if (VEG_EQUIL) then
          FORW = 1.0
          GAMMA = GAMMA_EQ
          KITER = ITER_EQ
        else
          FORW = 0.0

          if (DAY_TRIF .gt. 1) then
            GAMMA = DAY_YEAR/(segtim*real(DAY_TRIF))
          else
            GAMMA = DAY_YEAR/real(DAY_TRIF)
          endif

          KITER = 1
        endif

!       diagnose the total land carbon before triffid
        do j=2,jmt-1
          do i=2,imt-1
            L = land_map(i,j)
            if (L .ne. 0) then

              sbc(i,j,iburn) = CS(L) - RESP_S_DR(L)/GAMMA

              do N=1,NPFT
                sbc(i,j,iburn) = sbc(i,j,iburn) + FRAC(L,N)*C_VEG(L,N)
     &                         + FRAC(L,N)*NPP_DR(L,N)/GAMMA
              enddo
            endif
          enddo
        enddo

        do I=1,LAND_PTS
          L = LAND_INDEX(I)
	  FZN_SL(L)=0.
          do N=1,SCHDMAX
	    if(TGND(L,N).gt.274.15)then
	      FZN_SL(L)=REAL(N)
	    endif
          enddo
        enddo

        do I=1,KITER
          call TRIFFID (LAND_PTS, LAND_INDEX, FORW, GAMMA, FRAC_VS
     &,                 FRACA, FRAC_MIN, FRAC_SEED, DENOM_MIN, BF
     &,                 G_LEAF_DR, NPP_DR, RESP_S_DR, RESP_W_DR, CS

     &,                 RESP_S_DD, 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_S_DP,V_SAT,KOVSAT
     &,                 AFp, AF_DTa, iAF, PREG_C

     &,                 FRAC, HT, LAI, C_VEG, CV, LIT_C, LIT_C_T)
        enddo
        LAND_COUNTER = 0.

!----------------------------------------------------------------------
! Diagnose the amount of vegetation burnt (land use change emissions)
!----------------------------------------------------------------------

        RTSB = 1./(DAY_TRIF*segtim*SEC_DAY)

        do j=2,jmt-1
          do i=2,imt-1
            L = land_map(i,j)
            if (L .ne. 0) then

              sbc(i,j,iburn) = sbc(i,j,iburn) - CS(L)

              do N=1,NPFT
                sbc(i,j,iburn) = sbc(i,j,iburn) - FRAC(L,N)*C_VEG(L,N)
              enddo

              sbc(i,j,iburn) =  sbc(i,j,iburn)*RTSB
            else
              sbc(i,j,iburn) = 0.

            endif
          enddo
        enddo

!----------------------------------------------------------------------
! Zero the accumulated driving variables
!----------------------------------------------------------------------
        do I=1,LAND_PTS
          L = LAND_INDEX(I)
          RESP_S_DR(L) = 0.0

	  RESP_S_DD(L,:) = 0.0

	  RESP_S_DP(L,:) = 0.0
	  AF_DTa(L,:) = 0.0

          do N=1,NPFT

            NPP_DR(L,N) = 0.0
            G_LEAF_DR(L,N) = 0.0
            RESP_W_DR(L,N) = 0.0
          enddo
        enddo

!----------------------------------------------------------------------
! Derive vegetation parameters from the new areal fractions and
! structural properties.
!----------------------------------------------------------------------
        do N=1,NPFT

          call PFT_SPARM (LAND_PTS, LAND_INDEX, N, ALBSOIL
     &,                   HT(1,N), LAI(1,N), ROOTDEP(1,N)
     &,                   INFILFAC(1,N), Z0_SOIL, SNOW_HT
     &,                   ALBSNC(1,N), ALBSNF(1,N), CATCH(1,N)
     &,                   Z0(1,N))

        enddo

!----------------------------------------------------------------------
! Define other vegetation parameters
!----------------------------------------------------------------------
        do I=1,LAND_PTS
          L = LAND_INDEX(I)
          VEG_FRAC(L) = 0.0
          do N=1,NPFT
            VEG_FRAC(L) = VEG_FRAC(L) + FRAC(L,N)
          enddo

          FRAC(L,SOIL) = 1 - VEG_FRAC(L)

          FRAC_VS(L) = VEG_FRAC(L) + FRAC(L,SOIL)
        enddo

      endif

!-----------------------------------------------------------------------
!     accumulate time averages
!-----------------------------------------------------------------------
      if (timavgperts) call ta_mtlm_tavg (is, ie, js, je, 1)

!-----------------------------------------------------------------------
!     accumulate time step integrals
!-----------------------------------------------------------------------
      if (tsiperts) call ta_mtlm_tsi (is, ie, js, je, 1)

      return
      end