! *****************************COPYRIGHT*******************************
! (C) Crown copyright Met Office. All rights reserved.
! For further details please refer to the file COPYRIGHT.txt
! which you should have received as part of this distribution.
! *****************************COPYRIGHT*******************************
! Version 2A of vegetation section: models leaf phenology and vegetation
! competition

! Subroutine Interface:
SUBROUTINE veg2(                                                              &
               land_pts,land_index,nsurft,can_model                           &
,              a_step,asteps_since_triffid                                    &
,              phenol_period,triffid_period                                   &
,              l_phenol,l_triffid,l_trif_eq                                   &
,              atimestep,fraca,fracp,satcon                                   &
,              g_leaf_ac,g_leaf_phen_ac,npp_ac                                &
,              resp_s_ac,resp_w_ac                                            &
,              cs, frac, lai, clay_gb, z0m_soil, ht                         &
,              catch_s,catch_t,infil_t,z0_t, z0h_t                            &
,              c_veg,cv,lit_c,lit_c_mn,g_leaf_day,g_leaf_phen                 &
,              lai_phen,g_leaf_dr_out,npp_dr_out,resp_w_dr_out                &
,              resp_s_dr_out                                                  &
               )

USE atm_fields_bounds_mod
USE theta_field_sizes, ONLY : t_i_length
USE jules_surface_types_mod, ONLY : npft, nnpft, ncpft, ntype, soil

USE jules_vegetation_mod, ONLY : frac_min, l_crop, l_nitrogen,                &
                                 l_trif_crop, l_trait_phys

USE descent
USE veg_param, ONLY : agric
USE missing_data_mod, ONLY : rmdi

use jules_surface_mod, ONLY : l_aggregate, i_aggregate_opt, cmass

USE parkind1, ONLY: jprb, jpim
USE yomhook, ONLY: lhook, dr_hook

USE pftparm, ONLY: a_ws, eta_sl, a_wl, b_wl, lma, sigl

USE trif_vars_mod, ONLY: cnsrv_carbon_veg2_gb

USE prognostics, ONLY: wood_prod_fast_gb, wood_prod_med_gb,                   &
                       wood_prod_slow_gb

USE crop_vars_mod, ONLY : rootc_cpft, harvc_cpft, reservec_cpft,              &
                          stemc_diag_cpft, leafc_diag_cpft

USE jules_print_mgr, ONLY :                                                   &
    jules_message,                                                            &
    jules_print
IMPLICIT NONE

! Description:
!   Updates Leaf Area Index for Plant Functional Types (PFTs) and uses
!   this to derive new vegetation parameters for PFTs along with gridbox
!   mean values where appropriate.

! Method:
!   Calls PHENOL which models phenolgy and updates Leaf Area Index
!   (LAI), then calls TRIFFID to update vegetation and soil fractions,
!   LAI, canopy height, veg and soil carbon and carbon fluxes.  Passes
!   fractions, LAI and canopy height to SPARM which derives the
!   vegetation parameters for each PFT and also the gridbox means where
!   this is required.

INTEGER                                                                       &
 land_pts                                                                     &
                       ! IN Number of land points.
,land_index(land_pts)                                                         &
                       ! IN Index of land points
,nsurft                                                                       &
                       ! IN Number of land-surface tiles.
,can_model                                                                    &
                              ! IN Swith for thermal vegetation
,a_step                                                                       &
                       ! IN Atmospheric timestep number.
,asteps_since_triffid                                                         &
                       ! INOUT Number of atmosphere
!                                    timesteps since last call
!                                           to TRIFFID.
,phenol_period                                                                &
                       ! IN Phenology period (days).
,triffid_period        ! IN TRIFFID period (days).

INTEGER                                                                       &
 i,j,k,l,n                                                                    &
                              ! WORK loop counters.
,kiter                        ! WORK Number of TRIFFID iterations.

LOGICAL                                                                       &
 l_phenol                                                                     &
                              ! IN .T. for interactive leaf
!                                   !    phenology.
,l_triffid                                                                    &
                              ! IN .T. for interactive vegetation.
,l_trif_eq                    ! IN .T. for vegetation equilibrium.

REAL                                                                          &
 atimestep                                                                    &
                              ! IN Atmospheric timestep (s).
