Jules red veg3 couple soil carbon

CONTRIBUTORS.md

@@ -8,3 +8,4 @@
 | maggiehendry    | Maggie Hendry    | Met Office  | 2026-08-21 |
 | DanCopsey       | Dan Copsey       | Met Office  | 2026-08-24 |
 | Pierre-siddall  | Pierre Siddall   | Met Office  | 2026-08-26 |
+| mo-arthurargles | Arthur Argles    | Met Office  | 2026-09-04 |

src/control/shared/surf_couple_extra_mod.F90

@@ -8,6 +8,9 @@ MODULE surf_couple_extra_mod
 !
 ! [Met Office Ref SC0237]
 ! *****************************COPYRIGHT****************************************
+!
+! Some of the content of this file has been produced with the assistance of
+! Met Office Github Copilot Enterprise.
 USE um_types, ONLY: real_jlslsm
 
 IMPLICIT NONE
@@ -236,9 +239,9 @@ SUBROUTINE surf_couple_extra(                                                  &
 
 USE timestep_mod,             ONLY: timestep
 
-USE veg3_parm_mod,            ONLY: veg3_ctrl,litter_parms,red_parms
+USE veg3_parm_mod,            ONLY: veg3_ctrl,litter_parms,red_parms,soil_parms
 
-USE veg3_field_mod,           ONLY: veg_state,red_state
+USE veg3_field_mod,           ONLY: veg_state,red_state,soil_state
 
 USE water_constants_mod,      ONLY: rho_water
 
@@ -1073,12 +1076,12 @@ SUBROUTINE surf_couple_extra(                                                  &
 
         CALL next_gen_biogeochem(                                              &
           !IN control vars
-            asteps_since_triffid,land_pts,nnpft,nmasst,veg3_ctrl,              &
+            asteps_since_triffid,a_step,land_pts,nnpft,nmasst,veg3_ctrl,       &
             ainfo,                                                             &
           !IN parms
-            litter_parms,red_parms,                                            &
+            litter_parms,red_parms,soil_parms,                                 &
           !INOUT data structures
-            veg_state,red_state                                                &
+            veg_state,red_state,soil_state                                     &
           !OUT diagnostics
           )
 
@@ -1087,10 +1090,9 @@ SUBROUTINE surf_couple_extra(                                                  &
         ! - Uses non-veg tilepts
         ! - Accesses fields via USE statements
 
-        ! Note use lai_bal for now as we don't yet have phenology
         CALL sparm (land_pts,nsurft,surft_pts,ainfo%surft_index,               &
                     veg_state%frac,veg_state%canht,                            &
-                    veg_state%lai_bal,psparms%z0m_soil_gb,                     &
+                    veg_state%lai,psparms%z0m_soil_gb,                         &
                     psparms%catch_snow_surft,psparms%catch_surft,              &
                     psparms%z0_surft,psparms%z0h_bare_surft,urban_param%ztm_gb)
 

src/initialisation/standalone/init_mod.F90

@@ -150,9 +150,8 @@ SUBROUTINE init(nml_dir, crop_vars_data, crop_vars,                            &
 ! USE veg3_parm_mod,              ONLY: in_dev
                                         !
                                         !
-! USE veg3_field_mod,             ONLY: in_dev
-                                        !
-                                        !
+USE veg3_field_mod,               ONLY: veg3_field_assoc
+
 USE jules_chemvars_mod,           ONLY: chemvars_data_type,                    &
                                         chemvars_type,                         &
                                         chemvars_assoc
@@ -352,7 +351,7 @@ SUBROUTINE init(nml_dir, crop_vars_data, crop_vars,                            &
 CALL imgn_vars_assoc(imgn_vars,imgn_vars_data)
 CALL rivers_assoc(rivers,rivers_data)
 !CALL veg3_parm_assoc(in_dev)
-!CALL veg3_field_assoc(in_dev)
+CALL veg3_field_assoc(progs,ainfo,psparms,trifctl_data,trif_vars_data)
 CALL chemvars_assoc(chemvars,chemvars_data)
 CALL water_resources_assoc(water_resources,water_resources_data)
 CALL wtrac_jls_assoc(wtrac_jls,wtrac_jls_data)

src/initialisation/standalone/init_vars_tmp_mod.F90

@@ -20,7 +20,7 @@ SUBROUTINE init_vars_tmp(crop_vars,psparms,toppdm,ainfo,trif_vars, aerotype,   &
 USE Ancil_info, ONLY: lice_pts
 USE jules_deposition_mod, ONLY: dzl_const, l_deposition
 USE jules_vegetation_mod, ONLY: l_triffid, l_use_pft_psi, fsmc_shape,          &
-                                l_phenol
+                                l_phenol, l_red
 
 USE trif, ONLY: lai_min
 
@@ -207,7 +207,7 @@ SUBROUTINE init_vars_tmp(crop_vars,psparms,toppdm,ainfo,trif_vars, aerotype,   &
 ! This is necessary to ensure time average diagnostics are produced correctly.
 !-----------------------------------------------------------------------------
 
-IF (l_triffid) THEN
+IF (l_triffid .AND. .NOT. l_red) THEN
   trif_vars%rootC_pft       = 0.0
   trif_vars%woodC_pft       = 0.0
   trif_vars%leafC_pft       = 0.0

src/science/params/veg3_param_mod.F90

@@ -8,6 +8,9 @@
 ! Code Owner: Please refer to ModuleLeaders.txt
 ! This file belongs in Veg3 Ecosystem Demography
 ! *****************************COPYRIGHT****************************************
+!
+! Some of the content of this file has been produced with the assistance of
+! Met Office Github Copilot Enterprise.
 
 MODULE veg3_parm_mod
 
@@ -18,8 +21,9 @@ MODULE veg3_parm_mod
 !Set up object for veg3 control
 
 TYPE :: veg3_ctrl_type
-  INTEGER :: land_pts,nsurft,npft,nnpft,soil,triffid_period,nstep_trif,nmasst
-  REAL    :: timestep,frac_min
+  INTEGER :: land_pts,nsurft,npft,nnpft,soil,triffid_period,nstep_trif,nmasst, &
+             phenol_period,nstep_phen
+  REAL    :: timestep,frac_min,dt_red,dt_phen_360d
 END TYPE veg3_ctrl_type
 
 !Set up objects containing everything we need for litter calculation
@@ -56,20 +60,47 @@ MODULE veg3_parm_mod
 TYPE(litter_parm_type) :: litter_parms
 TYPE(red_parm_type)    :: red_parms
 
+!Set up object containing everything we need for soil carbon coupling
+
+TYPE :: soil_parm_type
+  LOGICAL :: l_layeredc
+              ! Layered (.TRUE.) or single-layer (.FALSE.) soil carbon.
+  INTEGER :: soil_bgc_model
+              ! Soil biogeochemistry model in use.
+  INTEGER :: dim_cslayer
+              ! Number of soil carbon layers.
+  INTEGER :: dim_cs1
+              ! Number of soil carbon pools.
+  REAL    :: tau_lit
+              ! Litter decomposition rate exponent for the vertical litter
+              ! profile (m-1).
+  REAL    :: litc_norm
+              ! Normalisation for the vertical litter profile.
+  REAL    :: resp_frac_a
+  REAL    :: resp_frac_b
+  REAL    :: resp_frac_c
+              ! Coefficients relating clay content to the fraction of soil
+              ! respiration that forms new soil C.
+  REAL, ALLOCATABLE :: dzsoil(:)
+              ! Soil layer thicknesses (m).
+END TYPE soil_parm_type
+
+TYPE(soil_parm_type)   :: soil_parms
+
 !Private by default
 PRIVATE
 
 !Expose routines
 PUBLIC :: veg3_parm_init, veg3_parm_allocate, check_jules_red_parms
 
 !Expose data
-PUBLIC :: veg3_ctrl, litter_parms, red_parms, l_red
+PUBLIC :: veg3_ctrl, litter_parms, red_parms, soil_parms, l_red
 
 !Expose data structures
-PUBLIC :: veg3_ctrl_type, litter_parm_type, red_parm_type
+PUBLIC :: veg3_ctrl_type, litter_parm_type, red_parm_type, soil_parm_type
 
 !Allow external code to read but not write
-PROTECTED :: litter_parms, veg3_ctrl, red_parms
+PROTECTED :: litter_parms, veg3_ctrl, red_parms, soil_parms
 
 CHARACTER(LEN=*), PARAMETER, PRIVATE :: ModuleName='VEG3_PARM_MOD'
 
@@ -79,6 +110,7 @@ MODULE veg3_parm_mod
 SUBROUTINE veg3_parm_allocate(land_pts,nsurft,nnpft,npft)
 
 USE missing_data_mod, ONLY: rmdi, imdi
+USE ancil_info,       ONLY: dim_cslayer
 
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: land_pts, nsurft, nnpft, npft
@@ -109,6 +141,9 @@ SUBROUTINE veg3_parm_allocate(land_pts,nsurft,nnpft,npft)
 ALLOCATE(red_parms%comp_coef          (nnpft,nnpft))
 ALLOCATE(red_parms%mclass_geom_mult   (nnpft))
 
+! Allocate soil_parm_type
+ALLOCATE(soil_parms%dzsoil(dim_cslayer))
+
 litter_parms%g_wood          = rmdi
 litter_parms%g_leaf          = rmdi
 litter_parms%g_root          = rmdi
@@ -130,6 +165,8 @@ SUBROUTINE veg3_parm_allocate(land_pts,nsurft,nnpft,npft)
 red_parms%comp_coef          = rmdi
 red_parms%mclass_geom_mult   = rmdi
 
+soil_parms%dzsoil             = rmdi
+
 RETURN
 END SUBROUTINE veg3_parm_allocate
 
@@ -141,7 +178,8 @@ SUBROUTINE veg3_set_parms(land_pts,nsurft,nnpft,npft,nmasst)
 USE pftparm,                  ONLY: g_leaf_0
 USE trif,                     ONLY: g_root, g_wood
 ! Above only allocated if triffid on - needs to be addressed
-USE jules_vegetation_mod,     ONLY: l_triffid, triffid_period,frac_min
+USE jules_vegetation_mod,     ONLY: l_triffid, triffid_period,frac_min,        &
+                                    phenol_period
 
 USE red_io,                   ONLY: alpha_recrt, crwn_area0, dom_order,        &
                                     height0, lai_bal0, mass0, massi, mclass,   &
@@ -152,12 +190,26 @@ SUBROUTINE veg3_set_parms(land_pts,nsurft,nnpft,npft,nmasst)
 
 USE jules_surface_types_mod,  ONLY: soil
 
+!Soil carbon coupling parameters
+USE jules_soil_biogeochem_mod, ONLY: l_layeredC, soil_bgc_model, tau_lit
+USE ancil_info,                ONLY: dim_cslayer, dim_cs1
+#if !defined(UM_JULES)
+USE jules_soil_mod,            ONLY: dzsoil
+USE veg_param,                 ONLY: litc_norm
+#endif
+
 IMPLICIT NONE
 
 INTEGER, INTENT(IN) :: land_pts, nsurft, nnpft, npft, nmasst
 
 INTEGER :: n,k
 
+! Coefficients relating clay content to the fraction of soil respiration
+! that forms new soil C (i.e. is NOT released to the atmosphere).
+REAL, PARAMETER :: resp_frac_a_local = 4.0895
+REAL, PARAMETER :: resp_frac_b_local = 2.672
+REAL, PARAMETER :: resp_frac_c_local = -0.0786
+
 !End of header
 
 IF (l_red .AND. l_triffid) THEN
@@ -174,9 +226,15 @@ SUBROUTINE veg3_set_parms(land_pts,nsurft,nnpft,npft,nmasst)
 
   veg3_ctrl%timestep = REAL(timestep)
   veg3_ctrl%triffid_period = triffid_period
+  veg3_ctrl%dt_red = rsec_per_day * REAL(veg3_ctrl%triffid_period)
   veg3_ctrl%nstep_trif = INT(rsec_per_day * veg3_ctrl%triffid_period           &
     / veg3_ctrl%timestep)
 
+  veg3_ctrl%phenol_period = phenol_period
+  veg3_ctrl%dt_phen_360d = REAL(veg3_ctrl%phenol_period) / 360.0
+  veg3_ctrl%nstep_phen = INT(rsec_per_day * veg3_ctrl%phenol_period            &
+    / veg3_ctrl%timestep)
+
   veg3_ctrl%land_pts = land_pts
   veg3_ctrl%nsurft   = nsurft
   veg3_ctrl%nnpft    = nnpft
@@ -199,18 +257,25 @@ SUBROUTINE veg3_set_parms(land_pts,nsurft,nnpft,npft,nmasst)
   red_parms%mass0(:)            = mass0(1:nnpft)
   red_parms%massi(:)            = massi(1:nnpft)
   red_parms%mclass(:)           = mclass(1:nnpft)
-  red_parms%mort_base(:)        = mort_base(1:nnpft)
+  red_parms%mort_base(:)        = mort_base(1:nnpft) / rsec_per_day / 360.0
   red_parms%phi_a(:)            = phi_a(1:nnpft)
   red_parms%phi_g(:)            = phi_g(1:nnpft)
   red_parms%phi_h(:)            = phi_h(1:nnpft)
   red_parms%phi_l(:)            = phi_l(1:nnpft)
-  red_parms%mclass_geom_mult(:) = 0.0
+  red_parms%mclass_geom_mult(:) = 1.0 ! Default assumes 1 mass class
   red_parms%frac_min(:)         = frac_min
 
   red_parms%comp_coef(:,:)      = 0.0
 
   DO n = 1,nnpft
     ! Cycle through the PFTs
+    ! Update mclass_geom_mult for each PFT
+    IF (red_parms%mclass(n) > 1) THEN
+      red_parms%mclass_geom_mult(n) =                                          &
+        (red_parms%massi(n) / red_parms%mass0(n))**                            &
+        (1.0 / REAL(red_parms%mclass(n)-1))
+    END IF
+
     DO k=1,nnpft
       ! If the n'th PFT is less dominant than k'th PFT, then k shades n.
       IF (dom_order(n)  <=  dom_order(k)) THEN
@@ -220,6 +285,23 @@ SUBROUTINE veg3_set_parms(land_pts,nsurft,nnpft,npft,nmasst)
     END DO
   END DO
 
+  ! Soil carbon coupling parameters
+  soil_parms%l_layeredc    = l_layeredC
+  soil_parms%soil_bgc_model = soil_bgc_model
+  soil_parms%dim_cslayer   = dim_cslayer
+  soil_parms%dim_cs1       = dim_cs1
+  soil_parms%tau_lit       = tau_lit
+  soil_parms%resp_frac_a   = resp_frac_a_local
+  soil_parms%resp_frac_b   = resp_frac_b_local
+  soil_parms%resp_frac_c   = resp_frac_c_local
+#if !defined(UM_JULES)
+  soil_parms%litc_norm     = litc_norm
+  soil_parms%dzsoil(:)     = dzsoil(1:dim_cslayer)
+#else
+  soil_parms%litc_norm     = 1.0
+  soil_parms%dzsoil(:)     = 0.0
+#endif
+
 END IF
 
 RETURN
@@ -246,6 +328,7 @@ SUBROUTINE check_jules_red_parms()
 
 USE ereport_mod,     ONLY: ereport
 USE jules_print_mgr, ONLY: jules_print, jules_message
+USE jules_soil_biogeochem_mod, ONLY: soil_model_4pool, soil_bgc_model
 
 IMPLICIT NONE
 
@@ -259,6 +342,15 @@ SUBROUTINE check_jules_red_parms()
 !-----------------------------------------------------------------------------
 error_sum = 0
 IF ( l_red ) THEN
+
+  ! veg3_soil_couple only supports the 4-pool soil carbon model (layered
+  ! or single-layer).
+  IF ( soil_bgc_model /= soil_model_4pool ) THEN
+    error_sum = error_sum + 1
+    CALL jules_print(RoutineName, "l_red requires soil_bgc_model=" //          &
+      "soil_model_4pool")
+  END IF
+
   IF ( ANY( red_parms%alpha_recrt(:) < 0 ) ) THEN
     error_sum = error_sum + 1
     CALL jules_print(RoutineName, "No value for alpha_recrt")

src/science/surface/cancap_jls.F90

@@ -127,6 +127,11 @@ SUBROUTINE cancap (land_pts,veg_pts,veg_index,can_model,ft,                    &
     IF ( ft > nnpft ) THEN
       leaf = leafc_from_prognostics(ft - nnpft, dvi_cpft(l,ft - nnpft), lai(l))
       wood = stemc_from_prognostics(ft - nnpft, ht(l) )
+
+    ELSE IF (l_red) THEN
+      leaf = veg_state%leafC(l,ft)
+      wood = veg_state%woodC(l,ft)
+
     ELSE
       lai_bal(l) = ( a_ws(ft) * eta_sl(ft) * ht(l) /                           &
                      a_wl(ft) )**(1.0 / (b_wl(ft) - 1))
@@ -140,11 +145,7 @@ SUBROUTINE cancap (land_pts,veg_pts,veg_index,can_model,ft,                    &
       wood = a_wl(ft) * (lai_bal(l)**b_wl(ft))
     END IF
 
-    IF (l_red) THEN
-      canhc(l) = hleaf * veg_state%leafC(l,ft) + hwood * veg_state%woodC(l,ft)
-    ELSE
-      canhc(l) = hleaf * leaf + hwood * wood
-    END IF
+    canhc(l) = hleaf * leaf + hwood * wood
 
     ! If using tree heights from alternative sources (such as SIMARD dataset),
     ! then the heat capacicy can become unrealisitcaly large for the tallest trees.

src/science/surface/sf_stom_jls_mod.F90

@@ -1473,10 +1473,14 @@ SUBROUTINE sf_stom  (land_pts,land_index                                       &
   ! calculate the total nitrogen content of the leaf, root and stem
   ! Assume that root biomass is equal to balanced growth leaf biomass
   !---------------------------------------------------------------------------
-  lai_bal(l) = (a_ws(ft) * eta_sl(ft) * ht(l) / a_wl(ft))                      &
+  IF (l_red) THEN
+    lai_bal(l) = veg_state%lai_bal(l,ft)
+
+  ELSE
+    lai_bal(l) = (a_ws(ft) * eta_sl(ft) * ht(l) / a_wl(ft))                    &
                **(1.0 / (b_wl(ft) - 1.0))
 
-  IF (l_red) lai_bal(l) = veg_state%lai_bal(l,ft)
+  END IF
 
   !---------------------------------------------------------------------------
   ! Calculate the total nitrogen content of the leaf, root and stem
@@ -1521,9 +1525,11 @@ SUBROUTINE sf_stom  (land_pts,land_index                                       &
       n_root(l) = nr(ft) * root(l) * cmass
 
       !Initial calculation of N content in respiring stem wood
-      n_stem(l) = eta_sl(ft) * ht(l) * lai_bal(l) * nsw(ft)
-      !USE veg3 allometry
-      IF (l_red) n_stem(l) = veg_state%woodC(l,ft) / a_ws(ft) * nsw(ft)
+      IF (l_red) THEN
+        n_stem(l) = veg_state%woodC(l,ft) / a_ws(ft) * nsw(ft)
+      ELSE
+        n_stem(l) = eta_sl(ft) * ht(l) * lai_bal(l) * nsw(ft)
+      END IF
 
       !Reduce n_stem for consistency with non-trait n_stem for now
       !This must be done to achieve realistic respiration rates.