#if defined(L19_2A)
! *****************************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,ntiles,can_model              &
,              row_length,rows                                    &
,              a_step,asteps_since_triffid                        &
,              phenol_period,triffid_period                       &
,              l_phenol,l_triffid,l_trif_eq                       &
,              atimestep,fraca,satcon                             &
,              g_leaf_ac,g_leaf_phen_ac,npp_ac                    &
,              resp_s_ac,resp_w_ac                                &
,              cs,frac,lai,clay_frac,ht                           &
,              catch_s,catch_t,infil_t,z0_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 nstypes, ONLY :                                               &
!      imported scalars with intent(in)
  npft,ntype,soil

USE descent
USE seed
USE veg_param, ONLY : agric
USE c_mdi, ONLY : rmdi

USE switches, ONLY : l_aggregate

USE parkind1, ONLY: jprb, jpim
USE yomhook, ONLY: lhook, dr_hook
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
,ntiles                                                           &
                       ! IN Number of land-surface tiles.
,can_model                                                        &
                              ! IN Swith for thermal vegetation
,row_length                                                       &
                  ! number of points on a row
,rows                                                             &
                  ! number of rows in a theta field
,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.
,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,ntiles)                                         &
                              ! OUT Snow capacity for tiles
!                                   !     (kg/m2).
,catch_t(land_pts,ntiles)                                         &
                              ! OUT Canopy capacity for tiles
!                                   !     (kg/m2).
,infil_t(land_pts,ntiles)                                         &
                              ! OUT Maximum surface infiltration
!                                   !     rate for tiles (kg/m2/s).
,z0_t(land_pts,ntiles)                                            &
                              ! OUT 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_frac(row_length, rows)                                      &
                                  ! IN  Clay fraction of soil
,clay_land(land_pts)                                              &
                                  ! IN  clay on land points
,resp_frac(land_pts)            ! respired fraction of RESP_S
!                                   !     driving TRIFFID
!                                   !     (kg C/m2/360days).

INTEGER                                                           &
 nstep_phen                                                       &
                              ! WORK Number of atmospheric
!                                   !      timesteps between calls to
!                                   !      PHENOL.
,nstep_trif                                                       &
                              ! WORK Number of atmospheric
!                                   !      timesteps between calls to
!                                   !      TRIFFID.
,tile_pts(ntype)                                                  &
                              ! WORK Number of land points which
!                                   !      include the nth surface type.
,tile_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

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).
,frac_agric(land_pts)                                             &
                              ! WORK Fraction of agriculture 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.

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

IF (lhook) CALL dr_hook('VEG2',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,ntiles
  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 (can_model == 4) THEN
  DO n=1,ntiles
    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
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,npft
    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) >= (npft*frac_min)) THEN
    trif_pts = trif_pts + 1
    trif_index(trif_pts) = l
  END IF
END DO

!-----------------------------------------------------------------------
! Create the TILE_INDEX array of land points with each surface type
!-----------------------------------------------------------------------
! DEPENDS ON: tilepts
CALL tilepts(land_pts,frac,tile_pts,tile_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,npft

!-----------------------------------------------------------------------
! 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,tile_pts(n)
      l=tile_index(j,n)
      g_leaf_day(l,n)=g_leaf_ac(l,n)/dtime_phen
    END DO

    WRITE(6,*) 'Calling phenology'

! DEPENDS ON: phenol
    CALL phenol (land_pts,tile_pts(n),tile_index(:,n),n,          &
                 g_leaf_day(:,n),ht(:,n),dtime_phen,              &
                 g_leaf_phen(:,n),lai(:,n))

    WRITE(6,*) 'Phenology completed normally'

    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,tile_pts(n)
      l=tile_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,npft
    DO j=1,tile_pts(n)
      l=tile_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,npft
      DO j=1,tile_pts(n)
        l=tile_index(j,n)
        g_leaf_dr(l,n)=g_leaf_phen_ac(l,n)*gam_trif
      END DO
    END DO
  ELSE
    DO n=1,npft
      DO j=1,tile_pts(n)
        l=tile_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)=fraca(l)
  END DO
  END IF

!-----------------------------------------------------------------------
! Take copies of TRIFFID input variables for output as diagnostics.
!-----------------------------------------------------------------------
  DO n=1,npft
    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)
    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
  j=(land_index(l)-1)/row_length + 1
  i=land_index(l) - (j-1)*row_length
  clay_land(l) = clay_frac(i,j)
! Set default value for missing data points
  IF (clay_land(l)  == rmdi ) clay_land(l)=0.23
  resp_frac(l) = 1.0 / (4.0895 + 2.672*                           &
                   EXP(-0.0786 * 100.0*clay_land(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(6,*) 'Calling TRIFFID'

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

    WRITE(6,*) 'TRIFFID completed normally'

  END DO

!-----------------------------------------------------------------------
! Update TILE_INDEX for new surface type fractions.
!-----------------------------------------------------------------------
! DEPENDS ON: tilepts
  CALL tilepts(land_pts,frac,tile_pts,tile_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,npft
      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,npft
      DO j=1,tile_pts(n)
        l=tile_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,npft
      DO l=1,land_pts
        g_leaf_phen_ac(l,n)=0.0
      END DO
    END DO
  ELSE
    DO n=1,npft
      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,ntiles,can_model,l_aggregate                 &
,           tile_pts,tile_index                                   &
,           frac,ht,lai,satcon,catch_s,catch_t,infil_t,z0_t)

IF (lhook) CALL dr_hook('VEG2',zhook_out,zhook_handle)
RETURN
END SUBROUTINE veg2
#endif