,fraca(land_pts)                                                              &
                              ! IN Fraction of agriculture.
,fracp(land_pts)                                                              &
                              ! IN Fraction of pasture.
,satcon(land_pts)                                                             &
                              ! IN Saturated hydraulic
!                                   !    conductivity (kg/m2/s).
,g_leaf_ac(land_pts,npft)                                                     &
                              ! INOUT Accumulated leaf turnover
!                                   !       rate.
,g_leaf_phen_ac(land_pts,npft)                                                &
                              ! INOUT Accumulated leaf turnover
!                                   !       rate including phenology.
,npp_ac(land_pts,npft)                                                        &
                              ! INOUT Accumulated NPP (kg C/m2).
,resp_w_ac(land_pts,npft)                                                     &
                              ! INOUT Accumulated wood respiration
!                                   !       (kg C/m2).
,resp_s_ac(land_pts,4)                                                        &
                           ! INOUT Accumulated soil respiration
!                                       !       (kg C/m2).
,cs(land_pts,4)                                                               &
                           ! INOUT Soil carbon content
,cs_tot(land_pts)                                                             &
                           ! INOUT Soil carbon content
!                                   !       (kg C/m2).
,frac(land_pts,ntype)                                                         &
                              ! INOUT Fractions of surface types.
,lai(land_pts,npft)                                                           &
                              ! INOUT LAI of plant functional
!                                   !       types.
,ht(land_pts,npft)                                                            &
                              ! INOUT Height of plant functional
!                                   !       types (m).
,catch_s(land_pts,nsurft)                                                     &
                              ! OUT Snow capacity for tiles
!                                   !     (kg/m2).
,catch_t(land_pts,nsurft)                                                     &
                              ! OUT Canopy capacity for tiles
!                                   !     (kg/m2).
,infil_t(land_pts,nsurft)                                                     &
                              ! OUT Maximum surface infiltration
!                                   !     rate for tiles (kg/m2/s).
,z0_t(land_pts,nsurft)                                                        &
                              ! OUT Roughness length for tiles (m)
,z0h_t(land_pts,nsurft)                                                       &
                              ! OUT Thermal roughness length for
                              !     tiles (m)
,c_veg(land_pts,npft)                                                         &
                              ! OUT Total carbon content of
!                                   !     the vegetation (kg C/m2).
,cv(land_pts)                                                                 &
                              ! OUT Gridbox mean vegetation
!                                   !     carbon (kg C/m2).
,g_leaf_day(land_pts,npft)                                                    &
                              ! OUT Mean leaf turnover rate for
!                                   !      input to PHENOL (/360days).
,g_leaf_dr_out(land_pts,npft)                                                 &
                              ! OUT Mean leaf turnover rate for
!                                   !       driving TRIFFID (/360days).
,lai_phen(land_pts,npft)                                                      &
                              ! OUT LAI of PFTs after phenology.
,lit_c(land_pts,npft)                                                         &
                              ! OUT Carbon Litter
!                                   !     (kg C/m2/360days).
,lit_c_mn(land_pts)                                                           &
                              ! OUT Gridbox mean carbon litter
!                                   !     (kg C/m2/360days).
,npp_dr_out(land_pts,npft)                                                    &
                              ! OUT Mean NPP for driving TRIFFID
!                                   !     (kg C/m2/360days).
,resp_w_dr_out(land_pts,npft)                                                 &
                              ! OUT Mean wood respiration for
!                                   !       driving TRIFFID
!                                   !       (kg C/m2/360days).
,resp_s_dr_out(land_pts,5)                                                    &
                           ! OUT Mean soil resp for
!                                   ! NB 5=dim_cs1+1. The 5th element
!                                   ! is used as workspace.
,clay_gb(land_pts)                                                            &
                                  ! IN  Clay fraction of soil

,resp_frac(land_pts)                                                          &
                                   ! respired fraction of RESP_S
                                   ! driving TRIFFID
                                   ! (kg C/m2/360days).
,z0m_soil(land_pts)

INTEGER                                                                       &
 nstep_phen                                                                   &
                              ! WORK Number of atmospheric
!                                   !      timesteps between calls to
!                                   !      PHENOL.
,nstep_trif                                                                   &
                              ! WORK Number of atmospheric
!                                   !      timesteps between calls to
!                                   !      TRIFFID.
,surft_pts(ntype)                                                             &
                              ! WORK Number of land points which
!                                   !      include the nth surface type.
,surft_index(land_pts,ntype)                                                  &
                              ! WORK Indices of land points which
!                                   !      include the nth surface type.
,trif_pts                                                                     &
                              ! WORK Number of points on which
!                                   !      TRIFFID may operate
,trif_index(land_pts)                                                         &
                              ! WORK Indices of land points on
!                                   !      which TRIFFID may operate
,t
                              ! WORK loop counter

REAL                                                                          &
 dtime_phen                                                                   &
                              ! WORK The phenology timestep (yr).
,forw                                                                         &
                              ! WORK Forward timestep weighting
!                                   !      for TRIFFID.
,GAMMA                                                                        &
                              ! WORK Inverse TRIFFID timestep
!                                   !      (/360days).
,gam_trif                                                                     &
                              ! WORK Inverse TRIFFID coupling
!                                   !      timestep (/360days).
,ratio                                                                        &
                              ! WORK Ratio of fractional
!                                   !      coverage before to that
!                                   !      after TRIFFID.
,dcs(land_pts)                                                                &
                              ! WORK Change in soil carbon
!                                   !      (kg C/m2).
,dcs_rothc(land_pts,5)                                                        &
                              ! WORK Change in soil carbon
!                                   !      (kg C/m2).
,frac_agric(land_pts)                                                         &
                              ! WORK Fraction of agriculture as
!                                   !      seen by TRIFFID.
,frac_past(land_pts)                                                          &
                              ! WORK Fraction of pasture as
!                                   !      seen by TRIFFID.
,frac_old(land_pts,ntype)                                                     &
                              ! WORK Fractions of surface types
!                                   !      before the call to TRIFFID.
,frac_vs(land_pts)                                                            &
                              ! WORK Total fraction of gridbox
!                                   !      covered by veg or soil.
,g_leaf_phen(land_pts,npft)                                                   &
                              ! WORK Mean leaf turnover rate over
!                                   !      phenology period (/360days).
,g_leaf_dr(land_pts,npft)                                                     &
                              ! WORK Mean leaf turnover rate
!                                   !      for driving TRIFFID
!                                   !      (/360days).
,npp_dr(land_pts,npft)                                                        &
                              ! WORK Mean NPP for driving
!                                   !      TRIFFID (kg C/m2/360days).
,resp_w_dr(land_pts,npft)                                                     &
                              ! WORK Mean wood respiration for
!                                   !      driving TRIFFID
!                                   !      (kg C/m2/360days).
,resp_s_dr(land_pts,5)     ! WORK Mean soil resp for
!                                   !      driving TRIFFID
!                                   !      (kg C/m2/360days).
!                                   ! NB 5=dim_cs1+1. The 5th element
!                                   ! is used as workspace. Not sure it
!                                   ! is even needed.

REAL                                                                          &
 cv_trif(land_pts)
                              ! WORK Gridbox mean natural pft vegetation
!                             !      carbon (kg C/m2).
!                             !      Not currently used

REAL                                                                          &
lai_work                                                                      &
,leaf_work                                                                    &
                  ! WORK used to calculate initial veg carbon for
                  !      cnsrv_carbon_veg2 diagnostic
,wood_prod_total(land_pts)                                                    &
                  ! WORK total of the wood product pools
,wp_total_out(land_pts)
                  ! WORK total flux out of wood products

REAL cnsrv_carbon_flux(land_pts)
                  ! WORK net carbon flux into land
                  !      (vegetation, soil, wood products)
                  !      used to assess carbon conservation

REAL,PARAMETER  :: resp_frac_a=4.0895
REAL,PARAMETER  :: resp_frac_b=2.672
REAL,PARAMETER  :: resp_frac_c=-0.0786
                  ! parameters used to calculate resp_frac from clay_gb

INTEGER(KIND=jpim), PARAMETER :: zhook_in  = 0
INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
REAL(KIND=jprb)               :: zhook_handle

CHARACTER(LEN=*), PARAMETER :: RoutineName='VEG2'

IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)

!-----------------------------------------------------------------------
! Initialisations
!-----------------------------------------------------------------------
DO n=1,npft
  DO l=1,land_pts
    g_leaf_phen(l,n)=0.0
    g_leaf_day(l,n)=0.0
    g_leaf_dr(l,n)=0.0
    npp_dr(l,n)=0.0
    resp_w_dr(l,n)=0.0
    c_veg(l,n)=0.0
    lit_c(l,n)=0.0
  END DO
END DO

DO n=1,nsurft
  DO l=1,land_pts
    catch_t(l,n)=0.0
    infil_t(l,n)=0.0
    z0_t(l,n)=0.0
  END DO
END DO

IF (l_aggregate .AND. i_aggregate_opt==1) z0h_t(:,:)=0.0

IF (can_model == 4) THEN
  DO n=1,nsurft
    DO l=1,land_pts
      catch_s(l,n)=0.0
    END DO
  END DO
END IF

DO l=1,land_pts
  cv(l)=0.0
  lit_c_mn(l)=0.0
  frac_vs(l) = 0.0
  frac_agric(l) = 0.0
  frac_past(l) = 0.0
END DO
DO n=1,5
  DO l=1,land_pts
    resp_s_dr(l,n)=0.0
  END DO
END DO

!-----------------------------------------------------------------------
! Calculate the number of atmospheric timesteps between calls to PHENOL
! and TRIFFID.
!-----------------------------------------------------------------------
nstep_phen=INT(86400.0*phenol_period/atimestep)
nstep_trif=INT(86400.0*triffid_period/atimestep)


!-----------------------------------------------------------------------
! Find total fraction of gridbox covered by vegetation and soil, and use
! this to set indices of land points on which TRIFFID may operate
!-----------------------------------------------------------------------
trif_pts = 0
DO l=1,land_pts
  DO n=1,nnpft
    frac_vs(l) = frac_vs(l) + frac(l,n)
  END DO
  n=soil
  frac_vs(l) = frac_vs(l) + frac(l,n)
  IF (frac_vs(l) >= (nnpft*frac_min)) THEN
    trif_pts = trif_pts + 1
    trif_index(trif_pts) = l
  END IF
END DO

!----------------------------------------------------------------------
! Check carbon conservation (1/3)
! Calculate total land carbon store and net carbon flux at start of
! routine. The wood product decay flux has not yet been calculated and
! will be added later.
!----------------------------------------------------------------------
DO t=1,trif_pts
  l=trif_index(t)
  cnsrv_carbon_veg2_gb(l)=0.0
  DO n=1,nnpft
     lai_work = (a_ws(n)*eta_sl(n)*ht(l,n)                                    &
                /a_wl(n))**(1.0/(b_wl(n)-1))
    IF(l_trait_phys) THEN
       leaf_work = cmass*lma(n)*lai_work
    ELSE
       leaf_work = sigl(n)*lai_work
    END IF
    cnsrv_carbon_veg2_gb(l) = cnsrv_carbon_veg2_gb(l) +                       &
       frac(l,n)*(2.0*leaf_work + a_wl(n)*(lai_work**b_wl(n)))
  END DO

  cnsrv_carbon_veg2_gb(l)= cnsrv_carbon_veg2_gb(l) +                          &
        SUM(cs(l,:))+                                                         &
        wood_prod_fast_gb(l)+wood_prod_med_gb(l)+wood_prod_slow_gb(l)

  resp_frac(l) = 1.0 / (resp_frac_a + resp_frac_b*                            &
                   EXP(resp_frac_c * 100.0*clay_gb(l)))
  cnsrv_carbon_flux(l) =                                                      &
        SUM(npp_ac(l,1:npft)*frac(l,1:npft))-                                 &
        (1.0-resp_frac(l))*SUM(resp_s_ac(l,1:4))
END DO

!-----------------------------------------------------------------------
! Create the surft_index array of land points with each surface type
!-----------------------------------------------------------------------
! DEPENDS ON: tilepts
CALL tilepts(land_pts,frac,surft_pts,surft_index)

IF (l_phenol .AND. MOD(a_step,nstep_phen) == 0) THEN

!-----------------------------------------------------------------------
! Calculate the phenology timestep in years.
!-----------------------------------------------------------------------
  dtime_phen=FLOAT(phenol_period)/360.0

  DO n=1,nnpft

!-----------------------------------------------------------------------
! Calculate the mean turnover rate and update the leaf phenological
! state, and take copy of updated LAI field for output as diagnostic.
!-----------------------------------------------------------------------
    DO j=1,surft_pts(n)
      l=surft_index(j,n)
      g_leaf_day(l,n)=g_leaf_ac(l,n)/dtime_phen
    END DO

    WRITE(jules_message,*) 'Calling phenology'
    CALL jules_print('veg-veg2a_jls',jules_message)
! DEPENDS ON: phenol
    CALL phenol (land_pts,surft_pts(n),surft_index(:,n),n,                    &
                 g_leaf_day(:,n),ht(:,n),dtime_phen,                          &
                 g_leaf_phen(:,n),lai(:,n))

    WRITE(jules_message,*) 'Phenology completed normally'
    CALL jules_print('veg-veg2a_jls',jules_message)

    DO l=1,land_pts
      lai_phen(l,n)=lai(l,n)
    END DO

!-----------------------------------------------------------------------
! Increment the leaf turnover rate for driving TRIFFID and reset the
! accumulation over atmospheric model timesteps to zero.
!-----------------------------------------------------------------------
    DO j=1,surft_pts(n)
      l=surft_index(j,n)
      g_leaf_phen_ac(l,n)=g_leaf_phen_ac(l,n)                                 &
                    +g_leaf_phen(l,n)*dtime_phen
    END DO

    DO l=1,land_pts
      g_leaf_ac(l,n)=0.0
    END DO

  END DO
END IF

!-----------------------------------------------------------------------
! Call TRIFFID vegetation model to update vegetation and terrestrial
! carbon storage.
!-----------------------------------------------------------------------
IF (l_triffid .AND.                                                           &
   (asteps_since_triffid == nstep_trif)) THEN

!-----------------------------------------------------------------------
! Calculate the TRIFFID inverse coupling timestep.
!-----------------------------------------------------------------------
  gam_trif=360.0/FLOAT(triffid_period)

!-----------------------------------------------------------------------
! Diagnose the mean fluxes over the coupling period.
!-----------------------------------------------------------------------
  DO l=1,land_pts
    resp_s_dr(l,1)=resp_s_ac(l,1)*gam_trif
    resp_s_dr(l,2)=resp_s_ac(l,2)*gam_trif
    resp_s_dr(l,3)=resp_s_ac(l,3)*gam_trif
    resp_s_dr(l,4)=resp_s_ac(l,4)*gam_trif
  END DO

  DO n=1,nnpft
    DO j=1,surft_pts(n)
      l=surft_index(j,n)
      g_leaf_dr(l,n)=g_leaf_phen_ac(l,n)*gam_trif
      npp_dr(l,n)=npp_ac(l,n)*gam_trif
      resp_w_dr(l,n)=resp_w_ac(l,n)*gam_trif
    END DO
  END DO

!-----------------------------------------------------------------------
! Diagnose the mean leaf turnover rates over the coupling period.
!-----------------------------------------------------------------------
  IF (l_phenol) THEN
    DO n=1,nnpft
      DO j=1,surft_pts(n)
        l=surft_index(j,n)
        g_leaf_dr(l,n)=g_leaf_phen_ac(l,n)*gam_trif
      END DO
    END DO
  ELSE
    DO n=1,nnpft
      DO j=1,surft_pts(n)
        l=surft_index(j,n)
        g_leaf_dr(l,n)=g_leaf_ac(l,n)*gam_trif
      END DO
    END DO
  END IF

!-----------------------------------------------------------------------
! Define the agricultural regions.
!-----------------------------------------------------------------------
  IF (agric) THEN
    DO l=1,land_pts
        frac_agric(l)=MIN(fraca(l),frac_vs(l)-((nnpft+1)*frac_min))
        frac_agric(l)=MAX(frac_agric(l),0.0)
    END DO
  END IF
  IF (l_trif_crop) THEN
    DO l=1,land_pts
        frac_past(l)=MIN(fracp(l),                                            &
                            frac_vs(l)-((nnpft+1)*frac_min)-frac_agric(l))
        frac_past(l)=MAX(frac_past(l),0.0)
    END DO
  END IF

!-----------------------------------------------------------------------
! Take copies of TRIFFID input variables for output as diagnostics.
!-----------------------------------------------------------------------
  DO n=1,nnpft
    DO l=1,land_pts
      g_leaf_dr_out(l,n)=g_leaf_dr(l,n)
      npp_dr_out(l,n)=npp_dr(l,n)
      resp_w_dr_out(l,n)=resp_w_dr(l,n)
      frac_old(l,n)=frac(l,n)
    END DO
  END DO
  DO n=1,5
    DO l=1,land_pts
      resp_s_dr_out(l,n)=resp_s_dr(l,n)
    END DO
  END DO
  DO l=1,land_pts
    dcs(l)=cs(l,1)+cs(l,2)+cs(l,3)+cs(l,4)

    dcs_rothc(l,1)=cs(l,1)
    dcs_rothc(l,2)=cs(l,2)
    dcs_rothc(l,3)=cs(l,3)
    dcs_rothc(l,4)=cs(l,4)

    resp_s_dr_out(l,5) = resp_s_dr_out(l,1) + resp_s_dr_out(l,2)              &
                       + resp_s_dr_out(l,3) + resp_s_dr_out(l,4)
  END DO

!-----------------------------------------------------------------------
! Select timestep and forward timestep weighting parameters for
! equilibrium or dynamic vegetation and call TRIFFID.
!-----------------------------------------------------------------------

! calculate RESP_FRAC.
! NOTE - here we want to use (1-X) - i.e. the fraction of
!        respiration NOT released to the atmos, so RESP_FRAC
!        here is 1-RESP_FRAC as calc'd in BL_CTL.

DO l=1,land_pts
  resp_frac(l) = 1.0 / (resp_frac_a + resp_frac_b*                            &
                   EXP(resp_frac_c * 100.0*clay_gb(l)))
END DO

  IF (l_trif_eq) THEN
    forw=1.0
    GAMMA=gamma_eq
    kiter=iter_eq
  ELSE
    forw=0.0
    GAMMA=gam_trif
    kiter=1
  END IF

  DO k=1,kiter

    WRITE(jules_message,*) 'Calling TRIFFID'
    CALL jules_print('veg-veg2a_jls',jules_message)

! DEPENDS ON: triffid
    CALL triffid (land_pts,trif_pts,trif_index,forw,GAMMA                     &
,                 frac_agric,frac_past                                        &
,                 g_leaf_dr,npp_dr,resp_s_dr                                  &
,                 resp_w_dr,cs,frac,ht,lai,resp_frac                          &
,                 c_veg,cv_trif,lit_c,lit_c_mn                                &
,                 wood_prod_total,wp_total_out)

    WRITE(jules_message,*) 'TRIFFID completed normally'
    CALL jules_print('veg-veg2a_jls',jules_message)

  END DO

!----------------------------------------------------------------------
! Check carbon conservation (2/3)
! Add the wood product decay flux to the net carbon flux
!----------------------------------------------------------------------
  DO t=1,trif_pts
    l=trif_index(t)
    cnsrv_carbon_flux(l)=cnsrv_carbon_flux(l)-wp_total_out(l)
  END DO

!-----------------------------------------------------------------------
! Update surft_index for new surface type fractions.
!-----------------------------------------------------------------------
! DEPENDS ON: tilepts
  CALL tilepts(land_pts,frac,surft_pts,surft_index)


!-----------------------------------------------------------------------
! Reset the accumulation fluxes to zero in equilibrium mode.
!-----------------------------------------------------------------------
  IF (l_trif_eq) THEN

    DO n=1,4
      DO l=1,land_pts
        resp_s_ac(l,n)=0.0
      END DO
    END DO

    DO n=1,nnpft
      DO l=1,land_pts
        npp_ac(l,n)=0.0
        resp_w_ac(l,n)=0.0
      END DO
    END DO

!-----------------------------------------------------------------------
! Reset the accumulation fluxes to the TRIFFID-diagnosed corrections
! in dynamic mode. Such corrections will be typically associated with
! the total depletion of a carbon reservoir or with the maintenance
! of the plant seed fraction. Any residual for the 4 component pools is
! assumed to be in proportion to the pool size itself.

! (in the case of zero total soil carbon, use 1e-10 to prevent zero
!  divide. This will have no impact because it is multiplied by the
!  constituent stores which are also zero...)
!-----------------------------------------------------------------------
  ELSE

    DO l=1,land_pts
      cs_tot(l)=MAX(1e-10,cs(l,1)+cs(l,2)+cs(l,3)+cs(l,4))
      dcs(l)=cs_tot(l)-dcs(l)
      resp_s_dr(l,1)=(1-resp_frac(l)) * (resp_s_dr_out(l,5))
      resp_s_ac(l,1)=(lit_c_mn(l)-GAMMA*dcs(l)-resp_s_dr(l,1)) *              &
                       cs(l,1) / (GAMMA*cs_tot(l))
      resp_s_ac(l,2)=(lit_c_mn(l)-GAMMA*dcs(l)-resp_s_dr(l,1)) *              &
                       cs(l,2) / (GAMMA*cs_tot(l))
      resp_s_ac(l,3)=(lit_c_mn(l)-GAMMA*dcs(l)-resp_s_dr(l,1)) *              &
                       cs(l,3) / (GAMMA*cs_tot(l))
      resp_s_ac(l,4)=(lit_c_mn(l)-GAMMA*dcs(l)-resp_s_dr(l,1)) *              &
                       cs(l,4) / (GAMMA*cs_tot(l))
    END DO

    DO n=1,nnpft
      DO j=1,surft_pts(n)
        l=surft_index(j,n)
        ratio=frac_old(l,n)/frac(l,n)
        npp_ac(l,n)=ratio*(npp_dr(l,n)-npp_dr_out(l,n))/gam_trif
        resp_w_ac(l,n)=ratio*(resp_w_dr(l,n)-resp_w_dr_out(l,n))              &
                             /gam_trif
      END DO
    END DO

  END IF

!-----------------------------------------------------------------------
! Reset the accumulated leaf turnover rates to zero.
!-----------------------------------------------------------------------
  IF (l_phenol) THEN
    DO n=1,nnpft
      DO l=1,land_pts
        g_leaf_phen_ac(l,n)=0.0
      END DO
    END DO
  ELSE
    DO n=1,nnpft
      DO l=1,land_pts
        g_leaf_ac(l,n)=0.0
      END DO
    END DO
  END IF

  asteps_since_triffid=0

END IF

!-----------------------------------------------------------------------
! Calculate gridbox mean vegetation parameters from fractions of
! surface functional types
!-----------------------------------------------------------------------
! DEPENDS ON: sparm
CALL sparm (land_pts, nsurft, can_model, l_aggregate,                         &
            surft_pts, surft_index,                                           &
            frac, ht, lai, satcon, z0m_soil, catch_s, catch_t,                &
            infil_t, z0_t, z0h_t)

DO t = 1, trif_pts
  l = trif_index(t)

  cv(l) = 0.0

  DO n = 1, nnpft
    cv(l) = cv(l) + frac(l,n) * c_veg(l,n)
  END DO

  IF ( l_crop ) THEN
    DO n = 1, ncpft
      c_veg(l, n + nnpft) = rootc_cpft(l,n) + harvc_cpft(l,n) +               &
                            reservec_cpft(l,n) + stemc_diag_cpft(l,n) +       &
                            leafc_diag_cpft(l,n)
      cv(l) = cv(l) + frac(l, n + nnpft) * c_veg(l, n + nnpft)
    END DO
  END IF
END DO

!----------------------------------------------------------------------
! Check carbon conservation (3/3)
! Calculate the error in carbon conservation as:
! error=(final store)-(initial store)-(net flux)
!      ="extra carbon in store not accounted for by fluxes"
!----------------------------------------------------------------------
  DO t=1,trif_pts
    l=trif_index(t)
    cnsrv_carbon_veg2_gb(l)=(SUM(frac(l,1:npft)*c_veg(l,1:npft))+             &
                           SUM(cs(l,:))+                                      &
                           wood_prod_total(l))-                               &
                         cnsrv_carbon_veg2_gb(l)-                             &
                         cnsrv_carbon_flux(l)-                                &
                         (SUM(npp_ac(l,1:npft)*frac(l,1:npft))-               &
                          SUM(resp_s_ac(l,1:4)))
  END DO

IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
RETURN

END SUBROUTINE veg2

