From c193294a65579cbce3113030ceebefe3440daa76 Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Fri, 21 Aug 2026 16:46:09 +0100 Subject: [PATCH 01/14] Couple veg3 to JULES phenology scheme and add RED diagnostics (#90) - Add veg3_phenol_couple to run phenology on its own timestep and diagnose the mean phenology-driven leaf turnover rate driving litter/RED dynamics, replacing the lai_bal-only estimate - Thread a_step and new phenology period controls (nstep_phen, dt_phen_360d) through veg3_ctrl/next_gen_biogeochem - Use veg3 lai/canht/leafC/woodC/lai_bal in sparm, cancap and sf_stom instead of TRIFFID lai_bal estimates when l_red is set - Add pointer associations for new veg3/RED diagnostic outputs (npp_dr_out, g_leaf_phen, g_leaf_phen_acc, lai_phen, litter fluxes) into trifctl/trif_vars - Refactor veg3_red_dynamic_mod, splitting growth_onto_mass_class and update_pft_size_structure out of veg3_red_dynamic - Fix mclass_geom_mult default and mort_base unit conversion in veg3_param_mod Refs #90 --- .../shared/next_gen_biogeochem_mod.F90 | 236 +++++++- src/control/shared/surf_couple_extra_mod.F90 | 8 +- src/control/shared/veg3_field_mod.F90 | 360 ++++++++++-- src/initialisation/standalone/init_mod.F90 | 7 +- .../standalone/init_vars_tmp_mod.F90 | 4 +- src/science/params/veg3_param_mod.F90 | 28 +- src/science/surface/cancap_jls.F90 | 11 +- src/science/surface/sf_stom_jls_mod.F90 | 16 +- src/science/vegetation/phenol_jls.F90 | 29 +- src/science/vegetation/veg3_litter_mod.F90 | 28 +- .../vegetation/veg3_red_dynamic_mod.F90 | 522 +++++++++++++----- 11 files changed, 987 insertions(+), 262 deletions(-) diff --git a/src/control/shared/next_gen_biogeochem_mod.F90 b/src/control/shared/next_gen_biogeochem_mod.F90 index b5619d3a..584d1153 100644 --- a/src/control/shared/next_gen_biogeochem_mod.F90 +++ b/src/control/shared/next_gen_biogeochem_mod.F90 @@ -23,7 +23,7 @@ MODULE next_gen_biogeochem_mod SUBROUTINE 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, & @@ -62,6 +62,9 @@ SUBROUTINE next_gen_biogeochem( & ! INTEGERS with INTENT in !---------------------------------------------------------------------------- INTEGER, INTENT(IN) :: land_pts, nnpft, nmasst +INTEGER, INTENT(IN) :: a_step + ! Current atmospheric timestep number, used to + ! determine when the phenology timestep falls due. !---------------------------------------------------------------------------- ! Variables with INTENT inout @@ -113,7 +116,7 @@ SUBROUTINE next_gen_biogeochem( & ! Call the Vegetation Biogeochemistry model IF (veg_index_pts > 0) CALL veg3_run_ctrl( & !IN Control vars - asteps_since_triffid,land_pts,nnpft,nmasst,veg_index_pts, & + asteps_since_triffid,a_step,land_pts,nnpft,nmasst,veg_index_pts, & veg_index,veg3_ctrl,ainfo, & !IN parms litter_parms,red_parms, & @@ -130,8 +133,8 @@ END SUBROUTINE next_gen_biogeochem !------------------------------------------------------------------------------ SUBROUTINE veg3_run_ctrl( & !IN Control vars - asteps_since_triffid,land_pts,nnpft,nmasst,veg_index_pts, & - veg_index,veg3_ctrl,ainfo, & + asteps_since_triffid,a_step,land_pts,nnpft,nmasst, & + veg_index_pts,veg_index,veg3_ctrl,ainfo, & !IN parms litter_parms,red_parms, & !IN state @@ -140,7 +143,6 @@ SUBROUTINE veg3_run_ctrl( & ) !Only get the data structures - the data comes through the calling tree - USE veg3_parm_mod, ONLY: veg3_ctrl_type,litter_parm_type,red_parm_type USE veg3_field_mod, ONLY: veg_state_type,red_state_type,red_veg3_couple USE veg3_litter_mod, ONLY: veg3_litter @@ -162,6 +164,10 @@ SUBROUTINE veg3_run_ctrl( & ! IN Number of atmospheric timesteps since last call ! to TRIFFID. +INTEGER, INTENT(IN) :: a_step + ! IN Current atmospheric timestep number, used to + ! determine when the phenology timestep falls due. + !----------------------------------------------------------------------------- ! Objects with INTENT IN !----------------------------------------------------------------------------- @@ -194,26 +200,37 @@ SUBROUTINE veg3_run_ctrl( & ! growth (kg C/m2/s). mort_add(land_pts,nnpft,nmasst), & ! mortality above baseline (/m2) -demographic_lit(land_pts,nnpft) - ! demographic litter aggregated to PFT +g_leaf_phen_dr(land_pts,nnpft) + ! Mean phenology-driven leaf turnover rate for driving vegetation + ! dynamics (/s) !End of headers !Initialise Arrays -npp_dr(:,:) = 0.0 -mort_add(:,:,:) = 0.0 +npp_dr(:,:) = 0.0 +mort_add(:,:,:) = 0.0 +!----------------------------------------------------------------------------- +! Work out the phenology at its own timestep, appending to the accumulated +! leaf turnover rates, and (on the vegetation dynamics timestep) diagnose +! the mean phenology-driven leaf turnover rate for driving litterfall and +! vegetation dynamics. This mirrors the phenology/TRIFFID timestep split in +! veg-veg2a_jls_mod. +!----------------------------------------------------------------------------- +CALL veg3_phenol_couple(veg_index_pts,veg_index,veg3_ctrl,land_pts,nnpft, & + a_step,asteps_since_triffid,veg_state,g_leaf_phen_dr) ! Now call vegetation model IF (asteps_since_triffid == veg3_ctrl%nstep_trif) THEN - !Call veg3_phenol() - + !Call Litter CALL veg3_Litter( & !IN Control vars veg_index_pts,veg_index,veg3_ctrl,land_pts,nnpft, & !IN parms litter_parms, & + !IN fields + g_leaf_phen_dr, & !IN state veg_state, & ! OUT Fields @@ -229,13 +246,16 @@ SUBROUTINE veg3_run_ctrl( & ! Use the accumulated npp from sf_expl. Copy to new variable for driving ! veg dynamics - npp_dr = veg_state%npp_acc / (rsec_per_day * veg3_ctrl%triffid_period) - ! Reset accumulation to zero - note this can be used to pass a negative flux back to JULES. + npp_dr = veg_state%npp_acc / veg3_ctrl%dt_red + + ! Record driving npp_dr and convert to s-1 -> (360 days)-1 + veg_state%npp_dr_out(:,:) = npp_dr * rsec_per_day * 360 + + ! Reset accumulation to zero - note this can be used to pass a negative flux + ! back to JULES. veg_state%npp_acc(:,:)=0.0 growth = npp_dr - local_litter - !Call Litter - !Call Allocation/Nitrogen/NSC !Now on Mass classes @@ -251,16 +271,14 @@ SUBROUTINE veg3_run_ctrl( & !Call Veg Dynamics - in this case RED CALL veg3_red_dynamic( & !IN control vars - rsec_per_day,veg_index_pts,veg_index,veg3_ctrl,land_pts, & + veg3_ctrl%dt_red,veg_index_pts,veg_index,veg3_ctrl,land_pts, & nnpft,nmasst, & !IN RED_parms red_parms, & !IN fields growth,mort_add, & !IN state - veg_state,red_state, & - ! OUT Fields - demographic_lit & + veg_state,red_state & !OUT Diagnostics ) @@ -282,4 +300,184 @@ SUBROUTINE veg3_run_ctrl( & END SUBROUTINE veg3_run_ctrl +!------------------------------------------------------------------------------ +SUBROUTINE veg3_phenol_couple( & + !IN Control vars + veg_index_pts,veg_index,veg3_ctrl,land_pts,nnpft, & + a_step,asteps_since_triffid, & + !IN state + veg_state, & + !OUT fields + g_leaf_phen_dr & + ) + +! Diagnoses leaf phenology and the mean phenology-driven leaf turnover rate +! that drives litterfall/vegetation dynamics, in veg1/veg2. + +!Only get the data structures - the data comes through the calling tree +USE jules_vegetation_mod, ONLY: l_phenol +USE veg3_parm_mod, ONLY: veg3_ctrl_type +USE veg3_field_mod, ONLY: veg_state_type +USE phenol_mod, ONLY: phenol + +!Access some parameters direct from module +USE conversions_mod, ONLY: rsec_per_day + +IMPLICIT NONE + +!---------------------------------------------------------------------------- +! Integers with INTENT IN +!---------------------------------------------------------------------------- +INTEGER, INTENT(IN) :: land_pts,nnpft,veg_index_pts,veg_index(land_pts) + +INTEGER, INTENT(IN) :: a_step + ! Current atmospheric timestep number, used to + ! determine when the phenology timestep falls due. + +INTEGER, INTENT(IN) :: asteps_since_triffid + ! Number of atmospheric timesteps since last call to + ! vegetation dynamics. + +!----------------------------------------------------------------------------- +! Objects with INTENT IN +!----------------------------------------------------------------------------- +TYPE(veg3_ctrl_type),INTENT(IN) :: veg3_ctrl + +!----------------------------------------------------------------------------- +! Objects with INTENT INOUT +!----------------------------------------------------------------------------- +TYPE(veg_state_type),INTENT(IN OUT) :: veg_state + +!----------------------------------------------------------------------------- +! Reals with INTENT OUT +!----------------------------------------------------------------------------- +REAL, INTENT(OUT) :: & +g_leaf_phen_dr(land_pts,nnpft) + ! Mean phenology-driven leaf turnover rate for driving vegetation + ! dynamics (s-1). + +!----------------------------------------------------------------------------- +!Local Vars +!----------------------------------------------------------------------------- +REAL :: & +gam_trif, & + ! Inverse vegetation dynamics coupling timestep ((360d)-1). +lai_bal_dummy(land_pts), & + ! Dummy lai to pass into phenol routine, gets around bug where + ! FORTRAN does not accept veg_state%lai_bal as optional + ! argument. +g_leaf_day(land_pts) + ! Mean leaf turnover rate driving phenology, diagnosed from + ! the physiological leaf turnover accumulated since the + ! previous phenology call. ((360d)-1) + +INTEGER :: l,n,k + ! Loop counters. + +!End of headers + +g_leaf_phen_dr(:,:) = 0.0 + +!----------------------------------------------------------------------------- +! Work out the phenology at its own timestep, appending to the accumulated +! leaf turnover rates. This is called independently of the vegetation +! dynamics timestep, mirroring the phenology/TRIFFID timestep split in +! veg-veg2a_jls_mod. +!----------------------------------------------------------------------------- +IF (l_phenol .AND. MOD(a_step,veg3_ctrl%nstep_phen) == 0) THEN + + veg_state%phen(:,:) = 1. + lai_bal_dummy(:) = 0. + g_leaf_day(:) = 0. + + DO n = 1,nnpft + lai_bal_dummy(:) = veg_state%lai_bal(:,n) + + ! Diagnose the mean leaf turnover rate driving phenology over the + ! elapsed phenology period, mirroring g_leaf_day in veg-veg2a_jls_mod. + DO k = 1,veg_index_pts + l = veg_index(k) + g_leaf_day(l) = veg_state%g_leaf_acc(l,n) / veg3_ctrl%dt_phen_360d + END DO + + CALL phenol(land_pts,veg_index_pts,n,veg_index,veg3_ctrl%dt_phen_360d, & + g_leaf_day,veg_state%canht(:,n), & + veg_state%lai(:,n),veg_state%g_leaf_phen(:,n),lai_bal_dummy) + + DO k = 1,veg_index_pts + l = veg_index(k) + + ! Save the diagnosed LAI immediately following the phenology update, + veg_state%lai_phen(l,n) = veg_state%lai(l,n) + + ! Accumulate the mean phenological leaf turnover rate for driving + ! vegetation dynamics. + veg_state%g_leaf_phen_acc(l,n) = veg_state%g_leaf_phen_acc(l,n) & + + veg_state%g_leaf_phen(l,n) * & + veg3_ctrl%dt_phen_360d + + ! Reset the accumulated physiological leaf turnover ready for the + ! next phenology period. + veg_state%g_leaf_acc(l,n) = 0.0 + + IF (veg_state%lai_bal(l,n) > 0) veg_state%phen(l,n) = & + veg_state%lai(l,n)/veg_state%lai_bal(l,n) + + END DO + + END DO + +END IF + +!----------------------------------------------------------------------------- +! On the vegetation dynamics timestep, diagnose the mean phenology-driven +! leaf turnover rate that will drive litterfall and vegetation dynamics. +!----------------------------------------------------------------------------- +IF (asteps_since_triffid == veg3_ctrl%nstep_trif) THEN + + ! Calculate the inverse vegetation dynamics coupling timestep. + gam_trif = 360.0 / REAL(veg3_ctrl%triffid_period) + + DO n = 1,nnpft + DO k = 1,veg_index_pts + l = veg_index(k) + + IF (l_phenol) THEN + ! Diagnose the mean phenological leaf turnover rate over the + ! coupling period and convert to JULES-standard per-second units. + g_leaf_phen_dr(l,n) = veg_state%g_leaf_phen_acc(l,n) * gam_trif / & + (rsec_per_day * 360.0) + + ! Reset the accumulated phenological turnover ready for the next + ! coupling period. + veg_state%g_leaf_phen_acc(l,n) = 0.0 + ELSE + ! No phenology - fall back to the raw accumulated physiological + ! leaf turnover rate, as in veg-veg2a_jls_mod. + g_leaf_phen_dr(l,n) = veg_state%g_leaf_acc(l,n) * gam_trif / & + (rsec_per_day * 360.0) + + veg_state%g_leaf_acc(l,n) = 0.0 + END IF + + ! Ensure the turnover will not remove more leaf than is present over + ! the vegetation dynamics timestep. If it does, reduce the rate so + ! that the turnover does not exceed the current LAI. + IF (veg_state%lai(l,n) > 0.0) THEN + IF (g_leaf_phen_dr(l,n) * veg3_ctrl%dt_red > 1.0) THEN + g_leaf_phen_dr(l,n) = 1.0 / veg3_ctrl%dt_red + + END IF + ELSE + g_leaf_phen_dr(l,n) = 0.0 + + END IF + + END DO + END DO + +END IF + +END SUBROUTINE veg3_phenol_couple + END MODULE next_gen_biogeochem_mod diff --git a/src/control/shared/surf_couple_extra_mod.F90 b/src/control/shared/surf_couple_extra_mod.F90 index 4b52d3c1..0b7628af 100644 --- a/src/control/shared/surf_couple_extra_mod.F90 +++ b/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 @@ -1073,7 +1076,7 @@ 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, & @@ -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) diff --git a/src/control/shared/veg3_field_mod.F90 b/src/control/shared/veg3_field_mod.F90 index 4e710113..9a2245d4 100644 --- a/src/control/shared/veg3_field_mod.F90 +++ b/src/control/shared/veg3_field_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_field_mod @@ -23,19 +26,62 @@ MODULE veg3_field_mod ! Structure to keep the vegetation state variables and fields TYPE :: veg_state_type - REAL, ALLOCATABLE :: & + REAL, POINTER :: & leafC(:,:), & + ! PFT leaf carbon per PFT area fraction. (kg C m-2) rootC(:,:), & + ! PFT root carbon per PFT area fraction. (kg C m-2) woodC(:,:), & + ! PFT woody carbon per PFT area fraction. (kg C m-2) lai_bal(:,:), & + ! Balanced LAI. (m2 m-2) vegCpft(:,:), & + ! Total PFT carbon density per PFT area fraction. (kg C m-2) vegC(:), & - npp_acc(:,:) - REAL, POINTER :: & + ! Gridbox mean vegetation carbon. (kg C m-2) + npp_acc(:,:), & + ! Accumulated NPP. (kg C m-2 s-1) + npp_dr_out(:,:), & + ! A diagnostic NPP variable driving RED. (kg C m-2 (360d)-1) frac(:,:), & + ! Gridbox area fraction of each tile type including PFTs. phen(:,:), & + ! PFT phenology state variable, diagnosed from lai and lai_bal. lai(:,:), & - canht(:,:) + ! Leaf area index. Associated with trif_vars%lai_pft (m2 m-2). + canht(:,:), & + ! Canopy height. Associated with trif_vars%canht_pft. (m) + g_leaf_phen(:,:), & + ! PFT leaf carbon mass density turnover. ((360d)-1) + g_leaf(:,:), & + ! Pointer necessary to get the diagnosed g_leaf rate from the + ! rest of JULES. ((360d)-1) + lai_phen(:,:), & + ! Diagnostic LAI immediately following the phenology update. + ! Associated with trifctl%lai_phen_pft. (m2 m-2) + g_leaf_acc(:,:), & + ! Accumulated leaf turnover rate since the last phenology call. + ! ((360d)-1) + g_leaf_phen_acc(:,:), & + ! Accumulated mean phenological leaf turnover rate since the + ! last vegetation dynamics call. ((360d)-1) + leaf_litC(:,:), & + ! Leaf litter carbon flux per PFT fraction. (kg C m-2 (360d)-1) + root_litC(:,:), & + ! Root litter carbon flux per PFT fraction. (kg C m-2 (360d)-1) + wood_litC(:,:), & + ! Wood litter carbon flux per PFT fraction. (kg C m-2 (360d)-1) + litCpft(:,:), & + ! Total litter carbon flux per PFT fraction, also includes any + ! litter from dynamics (e.g., mortality). (kg C m-2 (360d)-1) + litC(:) + ! Total litter carbon flux per gridbox. (kg C m-2 (360d)-1) + + REAL, ALLOCATABLE :: & + mort_litC(:,:) + ! Mortality/demographic litter carbon flux from vegetation + ! dynamics, normalised per unit PFT canopy area (kg C m-2 s-1). + END TYPE veg_state_type ! Structure to keep the RED state variables and fields @@ -68,13 +114,14 @@ MODULE veg3_field_mod PRIVATE !Expose routines -PUBLIC :: veg3_field_init, veg3_field_allocate, red_veg3_couple +PUBLIC :: veg3_field_init, veg3_field_allocate, veg3_field_assoc, & + red_veg3_couple !Expose data PUBLIC :: veg_state, red_state !Expose data structures -PUBLIC :: veg_state_type,red_state_type +PUBLIC :: veg_state_type, red_state_type !Allow external code to read but not write !PROTECTED :: @@ -100,19 +147,45 @@ SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) ALLOCATE(veg_state%lai ( land_pts, nnpft) ) ALLOCATE(veg_state%phen ( land_pts, nnpft) ) ALLOCATE(veg_state%npp_acc ( land_pts, nnpft) ) - +ALLOCATE(veg_state%npp_dr_out ( land_pts, nnpft) ) +ALLOCATE(veg_state%g_leaf_phen ( land_pts, nnpft) ) +ALLOCATE(veg_state%g_leaf ( land_pts, nnpft) ) +ALLOCATE(veg_state%lai_phen ( land_pts, nnpft) ) +ALLOCATE(veg_state%g_leaf_acc ( land_pts, nnpft) ) +ALLOCATE(veg_state%g_leaf_phen_acc ( land_pts, nnpft) ) ALLOCATE(veg_state%frac ( land_pts, nsurft) ) ALLOCATE(veg_state%vegC ( land_pts) ) +ALLOCATE(veg_state%leaf_litC ( land_pts, nnpft) ) +ALLOCATE(veg_state%root_litC ( land_pts, nnpft) ) +ALLOCATE(veg_state%wood_litC ( land_pts, nnpft) ) +ALLOCATE(veg_state%litCpft ( land_pts, nnpft) ) +ALLOCATE(veg_state%litC ( land_pts) ) +ALLOCATE(veg_state%mort_litC ( land_pts, nnpft) ) !Initialise -veg_state%leafC(:,:) = 0.0 -veg_state%rootC(:,:) = 0.0 -veg_state%woodC(:,:) = 0.0 -veg_state%vegCpft(:,:) = 0.0 -veg_state%lai_bal(:,:) = 0.0 -veg_state%vegC(:) = 0.0 -veg_state%npp_acc(:,:) = 0.0 - +veg_state%leafC(:,:) = 0.0 +veg_state%rootC(:,:) = 0.0 +veg_state%woodC(:,:) = 0.0 +veg_state%vegCpft(:,:) = 0.0 +veg_state%lai_bal(:,:) = 0.0 +veg_state%canht(:,:) = 0.0 +veg_state%lai(:,:) = 0.0 +veg_state%phen(:,:) = 0.0 +veg_state%npp_acc(:,:) = 0.0 +veg_state%npp_dr_out(:,:) = 0.0 +veg_state%g_leaf_phen(:,:) = 0.0 +veg_state%g_leaf(:,:) = 0.0 +veg_state%lai_phen(:,:) = 0.0 +veg_state%g_leaf_acc(:,:) = 0.0 +veg_state%g_leaf_phen_acc(:,:) = 0.0 +veg_state%frac(:,:) = 0.0 +veg_state%vegC(:) = 0.0 +veg_state%leaf_litC(:,:) = 0.0 +veg_state%root_litC(:,:) = 0.0 +veg_state%wood_litC(:,:) = 0.0 +veg_state%litCpft(:,:) = 0.0 +veg_state%litC(:) = 0.0 +veg_state%mort_litC(:,:) = 0.0 ! RED @@ -126,33 +199,91 @@ SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) ALLOCATE(red_state%mort (land_pts, nnpft, nmasst )) ! Initialise red_data_type - red_state%mass_mass(:,:) = 0.0 red_state%ht_mass(:,:) = 0.0 -red_state%lai_bal_mass(:,: ) = 0.0 +red_state%lai_bal_mass(:,:) = 0.0 red_state%crwn_area_mass(:,:) = 0.0 red_state%g_mass_scale(:,:) = 0.0 +red_state%plantNumDensity(:,:,:) = 0.0 red_state%mort(:,:,:) = 0.0 - RETURN END SUBROUTINE veg3_field_allocate +!------------------------------------------------------------------------------- +SUBROUTINE veg3_field_assoc(progs, ainfo) + +! Initial code to associate the veg3 and red fields to the rest of JULES +! This new routine moves out the pointers from veg3_set_fields to here +! to keep clean and more in line with in init.F90. + +! We still need to run with l_triffid, but this will need to be addressed in a +! future revision to allow for veg3 to run fully independently of the switch. + +USE jules_vegetation_mod, ONLY: l_triffid +USE jules_fields_mod, ONLY: trifctl_data, trif_vars_data + +IMPLICIT NONE + +TYPE(progs_type), INTENT(IN) :: progs +TYPE(ainfo_type), INTENT(IN) :: ainfo +! End of header +!------------------------------------------------------------------------------- +IF (l_red .AND. l_triffid) THEN + + ! Set pointers to the prognostic fields + ! Note: veg_state%phen is not associated to an external target here - it is + ! not an independent prognostic. It is diagnosed each phenology call as + ! lai / lai_bal (see veg3_phenol_couple in next_gen_biogeochem_mod), so it + ! remains the array allocated locally in veg3_field_allocate. + red_state%plantNumDensity => progs%plantNumDensity + veg_state%frac => ainfo%frac_surft + veg_state%canht => progs%canht_pft + veg_state%lai => progs%lai_pft + + ! Additional TRIFFID fields mimicking the pointer associations set up in + ! trif_vars_assoc and trifctl_assoc, so that veg3/RED diagnostics feed + ! straight back into the standard JULES carbon diagnostics and outputs. + + ! Firstly trif_vars_data + veg_state%lai_bal => trif_vars_data%lai_bal_pft + veg_state%leafC => trif_vars_data%leafc_pft + veg_state%rootC => trif_vars_data%rootc_pft + veg_state%woodC => trif_vars_data%woodc_pft + veg_state%leaf_litC => trif_vars_data%leaf_litc_pft + veg_state%root_litC => trif_vars_data%root_litc_pft + veg_state%wood_litC => trif_vars_data%wood_litc_pft + + ! Next trifctl_data + veg_state%npp_dr_out => trifctl_data%npp_dr_out_pft + veg_state%g_leaf => trifctl_data%g_leaf_pft + veg_state%g_leaf_phen => trifctl_data%g_leaf_phen_pft + veg_state%lai_phen => trifctl_data%lai_phen_pft + veg_state%g_leaf_acc => trifctl_data%g_leaf_acc_pft + veg_state%g_leaf_phen_acc => trifctl_data%g_leaf_phen_acc_pft + veg_state%vegCpft => trifctl_data%c_veg_pft + veg_state%litCpft => trifctl_data%lit_c_pft + veg_state%vegC => trifctl_data%cv_gb + veg_state%litC => trifctl_data%lit_c_mn_gb + +END IF + +RETURN + +END SUBROUTINE veg3_field_assoc + + !------------------------------------------------------------------------------- SUBROUTINE veg3_set_fields(land_pts,nsurft,nnpft,nmasst,ainfo,progs) !Source parms, etc from io modules - these should be available on and offline ! 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 USE conversions_mod, ONLY: rsec_per_day USE jules_surface_types_mod, ONLY: soil -!Functions -USE calc_c_comps_triffid_mod, ONLY: calc_c_comps_triffid - IMPLICIT NONE INTEGER, INTENT(IN) :: land_pts, nsurft, nnpft, nmasst @@ -161,38 +292,52 @@ SUBROUTINE veg3_set_fields(land_pts,nsurft,nnpft,nmasst,ainfo,progs) TYPE(progs_type), INTENT(IN) :: progs INTEGER :: l,n,k +REAL :: frac_check(land_pts,nnpft) !End of header !----------------------------------------------------------------------------- IF (l_red .AND. l_triffid) THEN + veg_state%phen(:,:) = 1.0 + frac_check(:,:) = 0.0 - ! Set pointers to the prognostic fields - red_state%plantNumDensity => progs%plantNumDensity - veg_state%frac => ainfo%frac_surft - veg_state%canht => progs%canht_pft - veg_state%lai => progs%lai_pft + DO l = 1, land_pts + DO n = 1, nnpft + DO k = 1, nmasst - !----------------------------------------------------------------------------- - ! AJW The following code section temporarily sets up the veg3 model - ! In further revisions these fields will imported from the dump - !----------------------------------------------------------------------------- + IF (k > red_parms%mclass(n)) EXIT + frac_check(l,n) = frac_check(l,n) + red_state%plantNumDensity(l,n,k) & + * red_state%crwn_area_mass(n,k) - veg_state%frac = 0.0 - DO n = 1, nnpft - red_state%plantNumDensity(:,n,1:red_parms%mclass(n)) = & - 0.1 / SUM(red_state%crwn_area_mass(n,:)) - END DO + END DO - veg_state%phen(:,:) = 1.0 - ! phen is a prognostic to be added to model - hardwired here to 1 for now + ! Check to see if frac_check is less than the min fraction, if so, add + ! the necessary amount to the plant number density to balance in the + ! lowest mass class. + IF (frac_check(l,n) < red_parms%frac_min(n)) THEN + red_state%plantNumDensity(l,n,1) = red_state%plantNumDensity(l,n,1) & + + (red_parms%frac_min(n) - frac_check(l,n)) / & + red_state%crwn_area_mass(n,1) + + END IF + + ! Here we estimate the phenology diagnosed from the lai and lai_bal from + ! the dump/intialisation. + IF (veg_state%lai_bal(l,n) > 0.0 .AND. veg_state%lai_bal(l,n) >= & + veg_state%lai(l,n)) THEN + veg_state%phen(l,n) = veg_state%lai(l,n) / veg_state%lai_bal(l,n) + + ELSE + veg_state%phen(l,n) = TINY(0.0) + + END IF + + END DO + END DO ! Initialise veg state from red prognostic CALL red_veg3_couple(ainfo) - ! Temp set rest to bareground - veg_state%frac(:,soil) = 1 - SUM(veg_state%frac(:,:)) - !----------------------------------------------------------------------------- ! AJW < red_parms%mclass(n)) EXIT veg_state%frac(l,n) = veg_state%frac(l,n) & + red_state%plantNumDensity(l,n,k) * red_state%crwn_area_mass(n,k) veg_state%vegCpft(l,n) = veg_state%vegCpft(l,n) & + red_state%plantNumDensity(l,n,k) * red_state%mass_mass(n,k) + veg_state%lai_bal(l,n) = veg_state%lai_bal(l,n) & + + red_state%plantNumDensity(l,n,k) * red_state%lai_bal_mass(n,k) * & + red_state%crwn_area_mass(n,k) + veg_state%canht(l,n) = veg_state%canht(l,n) + & + red_state%plantNumDensity(l,n,k) * red_state%ht_mass(n,k) * & + red_state%crwn_area_mass(n,k) END DO ! Convert to per m2 plant - veg_state%vegCpft(l,n) = veg_state%vegCpft(l,n) / veg_state%frac(l,n) - veg_state%lai_bal(l,n) = veg_state%lai_bal(l,n) / veg_state%frac(l,n) - veg_state%canht(l,n) = veg_state%canht(l,n) / veg_state%frac(l,n) + IF (veg_state%frac(l,n) > 0.0) THEN + veg_state%vegCpft(l,n) = veg_state%vegCpft(l,n) / veg_state%frac(l,n) + veg_state%lai_bal(l,n) = veg_state%lai_bal(l,n) / veg_state%frac(l,n) + veg_state%canht(l,n) = veg_state%canht(l,n) / veg_state%frac(l,n) + ELSE + + ! If frac is zero, set mean to lowest mass class value + veg_state%vegCpft(l,n) = red_state%mass_mass(n,1) + veg_state%lai_bal(l,n) = red_state%lai_bal_mass(n,1) + veg_state%canht(l,n) = red_state%ht_mass(n,1) + END IF veg_state%lai(l,n) = veg_state%phen(l,n) * veg_state%lai_bal(l,n) @@ -331,14 +496,95 @@ SUBROUTINE red_veg3_couple(ainfo) ! Wood carbon balance of total minus leaf and root veg_state%woodC(l,n) = veg_state%vegCpft(l,n) - veg_state%leafC(l,n) & - veg_state%rootC(l,n) + + ! Mean fraction over the update, approximating the changing canopy area + ! as vegetation dynamics were applied. + frac_mid(l,n) = 0.5 * (frac_old(l,n) + veg_state%frac(l,n)) + + ! NPP and its derived litter fluxes are normalised per unit PFT canopy + ! area, but that area has since changed from frac_old. Rescale them + ! onto the frac_mid basis to limit the error this introduces. + CALL veg3_implicit_flux_frac_adj(frac_old(l,n), veg_state%frac(l,n), & + veg_state%npp_dr_out(l,n), & + veg_state%leaf_litC(l,n), & + veg_state%root_litC(l,n), & + veg_state%wood_litC(l,n), & + veg_state%mort_litC(l,n)) + + ! Aggregate the leaf/root/wood turnover litter and the mortality/ + ! demographic litter for the total litter flux per PFT fraction. + veg_state%litCpft(l,n) = (veg_state%leaf_litC(l,n) + & + veg_state%root_litC(l,n) + & + veg_state%wood_litC(l,n) + & + veg_state%mort_litC(l,n) * rsec_per_day * & + 360.0) + ! Update bare soil + veg_state%frac(l,soil) = MAX(0.0, 1.0 - SUM(veg_state%frac(l,1:nnpft))) + END DO + END DO !Final aggregation to gridbox for vegetation carbon for diagnostic purposes veg_state%vegC = pfttiles_to_gbm(veg_state%vegCpft,ainfo,frac_surft_in & = veg_state%frac) +! Aggregate the per-PFT litter contributions for the gridbox total +veg_state%litC(:) = pfttiles_to_gbm(veg_state%litCpft,ainfo,frac_surft_in & + = veg_state%frac) + RETURN END SUBROUTINE red_veg3_couple !----------------------------------------------------------------------------- + +SUBROUTINE veg3_implicit_flux_frac_adj(frac_old, frac_new, npp, leaf_litC, & + root_litC, wood_litC, mort_litC) +!----------------------------------------------------------------------------- +! npp/leaf_litC/root_litC/wood_litC/mort_litC are normalised per unit PFT +! canopy area, i.e. on frac_old. If the PFT fraction has since moved to +! frac_new, that normalisation is stale, so rescale each flux onto the mean +! of frac_old and frac_new - an implicit approximation to the fraction +! evolving smoothly over the timestep - to limit the resulting error. +!----------------------------------------------------------------------------- + +IMPLICIT NONE + +REAL, INTENT(IN) :: frac_old + ! PFT fraction that npp/leaf_litC/root_litC/wood_litC/mort_litC are + ! normalised on, prior to this adjustment. +REAL, INTENT(IN) :: frac_new + ! Updated PFT fraction, following vegetation dynamics. + +REAL, INTENT(IN OUT) :: npp +REAL, INTENT(IN OUT) :: leaf_litC +REAL, INTENT(IN OUT) :: root_litC +REAL, INTENT(IN OUT) :: wood_litC +REAL, INTENT(IN OUT) :: mort_litC + ! Fluxes normalised per unit PFT canopy area (on the frac_old basis on + ! entry), rescaled onto the frac_mid basis on exit. + +REAL :: frac_mid + ! Mean of frac_old and frac_new; the basis the fluxes are rescaled onto. +REAL :: flux_scale + ! Ratio used to rescale each flux from the frac_old basis onto the + ! frac_mid basis. + +!End of header + +frac_mid = 0.5 * (frac_old + frac_new) + +IF (frac_mid > 0.0) THEN + flux_scale = frac_old / frac_mid +ELSE + flux_scale = 0.0 +END IF + +npp = npp * flux_scale +leaf_litC = leaf_litC * flux_scale +root_litC = root_litC * flux_scale +wood_litC = wood_litC * flux_scale +mort_litC = mort_litC * flux_scale + +END SUBROUTINE veg3_implicit_flux_frac_adj +!----------------------------------------------------------------------------- END MODULE veg3_field_mod diff --git a/src/initialisation/standalone/init_mod.F90 b/src/initialisation/standalone/init_mod.F90 index f94cdc16..81be9f59 100644 --- a/src/initialisation/standalone/init_mod.F90 +++ b/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) CALL chemvars_assoc(chemvars,chemvars_data) CALL water_resources_assoc(water_resources,water_resources_data) CALL wtrac_jls_assoc(wtrac_jls,wtrac_jls_data) diff --git a/src/initialisation/standalone/init_vars_tmp_mod.F90 b/src/initialisation/standalone/init_vars_tmp_mod.F90 index d6a27cb8..e610deb9 100644 --- a/src/initialisation/standalone/init_vars_tmp_mod.F90 +++ b/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 diff --git a/src/science/params/veg3_param_mod.F90 b/src/science/params/veg3_param_mod.F90 index e7a44513..656874af 100644 --- a/src/science/params/veg3_param_mod.F90 +++ b/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 @@ -141,7 +145,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, & @@ -174,9 +179,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 +210,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 diff --git a/src/science/surface/cancap_jls.F90 b/src/science/surface/cancap_jls.F90 index a4f326ca..3d939488 100644 --- a/src/science/surface/cancap_jls.F90 +++ b/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. diff --git a/src/science/surface/sf_stom_jls_mod.F90 b/src/science/surface/sf_stom_jls_mod.F90 index 23b46a4d..2a0b8a6a 100644 --- a/src/science/surface/sf_stom_jls_mod.F90 +++ b/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. diff --git a/src/science/vegetation/phenol_jls.F90 b/src/science/vegetation/phenol_jls.F90 index 04ff00b6..b83fe003 100644 --- a/src/science/vegetation/phenol_jls.F90 +++ b/src/science/vegetation/phenol_jls.F90 @@ -6,10 +6,10 @@ MODULE phenol_mod -USE um_types, ONLY: real_jlslsm IMPLICIT NONE + CHARACTER(LEN=*), PARAMETER, PRIVATE :: ModuleName='PHENOL_MOD' CONTAINS @@ -23,12 +23,14 @@ MODULE phenol_mod ! ! --------------------------------------------------------------------------- SUBROUTINE phenol (land_pts, veg_pts, n, veg_index, dtime_phen, g_leaf, ht, & - lai, g_leaf_phen) + lai, g_leaf_phen,lai_bal_red) -USE jules_vegetation_mod, ONLY: l_nitrogen +USE um_types, ONLY: real_jlslsm +USE jules_vegetation_mod, ONLY: l_nitrogen, & + ! TEMP RED + l_red USE pftparm, ONLY: a_wl, a_ws, b_wl, eta_sl, g_leaf_0 USE trif, ONLY: g_grow - USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook @@ -58,6 +60,12 @@ SUBROUTINE phenol (land_pts, veg_pts, n, veg_index, dtime_phen, g_leaf, ht, & ! Rate of leaf turnover (/360days). ht(land_pts) ! Canopy height (m). +!----------------------------------------------------------------------------- +! Optional Arguments with INTENT(IN) +!----------------------------------------------------------------------------- +REAL, INTENT(IN), OPTIONAL :: & + lai_bal_red(land_pts) + ! Balanced growth LAI from veg3. !----------------------------------------------------------------------------- ! Arguments with INTENT(INOUT). @@ -73,6 +81,7 @@ SUBROUTINE phenol (land_pts, veg_pts, n, veg_index, dtime_phen, g_leaf, ht, & g_leaf_phen(land_pts) ! Rate of leaf turnover including leaf phenology (/360days). + !----------------------------------------------------------------------------- ! Local variables. !----------------------------------------------------------------------------- @@ -102,10 +111,17 @@ SUBROUTINE phenol (land_pts, veg_pts, n, veg_index, dtime_phen, g_leaf, ht, & !----------------------------------------------------------------------------- ! Diagnose the phenological state !----------------------------------------------------------------------------- + DO j = 1,veg_pts l = veg_index(j) - lai_bal(l) = (a_ws(n) * eta_sl(n) * ht(l) & - /a_wl(n))**(1.0 / (b_wl(n) - 1.0)) + ! If we're running l_red we use lai_bal as an input, otherwise estimate the + ! classic TRIFFID way. + IF (l_red .AND. PRESENT(lai_bal_red)) THEN + lai_bal(l) = lai_bal_red(l) + ELSE + lai_bal(l) = (a_ws(n) * eta_sl(n) * ht(l) & + /a_wl(n))**(1.0 / (b_wl(n) - 1.0)) + END IF phen(l) = lai(l) / lai_bal(l) END DO @@ -137,7 +153,6 @@ SUBROUTINE phenol (land_pts, veg_pts, n, veg_index, dtime_phen, g_leaf, ht, & !--------------------------------------------------------------------------- phen(l) = phen(l) + dphen lai(l) = phen(l) * lai_bal(l) - END DO IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) diff --git a/src/science/vegetation/veg3_litter_mod.F90 b/src/science/vegetation/veg3_litter_mod.F90 index ef518ff0..cd669bb3 100644 --- a/src/science/vegetation/veg3_litter_mod.F90 +++ b/src/science/vegetation/veg3_litter_mod.F90 @@ -9,7 +9,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_litter_mod IMPLICIT NONE @@ -25,6 +27,8 @@ SUBROUTINE veg3_Litter( & veg_index_pts,veg_index,veg3_ctrl,land_pts,nnpft, & !IN parms litter_parms, & + !IN fields + g_leaf_phen_dr, & !IN state veg_state, & ! OUT Fields @@ -33,6 +37,8 @@ SUBROUTINE veg3_Litter( & ) !Only get the data structures - the data comes through the calling tree +USE conversions_mod, ONLY: rsec_per_day + USE veg3_parm_mod, ONLY: veg3_ctrl_type,litter_parm_type USE veg3_field_mod, ONLY: veg_state_type @@ -50,6 +56,15 @@ SUBROUTINE veg3_Litter( & TYPE(veg3_ctrl_type),INTENT(IN) :: veg3_ctrl TYPE(litter_parm_type),INTENT(IN) :: litter_parms +!----------------------------------------------------------------------------- +! Reals with INTENT IN +!----------------------------------------------------------------------------- +REAL, INTENT(IN) :: & +g_leaf_phen_dr(land_pts,nnpft) + ! Mean phenology-driven leaf turnover rate for driving vegetation + ! dynamics, accumulated and averaged over the vegetation dynamics + ! timestep in veg3_run_ctrl. (s-1) + !----------------------------------------------------------------------------- ! Objects with INTENT INOUT !----------------------------------------------------------------------------- @@ -85,15 +100,22 @@ SUBROUTINE veg3_Litter( & !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) COLLAPSE(2) & !$OMP PRIVATE(l,n,k) & !$OMP SHARED(litter_parms,veg_state,leaf_litter,root_litter,wood_litter, & -!$OMP local_litter,veg_index,veg_index_pts,nnpft) +!$OMP local_litter,veg_index,veg_index_pts,nnpft,g_leaf_phen_dr) DO n = 1, nnpft DO k = 1, veg_index_pts l = veg_index(k) - leaf_litter(l,n) = litter_parms%g_leaf(n) * veg_state%leafC(l,n) + leaf_litter(l,n) = g_leaf_phen_dr(l,n) * veg_state%leafC(l,n) * & + veg_state%phen(l,n) root_litter(l,n) = litter_parms%g_root(n) * veg_state%rootC(l,n) wood_litter(l,n) = litter_parms%g_wood(n) * veg_state%woodC(l,n) + ! Copy across to veg_state%leaf_litC, veg_state%root_litC, veg_state%wood_litC + ! Convert to kg C m-2 (360d)-1 + veg_state%leaf_litC(l,n) = leaf_litter(l,n) * rsec_per_day * 360. + veg_state%root_litC(l,n) = root_litter(l,n) * rsec_per_day * 360. + veg_state%wood_litC(l,n) = wood_litter(l,n) * rsec_per_day * 360. + local_litter(l,n) = leaf_litter(l,n) + root_litter(l,n) + wood_litter(l,n) END DO diff --git a/src/science/vegetation/veg3_red_dynamic_mod.F90 b/src/science/vegetation/veg3_red_dynamic_mod.F90 index 84b43607..f355e2f0 100644 --- a/src/science/vegetation/veg3_red_dynamic_mod.F90 +++ b/src/science/vegetation/veg3_red_dynamic_mod.F90 @@ -9,6 +9,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_red_dynamic_mod @@ -31,9 +34,7 @@ SUBROUTINE veg3_red_dynamic( & !IN fields growth,mort_add, & !IN state - veg_state,red_state, & - ! OUT Fields - demographic_lit & + veg_state,red_state & !OUT Diagnostics ) @@ -67,184 +68,401 @@ SUBROUTINE veg3_red_dynamic( & dt ! Dynamic vegetation time-step (s) -!----------------------------------------------------------------------------- -! Reals with INTENT OUT -!----------------------------------------------------------------------------- -REAL, INTENT(OUT) :: & -demographic_lit(land_pts,nnpft) - ! Demographic litter across PFT area. (kgC/m2/s) - !----------------------------------------------------------------------------- !Local Vars !----------------------------------------------------------------------------- INTEGER ::l,n,k,j REAL :: & -p(land_pts,nnpft), & - ! The total PFT carbon assimilate across the gridbox. (kgC/m2/s) P_s(land_pts,nnpft), & - ! Total gridbox carbon assimilate devoted to recruitment. (kgC/m2/s) -g(land_pts,nnpft), & - ! Total gridbox carbon assimilate devoted to vegetation structural growth. (kgC/m2/s) + ! Total gridbox carbon assimilate devoted to recruitment. + ! (kgC/m2/s) g0(land_pts,nnpft), & - ! Boundary growth for an individual member of the smallest mass cohort. (kgC/year) -g_mass(land_pts,nnpft,nmasst), & - ! Individual growth across the mass cohorts. (kgC/s) -plantNumDensity_g_sum(land_pts,nnpft), & - ! Summation of the relative cohort contribution towards the total PFT assimilate (/m2) -frac_shade(land_pts,nnpft), & - ! Competitive shading of seedlings in each PFT -dplantNumDensity_dt(land_pts,nnpft,nmasst), & - ! Net rate of change of population density within each mass cohort. (m2/s) -flux_in(land_pts,nnpft,nmasst), & - ! Rate of change of population growing into a mass cohort. (m2/s) -flux_out(land_pts,nnpft,nmasst), & - ! Rate of change of population growing out of a mass cohort. (m2/s) -frac_check(land_pts,nnpft) - ! The difference between the minimum vegetation fraction and the updated fraction. (-) + ! Boundary growth for an individual member of the smallest mass + ! cohort. (kgC/year) +neg_growth(land_pts,nnpft), & + ! Uniform PFT-carbon-proportional loss rate for when growth is + ! negative (e.g., local litterfall exceeds NPP). Drives a + ! downward mass-class shrinkage flux, and is only applied as a + ! mortality rate for the lowest mass class. (/s) +frac_shade(land_pts,nnpft) + ! Competitive shading of seedlings in each PFT. !End of headers ! Initialise vars -demographic_lit(:,:) = 0.0 -p(:,:) = 0.0 +veg_state%mort_litC(:,:) = 0.0 P_s(:,:) = 0.0 -g(:,:) = 0.0 g0(:,:) = 0.0 -g_mass(:,:,:) = 0.0 -plantNumDensity_g_sum(:,:) = 0.0 +neg_growth(:,:) = 0.0 frac_shade(:,:) = 0.0 -frac_check(:,:) = 0.0 -dplantNumDensity_dt(:,:,:) = 0.0 -flux_in(:,:,:) = 0.0 -flux_out(:,:,:) = 0.0 - -! Loop finds the growth boundary condition, the fraction of shade for seedlings -DO k = 1,nmasst +! Dynamic demographic loop to update the number density of each PFT across the +! PFT mass classes. +DO l = 1,land_pts DO n = 1,nnpft - DO l = 1,land_pts - IF (k < red_parms%mclass(n)) THEN - ! Sum product of the number density and the allometric scaling - plantNumDensity_g_sum(l,n) = plantNumDensity_g_sum(l,n) & - + red_state%plantNumDensity(l,n,k) * red_state%g_mass_scale(n,k) - - ELSE IF (k == red_parms%mclass(n)) THEN - plantNumDensity_g_sum(l,n) = plantNumDensity_g_sum(l,n) & - + red_state%plantNumDensity(l,n,k) * red_state%g_mass_scale(n,k) - - ! Partition the growth - p(l,n) = veg_state%frac(l,n) * growth(l,n) - P_s(l,n) = red_parms%alpha_recrt(n) * p(l,n) - g(l,n) = (1.0 - red_parms%alpha_recrt(n)) * p(l,n) - - IF (plantNumDensity_g_sum(l,n) > 0) THEN - g0(l,n) = g(l,n) / plantNumDensity_g_sum(l,n) - END IF - - IF (growth(l,n)<0) THEN - p(l,n) = 0.0 - P_s(l,n) = 0.0 - g(l,n) = 0.0 - END IF - - !Extra loop through to calculate shading - DO j = 1, nnpft - frac_shade(l,n) = MIN(1.0,frac_shade(l,n) + red_parms%comp_coef(n,j) & - * veg_state%frac(l,j)) - END DO + ! Call to partition the PFT growth onto the mass class structure. + CALL growth_onto_mass_class( & + !IN sizing + red_parms%mclass(n), & + !IN PFT parameters + red_parms%alpha_recrt(n), & + !IN fields + veg_state%frac(l,n),growth(l,n),veg_state%vegCpft(l,n), & + !IN mass-cohort properties + red_state%plantNumDensity(l,n,1:red_parms%mclass(n)), & + red_state%g_mass_scale(n,1:red_parms%mclass(n)), & + !OUT fields + P_s(l,n),g0(l,n),neg_growth(l,n) & + ) + + !Estimate the inter-PFT competition. + DO j = 1, nnpft + frac_shade(l,n) = MIN(1.0,frac_shade(l,n) + red_parms%comp_coef(n,j) & + * veg_state%frac(l,j)) + END DO + ! Call to update the PFT number density. + CALL update_pft_size_structure( & + !IN sizing + red_parms%mclass(n), & + !IN Control vars + dt, & + !IN PFT parameters + red_parms%mort_base(n),red_parms%frac_min(n), & + !IN fields + mort_add(l,n,1:red_parms%mclass(n)),growth(l,n),neg_growth(l,n), & + P_s(l,n),g0(l,n),frac_shade(l,n), & + !IN mass-cohort properties + red_state%g_mass_scale(n,1:red_parms%mclass(n)), & + red_state%mass_mass(n,1:red_parms%mclass(n)), & + red_state%crwn_area_mass(n,1:red_parms%mclass(n)), & + !INOUT state + red_state%plantNumDensity(l,n,1:red_parms%mclass(n)), & + red_state%mort(l,n,1:red_parms%mclass(n)), & + !OUT diagnostics + veg_state%mort_litC(l,n) & + ) + + ! Divide mort_litC by the PFT fraction (not restimated yet) + IF (veg_state%frac(l,n) > 0.0) THEN + veg_state%mort_litC(l,n) = veg_state%mort_litC(l,n) / & + veg_state%frac(l,n) + ELSE + veg_state%mort_litC(l,n) = 0.0 END IF - END DO + END DO END DO -! Loop updates the demographic distribution. -DO k = 1, nmasst - DO n = 1, nnpft - DO l = 1, land_pts +END SUBROUTINE veg3_red_dynamic +!----------------------------------------------------------------------------- - g_mass(l,n,k) = g0(l,n) * red_state%g_mass_scale(n,k) +!----------------------------------------------------------------------------- +SUBROUTINE growth_onto_mass_class( & + !IN sizing + mclass, & + !IN PFT parameters + alpha_recrt, & + !IN fields + frac,growth,vegCpft, & + !IN mass-cohort properties + plantNumDensity,g_mass_scale, & + !OUT fields + P_s,g0,neg_growth & + ) - IF (growth(l,n)<0) THEN - ! For negative growth we add to the mortality rate - red_state%mort(l,n,k) = red_parms%mort_base(n) + mort_add(l,n,k) & - - g_mass(l,n,k) / red_state%mass_mass(n,k) - g_mass(l,n,k) = 0.0 - ELSE - red_state%mort(l,n,k) = red_parms%mort_base(n) + mort_add(l,n,k) - END IF +IMPLICIT NONE - IF (k == 1) THEN - ! Seedling flux - flux_in(l,n,k) = P_s(l,n) / red_state%mass_mass(n,k) & - * (1.0 - frac_shade(l,n)) - demographic_lit(l,n) = demographic_lit(l,n) + frac_shade(l,n) * P_s(l,n) - ELSE - ! Flux into mass class - flux_in(l,n,k) = flux_out(l,n,k-1) - END IF +!----------------------------------------------------------------------------- +! Integers with INTENT IN +!----------------------------------------------------------------------------- +INTEGER, INTENT(IN) :: mclass + ! Number of mass classes for this PFT. - IF (k == red_parms%mclass(n)) THEN - ! Truncate growth at the top mass class - flux_out(l,n,k) = 0.0 - demographic_lit(l,n) = demographic_lit(l,n) & - + red_state%plantNumDensity(l,n,k) * g_mass(l,n,k) - ELSE IF (k < red_parms%mclass(n)) THEN - flux_out(l,n,k) = red_state%plantNumDensity(l,n,k) & - * g_mass(l,n,k) / (red_state%mass_mass(n,k+1) & - - red_state%mass_mass(n,k)) - END IF - ! Update the number density - - !Only update over the given mass classes - IF (k <= red_parms%mclass(n)) THEN - dplantNumDensity_dt(l,n,k) = flux_in(l,n,k) - flux_out(l,n,k) & - - red_state%mort(l,n,k) * red_state%plantNumDensity(l,n,k) - - !Prevent the mass class from being exhausted over a timestep - IF (red_state%plantNumDensity(l,n,k) & - + (dplantNumDensity_dt(l,n,k) * dt) < 0.0 ) THEN - dplantNumDensity_dt(l,n,k) = -red_state%plantNumDensity(l,n,k) / dt - demographic_lit(l,n) = demographic_lit(l,n) & - + (dplantNumDensity_dt(l,n,k) - flux_out(l,n,k)) & - * red_state%mass_mass(n,k) - ELSE - demographic_lit(l,n) = demographic_lit(l,n) + red_state%mort(l,n,k) & - * red_state%plantNumDensity(l,n,k) * red_state%mass_mass(n,k) - - END IF - red_state%plantNumDensity(l,n,k) = red_state%plantNumDensity(l,n,k) & - + dplantNumDensity_dt(l,n,k) * dt - frac_check(l,n) = frac_check(l,n) + red_state%plantNumDensity(l,n,k) & - * red_state%crwn_area_mass(n,k) - - ! If the resultant vegetation fraction is less than the minimum - ! fraction, add trees to the lowest mass class to make up the difference - IF (k == red_parms%mclass(n) .AND. frac_check(l,n) & - < red_parms%frac_min(n)) THEN - red_state%plantNumDensity(l,n,1) = red_state%plantNumDensity(l,n,1) & - +(red_parms%frac_min(n) - frac_check(l,n)) & - / red_state%crwn_area_mass(n,1) - demographic_lit(l,n) = demographic_lit(l,n) & - - (red_parms%frac_min(n) - frac_check(l,n)) & - * red_state%mass_mass(n,1) / red_state%crwn_area_mass(n,1) / dt - frac_check(l,n) = red_parms%frac_min(n) - - END IF - - IF (k == red_parms%mclass(n)) THEN - ! Convert demographic litter into correct dimensions (per PFT area - ! rather than per grid-box area) - demographic_lit(l,n) = demographic_lit(l,n) / frac_check(l,n) - END IF - END IF - END DO - END DO +!----------------------------------------------------------------------------- +! Reals with INTENT IN +!----------------------------------------------------------------------------- +REAL, INTENT(IN) :: & +alpha_recrt, & + ! Fraction of PFT growth devoted to recruitment. (-) +frac, & + ! PFT fraction across the gridbox. (-) +growth, & + ! The total carbon assimilate across the PFT area. (kgC/m2/s) +vegCpft, & + ! The total PFT vegetation carbon. (kgC/m2) +plantNumDensity(mclass), & + ! Population density within each mass cohort. (/m2) +g_mass_scale(mclass) + ! Allometric scaling of growth across the mass cohorts. + +!----------------------------------------------------------------------------- +! Reals with INTENT OUT +!----------------------------------------------------------------------------- +REAL, INTENT(OUT) :: & +P_s, & + ! Gridbox carbon assimilate devoted to recruitment. (kgC/m2/s) +g0, & + ! Boundary growth for an individual member of the smallest mass + ! cohort. (kgC/year) +neg_growth + ! Uniform PFT-carbon-proportional loss rate for when growth is + ! negative (e.g., local litterfall exceeds NPP). Drives a + ! downward mass-class shrinkage flux, and is only applied as a + ! mortality rate for the lowest mass class. (/s) + +!----------------------------------------------------------------------------- +!Local Vars +!----------------------------------------------------------------------------- +INTEGER :: k + +REAL :: & +p, & + ! The total PFT carbon assimilate across the gridbox. (kgC/m2/s) +g, & + ! Total gridbox carbon assimilate devoted to vegetation + ! structural growth. (kgC/m2/s) +plantNumDensity_g_sum + ! Summation of the relative cohort contribution towards the + ! total PFT assimilate (/m2) + +!End of headers + +! Initialise vars +P_s = 0.0 +g0 = 0.0 +neg_growth = 0.0 +plantNumDensity_g_sum = 0.0 + +! Sum product of the number density and the allometric scaling +DO k = 1, mclass + plantNumDensity_g_sum = plantNumDensity_g_sum & + + plantNumDensity(k) * g_mass_scale(k) END DO -END SUBROUTINE veg3_red_dynamic +! Partition the growth into recruitment and structural growth +p = frac * growth +P_s = alpha_recrt * p +g = (1.0 - alpha_recrt) * p + +IF (plantNumDensity_g_sum > 0) THEN + g0 = g / plantNumDensity_g_sum +END IF + +IF (growth < 0.0) THEN + P_s = 0.0 + + ! Estimate the uniform PFT-carbon-proportional loss rate required to + ! account for the negative growth (local litterfall exceeding NPP). + IF (vegCpft > 0) THEN + neg_growth = - growth / vegCpft + + END IF + +END IF + +END SUBROUTINE growth_onto_mass_class +!----------------------------------------------------------------------------- + +!----------------------------------------------------------------------------- +SUBROUTINE update_pft_size_structure( & + !IN sizing + mclass, & + !IN Control vars + dt, & + !IN PFT parameters + mort_base,frac_min, & + !IN fields + mort_add,growth,neg_growth,P_s,g0,frac_shade, & + !IN mass-cohort properties + g_mass_scale,mass_mass,crwn_area_mass, & + !INOUT state + plantNumDensity,mort, & + !OUT diagnostics + mort_litC & + ) + +IMPLICIT NONE + +!----------------------------------------------------------------------------- +! Integers with INTENT IN +!----------------------------------------------------------------------------- +INTEGER, INTENT(IN) :: mclass + ! Number of mass classes for this PFT. + +!----------------------------------------------------------------------------- +! Reals with INTENT IN +!----------------------------------------------------------------------------- +REAL, INTENT(IN) :: & +dt, & + ! Dynamic vegetation time-step (s) +mort_base, & + ! Background mortality rate for this PFT. (/s) +frac_min, & + ! Minimum vegetation fraction for this PFT. (-) +mort_add(mclass), & + ! Additional plant mortality across plant mass (/s) +growth, & + ! The total carbon assimilate across the PFT area. (kgC/m2/s) +neg_growth, & + ! Uniform PFT-carbon-proportional loss rate for when growth is + ! negative (e.g., local litterfall exceeds NPP). Drives a + ! downward mass-class shrinkage flux, and is only applied as a + ! mortality rate for the lowest mass class. (/s) +P_s, & + ! Gridbox carbon assimilate devoted to recruitment. (kgC/m2/s) +g0, & + ! Boundary growth for an individual member of the smallest mass + ! cohort. (kgC/year) +frac_shade, & + ! Competitive shading of seedlings in this PFT. +g_mass_scale(mclass), & + ! Allometric scaling of growth across the mass cohorts. +mass_mass(mclass), & + ! Mass of an individual member of each mass cohort. (kgC) +crwn_area_mass(mclass) + ! Crown area of an individual member of each mass cohort. (m2) + +!----------------------------------------------------------------------------- +! Reals with INTENT INOUT +!----------------------------------------------------------------------------- +REAL, INTENT(INOUT) :: & +plantNumDensity(mclass), & + ! Population density within each mass cohort. (/m2) +mort(mclass) + ! Mortality rate within each mass cohort. (/s) + +!----------------------------------------------------------------------------- +! Reals with INTENT OUT +!----------------------------------------------------------------------------- +REAL, INTENT(OUT) :: mort_litC + ! Mortality/demographic litter for this PFT, normalised per + ! unit PFT canopy area. (kgC/m2/s) + +!----------------------------------------------------------------------------- +!Local Vars +!----------------------------------------------------------------------------- +INTEGER :: k + +REAL :: & +g_mass(mclass), & + ! Individual growth across the mass cohorts. (kgC/s) +dplantNumDensity_dt(mclass), & + ! Net rate of change of population density within each mass + ! cohort. (m2/s) +flux_in(mclass), & + ! Rate of change of population growing into a mass cohort. (m2/s) +flux_out(mclass), & + ! Rate of change of population growing out of a mass cohort. + ! (m2/s) +frac_check + ! The difference between the minimum vegetation fraction and the + ! updated fraction. (-) + +!End of headers + +! Initialise vars +mort_litC = 0.0 +frac_check = 0.0 +g_mass(:) = 0.0 +dplantNumDensity_dt(:) = 0.0 +flux_in(:) = 0.0 +flux_out(:) = 0.0 + +DO k = 1, mclass + + IF (growth < 0.0) THEN + ! Negative growth is represented as a downward shrinkage flux of + ! individuals through the mass classes, the lowest mass class has + ! nowhere lower to shrink into, so its share of the loss is instead matched + ! by an equivalent mortality rate. + g_mass(k) = 0.0 + + IF (k == 1) THEN + mort(k) = mort_base + mort_add(k) + neg_growth + ELSE + mort(k) = mort_base + mort_add(k) + END IF + + IF (k < mclass) THEN + ! Flux shrinking down into this class from the class above + flux_in(k) = plantNumDensity(k+1) * neg_growth * mass_mass(k+1) & + / (mass_mass(k+1) - mass_mass(k)) + ELSE + flux_in(k) = 0.0 + END IF + + IF (k > 1) THEN + ! Flux shrinking out of this class into the class below + flux_out(k) = plantNumDensity(k) * neg_growth * mass_mass(k) & + / (mass_mass(k) - mass_mass(k-1)) + ELSE + flux_out(k) = 0.0 + END IF + + ELSE + g_mass(k) = g0 * g_mass_scale(k) + mort(k) = mort_base + mort_add(k) + + IF (k == 1) THEN + ! Seedling flux + flux_in(k) = P_s / mass_mass(k) * (1.0 - frac_shade) + mort_litC = mort_litC + frac_shade * P_s + + ELSE + ! Flux into mass class + flux_in(k) = flux_out(k-1) + + END IF + + IF (k == mclass) THEN + ! Truncate growth at the top mass class + flux_out(k) = 0.0 + mort_litC = mort_litC + plantNumDensity(k) * g_mass(k) + + ELSE + flux_out(k) = plantNumDensity(k) * g_mass(k) / (mass_mass(k+1) & + - mass_mass(k)) + + END IF + + END IF + + ! Update the number density + dplantNumDensity_dt(k) = flux_in(k) - flux_out(k) & + - mort(k) * plantNumDensity(k) + + !Prevent the mass class from being exhausted over a timestep + IF (plantNumDensity(k) & + + (dplantNumDensity_dt(k) * dt) < 0.0 ) THEN + dplantNumDensity_dt(k) = -plantNumDensity(k) / dt + mort_litC = mort_litC & + + (dplantNumDensity_dt(k) - flux_out(k)) & + * mass_mass(k) + + ELSE + mort_litC = mort_litC + mort(k) & + * plantNumDensity(k) * mass_mass(k) + + END IF + + plantNumDensity(k) = plantNumDensity(k) & + + dplantNumDensity_dt(k) * dt + frac_check = frac_check + plantNumDensity(k) & + * crwn_area_mass(k) + +END DO + +! If the resultant vegetation fraction is less than the minimum +! fraction, add trees to the lowest mass class to make up the difference +IF (frac_check < frac_min) THEN + plantNumDensity(1) = plantNumDensity(1) & + +(frac_min - frac_check) / crwn_area_mass(1) +END IF + +END SUBROUTINE update_pft_size_structure !----------------------------------------------------------------------------- END MODULE veg3_red_dynamic_mod From 7af5d8c790768b975ded8d01090d4e9db1758f3c Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Thu, 3 Sep 2026 16:26:55 +0100 Subject: [PATCH 02/14] For #90, fix bug with the re-initillisation of phenological state. Deleted redundant bits of code we don't need to get phenology correctly estimated in veg3_field_mod.F90. Tests now pass umpd3_checker and rose stem tests. --- .../shared/next_gen_biogeochem_mod.F90 | 8 +- src/control/shared/veg3_field_mod.F90 | 256 ++++++++++-------- src/science/vegetation/veg3_litter_mod.F90 | 6 +- .../vegetation/veg3_red_dynamic_mod.F90 | 185 ++++++------- 4 files changed, 236 insertions(+), 219 deletions(-) diff --git a/src/control/shared/next_gen_biogeochem_mod.F90 b/src/control/shared/next_gen_biogeochem_mod.F90 index 584d1153..3c6e75ae 100644 --- a/src/control/shared/next_gen_biogeochem_mod.F90 +++ b/src/control/shared/next_gen_biogeochem_mod.F90 @@ -248,7 +248,7 @@ SUBROUTINE veg3_run_ctrl( & npp_dr = veg_state%npp_acc / veg3_ctrl%dt_red - ! Record driving npp_dr and convert to s-1 -> (360 days)-1 + ! Record driving npp_dr and convert to s-1 -> (360 days)-1 veg_state%npp_dr_out(:,:) = npp_dr * rsec_per_day * 360 ! Reset accumulation to zero - note this can be used to pass a negative flux @@ -386,9 +386,9 @@ SUBROUTINE veg3_phenol_couple( & !----------------------------------------------------------------------------- IF (l_phenol .AND. MOD(a_step,veg3_ctrl%nstep_phen) == 0) THEN - veg_state%phen(:,:) = 1. - lai_bal_dummy(:) = 0. - g_leaf_day(:) = 0. + veg_state%phen(:,:) = 1.0 + lai_bal_dummy(:) = 0.0 + g_leaf_day(:) = 0.0 DO n = 1,nnpft lai_bal_dummy(:) = veg_state%lai_bal(:,n) diff --git a/src/control/shared/veg3_field_mod.F90 b/src/control/shared/veg3_field_mod.F90 index 9a2245d4..9745320f 100644 --- a/src/control/shared/veg3_field_mod.F90 +++ b/src/control/shared/veg3_field_mod.F90 @@ -42,7 +42,7 @@ MODULE veg3_field_mod npp_acc(:,:), & ! Accumulated NPP. (kg C m-2 s-1) npp_dr_out(:,:), & - ! A diagnostic NPP variable driving RED. (kg C m-2 (360d)-1) + ! A diagnostic NPP variable driving RED. (kg C m-2 (360d)-1) frac(:,:), & ! Gridbox area fraction of each tile type including PFTs. phen(:,:), & @@ -54,7 +54,7 @@ MODULE veg3_field_mod g_leaf_phen(:,:), & ! PFT leaf carbon mass density turnover. ((360d)-1) g_leaf(:,:), & - ! Pointer necessary to get the diagnosed g_leaf rate from the + ! Pointer necessary to get the diagnosed g_leaf rate from the ! rest of JULES. ((360d)-1) lai_phen(:,:), & ! Diagnostic LAI immediately following the phenology update. @@ -72,7 +72,7 @@ MODULE veg3_field_mod wood_litC(:,:), & ! Wood litter carbon flux per PFT fraction. (kg C m-2 (360d)-1) litCpft(:,:), & - ! Total litter carbon flux per PFT fraction, also includes any + ! Total litter carbon flux per PFT fraction, also includes any ! litter from dynamics (e.g., mortality). (kg C m-2 (360d)-1) litC(:) ! Total litter carbon flux per gridbox. (kg C m-2 (360d)-1) @@ -115,7 +115,7 @@ MODULE veg3_field_mod !Expose routines PUBLIC :: veg3_field_init, veg3_field_allocate, veg3_field_assoc, & - red_veg3_couple + red_veg3_couple !Expose data PUBLIC :: veg_state, red_state @@ -214,10 +214,10 @@ END SUBROUTINE veg3_field_allocate SUBROUTINE veg3_field_assoc(progs, ainfo) ! Initial code to associate the veg3 and red fields to the rest of JULES -! This new routine moves out the pointers from veg3_set_fields to here +! This new routine moves out the pointers from veg3_set_fields to here ! to keep clean and more in line with in init.F90. -! We still need to run with l_triffid, but this will need to be addressed in a +! We still need to run with l_triffid, but this will need to be addressed in a ! future revision to allow for veg3 to run fully independently of the switch. USE jules_vegetation_mod, ONLY: l_triffid @@ -291,37 +291,55 @@ SUBROUTINE veg3_set_fields(land_pts,nsurft,nnpft,nmasst,ainfo,progs) TYPE(ainfo_type), INTENT(IN OUT) :: ainfo TYPE(progs_type), INTENT(IN) :: progs -INTEGER :: l,n,k -REAL :: frac_check(land_pts,nnpft) +INTEGER :: l,n !End of header !----------------------------------------------------------------------------- IF (l_red .AND. l_triffid) THEN veg_state%phen(:,:) = 1.0 - frac_check(:,:) = 0.0 DO l = 1, land_pts DO n = 1, nnpft - DO k = 1, nmasst - IF (k > red_parms%mclass(n)) EXIT - frac_check(l,n) = frac_check(l,n) + red_state%plantNumDensity(l,n,k) & - * red_state%crwn_area_mass(n,k) - - END DO - - ! Check to see if frac_check is less than the min fraction, if so, add + ! We estimate the PFT physical properties from the prognostic field + ! plantNumDensity for the following reasons: + ! + ! i.) Check to see if plantNumDensity derived frac is less than the + ! minimum fraction, this should introduce a small error in carbon + ! flux but probably only once. + + ! ii.) red_veg3_couple needs to know the phen state variable, which + ! is diagnosed from lai_bal. Both are not prognostic variables, + ! but lai is, so we can infer then phen state from this and + ! plantNumDensity for lai_bal. + + CALL pft_mean_from_mass_class( & + !IN sizing + red_parms%mclass(n), & + !IN mass-cohort properties + red_state%plantNumDensity(l,n,1:red_parms%mclass(n)), & + red_state%mass_mass(n,1:red_parms%mclass(n)), & + red_state%lai_bal_mass(n,1:red_parms%mclass(n)), & + red_state%ht_mass(n,1:red_parms%mclass(n)), & + red_state%crwn_area_mass(n,1:red_parms%mclass(n)), & + !OUT fields + veg_state%frac(l,n),veg_state%vegCpft(l,n),veg_state%lai_bal(l,n), & + veg_state%canht(l,n) & + ) + + ! Check to see if frac is less than the min fraction, if so, add ! the necessary amount to the plant number density to balance in the - ! lowest mass class. - IF (frac_check(l,n) < red_parms%frac_min(n)) THEN + ! lowest mass class, and set the fraction to the minimum directly. + IF (veg_state%frac(l,n) < red_parms%frac_min(n)) THEN red_state%plantNumDensity(l,n,1) = red_state%plantNumDensity(l,n,1) & - + (red_parms%frac_min(n) - frac_check(l,n)) / & + + (red_parms%frac_min(n) - veg_state%frac(l,n)) / & red_state%crwn_area_mass(n,1) + veg_state%frac(l,n) = red_parms%frac_min(n) END IF - ! Here we estimate the phenology diagnosed from the lai and lai_bal from + ! Here we estimate the phenology diagnosed from the lai and lai_bal from ! the dump/intialisation. IF (veg_state%lai_bal(l,n) > 0.0 .AND. veg_state%lai_bal(l,n) >= & veg_state%lai(l,n)) THEN @@ -434,15 +452,11 @@ SUBROUTINE red_veg3_couple(ainfo) TYPE(ainfo_type),INTENT(IN OUT) :: ainfo !Local -INTEGER :: l,n,k ! Index variables. +INTEGER :: l,n ! Index variables. REAL :: frac_old(land_pts,nnpft) ! PFT fraction prior to updating veg_state%frac below, i.e. before the ! change in canopy area resulting from vegetation dynamics. -REAL :: frac_mid(land_pts,nnpft) - ! Mean of the fraction before and after vegetation dynamics, used to - ! convert the per-canopy-area litter fluxes to a gridbox-area flux, - ! approximating the changing canopy area over the timestep. !----------------------------------------------------------------------------- !end of header @@ -460,63 +474,40 @@ SUBROUTINE red_veg3_couple(ainfo) DO n = 1,nnpft DO l = 1,land_pts - ! Convert from plant number on mass classes to area on PFTS - ! and biomass to carbon - DO k = 1,nmasst - IF (k > red_parms%mclass(n)) EXIT - veg_state%frac(l,n) = veg_state%frac(l,n) & - + red_state%plantNumDensity(l,n,k) * red_state%crwn_area_mass(n,k) - veg_state%vegCpft(l,n) = veg_state%vegCpft(l,n) & - + red_state%plantNumDensity(l,n,k) * red_state%mass_mass(n,k) - veg_state%lai_bal(l,n) = veg_state%lai_bal(l,n) & - + red_state%plantNumDensity(l,n,k) * red_state%lai_bal_mass(n,k) * & - red_state%crwn_area_mass(n,k) - veg_state%canht(l,n) = veg_state%canht(l,n) + & - red_state%plantNumDensity(l,n,k) * red_state%ht_mass(n,k) * & - red_state%crwn_area_mass(n,k) - END DO - - ! Convert to per m2 plant - IF (veg_state%frac(l,n) > 0.0) THEN - veg_state%vegCpft(l,n) = veg_state%vegCpft(l,n) / veg_state%frac(l,n) - veg_state%lai_bal(l,n) = veg_state%lai_bal(l,n) / veg_state%frac(l,n) - veg_state%canht(l,n) = veg_state%canht(l,n) / veg_state%frac(l,n) - ELSE - - ! If frac is zero, set mean to lowest mass class value - veg_state%vegCpft(l,n) = red_state%mass_mass(n,1) - veg_state%lai_bal(l,n) = red_state%lai_bal_mass(n,1) - veg_state%canht(l,n) = red_state%ht_mass(n,1) - END IF - - veg_state%lai(l,n) = veg_state%phen(l,n) * veg_state%lai_bal(l,n) + ! Estimate the PFT mean physical properties by aggregating across the + ! mass class structure. + CALL pft_mean_from_mass_class( & + !IN sizing + red_parms%mclass(n), & + !IN mass-cohort properties + red_state%plantNumDensity(l,n,1:red_parms%mclass(n)), & + red_state%mass_mass(n,1:red_parms%mclass(n)), & + red_state%lai_bal_mass(n,1:red_parms%mclass(n)), & + red_state%ht_mass(n,1:red_parms%mclass(n)), & + red_state%crwn_area_mass(n,1:red_parms%mclass(n)), & + !OUT fields + veg_state%frac(l,n),veg_state%vegCpft(l,n),veg_state%lai_bal(l,n), & + veg_state%canht(l,n) & + ) + + ! Update the phenology + veg_state%lai(l,n) = veg_state%phen(l,n) * veg_state%lai_bal(l,n) + + ! Update the leaf the root pools veg_state%leafC(l,n) = cmass * lma(n) * veg_state%lai_bal(l,n) veg_state%rootC(l,n) = cmass * lma(n) * veg_state%lai_bal(l,n) ! Wood carbon balance of total minus leaf and root veg_state%woodC(l,n) = veg_state%vegCpft(l,n) - veg_state%leafC(l,n) & - veg_state%rootC(l,n) - ! Mean fraction over the update, approximating the changing canopy area - ! as vegetation dynamics were applied. - frac_mid(l,n) = 0.5 * (frac_old(l,n) + veg_state%frac(l,n)) - ! NPP and its derived litter fluxes are normalised per unit PFT canopy - ! area, but that area has since changed from frac_old. Rescale them - ! onto the frac_mid basis to limit the error this introduces. - CALL veg3_implicit_flux_frac_adj(frac_old(l,n), veg_state%frac(l,n), & - veg_state%npp_dr_out(l,n), & - veg_state%leaf_litC(l,n), & - veg_state%root_litC(l,n), & - veg_state%wood_litC(l,n), & - veg_state%mort_litC(l,n)) - - ! Aggregate the leaf/root/wood turnover litter and the mortality/ + ! Aggregate the leaf/root/wood turnover litter and the mortality litter ! demographic litter for the total litter flux per PFT fraction. - veg_state%litCpft(l,n) = (veg_state%leaf_litC(l,n) + & - veg_state%root_litC(l,n) + & - veg_state%wood_litC(l,n) + & - veg_state%mort_litC(l,n) * rsec_per_day * & + veg_state%litCpft(l,n) = (veg_state%leaf_litC(l,n) + & + veg_state%root_litC(l,n) + & + veg_state%wood_litC(l,n) + & + veg_state%mort_litC(l,n) * rsec_per_day * & 360.0) ! Update bare soil veg_state%frac(l,soil) = MAX(0.0, 1.0 - SUM(veg_state%frac(l,1:nnpft))) @@ -535,56 +526,95 @@ SUBROUTINE red_veg3_couple(ainfo) RETURN END SUBROUTINE red_veg3_couple -!----------------------------------------------------------------------------- -SUBROUTINE veg3_implicit_flux_frac_adj(frac_old, frac_new, npp, leaf_litC, & - root_litC, wood_litC, mort_litC) !----------------------------------------------------------------------------- -! npp/leaf_litC/root_litC/wood_litC/mort_litC are normalised per unit PFT -! canopy area, i.e. on frac_old. If the PFT fraction has since moved to -! frac_new, that normalisation is stale, so rescale each flux onto the mean -! of frac_old and frac_new - an implicit approximation to the fraction -! evolving smoothly over the timestep - to limit the resulting error. +SUBROUTINE pft_mean_from_mass_class( & + !IN sizing + mclass, & + !IN mass-cohort properties + plantNumDensity,mass_mass,lai_bal_mass,ht_mass,crwn_area_mass, & + !OUT fields + frac,vegCpft,lai_bal,canht & + ) +!----------------------------------------------------------------------------- +! Aggregates the plant number density across the mass class structure of a +! single PFT at a single point into the PFT mean physical properties +! (fraction, carbon density, balanced LAI and canopy height). !----------------------------------------------------------------------------- IMPLICIT NONE -REAL, INTENT(IN) :: frac_old - ! PFT fraction that npp/leaf_litC/root_litC/wood_litC/mort_litC are - ! normalised on, prior to this adjustment. -REAL, INTENT(IN) :: frac_new - ! Updated PFT fraction, following vegetation dynamics. - -REAL, INTENT(IN OUT) :: npp -REAL, INTENT(IN OUT) :: leaf_litC -REAL, INTENT(IN OUT) :: root_litC -REAL, INTENT(IN OUT) :: wood_litC -REAL, INTENT(IN OUT) :: mort_litC - ! Fluxes normalised per unit PFT canopy area (on the frac_old basis on - ! entry), rescaled onto the frac_mid basis on exit. - -REAL :: frac_mid - ! Mean of frac_old and frac_new; the basis the fluxes are rescaled onto. -REAL :: flux_scale - ! Ratio used to rescale each flux from the frac_old basis onto the - ! frac_mid basis. +!----------------------------------------------------------------------------- +! Integers with INTENT IN +!----------------------------------------------------------------------------- +INTEGER, INTENT(IN) :: mclass + ! Number of mass classes for this PFT. -!End of header +!----------------------------------------------------------------------------- +! Reals with INTENT IN +!----------------------------------------------------------------------------- +REAL, INTENT(IN) :: & +plantNumDensity(mclass), & + ! Population density within each mass cohort. (m-2) +mass_mass(mclass), & + ! Mass of an individual member of each mass cohort. (kgC) +lai_bal_mass(mclass), & + ! Balanced leaf area index of an individual member of each mass + ! cohort. (m2 m-2) +ht_mass(mclass), & + ! Height of an individual member of each mass cohort. (m) +crwn_area_mass(mclass) + ! Crown area of an individual member of each mass cohort. (m2) -frac_mid = 0.5 * (frac_old + frac_new) +!----------------------------------------------------------------------------- +! Reals with INTENT OUT +!----------------------------------------------------------------------------- +REAL, INTENT(OUT) :: & +frac, & + ! PFT fraction across the gridbox. (-) +vegCpft, & + ! Total PFT carbon density per PFT area fraction. (kg C m-2) +lai_bal, & + ! Balanced LAI per PFT area fraction. (m2 m-2) +canht + ! Canopy height per PFT area fraction. (m) -IF (frac_mid > 0.0) THEN - flux_scale = frac_old / frac_mid +!----------------------------------------------------------------------------- +!Local Vars +!----------------------------------------------------------------------------- +INTEGER :: k + +!End of headers + +! Initialise vars +frac = 0.0 +vegCpft = 0.0 +lai_bal = 0.0 +canht = 0.0 + +! Convert from plant number on mass classes to area on PFTs +! and biomass to carbon +DO k = 1, mclass + frac = frac + plantNumDensity(k) * crwn_area_mass(k) + vegCpft = vegCpft + plantNumDensity(k) * mass_mass(k) + lai_bal = lai_bal + plantNumDensity(k) * lai_bal_mass(k) * crwn_area_mass(k) + canht = canht + plantNumDensity(k) * ht_mass(k) * crwn_area_mass(k) +END DO + +! Convert to per m2 plant +IF (frac > 0.0) THEN + vegCpft = vegCpft / frac + lai_bal = lai_bal / frac + canht = canht / frac ELSE - flux_scale = 0.0 -END IF -npp = npp * flux_scale -leaf_litC = leaf_litC * flux_scale -root_litC = root_litC * flux_scale -wood_litC = wood_litC * flux_scale -mort_litC = mort_litC * flux_scale + ! If frac is zero, set mean to lowest mass class value + vegCpft = mass_mass(1) + lai_bal = lai_bal_mass(1) + canht = ht_mass(1) +END IF -END SUBROUTINE veg3_implicit_flux_frac_adj +END SUBROUTINE pft_mean_from_mass_class !----------------------------------------------------------------------------- + END MODULE veg3_field_mod diff --git a/src/science/vegetation/veg3_litter_mod.F90 b/src/science/vegetation/veg3_litter_mod.F90 index cd669bb3..bb7b285f 100644 --- a/src/science/vegetation/veg3_litter_mod.F90 +++ b/src/science/vegetation/veg3_litter_mod.F90 @@ -112,9 +112,9 @@ SUBROUTINE veg3_Litter( & ! Copy across to veg_state%leaf_litC, veg_state%root_litC, veg_state%wood_litC ! Convert to kg C m-2 (360d)-1 - veg_state%leaf_litC(l,n) = leaf_litter(l,n) * rsec_per_day * 360. - veg_state%root_litC(l,n) = root_litter(l,n) * rsec_per_day * 360. - veg_state%wood_litC(l,n) = wood_litter(l,n) * rsec_per_day * 360. + veg_state%leaf_litC(l,n) = leaf_litter(l,n) * rsec_per_day * 360.0 + veg_state%root_litC(l,n) = root_litter(l,n) * rsec_per_day * 360.0 + veg_state%wood_litC(l,n) = wood_litter(l,n) * rsec_per_day * 360.0 local_litter(l,n) = leaf_litter(l,n) + root_litter(l,n) + wood_litter(l,n) diff --git a/src/science/vegetation/veg3_red_dynamic_mod.F90 b/src/science/vegetation/veg3_red_dynamic_mod.F90 index f355e2f0..3dbd25d6 100644 --- a/src/science/vegetation/veg3_red_dynamic_mod.F90 +++ b/src/science/vegetation/veg3_red_dynamic_mod.F90 @@ -62,9 +62,9 @@ SUBROUTINE veg3_red_dynamic( & !---------------------------------------------------------------------------- REAL, INTENT(IN) :: & growth(land_pts,nnpft), & - ! The total carbon assimilate across the PFT area. (kgC/m2/s) + ! The total carbon assimilate across the PFT area. (kgC m-2 s-1) mort_add(land_pts,nnpft,nmasst), & - ! Additional plant mortality across plant mass (/s) + ! Additional plant mortality across plant mass (s-1) dt ! Dynamic vegetation time-step (s) @@ -75,18 +75,13 @@ SUBROUTINE veg3_red_dynamic( & REAL :: & P_s(land_pts,nnpft), & - ! Total gridbox carbon assimilate devoted to recruitment. - ! (kgC/m2/s) + ! Total gridbox carbon assimilate devoted to recruitment. + ! (kgC m-2 s-1) g0(land_pts,nnpft), & ! Boundary growth for an individual member of the smallest mass - ! cohort. (kgC/year) -neg_growth(land_pts,nnpft), & - ! Uniform PFT-carbon-proportional loss rate for when growth is - ! negative (e.g., local litterfall exceeds NPP). Drives a - ! downward mass-class shrinkage flux, and is only applied as a - ! mortality rate for the lowest mass class. (/s) + ! cohort. (kgC s-1) frac_shade(land_pts,nnpft) - ! Competitive shading of seedlings in each PFT. + ! Competitive shading of seedlings in each PFT. (-) !End of headers @@ -95,7 +90,6 @@ SUBROUTINE veg3_red_dynamic( & veg_state%mort_litC(:,:) = 0.0 P_s(:,:) = 0.0 g0(:,:) = 0.0 -neg_growth(:,:) = 0.0 frac_shade(:,:) = 0.0 ! Dynamic demographic loop to update the number density of each PFT across the @@ -109,49 +103,49 @@ SUBROUTINE veg3_red_dynamic( & !IN PFT parameters red_parms%alpha_recrt(n), & !IN fields - veg_state%frac(l,n),growth(l,n),veg_state%vegCpft(l,n), & + veg_state%frac(l,n),growth(l,n), & !IN mass-cohort properties red_state%plantNumDensity(l,n,1:red_parms%mclass(n)), & red_state%g_mass_scale(n,1:red_parms%mclass(n)), & !OUT fields - P_s(l,n),g0(l,n),neg_growth(l,n) & + P_s(l,n),g0(l,n) & ) !Estimate the inter-PFT competition. DO j = 1, nnpft - frac_shade(l,n) = MIN(1.0,frac_shade(l,n) + red_parms%comp_coef(n,j) & + frac_shade(l,n) = MIN(1.0,frac_shade(l,n) + red_parms%comp_coef(n,j) & * veg_state%frac(l,j)) END DO ! Call to update the PFT number density. - CALL update_pft_size_structure( & + CALL update_pft_size_structure( & !IN sizing - red_parms%mclass(n), & + red_parms%mclass(n), & !IN Control vars - dt, & + dt, & !IN PFT parameters - red_parms%mort_base(n),red_parms%frac_min(n), & + red_parms%mort_base(n),red_parms%frac_min(n), & !IN fields - mort_add(l,n,1:red_parms%mclass(n)),growth(l,n),neg_growth(l,n), & - P_s(l,n),g0(l,n),frac_shade(l,n), & + mort_add(l,n,1:red_parms%mclass(n)),growth(l,n), & + P_s(l,n),g0(l,n),frac_shade(l,n), & !IN mass-cohort properties - red_state%g_mass_scale(n,1:red_parms%mclass(n)), & - red_state%mass_mass(n,1:red_parms%mclass(n)), & - red_state%crwn_area_mass(n,1:red_parms%mclass(n)), & + red_state%g_mass_scale(n,1:red_parms%mclass(n)), & + red_state%mass_mass(n,1:red_parms%mclass(n)), & + red_state%crwn_area_mass(n,1:red_parms%mclass(n)), & !INOUT state - red_state%plantNumDensity(l,n,1:red_parms%mclass(n)), & - red_state%mort(l,n,1:red_parms%mclass(n)), & + red_state%plantNumDensity(l,n,1:red_parms%mclass(n)), & + red_state%mort(l,n,1:red_parms%mclass(n)), & !OUT diagnostics - veg_state%mort_litC(l,n) & + veg_state%mort_litC(l,n) & ) ! Divide mort_litC by the PFT fraction (not restimated yet) - IF (veg_state%frac(l,n) > 0.0) THEN - veg_state%mort_litC(l,n) = veg_state%mort_litC(l,n) / & - veg_state%frac(l,n) - ELSE - veg_state%mort_litC(l,n) = 0.0 - END IF + IF (veg_state%frac(l,n) > 0.0) THEN + veg_state%mort_litC(l,n) = veg_state%mort_litC(l,n) / & + veg_state%frac(l,n) + ELSE + veg_state%mort_litC(l,n) = 0.0 + END IF END DO END DO @@ -166,11 +160,11 @@ SUBROUTINE growth_onto_mass_class( & !IN PFT parameters alpha_recrt, & !IN fields - frac,growth,vegCpft, & + frac,growth, & !IN mass-cohort properties plantNumDensity,g_mass_scale, & !OUT fields - P_s,g0,neg_growth & + P_s, g0 & ) IMPLICIT NONE @@ -190,11 +184,9 @@ SUBROUTINE growth_onto_mass_class( & frac, & ! PFT fraction across the gridbox. (-) growth, & - ! The total carbon assimilate across the PFT area. (kgC/m2/s) -vegCpft, & - ! The total PFT vegetation carbon. (kgC/m2) + ! The total carbon assimilate across the PFT area. (kgC m-2 s-1) plantNumDensity(mclass), & - ! Population density within each mass cohort. (/m2) + ! Population density within each mass cohort. (m-2) g_mass_scale(mclass) ! Allometric scaling of growth across the mass cohorts. @@ -203,15 +195,10 @@ SUBROUTINE growth_onto_mass_class( & !----------------------------------------------------------------------------- REAL, INTENT(OUT) :: & P_s, & - ! Gridbox carbon assimilate devoted to recruitment. (kgC/m2/s) -g0, & + ! Gridbox carbon assimilate devoted to recruitment. (kgC m-2 s-1) +g0 ! Boundary growth for an individual member of the smallest mass - ! cohort. (kgC/year) -neg_growth - ! Uniform PFT-carbon-proportional loss rate for when growth is - ! negative (e.g., local litterfall exceeds NPP). Drives a - ! downward mass-class shrinkage flux, and is only applied as a - ! mortality rate for the lowest mass class. (/s) + ! cohort. (kgC s-1) !----------------------------------------------------------------------------- !Local Vars @@ -220,20 +207,20 @@ SUBROUTINE growth_onto_mass_class( & REAL :: & p, & - ! The total PFT carbon assimilate across the gridbox. (kgC/m2/s) + ! The total PFT carbon assimilate across the gridbox. + ! (kgC m-2 s-1) g, & ! Total gridbox carbon assimilate devoted to vegetation - ! structural growth. (kgC/m2/s) + ! structural growth. (kgC m-2 s-1) plantNumDensity_g_sum ! Summation of the relative cohort contribution towards the - ! total PFT assimilate (/m2) + ! total PFT assimilate (m-2) !End of headers ! Initialise vars P_s = 0.0 g0 = 0.0 -neg_growth = 0.0 plantNumDensity_g_sum = 0.0 ! Sum product of the number density and the allometric scaling @@ -243,7 +230,7 @@ SUBROUTINE growth_onto_mass_class( & END DO ! Partition the growth into recruitment and structural growth -p = frac * growth +p= frac * growth P_s = alpha_recrt * p g = (1.0 - alpha_recrt) * p @@ -252,13 +239,14 @@ SUBROUTINE growth_onto_mass_class( & END IF IF (growth < 0.0) THEN + ! Recruitment does not apply when growth is negative. The associated + ! shrinkage is instead applied via g0/g_mass_scale metabolic (phi_g) + ! scaling in update_pft_size_structure. P_s = 0.0 - ! Estimate the uniform PFT-carbon-proportional loss rate required to - ! account for the negative growth (local litterfall exceeding NPP). - IF (vegCpft > 0) THEN - neg_growth = - growth / vegCpft - + ! Adjust g0 to account for the loss of recruitment + IF (plantNumDensity_g_sum > 0) THEN + g0 = p / plantNumDensity_g_sum END IF END IF @@ -275,7 +263,7 @@ SUBROUTINE update_pft_size_structure( & !IN PFT parameters mort_base,frac_min, & !IN fields - mort_add,growth,neg_growth,P_s,g0,frac_shade, & + mort_add,growth,P_s,g0,frac_shade, & !IN mass-cohort properties g_mass_scale,mass_mass,crwn_area_mass, & !INOUT state @@ -295,32 +283,27 @@ SUBROUTINE update_pft_size_structure( & !----------------------------------------------------------------------------- ! Reals with INTENT IN !----------------------------------------------------------------------------- -REAL, INTENT(IN) :: & -dt, & +REAL, INTENT(IN) :: & +dt, & ! Dynamic vegetation time-step (s) -mort_base, & - ! Background mortality rate for this PFT. (/s) -frac_min, & +mort_base, & + ! Background mortality rate for this PFT. (s-1) +frac_min, & ! Minimum vegetation fraction for this PFT. (-) -mort_add(mclass), & - ! Additional plant mortality across plant mass (/s) -growth, & - ! The total carbon assimilate across the PFT area. (kgC/m2/s) -neg_growth, & - ! Uniform PFT-carbon-proportional loss rate for when growth is - ! negative (e.g., local litterfall exceeds NPP). Drives a - ! downward mass-class shrinkage flux, and is only applied as a - ! mortality rate for the lowest mass class. (/s) -P_s, & - ! Gridbox carbon assimilate devoted to recruitment. (kgC/m2/s) -g0, & +mort_add(mclass), & + ! Additional plant mortality across plant mass (s-1) +growth, & + ! The total carbon assimilate across the PFT area. (kgC m-2 s-1) +P_s, & + ! Gridbox carbon assimilate devoted to recruitment. (kgC m-2 s-1) +g0, & ! Boundary growth for an individual member of the smallest mass - ! cohort. (kgC/year) -frac_shade, & + ! cohort. (kgC s-1) +frac_shade, & ! Competitive shading of seedlings in this PFT. -g_mass_scale(mclass), & +g_mass_scale(mclass), & ! Allometric scaling of growth across the mass cohorts. -mass_mass(mclass), & +mass_mass(mclass), & ! Mass of an individual member of each mass cohort. (kgC) crwn_area_mass(mclass) ! Crown area of an individual member of each mass cohort. (m2) @@ -328,18 +311,18 @@ SUBROUTINE update_pft_size_structure( & !----------------------------------------------------------------------------- ! Reals with INTENT INOUT !----------------------------------------------------------------------------- -REAL, INTENT(INOUT) :: & +REAL, INTENT(IN OUT) :: & plantNumDensity(mclass), & - ! Population density within each mass cohort. (/m2) + ! Population density within each mass cohort. (m-2) mort(mclass) - ! Mortality rate within each mass cohort. (/s) + ! Mortality rate within each mass cohort. (s-1) !----------------------------------------------------------------------------- ! Reals with INTENT OUT !----------------------------------------------------------------------------- REAL, INTENT(OUT) :: mort_litC ! Mortality/demographic litter for this PFT, normalised per - ! unit PFT canopy area. (kgC/m2/s) + ! unit PFT canopy area. (kgC m-2 s-1) !----------------------------------------------------------------------------- !Local Vars @@ -348,18 +331,19 @@ SUBROUTINE update_pft_size_structure( & REAL :: & g_mass(mclass), & - ! Individual growth across the mass cohorts. (kgC/s) + ! Individual growth across the mass cohorts. (kgC s-1) dplantNumDensity_dt(mclass), & ! Net rate of change of population density within each mass - ! cohort. (m2/s) + ! cohort. (m-2 s-1) flux_in(mclass), & - ! Rate of change of population growing into a mass cohort. (m2/s) + ! Rate of change of population growing into a mass cohort. + ! (m-2 s-1) flux_out(mclass), & - ! Rate of change of population growing out of a mass cohort. - ! (m2/s) + ! Rate of change of population growing out of a mass cohort. + ! (m-2 s-1) frac_check - ! The difference between the minimum vegetation fraction and the - ! updated fraction. (-) + ! The difference between the minimum vegetation fraction and the + ! updated fraction. (-) !End of headers @@ -373,22 +357,26 @@ SUBROUTINE update_pft_size_structure( & DO k = 1, mclass + ! Metabolic (phi_g) scaling of growth across the mass classes. When growth + ! is negative, this represents the equivalent scaling of litterfall losses: + ! litter is dominated by leaf and fine root turnover, which scale with + ! metabolic rate (mass**phi_g) rather than with mass itself. + g_mass(k) = g0 * g_mass_scale(k) + IF (growth < 0.0) THEN ! Negative growth is represented as a downward shrinkage flux of - ! individuals through the mass classes, the lowest mass class has + ! individuals through the mass classes, the lowest mass class has ! nowhere lower to shrink into, so its share of the loss is instead matched ! by an equivalent mortality rate. - g_mass(k) = 0.0 + mort(k) = mort_base + mort_add(k) IF (k == 1) THEN - mort(k) = mort_base + mort_add(k) + neg_growth - ELSE - mort(k) = mort_base + mort_add(k) + mort(k) = mort(k) - g_mass(k) / mass_mass(k) END IF IF (k < mclass) THEN ! Flux shrinking down into this class from the class above - flux_in(k) = plantNumDensity(k+1) * neg_growth * mass_mass(k+1) & + flux_in(k) = - plantNumDensity(k+1) * g_mass(k+1) & / (mass_mass(k+1) - mass_mass(k)) ELSE flux_in(k) = 0.0 @@ -396,14 +384,13 @@ SUBROUTINE update_pft_size_structure( & IF (k > 1) THEN ! Flux shrinking out of this class into the class below - flux_out(k) = plantNumDensity(k) * neg_growth * mass_mass(k) & + flux_out(k) = - plantNumDensity(k) * g_mass(k) & / (mass_mass(k) - mass_mass(k-1)) ELSE flux_out(k) = 0.0 END IF ELSE - g_mass(k) = g0 * g_mass_scale(k) mort(k) = mort_base + mort_add(k) IF (k == 1) THEN @@ -438,7 +425,7 @@ SUBROUTINE update_pft_size_structure( & IF (plantNumDensity(k) & + (dplantNumDensity_dt(k) * dt) < 0.0 ) THEN dplantNumDensity_dt(k) = -plantNumDensity(k) / dt - mort_litC = mort_litC & + mort_litC = mort_litC & + (dplantNumDensity_dt(k) - flux_out(k)) & * mass_mass(k) From b338a2f15cc0bda6293b6a2785776c905a314270 Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Fri, 4 Sep 2026 10:25:09 +0100 Subject: [PATCH 03/14] Adding my name to contributors. --- CONTRIBUTORS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index db1cc110..f281d6d3 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -5,3 +5,4 @@ | yaswant | Yaswant Pradhan | Met Office | 2026-06-30 | | doucla | Douglas Clark | UKCEH | 2026-07-03 | | james-bruten-mo | James Bruten | Met Office | 2026-07-09 | +| mo-athurargles | Arthur Argles | Met Office | 2026-09-04 | \ No newline at end of file From 916cd6eaa1111b04956b9f5007e8948f38fc54ab Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Fri, 4 Sep 2026 10:35:45 +0100 Subject: [PATCH 04/14] Fxing typo in GitHub username. --- CONTRIBUTORS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index f281d6d3..194adebc 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -5,4 +5,4 @@ | yaswant | Yaswant Pradhan | Met Office | 2026-06-30 | | doucla | Douglas Clark | UKCEH | 2026-07-03 | | james-bruten-mo | James Bruten | Met Office | 2026-07-09 | -| mo-athurargles | Arthur Argles | Met Office | 2026-09-04 | \ No newline at end of file +| mo-arthurargles | Arthur Argles | Met Office | 2026-09-04 | \ No newline at end of file From 58c3039ddf9df259cd2d274cdc90f0fcdfdf166a Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Sun, 6 Sep 2026 14:08:49 +0100 Subject: [PATCH 05/14] Changes address comments from PR GitHub co-pilot for issue #90: 1. In phenol_jls.F90 set in variable lai_bal_red to be be REAL(KIND=real_jlslsm) to ensure consistency between other REAL inputs. 2. From veg3_field_mod.F90, remove frac_old from red_veg3_couple_mod.F90 not necessary for this change and not used in the script. --- src/control/shared/veg3_field_mod.F90 | 7 ------- src/science/vegetation/phenol_jls.F90 | 2 +- 2 files changed, 1 insertion(+), 8 deletions(-) diff --git a/src/control/shared/veg3_field_mod.F90 b/src/control/shared/veg3_field_mod.F90 index 9745320f..c242b6da 100644 --- a/src/control/shared/veg3_field_mod.F90 +++ b/src/control/shared/veg3_field_mod.F90 @@ -454,16 +454,9 @@ SUBROUTINE red_veg3_couple(ainfo) !Local INTEGER :: l,n ! Index variables. -REAL :: frac_old(land_pts,nnpft) - ! PFT fraction prior to updating veg_state%frac below, i.e. before the - ! change in canopy area resulting from vegetation dynamics. - !----------------------------------------------------------------------------- !end of header -! Record the vegetation fraction prior to updating it below, for use in -! weighting litter fluxes that occurred over the previous fraction. -frac_old(:,:) = veg_state%frac(:,1:nnpft) veg_state%vegCpft(:,:) = 0.0 veg_state%lai_bal(:,:) = 0.0 diff --git a/src/science/vegetation/phenol_jls.F90 b/src/science/vegetation/phenol_jls.F90 index b83fe003..9f821f6d 100644 --- a/src/science/vegetation/phenol_jls.F90 +++ b/src/science/vegetation/phenol_jls.F90 @@ -63,7 +63,7 @@ SUBROUTINE phenol (land_pts, veg_pts, n, veg_index, dtime_phen, g_leaf, ht, & !----------------------------------------------------------------------------- ! Optional Arguments with INTENT(IN) !----------------------------------------------------------------------------- -REAL, INTENT(IN), OPTIONAL :: & +REAL(KIND=real_jlslsm), INTENT(IN), OPTIONAL :: & lai_bal_red(land_pts) ! Balanced growth LAI from veg3. From a0fcfe594b24915a163c4dcfa053e57682c692ef Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Sun, 6 Sep 2026 15:36:17 +0100 Subject: [PATCH 06/14] Uploading changes to hopefully fix the coupling with UM in JULES shared routines related to veg3_field_mod.F90. --- src/control/shared/veg3_field_mod.F90 | 12 ++++++++++-- src/initialisation/standalone/init_mod.F90 | 2 +- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/control/shared/veg3_field_mod.F90 b/src/control/shared/veg3_field_mod.F90 index c242b6da..caa8b9da 100644 --- a/src/control/shared/veg3_field_mod.F90 +++ b/src/control/shared/veg3_field_mod.F90 @@ -211,7 +211,7 @@ SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) END SUBROUTINE veg3_field_allocate !------------------------------------------------------------------------------- -SUBROUTINE veg3_field_assoc(progs, ainfo) +SUBROUTINE veg3_field_assoc(progs, ainfo, trifctl_data, trif_vars_data) ! Initial code to associate the veg3 and red fields to the rest of JULES ! This new routine moves out the pointers from veg3_set_fields to here @@ -220,13 +220,21 @@ SUBROUTINE veg3_field_assoc(progs, ainfo) ! We still need to run with l_triffid, but this will need to be addressed in a ! future revision to allow for veg3 to run fully independently of the switch. +! Note: trifctl_data and trif_vars_data are passed in as arguments (rather than +! USE-associated from jules_fields_mod) because this module lives in +! src/control/shared and must also build for the UM, where jules_fields_mod +! (a standalone-only module) is not available. + USE jules_vegetation_mod, ONLY: l_triffid -USE jules_fields_mod, ONLY: trifctl_data, trif_vars_data +USE trifctl, ONLY: trifctl_data_type +USE trif_vars_mod, ONLY: trif_vars_data_type IMPLICIT NONE TYPE(progs_type), INTENT(IN) :: progs TYPE(ainfo_type), INTENT(IN) :: ainfo +TYPE(trifctl_data_type), INTENT(IN), TARGET :: trifctl_data +TYPE(trif_vars_data_type), INTENT(IN), TARGET :: trif_vars_data ! End of header !------------------------------------------------------------------------------- IF (l_red .AND. l_triffid) THEN diff --git a/src/initialisation/standalone/init_mod.F90 b/src/initialisation/standalone/init_mod.F90 index abd8f711..d4dd9879 100644 --- a/src/initialisation/standalone/init_mod.F90 +++ b/src/initialisation/standalone/init_mod.F90 @@ -351,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(progs,ainfo) +CALL veg3_field_assoc(progs,ainfo,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) From 52258c02e2756799a4eb1b54f961624a96b65c68 Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Mon, 7 Sep 2026 10:31:23 +0100 Subject: [PATCH 07/14] Fixing a few issues after Copilot review in MetOffice/jules/pull/160 (#90). Added veg3_field_deallocate during veg3_field_assoc, such that we avoid potential memory issues with pointer when associating with prognostic & diagnositic variables. Fixed typos and inconsistencies with mort_litC description. --- src/control/shared/veg3_field_mod.F90 | 54 ++++++++++++++++++- .../vegetation/veg3_red_dynamic_mod.F90 | 5 +- 2 files changed, 55 insertions(+), 4 deletions(-) diff --git a/src/control/shared/veg3_field_mod.F90 b/src/control/shared/veg3_field_mod.F90 index caa8b9da..ba89862b 100644 --- a/src/control/shared/veg3_field_mod.F90 +++ b/src/control/shared/veg3_field_mod.F90 @@ -114,8 +114,8 @@ MODULE veg3_field_mod PRIVATE !Expose routines -PUBLIC :: veg3_field_init, veg3_field_allocate, veg3_field_assoc, & - red_veg3_couple +PUBLIC :: veg3_field_init, veg3_field_allocate, veg3_field_deallocate, & + veg3_field_assoc, red_veg3_couple !Expose data PUBLIC :: veg_state, red_state @@ -210,6 +210,51 @@ SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) RETURN END SUBROUTINE veg3_field_allocate +!------------------------------------------------------------------------------- + +SUBROUTINE veg3_field_deallocate() + +! Deallocates the veg_state and red_state POINTER components that were +! allocated locally by veg3_field_allocate but which veg3_field_assoc later +! re-associates onto external targets (e.g. progs, trif_vars_data and +! trifctl_data fields). This must be called before that re-association is +! done, otherwise the original local allocations become orphaned (leaked) - +! this can be a significant leak for red_state%plantNumDensity in particular. +! +! Note: veg_state%phen and veg_state%npp_acc are deliberately NOT included +! here. They are never re-associated in veg3_field_assoc and remain the +! arrays allocated in veg3_field_allocate for the lifetime of the run. + +IMPLICIT NONE + +!End of header + +DEALLOCATE(veg_state%leafC) +DEALLOCATE(veg_state%rootC) +DEALLOCATE(veg_state%woodC) +DEALLOCATE(veg_state%vegCpft) +DEALLOCATE(veg_state%lai_bal) +DEALLOCATE(veg_state%canht) +DEALLOCATE(veg_state%lai) +DEALLOCATE(veg_state%npp_dr_out) +DEALLOCATE(veg_state%frac) +DEALLOCATE(veg_state%g_leaf_phen) +DEALLOCATE(veg_state%g_leaf) +DEALLOCATE(veg_state%lai_phen) +DEALLOCATE(veg_state%g_leaf_acc) +DEALLOCATE(veg_state%g_leaf_phen_acc) +DEALLOCATE(veg_state%leaf_litC) +DEALLOCATE(veg_state%root_litC) +DEALLOCATE(veg_state%wood_litC) +DEALLOCATE(veg_state%litCpft) +DEALLOCATE(veg_state%litC) +DEALLOCATE(veg_state%vegC) + +DEALLOCATE(red_state%plantNumDensity) + +RETURN +END SUBROUTINE veg3_field_deallocate + !------------------------------------------------------------------------------- SUBROUTINE veg3_field_assoc(progs, ainfo, trifctl_data, trif_vars_data) @@ -239,6 +284,11 @@ SUBROUTINE veg3_field_assoc(progs, ainfo, trifctl_data, trif_vars_data) !------------------------------------------------------------------------------- IF (l_red .AND. l_triffid) THEN + ! Deallocate the local pointer targets set up in veg3_field_allocate before + ! they are re-associated onto external targets below, to avoid leaking the + ! original allocations (see veg3_field_deallocate for details). + CALL veg3_field_deallocate() + ! Set pointers to the prognostic fields ! Note: veg_state%phen is not associated to an external target here - it is ! not an independent prognostic. It is diagnosed each phenology call as diff --git a/src/science/vegetation/veg3_red_dynamic_mod.F90 b/src/science/vegetation/veg3_red_dynamic_mod.F90 index 3dbd25d6..3a9df36b 100644 --- a/src/science/vegetation/veg3_red_dynamic_mod.F90 +++ b/src/science/vegetation/veg3_red_dynamic_mod.F90 @@ -139,7 +139,7 @@ SUBROUTINE veg3_red_dynamic( & veg_state%mort_litC(l,n) & ) - ! Divide mort_litC by the PFT fraction (not restimated yet) + ! Divide mort_litC by the PFT fraction (not re-estimated yet) IF (veg_state%frac(l,n) > 0.0) THEN veg_state%mort_litC(l,n) = veg_state%mort_litC(l,n) / & veg_state%frac(l,n) @@ -322,7 +322,8 @@ SUBROUTINE update_pft_size_structure( & !----------------------------------------------------------------------------- REAL, INTENT(OUT) :: mort_litC ! Mortality/demographic litter for this PFT, normalised per - ! unit PFT canopy area. (kgC m-2 s-1) + ! unit gridbox area, veg3_red_dynamics renormalises to + ! PFT canopy area. (kgC m-2 s-1) !----------------------------------------------------------------------------- !Local Vars From 44006feee334f556deafaf2982486cf6a4d1ab2a Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Mon, 7 Sep 2026 10:38:29 +0100 Subject: [PATCH 08/14] Using umdp3_fixer.py to fix veg3_red_dynamic_mod.F90. --- src/science/vegetation/veg3_red_dynamic_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/science/vegetation/veg3_red_dynamic_mod.F90 b/src/science/vegetation/veg3_red_dynamic_mod.F90 index 3a9df36b..6896138d 100644 --- a/src/science/vegetation/veg3_red_dynamic_mod.F90 +++ b/src/science/vegetation/veg3_red_dynamic_mod.F90 @@ -322,7 +322,7 @@ SUBROUTINE update_pft_size_structure( & !----------------------------------------------------------------------------- REAL, INTENT(OUT) :: mort_litC ! Mortality/demographic litter for this PFT, normalised per - ! unit gridbox area, veg3_red_dynamics renormalises to + ! unit gridbox area, veg3_red_dynamics renormalises to ! PFT canopy area. (kgC m-2 s-1) !----------------------------------------------------------------------------- From b9a659af125ca6ef76c83bc48d8091000f76f864 Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Mon, 14 Sep 2026 15:12:25 +0100 Subject: [PATCH 09/14] Updating code so we could output more diagnostics related to leaf turnover (g_leaf, g_leaf_day, g_leaf_dr_out). Related issue #90. --- .../shared/next_gen_biogeochem_mod.F90 | 49 ++++++------------- src/control/shared/veg3_field_mod.F90 | 13 +++++ 2 files changed, 28 insertions(+), 34 deletions(-) diff --git a/src/control/shared/next_gen_biogeochem_mod.F90 b/src/control/shared/next_gen_biogeochem_mod.F90 index 3c6e75ae..91ec43c0 100644 --- a/src/control/shared/next_gen_biogeochem_mod.F90 +++ b/src/control/shared/next_gen_biogeochem_mod.F90 @@ -198,11 +198,8 @@ SUBROUTINE veg3_run_ctrl( & ! Litter production (kg C/m2/s). growth(land_pts,nnpft), & ! growth (kg C/m2/s). -mort_add(land_pts,nnpft,nmasst), & +mort_add(land_pts,nnpft,nmasst) ! mortality above baseline (/m2) -g_leaf_phen_dr(land_pts,nnpft) - ! Mean phenology-driven leaf turnover rate for driving vegetation - ! dynamics (/s) !End of headers @@ -218,7 +215,7 @@ SUBROUTINE veg3_run_ctrl( & ! veg-veg2a_jls_mod. !----------------------------------------------------------------------------- CALL veg3_phenol_couple(veg_index_pts,veg_index,veg3_ctrl,land_pts,nnpft, & - a_step,asteps_since_triffid,veg_state,g_leaf_phen_dr) + a_step,asteps_since_triffid,veg_state) ! Now call vegetation model IF (asteps_since_triffid == veg3_ctrl%nstep_trif) THEN @@ -230,7 +227,7 @@ SUBROUTINE veg3_run_ctrl( & !IN parms litter_parms, & !IN fields - g_leaf_phen_dr, & + veg_state%g_leaf_dr_out, & !IN state veg_state, & ! OUT Fields @@ -306,9 +303,7 @@ SUBROUTINE veg3_phenol_couple( & veg_index_pts,veg_index,veg3_ctrl,land_pts,nnpft, & a_step,asteps_since_triffid, & !IN state - veg_state, & - !OUT fields - g_leaf_phen_dr & + veg_state & ) ! Diagnoses leaf phenology and the mean phenology-driven leaf turnover rate @@ -348,36 +343,22 @@ SUBROUTINE veg3_phenol_couple( & !----------------------------------------------------------------------------- TYPE(veg_state_type),INTENT(IN OUT) :: veg_state -!----------------------------------------------------------------------------- -! Reals with INTENT OUT -!----------------------------------------------------------------------------- -REAL, INTENT(OUT) :: & -g_leaf_phen_dr(land_pts,nnpft) - ! Mean phenology-driven leaf turnover rate for driving vegetation - ! dynamics (s-1). - !----------------------------------------------------------------------------- !Local Vars !----------------------------------------------------------------------------- REAL :: & gam_trif, & ! Inverse vegetation dynamics coupling timestep ((360d)-1). -lai_bal_dummy(land_pts), & +lai_bal_dummy(land_pts) ! Dummy lai to pass into phenol routine, gets around bug where ! FORTRAN does not accept veg_state%lai_bal as optional ! argument. -g_leaf_day(land_pts) - ! Mean leaf turnover rate driving phenology, diagnosed from - ! the physiological leaf turnover accumulated since the - ! previous phenology call. ((360d)-1) INTEGER :: l,n,k ! Loop counters. !End of headers -g_leaf_phen_dr(:,:) = 0.0 - !----------------------------------------------------------------------------- ! Work out the phenology at its own timestep, appending to the accumulated ! leaf turnover rates. This is called independently of the vegetation @@ -388,7 +369,6 @@ SUBROUTINE veg3_phenol_couple( & veg_state%phen(:,:) = 1.0 lai_bal_dummy(:) = 0.0 - g_leaf_day(:) = 0.0 DO n = 1,nnpft lai_bal_dummy(:) = veg_state%lai_bal(:,n) @@ -397,11 +377,12 @@ SUBROUTINE veg3_phenol_couple( & ! elapsed phenology period, mirroring g_leaf_day in veg-veg2a_jls_mod. DO k = 1,veg_index_pts l = veg_index(k) - g_leaf_day(l) = veg_state%g_leaf_acc(l,n) / veg3_ctrl%dt_phen_360d + veg_state%g_leaf_day(l,n) = veg_state%g_leaf_acc(l,n) / & + veg3_ctrl%dt_phen_360d END DO CALL phenol(land_pts,veg_index_pts,n,veg_index,veg3_ctrl%dt_phen_360d, & - g_leaf_day,veg_state%canht(:,n), & + veg_state%g_leaf_day(:,n),veg_state%canht(:,n), & veg_state%lai(:,n),veg_state%g_leaf_phen(:,n),lai_bal_dummy) DO k = 1,veg_index_pts @@ -445,8 +426,8 @@ SUBROUTINE veg3_phenol_couple( & IF (l_phenol) THEN ! Diagnose the mean phenological leaf turnover rate over the ! coupling period and convert to JULES-standard per-second units. - g_leaf_phen_dr(l,n) = veg_state%g_leaf_phen_acc(l,n) * gam_trif / & - (rsec_per_day * 360.0) + veg_state%g_leaf_dr_out(l,n) = veg_state%g_leaf_phen_acc(l,n) * & + gam_trif / (rsec_per_day * 360.0) ! Reset the accumulated phenological turnover ready for the next ! coupling period. @@ -454,8 +435,8 @@ SUBROUTINE veg3_phenol_couple( & ELSE ! No phenology - fall back to the raw accumulated physiological ! leaf turnover rate, as in veg-veg2a_jls_mod. - g_leaf_phen_dr(l,n) = veg_state%g_leaf_acc(l,n) * gam_trif / & - (rsec_per_day * 360.0) + veg_state%g_leaf_dr_out(l,n) = veg_state%g_leaf_acc(l,n) * gam_trif / & + (rsec_per_day * 360.0) veg_state%g_leaf_acc(l,n) = 0.0 END IF @@ -464,12 +445,12 @@ SUBROUTINE veg3_phenol_couple( & ! the vegetation dynamics timestep. If it does, reduce the rate so ! that the turnover does not exceed the current LAI. IF (veg_state%lai(l,n) > 0.0) THEN - IF (g_leaf_phen_dr(l,n) * veg3_ctrl%dt_red > 1.0) THEN - g_leaf_phen_dr(l,n) = 1.0 / veg3_ctrl%dt_red + IF (veg_state%g_leaf_dr_out(l,n) * veg3_ctrl%dt_red > 1.0) THEN + veg_state%g_leaf_dr_out(l,n) = 1.0 / veg3_ctrl%dt_red END IF ELSE - g_leaf_phen_dr(l,n) = 0.0 + veg_state%g_leaf_dr_out(l,n) = 0.0 END IF diff --git a/src/control/shared/veg3_field_mod.F90 b/src/control/shared/veg3_field_mod.F90 index ba89862b..bc44a065 100644 --- a/src/control/shared/veg3_field_mod.F90 +++ b/src/control/shared/veg3_field_mod.F90 @@ -56,6 +56,11 @@ MODULE veg3_field_mod g_leaf(:,:), & ! Pointer necessary to get the diagnosed g_leaf rate from the ! rest of JULES. ((360d)-1) + g_leaf_day(:,:), & + ! PFT mean leaf turnover rate for input to phenology. ((360d)-1) + g_leaf_dr_out(:,:), & + ! PFT mean leaf turnover rate for driving vegetation dynamics. + ! ((360d)-1) lai_phen(:,:), & ! Diagnostic LAI immediately following the phenology update. ! Associated with trifctl%lai_phen_pft. (m2 m-2) @@ -150,6 +155,8 @@ SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) ALLOCATE(veg_state%npp_dr_out ( land_pts, nnpft) ) ALLOCATE(veg_state%g_leaf_phen ( land_pts, nnpft) ) ALLOCATE(veg_state%g_leaf ( land_pts, nnpft) ) +ALLOCATE(veg_state%g_leaf_day ( land_pts, nnpft) ) +ALLOCATE(veg_state%g_leaf_dr_out( land_pts, nnpft) ) ALLOCATE(veg_state%lai_phen ( land_pts, nnpft) ) ALLOCATE(veg_state%g_leaf_acc ( land_pts, nnpft) ) ALLOCATE(veg_state%g_leaf_phen_acc ( land_pts, nnpft) ) @@ -175,6 +182,8 @@ SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) veg_state%npp_dr_out(:,:) = 0.0 veg_state%g_leaf_phen(:,:) = 0.0 veg_state%g_leaf(:,:) = 0.0 +veg_state%g_leaf_day(:,:) = 0.0 +veg_state%g_leaf_dr_out(:,:) = 0.0 veg_state%lai_phen(:,:) = 0.0 veg_state%g_leaf_acc(:,:) = 0.0 veg_state%g_leaf_phen_acc(:,:) = 0.0 @@ -240,6 +249,8 @@ SUBROUTINE veg3_field_deallocate() DEALLOCATE(veg_state%frac) DEALLOCATE(veg_state%g_leaf_phen) DEALLOCATE(veg_state%g_leaf) +DEALLOCATE(veg_state%g_leaf_day) +DEALLOCATE(veg_state%g_leaf_dr_out) DEALLOCATE(veg_state%lai_phen) DEALLOCATE(veg_state%g_leaf_acc) DEALLOCATE(veg_state%g_leaf_phen_acc) @@ -316,6 +327,8 @@ SUBROUTINE veg3_field_assoc(progs, ainfo, trifctl_data, trif_vars_data) veg_state%npp_dr_out => trifctl_data%npp_dr_out_pft veg_state%g_leaf => trifctl_data%g_leaf_pft veg_state%g_leaf_phen => trifctl_data%g_leaf_phen_pft + veg_state%g_leaf_day => trifctl_data%g_leaf_day_pft + veg_state%g_leaf_dr_out => trifctl_data%g_leaf_dr_out_pft veg_state%lai_phen => trifctl_data%lai_phen_pft veg_state%g_leaf_acc => trifctl_data%g_leaf_acc_pft veg_state%g_leaf_phen_acc => trifctl_data%g_leaf_phen_acc_pft From 51393ae3b6e0e1d307f53390c33e9c371c6dffc5 Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Mon, 14 Sep 2026 18:06:10 +0100 Subject: [PATCH 10/14] Updating to ensure that "g_leaf_dr_out" is in consistent units for diagnostic purposes, while having a seperate "g_leaf_dr" in turnover per second. --- .../shared/next_gen_biogeochem_mod.F90 | 40 ++++++++++++++----- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/src/control/shared/next_gen_biogeochem_mod.F90 b/src/control/shared/next_gen_biogeochem_mod.F90 index 91ec43c0..34096a36 100644 --- a/src/control/shared/next_gen_biogeochem_mod.F90 +++ b/src/control/shared/next_gen_biogeochem_mod.F90 @@ -194,6 +194,9 @@ SUBROUTINE veg3_run_ctrl( & REAL:: & npp_dr(land_pts,nnpft), & ! Mean NPP for driving vegetation (kg C/m2/s). +g_leaf_dr(land_pts,nnpft), & + ! Mean phenology-driven leaf turnover rate for driving litterfall and + ! vegetation dynamics (s-1). local_litter(land_pts,nnpft), & ! Litter production (kg C/m2/s). growth(land_pts,nnpft), & @@ -205,6 +208,7 @@ SUBROUTINE veg3_run_ctrl( & !Initialise Arrays npp_dr(:,:) = 0.0 +g_leaf_dr(:,:) = 0.0 mort_add(:,:,:) = 0.0 !----------------------------------------------------------------------------- @@ -215,7 +219,7 @@ SUBROUTINE veg3_run_ctrl( & ! veg-veg2a_jls_mod. !----------------------------------------------------------------------------- CALL veg3_phenol_couple(veg_index_pts,veg_index,veg3_ctrl,land_pts,nnpft, & - a_step,asteps_since_triffid,veg_state) + a_step,asteps_since_triffid,veg_state,g_leaf_dr) ! Now call vegetation model IF (asteps_since_triffid == veg3_ctrl%nstep_trif) THEN @@ -227,7 +231,7 @@ SUBROUTINE veg3_run_ctrl( & !IN parms litter_parms, & !IN fields - veg_state%g_leaf_dr_out, & + g_leaf_dr, & !IN state veg_state, & ! OUT Fields @@ -235,6 +239,9 @@ SUBROUTINE veg3_run_ctrl( & !OUT Diagnostics ) + ! Record driving g_leaf_dr and convert to s-1 -> (360 days)-1 + veg_state%g_leaf_dr_out(:,:) = g_leaf_dr * rsec_per_day * 360.0 + !CALL Allocation/Nitrogen/NSC @@ -303,7 +310,9 @@ SUBROUTINE veg3_phenol_couple( & veg_index_pts,veg_index,veg3_ctrl,land_pts,nnpft, & a_step,asteps_since_triffid, & !IN state - veg_state & + veg_state, & + !OUT Diagnostics + g_leaf_dr & ) ! Diagnoses leaf phenology and the mean phenology-driven leaf turnover rate @@ -343,6 +352,13 @@ SUBROUTINE veg3_phenol_couple( & !----------------------------------------------------------------------------- TYPE(veg_state_type),INTENT(IN OUT) :: veg_state +!----------------------------------------------------------------------------- +! Reals with INTENT OUT +!----------------------------------------------------------------------------- +REAL, INTENT(OUT) :: g_leaf_dr(land_pts,nnpft) + ! Mean phenology-driven leaf turnover rate for driving + ! litterfall and vegetation dynamics (s-1). + !----------------------------------------------------------------------------- !Local Vars !----------------------------------------------------------------------------- @@ -359,6 +375,8 @@ SUBROUTINE veg3_phenol_couple( & !End of headers +g_leaf_dr(:,:) = 0.0 + !----------------------------------------------------------------------------- ! Work out the phenology at its own timestep, appending to the accumulated ! leaf turnover rates. This is called independently of the vegetation @@ -425,9 +443,9 @@ SUBROUTINE veg3_phenol_couple( & IF (l_phenol) THEN ! Diagnose the mean phenological leaf turnover rate over the - ! coupling period and convert to JULES-standard per-second units. - veg_state%g_leaf_dr_out(l,n) = veg_state%g_leaf_phen_acc(l,n) * & - gam_trif / (rsec_per_day * 360.0) + ! coupling period, in JULES-standard per-second units. + g_leaf_dr(l,n) = veg_state%g_leaf_phen_acc(l,n) * & + gam_trif / (rsec_per_day * 360.0) ! Reset the accumulated phenological turnover ready for the next ! coupling period. @@ -435,8 +453,8 @@ SUBROUTINE veg3_phenol_couple( & ELSE ! No phenology - fall back to the raw accumulated physiological ! leaf turnover rate, as in veg-veg2a_jls_mod. - veg_state%g_leaf_dr_out(l,n) = veg_state%g_leaf_acc(l,n) * gam_trif / & - (rsec_per_day * 360.0) + g_leaf_dr(l,n) = veg_state%g_leaf_acc(l,n) * gam_trif / & + (rsec_per_day * 360.0) veg_state%g_leaf_acc(l,n) = 0.0 END IF @@ -445,12 +463,12 @@ SUBROUTINE veg3_phenol_couple( & ! the vegetation dynamics timestep. If it does, reduce the rate so ! that the turnover does not exceed the current LAI. IF (veg_state%lai(l,n) > 0.0) THEN - IF (veg_state%g_leaf_dr_out(l,n) * veg3_ctrl%dt_red > 1.0) THEN - veg_state%g_leaf_dr_out(l,n) = 1.0 / veg3_ctrl%dt_red + IF (g_leaf_dr(l,n) * veg3_ctrl%dt_red > 1.0) THEN + g_leaf_dr(l,n) = 1.0 / veg3_ctrl%dt_red END IF ELSE - veg_state%g_leaf_dr_out(l,n) = 0.0 + g_leaf_dr(l,n) = 0.0 END IF From 5c806bdc1508916cb70fe447a5854e1cff1df148 Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Tue, 15 Sep 2026 15:58:34 +0100 Subject: [PATCH 11/14] Merging jules_red_veg3_couple_phenology into jules_red_veg3_couple_soil (#91). --- .../shared/next_gen_biogeochem_mod.F90 | 292 +++++++++++++++++- src/control/shared/surf_couple_extra_mod.F90 | 8 +- src/control/shared/veg3_field_mod.F90 | 186 ++++++++++- src/initialisation/standalone/init_mod.F90 | 2 +- src/science/params/veg3_param_mod.F90 | 80 ++++- .../vegetation/veg3_red_dynamic_mod.F90 | 19 +- 6 files changed, 561 insertions(+), 26 deletions(-) diff --git a/src/control/shared/next_gen_biogeochem_mod.F90 b/src/control/shared/next_gen_biogeochem_mod.F90 index 34096a36..25e29d99 100644 --- a/src/control/shared/next_gen_biogeochem_mod.F90 +++ b/src/control/shared/next_gen_biogeochem_mod.F90 @@ -26,19 +26,19 @@ SUBROUTINE next_gen_biogeochem( & 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 ) !Only get the data structures - the data comes through the calling tree USE veg3_parm_mod, ONLY: & - veg3_ctrl_type,litter_parm_type,red_parm_type + veg3_ctrl_type,litter_parm_type,red_parm_type,soil_parm_type USE veg3_field_mod, ONLY: & - veg_state_type,red_state_type + veg_state_type,red_state_type,soil_state_type USE ancil_info, ONLY: ainfo_type @@ -50,12 +50,14 @@ SUBROUTINE next_gen_biogeochem( & TYPE(veg3_ctrl_type), INTENT(IN) :: veg3_ctrl TYPE(litter_parm_type), INTENT(IN) :: litter_parms TYPE(red_parm_type), INTENT(IN) :: red_parms +TYPE(soil_parm_type), INTENT(IN) :: soil_parms !---------------------------------------------------------------------------- ! Objects with INTENT inout !---------------------------------------------------------------------------- TYPE(veg_state_type), INTENT(IN OUT) :: veg_state TYPE(red_state_type), INTENT(IN OUT) :: red_state +TYPE(soil_state_type), INTENT(IN OUT) :: soil_state TYPE(ainfo_type), INTENT(IN OUT) :: ainfo !---------------------------------------------------------------------------- @@ -94,6 +96,10 @@ SUBROUTINE next_gen_biogeochem( & frac_vs(land_pts) ! Veg/Soil Fractional coverage +LOGICAL :: & + l_veg_step + ! Flag to indicate whether this is a vegetation dynamics timestep. + ! End of header !--------------------------------------------------------------------- @@ -113,6 +119,9 @@ SUBROUTINE next_gen_biogeochem( & END IF END DO + +l_veg_step = (asteps_since_triffid == veg3_ctrl%nstep_trif) + ! Call the Vegetation Biogeochemistry model IF (veg_index_pts > 0) CALL veg3_run_ctrl( & !IN Control vars @@ -126,7 +135,12 @@ SUBROUTINE next_gen_biogeochem( & ) ! Call the soil Biogeochemistry model -! This is where the soil biogechemistry ctrl call will be made +IF (veg_index_pts > 0 .AND. l_veg_step) CALL veg3_soil_couple( & + !IN Control vars + veg_index_pts,veg_index,land_pts,nnpft,veg3_ctrl,soil_parms, & + !INOUT state + veg_state,soil_state & + ) END SUBROUTINE next_gen_biogeochem @@ -479,4 +493,272 @@ SUBROUTINE veg3_phenol_couple( & END SUBROUTINE veg3_phenol_couple +!------------------------------------------------------------------------------ +SUBROUTINE veg3_soil_couple( & + !IN Control vars + veg_index_pts,veg_index,land_pts,nnpft,veg3_ctrl,soil_parms, & + !INOUT state + veg_state,soil_state & + ) + +! Couples RED litter carbon to the soil carbon model, mirroring the soil +! carbon section of TRIFFID (see triffid_jls.F90 and +! soilcarb_layers_jls_mod.F90/soilcarb_jls.F90). Requires the 4-pool soil +! carbon model (soil_bgc_model=soil_model_4pool), either layered +! (l_layeredC=.TRUE.) or single-layer; nitrogen is not yet coupled to +! veg3/RED. + +!Only get the data structures - the data comes through the calling tree +USE veg3_parm_mod, ONLY: veg3_ctrl_type, soil_parm_type +USE veg3_field_mod, ONLY: veg_state_type, soil_state_type + +!Access subroutines +USE soilcarb_layers_mod, ONLY: soilcarb_layers +USE soilcarb_mod, ONLY: soilcarb +#if !defined(UM_JULES) +USE soilcarb_mix_mod, ONLY: soilcarb_mix +#endif + +!Access some parameters direct from module +USE conversions_mod, ONLY: rsec_per_day + +IMPLICIT NONE + +!---------------------------------------------------------------------------- +! Integers with INTENT IN +!---------------------------------------------------------------------------- +INTEGER, INTENT(IN) :: land_pts,nnpft,veg_index_pts,veg_index(land_pts) + +!----------------------------------------------------------------------------- +! Objects with INTENT IN +!----------------------------------------------------------------------------- +TYPE(veg3_ctrl_type),INTENT(IN) :: veg3_ctrl +TYPE(soil_parm_type),INTENT(IN) :: soil_parms + +!----------------------------------------------------------------------------- +! Objects with INTENT INOUT +!----------------------------------------------------------------------------- +TYPE(veg_state_type),INTENT(IN OUT) :: veg_state +TYPE(soil_state_type),INTENT(IN OUT) :: soil_state + +!----------------------------------------------------------------------------- +!Local Vars +!----------------------------------------------------------------------------- +INTEGER :: l,n,k + ! Loop counters. + +REAL :: & +inv_timestep, & + ! Inverse soil carbon coupling timestep ((360d)-1). +resp_frac(land_pts,soil_parms%dim_cslayer), & + ! The fraction of soil respiration that forms new soil C (i.e. is NOT + ! released to the atmosphere). +resp_frac_cspool(land_pts,soil_parms%dim_cslayer,soil_parms%dim_cs1), & + ! As resp_frac, but broken down by soil C pool. Only used for + ! l_layeredC. +resp_s_dr(land_pts,soil_parms%dim_cslayer,5), & + ! Mean soil respiration for driving the soil carbon update + ! (kg C/m2/360days). NB 5=dim_cs1+1; the 5th element is workspace. +lit_n_t_gb(land_pts), & + ! Total nitrogen litter (kg N/m2/360days). The nitrogen cycle is not + ! yet coupled to veg3/RED, so this is always zero. +cs_tot(land_pts,soil_parms%dim_cslayer), & + ! Soil carbon content (kg C/m2). +ns_gb(land_pts,soil_parms%dim_cslayer), & + ! Total soil N on layers (kg N/m2). Always zero (see lit_n_t_gb). +neg_n(land_pts), & + ! Negative N required to prevent ns<0 (kg N). Unused (l_nitrogen=F). +implicit_resp_correction(land_pts), & + ! Respiration carried to next coupling period to account for applying + ! the minimum soil carbon constraint (kg m-2). +isunfrozen(land_pts,soil_parms%dim_cslayer), & + ! Matrix to mask out frozen layers (inaccessible to plants). Assumed + ! unfrozen throughout, consistent with the nitrogen cycle being off. +burnt_soil(land_pts), & + ! Burnt C in RPM and DPM pools (kg m-2 360d-1). Fire is not yet coupled + ! to veg3/RED, so this is always zero. +lit_frac(soil_parms%dim_cslayer), & + ! Litter fraction into each soil layer. +dcs(land_pts,soil_parms%dim_cslayer), & + ! Change in soil carbon over the coupling period (kg C/m2). +denom_resp, & + ! Denominator for calculating resp_s_acc_soilt. +#if !defined(UM_JULES) +dcs_pools(land_pts,soil_parms%dim_cslayer,4), & + ! Soil carbon by pool at the start of the coupling period, used to + ! calculate the layer mixing term (kg C/m2). +mix_s(land_pts,soil_parms%dim_cslayer-1,4), & + ! Diffusion coefficient for soil C between soil layers (m^2/360days). + ! Equation 15 of Burke et al. (2017), + ! https://www.geosci-model-dev.net/10/959/2017/gmd-10-959-2017.pdf +#endif +mix_term(land_pts,soil_parms%dim_cslayer,4), & + ! Mixing term for calculating the respiration correction + ! (kg C/m2/360days). +lit_resp + ! Net litter carbon reaching the soil after allowing for the change in + ! soil carbon and respiration (kg C/m2/360days). + +!End of headers + +! Nitrogen and fire are not yet coupled to veg3/RED soil carbon - keep the +! associated inputs at zero/unfrozen so that the shared soilcarb/ +! soilcarb_layers routines behave as a pure carbon-only 4-pool model. +lit_n_t_gb(:) = 0.0 +ns_gb(:,:) = 0.0 +burnt_soil(:) = 0.0 +isunfrozen(:,:) = 1.0 + +! Inverse soil carbon coupling timestep (/360days). Note this must be based +! on the vegetation dynamics/soil carbon coupling period (dt_red), not the +! raw physics timestep - resp_s_acc_soilt is accumulated over a full +! coupling period between calls to this routine, mirroring r_gamma/gam_trif +! in TRIFFID (veg-veg2a_jls_mod.F90). +inv_timestep = 1.0 / (veg3_ctrl%dt_red / rsec_per_day / 360.0) + +resp_s_dr(:,:,:) = 0.0 + +DO k = 1,veg_index_pts + l = veg_index(k) + DO n = 1,soil_parms%dim_cslayer + resp_s_dr(l,n,1) = soil_state%resp_s_acc_soilt(l,1,n,1) * inv_timestep + resp_s_dr(l,n,2) = soil_state%resp_s_acc_soilt(l,1,n,2) * inv_timestep + resp_s_dr(l,n,3) = soil_state%resp_s_acc_soilt(l,1,n,3) * inv_timestep + resp_s_dr(l,n,4) = soil_state%resp_s_acc_soilt(l,1,n,4) * inv_timestep + + soil_state%resp_s_dr_out_gb(l,n,1) = resp_s_dr(l,n,1) + soil_state%resp_s_dr_out_gb(l,n,2) = resp_s_dr(l,n,2) + soil_state%resp_s_dr_out_gb(l,n,3) = resp_s_dr(l,n,3) + soil_state%resp_s_dr_out_gb(l,n,4) = resp_s_dr(l,n,4) + soil_state%resp_s_dr_out_gb(l,n,5) = & + SUM(soil_state%resp_s_dr_out_gb(l,n,1:4)) + + ! Save the soil carbon at the start of the coupling period, used below + ! to diagnose the change in soil carbon. + dcs(l,n) = SUM(soil_state%cs_pool_soilt(l,1,n,1:4)) + +#if !defined(UM_JULES) + dcs_pools(l,n,:) = soil_state%cs_pool_soilt(l,1,n,1:4) +#endif + + ! Fraction of soil respiration that forms new soil C (i.e. is NOT + ! released to the atmosphere), calculated from the clay content. + resp_frac(l,n) = 1.0 / (soil_parms%resp_frac_a + soil_parms%resp_frac_b * & + EXP(soil_parms%resp_frac_c * 100.0 * & + soil_state%clay_soilt(l,1,n))) + resp_frac_cspool(l,n,:) = resp_frac(l,n) + + END DO + + ! Convert the litter carbon from per-PFT-area to gridbox mean units, as + ! required by soilcarb/soilcarb_layers. + DO n = 1,nnpft + veg_state%litCpft(l,n) = veg_state%litCpft(l,n) * veg_state%frac(l,n) + END DO + +END DO + +! Layered vs single-layer 4-pool soil carbon, as per triffid_jls.F90. +IF (soil_parms%l_layeredc) THEN + CALL soilcarb_layers(land_pts, veg_index_pts, veg_index, 0.0, & + inv_timestep, veg_state%litCpft, veg_state%litC, & + lit_n_t_gb, resp_frac_cspool, resp_s_dr, & + soil_state%cs_pool_soilt(:,1,:,:), & + soil_state%frac_c_label_pool_soilt(:,1,:,:), & + ns_gb, neg_n, implicit_resp_correction, burnt_soil, & + isunfrozen, soil_state%ns_pool_gb, & + soil_state%n_inorg_soilt_lyrs, & + soil_state%n_inorg_avail_pft, & + soil_state%t_soil_soilt_acc, & + soil_state%burnt_carbon_dpm, soil_state%g_burn_gb, & + soil_state%burnt_carbon_rpm, soil_state%minl_n_gb, & + soil_state%minl_n_pot_gb, soil_state%immob_n_gb, & + soil_state%immob_n_pot_gb, soil_state%fn_gb, & + soil_state%resp_s_diag_gb, & + soil_state%resp_s_pot_diag_gb, & + soil_state%dpm_ratio_gb, soil_state%n_gas_gb, & + soil_state%resp_s_to_atmos_gb, soil_state%sthu_soilt) +ELSE + CALL soilcarb(land_pts, veg_index_pts, veg_index, 0.0, & + inv_timestep, veg_state%litCpft, veg_state%litC, & + lit_n_t_gb, resp_frac(:,1), resp_s_dr, & + soil_state%cs_pool_soilt(:,1,:,:), & + ns_gb, neg_n, implicit_resp_correction, burnt_soil, & + soil_state%ns_pool_gb, soil_state%n_inorg_soilt_lyrs, & + soil_state%burnt_carbon_dpm, soil_state%g_burn_gb, & + soil_state%burnt_carbon_rpm, soil_state%minl_n_gb, & + soil_state%minl_n_pot_gb, soil_state%immob_n_gb, & + soil_state%immob_n_pot_gb, soil_state%fn_gb, & + soil_state%resp_s_diag_gb, soil_state%resp_s_pot_diag_gb, & + soil_state%dpm_ratio_gb, soil_state%n_gas_gb, & + soil_state%resp_s_to_atmos_gb, soil_state%sthu_soilt) +END IF + +! Convert the litter carbon back to per-PFT-area units. +DO k = 1,veg_index_pts + l = veg_index(k) + DO n = 1,nnpft + IF (veg_state%frac(l,n) > 0.0) THEN + veg_state%litCpft(l,n) = veg_state%litCpft(l,n) / veg_state%frac(l,n) + ELSE + veg_state%litCpft(l,n) = 0.0 + END IF + END DO +END DO + +!----------------------------------------------------------------------------- +! Diagnose the mean soil respiration to drive the soil carbon update over +! the next coupling period, correcting for any minimum soil carbon +! constraint applied above. +!----------------------------------------------------------------------------- +mix_term(:,:,:) = 0.0 +lit_frac(:) = 1.0 + +#if !defined(UM_JULES) +! Layer profile/mixing only apply with more than one soil C layer; for a +! single layer all litter reaches the one layer (lit_frac=1.0) and there is +! nothing to mix. +IF (soil_parms%l_layeredc) THEN + ! Calculate vertical profile of litter inputs. + lit_frac(1) = soil_parms%dzsoil(1) * & + EXP( -soil_parms%tau_lit * 0.5 * soil_parms%dzsoil(1) ) / & + soil_parms%litc_norm + DO n = 2,soil_parms%dim_cslayer + lit_frac(n) = soil_parms%dzsoil(n) * EXP( -soil_parms%tau_lit * & + (SUM(soil_parms%dzsoil(1:n-1)) + 0.5 * & + soil_parms%dzsoil(n)) ) / soil_parms%litc_norm + END DO + + ! Calculate the mixing term between soil layers. + CALL soilcarb_mix(land_pts, veg_index_pts, veg_index, dcs_pools, & + soil_state%t_soil_soilt_acc, mix_term, mix_s) +END IF +#endif + +DO k = 1,veg_index_pts + l = veg_index(k) + DO n = 1,soil_parms%dim_cslayer + cs_tot(l,n) = MAX(1.0e-10, SUM(soil_state%cs_pool_soilt(l,1,n,1:4))) + denom_resp = 1.0 / (cs_tot(l,n) * inv_timestep) + dcs(l,n) = cs_tot(l,n) - dcs(l,n) + + resp_s_dr(l,n,1) = SUM((1.0 - resp_frac_cspool(l,n,1:4)) * & + resp_s_dr(l,n,1:4)) + + lit_resp = veg_state%litC(l) * lit_frac(n) - (inv_timestep * dcs(l,n)) & + - resp_s_dr(l,n,1) - burnt_soil(l) + SUM(mix_term(l,n,:)) + + soil_state%resp_s_acc_soilt(l,1,n,1) = lit_resp * & + soil_state%cs_pool_soilt(l,1,n,1) * denom_resp + soil_state%resp_s_acc_soilt(l,1,n,2) = lit_resp * & + soil_state%cs_pool_soilt(l,1,n,2) * denom_resp + soil_state%resp_s_acc_soilt(l,1,n,3) = lit_resp * & + soil_state%cs_pool_soilt(l,1,n,3) * denom_resp + soil_state%resp_s_acc_soilt(l,1,n,4) = lit_resp * & + soil_state%cs_pool_soilt(l,1,n,4) * denom_resp + END DO +END DO + +END SUBROUTINE veg3_soil_couple + END MODULE next_gen_biogeochem_mod diff --git a/src/control/shared/surf_couple_extra_mod.F90 b/src/control/shared/surf_couple_extra_mod.F90 index d3f0d439..dd7b08a2 100644 --- a/src/control/shared/surf_couple_extra_mod.F90 +++ b/src/control/shared/surf_couple_extra_mod.F90 @@ -239,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 @@ -1079,9 +1079,9 @@ SUBROUTINE surf_couple_extra( & 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 ) diff --git a/src/control/shared/veg3_field_mod.F90 b/src/control/shared/veg3_field_mod.F90 index bc44a065..3da8a16e 100644 --- a/src/control/shared/veg3_field_mod.F90 +++ b/src/control/shared/veg3_field_mod.F90 @@ -21,6 +21,7 @@ MODULE veg3_field_mod !Use at module level USE ancil_info, ONLY: ainfo_type USE prognostics, ONLY: progs_type +USE p_s_parms, ONLY: psparms_type IMPLICIT NONE @@ -112,8 +113,65 @@ MODULE veg3_field_mod ! PFT number density across plant mass. (/m2) END TYPE red_state_type +! Structure to keep the soil state variables and fields used in coupling +! veg3/RED to the soil carbon (and, in future, nitrogen) model. +TYPE :: soil_state_type + REAL, POINTER :: & + cs_pool_soilt(:,:,:,:), & + ! Soil carbon pools (DPM,RPM,bio,humus), by soil tile and + ! layer. (kg C m-2) + frac_c_label_pool_soilt(:,:,:,:), & + ! Fraction of each soil C pool from labelled carbon. + ns_pool_gb(:,:,:), & + ! Soil nitrogen pools, by layer. (kg N m-2) + n_inorg_soilt_lyrs(:,:,:), & + ! Gridbox inorganic N pool on soil levels. (kg N m-2) + n_inorg_avail_pft(:,:,:), & + ! Available inorganic N for PFTs. (kg N m-2) + t_soil_soilt_acc(:,:,:), & + ! Accumulated soil temperature, used for soil C mixing. (K) + clay_soilt(:,:,:), & + ! Soil clay fraction, by layer. + sthu_soilt(:,:,:), & + ! Unfrozen soil moisture content as a fraction of saturation. + resp_s_acc_soilt(:,:,:,:), & + ! Accumulated soil respiration since the last coupling call. + ! (kg C m-2) + resp_s_dr_out_gb(:,:,:), & + ! Mean soil respiration driving the soil C update. (kg C m-2 + ! (360d)-1) + burnt_carbon_dpm(:), & + ! Burnt DPM carbon. (kg C m-2 (360d)-1) + burnt_carbon_rpm(:), & + ! Burnt RPM carbon. (kg C m-2 (360d)-1) + g_burn_gb(:), & + ! Gridbox mean fire disturbance rate. ((360d)-1) + minl_n_gb(:,:,:), & + ! Gross mineralisation of N. (kg N m-2 (360d)-1) + minl_n_pot_gb(:,:,:), & + ! Potential gross mineralisation of N. (kg N m-2 (360d)-1) + immob_n_gb(:,:,:), & + ! Immobilisation of N. (kg N m-2 (360d)-1) + immob_n_pot_gb(:,:,:), & + ! Potential immobilisation of N. (kg N m-2 (360d)-1) + fn_gb(:,:), & + ! Nitrogen decomposition rate modifier. + resp_s_diag_gb(:,:,:), & + ! Diagnosed soil respiration by pool. (kg C m-2 (360d)-1) + resp_s_pot_diag_gb(:,:,:), & + ! Diagnosed potential soil respiration by pool. + ! (kg C m-2 (360d)-1) + dpm_ratio_gb(:), & + ! Ratio of DPM carbon to total litter carbon. + n_gas_gb(:,:), & + ! Gaseous N loss. (kg N m-2 (360d)-1) + resp_s_to_atmos_gb(:,:) + ! Soil-to-atmosphere respiration flux. (kg C m-2 (360d)-1) +END TYPE soil_state_type + TYPE(veg_state_type) :: veg_state TYPE(red_state_type) :: red_state +TYPE(soil_state_type) :: soil_state !Private by default PRIVATE @@ -123,10 +181,10 @@ MODULE veg3_field_mod veg3_field_assoc, red_veg3_couple !Expose data -PUBLIC :: veg_state, red_state +PUBLIC :: veg_state, red_state, soil_state !Expose data structures -PUBLIC :: veg_state_type, red_state_type +PUBLIC :: veg_state_type, red_state_type, soil_state_type !Allow external code to read but not write !PROTECTED :: @@ -138,6 +196,9 @@ MODULE veg3_field_mod SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) +USE ancil_info, ONLY: dim_cslayer, dim_cs1, nsoilt +USE jules_soil_mod, ONLY: sm_levels + IMPLICIT NONE INTEGER, INTENT(IN) :: land_pts, nsurft, nnpft, nmasst @@ -216,6 +277,58 @@ SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) red_state%plantNumDensity(:,:,:) = 0.0 red_state%mort(:,:,:) = 0.0 +! Soil state + +! Allocate soil_state_type +ALLOCATE(soil_state%cs_pool_soilt (land_pts,nsoilt,dim_cslayer,dim_cs1)) +ALLOCATE(soil_state%frac_c_label_pool_soilt(land_pts,nsoilt,dim_cslayer,dim_cs1)) +ALLOCATE(soil_state%ns_pool_gb (land_pts,dim_cslayer,dim_cs1)) +ALLOCATE(soil_state%n_inorg_soilt_lyrs (land_pts,nsoilt,dim_cslayer)) +ALLOCATE(soil_state%n_inorg_avail_pft (land_pts,nnpft,dim_cslayer)) +ALLOCATE(soil_state%t_soil_soilt_acc (land_pts,nsoilt,sm_levels)) +ALLOCATE(soil_state%clay_soilt (land_pts,nsoilt,dim_cslayer)) +ALLOCATE(soil_state%sthu_soilt (land_pts,nsoilt,sm_levels)) +ALLOCATE(soil_state%resp_s_acc_soilt (land_pts,nsoilt,dim_cslayer,dim_cs1)) +ALLOCATE(soil_state%resp_s_dr_out_gb (land_pts,dim_cslayer,dim_cs1+1)) +ALLOCATE(soil_state%burnt_carbon_dpm (land_pts)) +ALLOCATE(soil_state%burnt_carbon_rpm (land_pts)) +ALLOCATE(soil_state%g_burn_gb (land_pts)) +ALLOCATE(soil_state%minl_n_gb (land_pts,dim_cslayer,dim_cs1+1)) +ALLOCATE(soil_state%minl_n_pot_gb (land_pts,dim_cslayer,dim_cs1+1)) +ALLOCATE(soil_state%immob_n_gb (land_pts,dim_cslayer,dim_cs1+1)) +ALLOCATE(soil_state%immob_n_pot_gb (land_pts,dim_cslayer,dim_cs1+1)) +ALLOCATE(soil_state%fn_gb (land_pts,dim_cslayer)) +ALLOCATE(soil_state%resp_s_diag_gb (land_pts,dim_cslayer,dim_cs1+1)) +ALLOCATE(soil_state%resp_s_pot_diag_gb (land_pts,dim_cslayer,dim_cs1+1)) +ALLOCATE(soil_state%dpm_ratio_gb (land_pts)) +ALLOCATE(soil_state%n_gas_gb (land_pts,dim_cslayer)) +ALLOCATE(soil_state%resp_s_to_atmos_gb (land_pts,dim_cslayer)) + +! Initialise soil_state_type +soil_state%cs_pool_soilt(:,:,:,:) = 0.0 +soil_state%frac_c_label_pool_soilt(:,:,:,:) = 0.0 +soil_state%ns_pool_gb(:,:,:) = 0.0 +soil_state%n_inorg_soilt_lyrs(:,:,:) = 0.0 +soil_state%n_inorg_avail_pft(:,:,:) = 0.0 +soil_state%t_soil_soilt_acc(:,:,:) = 0.0 +soil_state%clay_soilt(:,:,:) = 0.0 +soil_state%sthu_soilt(:,:,:) = 0.0 +soil_state%resp_s_acc_soilt(:,:,:,:) = 0.0 +soil_state%resp_s_dr_out_gb(:,:,:) = 0.0 +soil_state%burnt_carbon_dpm(:) = 0.0 +soil_state%burnt_carbon_rpm(:) = 0.0 +soil_state%g_burn_gb(:) = 0.0 +soil_state%minl_n_gb(:,:,:) = 0.0 +soil_state%minl_n_pot_gb(:,:,:) = 0.0 +soil_state%immob_n_gb(:,:,:) = 0.0 +soil_state%immob_n_pot_gb(:,:,:) = 0.0 +soil_state%fn_gb(:,:) = 0.0 +soil_state%resp_s_diag_gb(:,:,:) = 0.0 +soil_state%resp_s_pot_diag_gb(:,:,:) = 0.0 +soil_state%dpm_ratio_gb(:) = 0.0 +soil_state%n_gas_gb(:,:) = 0.0 +soil_state%resp_s_to_atmos_gb(:,:) = 0.0 + RETURN END SUBROUTINE veg3_field_allocate @@ -263,11 +376,35 @@ SUBROUTINE veg3_field_deallocate() DEALLOCATE(red_state%plantNumDensity) +DEALLOCATE(soil_state%cs_pool_soilt) +DEALLOCATE(soil_state%frac_c_label_pool_soilt) +DEALLOCATE(soil_state%ns_pool_gb) +DEALLOCATE(soil_state%n_inorg_soilt_lyrs) +DEALLOCATE(soil_state%n_inorg_avail_pft) +DEALLOCATE(soil_state%t_soil_soilt_acc) +DEALLOCATE(soil_state%clay_soilt) +DEALLOCATE(soil_state%sthu_soilt) +DEALLOCATE(soil_state%resp_s_acc_soilt) +DEALLOCATE(soil_state%resp_s_dr_out_gb) +DEALLOCATE(soil_state%burnt_carbon_dpm) +DEALLOCATE(soil_state%burnt_carbon_rpm) +DEALLOCATE(soil_state%g_burn_gb) +DEALLOCATE(soil_state%minl_n_gb) +DEALLOCATE(soil_state%minl_n_pot_gb) +DEALLOCATE(soil_state%immob_n_gb) +DEALLOCATE(soil_state%immob_n_pot_gb) +DEALLOCATE(soil_state%fn_gb) +DEALLOCATE(soil_state%resp_s_diag_gb) +DEALLOCATE(soil_state%resp_s_pot_diag_gb) +DEALLOCATE(soil_state%dpm_ratio_gb) +DEALLOCATE(soil_state%n_gas_gb) +DEALLOCATE(soil_state%resp_s_to_atmos_gb) + RETURN END SUBROUTINE veg3_field_deallocate !------------------------------------------------------------------------------- -SUBROUTINE veg3_field_assoc(progs, ainfo, trifctl_data, trif_vars_data) +SUBROUTINE veg3_field_assoc(progs, ainfo, psparms, trifctl_data, trif_vars_data) ! Initial code to associate the veg3 and red fields to the rest of JULES ! This new routine moves out the pointers from veg3_set_fields to here @@ -289,6 +426,7 @@ SUBROUTINE veg3_field_assoc(progs, ainfo, trifctl_data, trif_vars_data) TYPE(progs_type), INTENT(IN) :: progs TYPE(ainfo_type), INTENT(IN) :: ainfo +TYPE(psparms_type), INTENT(IN) :: psparms TYPE(trifctl_data_type), INTENT(IN), TARGET :: trifctl_data TYPE(trif_vars_data_type), INTENT(IN), TARGET :: trif_vars_data ! End of header @@ -337,6 +475,38 @@ SUBROUTINE veg3_field_assoc(progs, ainfo, trifctl_data, trif_vars_data) veg_state%vegC => trifctl_data%cv_gb veg_state%litC => trifctl_data%lit_c_mn_gb + ! Soil state fields used in the veg3/RED soil carbon coupling. + ! Firstly progs + soil_state%cs_pool_soilt => progs%cs_pool_soilt + soil_state%frac_c_label_pool_soilt => progs%frac_c_label_pool_soilt + soil_state%ns_pool_gb => progs%ns_pool_gb + soil_state%n_inorg_soilt_lyrs => progs%n_inorg_soilt_lyrs + soil_state%n_inorg_avail_pft => progs%n_inorg_avail_pft + soil_state%t_soil_soilt_acc => progs%t_soil_soilt_acc + + ! Next psparms + soil_state%clay_soilt => psparms%clay_soilt + soil_state%sthu_soilt => psparms%sthu_soilt + + ! Next trifctl_data + soil_state%resp_s_acc_soilt => trifctl_data%resp_s_acc_soilt + soil_state%resp_s_dr_out_gb => trifctl_data%resp_s_dr_out_gb + + ! Finally trif_vars_data + soil_state%burnt_carbon_dpm => trif_vars_data%burnt_carbon_dpm + soil_state%burnt_carbon_rpm => trif_vars_data%burnt_carbon_rpm + soil_state%g_burn_gb => trif_vars_data%g_burn_gb + soil_state%minl_n_gb => trif_vars_data%minl_n_gb + soil_state%minl_n_pot_gb => trif_vars_data%minl_n_pot_gb + soil_state%immob_n_gb => trif_vars_data%immob_n_gb + soil_state%immob_n_pot_gb => trif_vars_data%immob_n_pot_gb + soil_state%fn_gb => trif_vars_data%fn_gb + soil_state%resp_s_diag_gb => trif_vars_data%resp_s_diag_gb + soil_state%resp_s_pot_diag_gb => trif_vars_data%resp_s_pot_diag_gb + soil_state%dpm_ratio_gb => trif_vars_data%dpm_ratio_gb + soil_state%n_gas_gb => trif_vars_data%n_gas_gb + soil_state%resp_s_to_atmos_gb => trif_vars_data%resp_s_to_atmos_gb + END IF RETURN @@ -580,6 +750,16 @@ SUBROUTINE red_veg3_couple(ainfo) END DO +! Litter must be positive for the 4-pool soil carbon model (soilcarb_layers +! calculates the DPM:RPM ratio as a fraction of the total litter, which +! would otherwise divide by zero when a point has no litter). If there is +! no litter, make it a very small positive value. +DO l = 1,land_pts + IF (SUM(veg_state%litCpft(l,:)) == 0.0) THEN + veg_state%litCpft(l,1) = TINY(0.0) + END IF +END DO + !Final aggregation to gridbox for vegetation carbon for diagnostic purposes veg_state%vegC = pfttiles_to_gbm(veg_state%vegCpft,ainfo,frac_surft_in & = veg_state%frac) diff --git a/src/initialisation/standalone/init_mod.F90 b/src/initialisation/standalone/init_mod.F90 index d4dd9879..fea0d31c 100644 --- a/src/initialisation/standalone/init_mod.F90 +++ b/src/initialisation/standalone/init_mod.F90 @@ -351,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(progs,ainfo,trifctl_data,trif_vars_data) +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) diff --git a/src/science/params/veg3_param_mod.F90 b/src/science/params/veg3_param_mod.F90 index 656874af..f215f24f 100644 --- a/src/science/params/veg3_param_mod.F90 +++ b/src/science/params/veg3_param_mod.F90 @@ -60,6 +60,33 @@ 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 @@ -67,13 +94,13 @@ MODULE veg3_parm_mod 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' @@ -83,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 @@ -113,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 @@ -134,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 @@ -157,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 @@ -238,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 @@ -264,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 @@ -277,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") diff --git a/src/science/vegetation/veg3_red_dynamic_mod.F90 b/src/science/vegetation/veg3_red_dynamic_mod.F90 index 6896138d..b6a37ef7 100644 --- a/src/science/vegetation/veg3_red_dynamic_mod.F90 +++ b/src/science/vegetation/veg3_red_dynamic_mod.F90 @@ -366,15 +366,9 @@ SUBROUTINE update_pft_size_structure( & IF (growth < 0.0) THEN ! Negative growth is represented as a downward shrinkage flux of - ! individuals through the mass classes, the lowest mass class has - ! nowhere lower to shrink into, so its share of the loss is instead matched - ! by an equivalent mortality rate. + ! individuals through the mass classes. mort(k) = mort_base + mort_add(k) - IF (k == 1) THEN - mort(k) = mort(k) - g_mass(k) / mass_mass(k) - END IF - IF (k < mclass) THEN ! Flux shrinking down into this class from the class above flux_in(k) = - plantNumDensity(k+1) * g_mass(k+1) & @@ -383,12 +377,17 @@ SUBROUTINE update_pft_size_structure( & flux_in(k) = 0.0 END IF - IF (k > 1) THEN + IF (k == 1) THEN + ! Truncate shrinkage at the lowest mass class - individuals cannot + ! shrink below the lowest mass class, so we reduce the litterfall flux + ! to account for the negative growth. + flux_out(k) = 0.0 + mort_litC = mort_litC + plantNumDensity(k) * g_mass(k) + + ELSE ! Flux shrinking out of this class into the class below flux_out(k) = - plantNumDensity(k) * g_mass(k) & / (mass_mass(k) - mass_mass(k-1)) - ELSE - flux_out(k) = 0.0 END IF ELSE From a478a7431fa01d89e5b69a7f5fe2b85df4eb8de1 Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Fri, 18 Sep 2026 11:51:55 +0100 Subject: [PATCH 12/14] Add NBP diagnostic and carbon-conserving litter flux for veg3/RED soil coupling (#91) Adds the net biosphere productivity (nbp_gb) diagnostic to veg3/RED, calculated in veg3_soil_couple as NPP minus soil respiration to the atmosphere, and exposes it via extract_var (falling back to the existing TRIFFID calculation when l_red is off). Gridbox mean NPP (npp_gb/npp_n_gb) is now aggregated from the PFT-level npp_dr_out diagnostic in red_veg3_couple. Replaces the explicit leaf/root/wood/mortality litter summation with an implicit calculation (litCpft_implicit) that derives litCpft from carbon conservation (NPP minus the change in standing PFT carbon), guaranteeing the litter flux is consistent with the actual change in vegetation carbon. NPP and litterfall fluxes are rescaled by frac_old/frac_flux so the total gridbox flux is preserved when PFT fraction changes. Also: capitalise P/G in growth_onto_mass_class for consistency with the P_s notation used in the RED growth partition equations, fix a stale unit comment (m2 360d-1), and apply UMDP3 styling. --- .../shared/next_gen_biogeochem_mod.F90 | 41 ++++-- src/control/shared/veg3_field_mod.F90 | 129 ++++++++++++++++-- src/io/model_interface/extract_var.inc | 27 ++-- .../vegetation/veg3_red_dynamic_mod.F90 | 22 +-- 4 files changed, 177 insertions(+), 42 deletions(-) diff --git a/src/control/shared/next_gen_biogeochem_mod.F90 b/src/control/shared/next_gen_biogeochem_mod.F90 index 25e29d99..c4534334 100644 --- a/src/control/shared/next_gen_biogeochem_mod.F90 +++ b/src/control/shared/next_gen_biogeochem_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 next_gen_biogeochem_mod @@ -28,7 +31,7 @@ SUBROUTINE next_gen_biogeochem( & !IN parms litter_parms,red_parms,soil_parms, & !INOUT data structures - veg_state,red_state,soil_state & + veg_state,red_state,soil_state & !OUT diagnostics ) @@ -97,7 +100,7 @@ SUBROUTINE next_gen_biogeochem( & ! Veg/Soil Fractional coverage LOGICAL :: & - l_veg_step + l_veg_step ! Flag to indicate whether this is a vegetation dynamics timestep. ! End of header @@ -135,7 +138,7 @@ SUBROUTINE next_gen_biogeochem( & ) ! Call the soil Biogeochemistry model -IF (veg_index_pts > 0 .AND. l_veg_step) CALL veg3_soil_couple( & +IF (veg_index_pts > 0 .AND. l_veg_step) CALL veg3_soil_couple( & !IN Control vars veg_index_pts,veg_index,land_pts,nnpft,veg3_ctrl,soil_parms, & !INOUT state @@ -498,7 +501,7 @@ SUBROUTINE veg3_soil_couple( & !IN Control vars veg_index_pts,veg_index,land_pts,nnpft,veg3_ctrl,soil_parms, & !INOUT state - veg_state,soil_state & + veg_state,soil_state & ) ! Couples RED litter carbon to the soil carbon model, mirroring the soil @@ -568,7 +571,7 @@ SUBROUTINE veg3_soil_couple( & ! Total soil N on layers (kg N/m2). Always zero (see lit_n_t_gb). neg_n(land_pts), & ! Negative N required to prevent ns<0 (kg N). Unused (l_nitrogen=F). -implicit_resp_correction(land_pts), & +implicit_resp_correction(land_pts), & ! Respiration carried to next coupling period to account for applying ! the minimum soil carbon constraint (kg m-2). isunfrozen(land_pts,soil_parms%dim_cslayer), & @@ -577,22 +580,22 @@ SUBROUTINE veg3_soil_couple( & burnt_soil(land_pts), & ! Burnt C in RPM and DPM pools (kg m-2 360d-1). Fire is not yet coupled ! to veg3/RED, so this is always zero. -lit_frac(soil_parms%dim_cslayer), & +lit_frac(soil_parms%dim_cslayer), & ! Litter fraction into each soil layer. -dcs(land_pts,soil_parms%dim_cslayer), & +dcs(land_pts,soil_parms%dim_cslayer), & ! Change in soil carbon over the coupling period (kg C/m2). -denom_resp, & +denom_resp, & ! Denominator for calculating resp_s_acc_soilt. #if !defined(UM_JULES) -dcs_pools(land_pts,soil_parms%dim_cslayer,4), & +dcs_pools(land_pts,soil_parms%dim_cslayer,4), & ! Soil carbon by pool at the start of the coupling period, used to ! calculate the layer mixing term (kg C/m2). -mix_s(land_pts,soil_parms%dim_cslayer-1,4), & - ! Diffusion coefficient for soil C between soil layers (m^2/360days). +mix_s(land_pts,soil_parms%dim_cslayer-1,4), & + ! Diffusion coefficient for soil C between soil layers (m2 360d-1). ! Equation 15 of Burke et al. (2017), ! https://www.geosci-model-dev.net/10/959/2017/gmd-10-959-2017.pdf #endif -mix_term(land_pts,soil_parms%dim_cslayer,4), & +mix_term(land_pts,soil_parms%dim_cslayer,4), & ! Mixing term for calculating the respiration correction ! (kg C/m2/360days). lit_resp @@ -720,12 +723,12 @@ SUBROUTINE veg3_soil_couple( & ! nothing to mix. IF (soil_parms%l_layeredc) THEN ! Calculate vertical profile of litter inputs. - lit_frac(1) = soil_parms%dzsoil(1) * & - EXP( -soil_parms%tau_lit * 0.5 * soil_parms%dzsoil(1) ) / & + lit_frac(1) = soil_parms%dzsoil(1) * & + EXP( -soil_parms%tau_lit * 0.5 * soil_parms%dzsoil(1) ) / & soil_parms%litc_norm DO n = 2,soil_parms%dim_cslayer lit_frac(n) = soil_parms%dzsoil(n) * EXP( -soil_parms%tau_lit * & - (SUM(soil_parms%dzsoil(1:n-1)) + 0.5 * & + (SUM(soil_parms%dzsoil(1:n-1)) + 0.5 * & soil_parms%dzsoil(n)) ) / soil_parms%litc_norm END DO @@ -759,6 +762,14 @@ SUBROUTINE veg3_soil_couple( & END DO END DO +! Diagnose the net biosphere productivity (NPP minus soil respiration to +! atmosphere) for each land point. +DO k = 1,veg_index_pts + l = veg_index(k) + veg_state%nbp_gb(l) = veg_state%npp_n_gb(l) - & + SUM(soil_state%resp_s_to_atmos_gb(l,:)) +END DO + END SUBROUTINE veg3_soil_couple END MODULE next_gen_biogeochem_mod diff --git a/src/control/shared/veg3_field_mod.F90 b/src/control/shared/veg3_field_mod.F90 index 3da8a16e..71b0ec8a 100644 --- a/src/control/shared/veg3_field_mod.F90 +++ b/src/control/shared/veg3_field_mod.F90 @@ -40,6 +40,13 @@ MODULE veg3_field_mod ! Total PFT carbon density per PFT area fraction. (kg C m-2) vegC(:), & ! Gridbox mean vegetation carbon. (kg C m-2) + npp_gb(:), & + ! Gridbox mean NPP driving RED. (kg C m-2 s-1) + npp_n_gb(:), & + ! Gridbox mean NPP after nitrogen limitation. Nitrogen is not + ! yet coupled to veg3/RED, so this is npp_gb converted to + ! (360d)-1 units, matching trif_vars_data%npp_n_gb. + ! (kg C m-2 (360d)-1) npp_acc(:,:), & ! Accumulated NPP. (kg C m-2 s-1) npp_dr_out(:,:), & @@ -84,9 +91,15 @@ MODULE veg3_field_mod ! Total litter carbon flux per gridbox. (kg C m-2 (360d)-1) REAL, ALLOCATABLE :: & - mort_litC(:,:) + mort_litC(:,:), & ! Mortality/demographic litter carbon flux from vegetation ! dynamics, normalised per unit PFT canopy area (kg C m-2 s-1). + nbp_gb(:) + ! Gridbox mean net biosphere productivity (NPP minus all + ! carbon fluxes out of land). Only soil respiration is + ! currently coupled to veg3/RED, so this is npp_n_gb minus the + ! soil-to-atmosphere respiration flux (see veg3_soil_couple). + ! (kg C m-2 (360d)-1) END TYPE veg_state_type @@ -223,12 +236,15 @@ SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) ALLOCATE(veg_state%g_leaf_phen_acc ( land_pts, nnpft) ) ALLOCATE(veg_state%frac ( land_pts, nsurft) ) ALLOCATE(veg_state%vegC ( land_pts) ) +ALLOCATE(veg_state%npp_gb ( land_pts) ) +ALLOCATE(veg_state%npp_n_gb ( land_pts) ) ALLOCATE(veg_state%leaf_litC ( land_pts, nnpft) ) ALLOCATE(veg_state%root_litC ( land_pts, nnpft) ) ALLOCATE(veg_state%wood_litC ( land_pts, nnpft) ) ALLOCATE(veg_state%litCpft ( land_pts, nnpft) ) ALLOCATE(veg_state%litC ( land_pts) ) ALLOCATE(veg_state%mort_litC ( land_pts, nnpft) ) +ALLOCATE(veg_state%nbp_gb ( land_pts) ) !Initialise veg_state%leafC(:,:) = 0.0 @@ -250,12 +266,15 @@ SUBROUTINE veg3_field_allocate(land_pts,nsurft,nnpft,nmasst) veg_state%g_leaf_phen_acc(:,:) = 0.0 veg_state%frac(:,:) = 0.0 veg_state%vegC(:) = 0.0 +veg_state%npp_gb(:) = 0.0 +veg_state%npp_n_gb(:) = 0.0 veg_state%leaf_litC(:,:) = 0.0 veg_state%root_litC(:,:) = 0.0 veg_state%wood_litC(:,:) = 0.0 veg_state%litCpft(:,:) = 0.0 veg_state%litC(:) = 0.0 veg_state%mort_litC(:,:) = 0.0 +veg_state%nbp_gb(:) = 0.0 ! RED @@ -373,6 +392,8 @@ SUBROUTINE veg3_field_deallocate() DEALLOCATE(veg_state%litCpft) DEALLOCATE(veg_state%litC) DEALLOCATE(veg_state%vegC) +DEALLOCATE(veg_state%npp_gb) +DEALLOCATE(veg_state%npp_n_gb) DEALLOCATE(red_state%plantNumDensity) @@ -474,6 +495,8 @@ SUBROUTINE veg3_field_assoc(progs, ainfo, psparms, trifctl_data, trif_vars_data) veg_state%litCpft => trifctl_data%lit_c_pft veg_state%vegC => trifctl_data%cv_gb veg_state%litC => trifctl_data%lit_c_mn_gb + veg_state%npp_gb => trifctl_data%npp_gb + veg_state%npp_n_gb => trif_vars_data%npp_n_gb ! Soil state fields used in the veg3/RED soil carbon coupling. ! Firstly progs @@ -684,6 +707,7 @@ SUBROUTINE red_veg3_couple(ainfo) USE gridbox_mean_mod, ONLY: pfttiles_to_gbm, & masstiles_to_pfttiles USE conversions_mod, ONLY: rsec_per_day +USE veg3_parm_mod, ONLY: veg3_ctrl IMPLICIT NONE @@ -695,9 +719,21 @@ SUBROUTINE red_veg3_couple(ainfo) !Local INTEGER :: l,n ! Index variables. +REAL :: frac_old(land_pts,nnpft) + ! PFT fraction before this call updates veg_state%frac. +REAL :: vegCpft_old(land_pts,nnpft) + ! PFT carbon density before this call updates veg_state%vegCpft. + ! Used to derive the implicit litter flux below. +REAL :: frac_flux(land_pts,nnpft) + ! Representative PFT fraction used to convert per-PFT-area fluxes + ! to/from the gridbox mean over this coupling step, taken as the + ! midpoint of the old and new PFT fraction. + !----------------------------------------------------------------------------- !end of header +frac_old(:,:) = veg_state%frac(:,1:nnpft) +vegCpft_old(:,:) = veg_state%vegCpft(:,:) veg_state%vegCpft(:,:) = 0.0 veg_state%lai_bal(:,:) = 0.0 @@ -735,14 +771,38 @@ SUBROUTINE red_veg3_couple(ainfo) veg_state%woodC(l,n) = veg_state%vegCpft(l,n) - veg_state%leafC(l,n) & - veg_state%rootC(l,n) + ! NPP and its derived litter fluxes are normalised per unit PFT canopy - ! Aggregate the leaf/root/wood turnover litter and the mortality litter - ! demographic litter for the total litter flux per PFT fraction. - veg_state%litCpft(l,n) = (veg_state%leaf_litC(l,n) + & - veg_state%root_litC(l,n) + & - veg_state%wood_litC(l,n) + & - veg_state%mort_litC(l,n) * rsec_per_day * & - 360.0) + ! area. frac_flux is the representative PFT fraction used to convert + ! between per-PFT-area and gridbox-mean quantities over this coupling + ! step, taken as the midpoint of the old and new PFT fraction (mirroring + ! TRIFFID's own frac_flux normalisation of lit_c). + frac_flux(l,n) = 0.5 * (frac_old(l,n) + veg_state%frac(l,n)) + + ! Derive litCpft implicitly from carbon conservation (NPP minus the + ! change in standing PFT carbon), which guarantees the litter flux is + ! consistent with the actual change in vegetation carbon. + CALL litCpft_implicit(veg_state%npp_dr_out(l,n), frac_old(l,n), & + veg_state%frac(l,n), frac_flux(l,n), & + vegCpft_old(l,n), veg_state%vegCpft(l,n), & + veg3_ctrl%dt_red, veg_state%litCpft(l,n)) + + ! Rescale NPP and litterfall fluxes by frac_old/frac_flux so the total + ! gridbox flux is retained under the new PFT area. Left unchanged if + ! frac_flux is zero. + IF (frac_flux(l,n) > 0.0) THEN + veg_state%npp_dr_out(l,n) = veg_state%npp_dr_out(l,n) * frac_old(l,n) & + / frac_flux(l,n) + veg_state%leaf_litC(l,n) = veg_state%leaf_litC(l,n) * frac_old(l,n) & + / frac_flux(l,n) + veg_state%root_litC(l,n) = veg_state%root_litC(l,n) * frac_old(l,n) & + / frac_flux(l,n) + veg_state%wood_litC(l,n) = veg_state%wood_litC(l,n) * frac_old(l,n) & + / frac_flux(l,n) + veg_state%mort_litC(l,n) = veg_state%mort_litC(l,n) * frac_old(l,n) & + / frac_flux(l,n) + END IF + ! Update bare soil veg_state%frac(l,soil) = MAX(0.0, 1.0 - SUM(veg_state%frac(l,1:nnpft))) @@ -768,6 +828,15 @@ SUBROUTINE red_veg3_couple(ainfo) veg_state%litC(:) = pfttiles_to_gbm(veg_state%litCpft,ainfo,frac_surft_in & = veg_state%frac) +! Aggregate the per-PFT NPP driving RED to gridbox mean diagnostics. +! npp_dr_out is in kg C m-2 (360d)-1, matching npp_n_gb's units directly. +! npp_gb instead uses kg C m-2 s-1, so is converted back from (360d)-1. +! Nitrogen is not yet coupled to veg3/RED, so npp_n_gb is currently just +! the (360d)-1 equivalent of npp_gb. +veg_state%npp_n_gb(:) = pfttiles_to_gbm(veg_state%npp_dr_out,ainfo, & + frac_surft_in = veg_state%frac) +veg_state%npp_gb(:) = veg_state%npp_n_gb(:) / (rsec_per_day * 360.0) + RETURN END SUBROUTINE red_veg3_couple @@ -861,4 +930,48 @@ SUBROUTINE pft_mean_from_mass_class( & END SUBROUTINE pft_mean_from_mass_class !----------------------------------------------------------------------------- +SUBROUTINE litCpft_implicit(npp_dr, frac_old, frac_new, frac_flux, & + vegCpft_old, vegCpft_new, dt, litCpft) +!----------------------------------------------------------------------------- +! Derives the PFT litter flux implicitly from carbon conservation (NPP minus +! the change in standing PFT carbon), rather than summing the explicit +! leaf/root/wood/mortality litter fluxes: +! litCpft = npp_dr - d(vegCpft*frac) / (frac_flux*dt) * rsec_per_day*360 +! frac_flux (the midpoint of frac_old and frac_new) normalises the fluxes, +! mirroring TRIFFID's own frac_flux normalisation of lit_c. +!----------------------------------------------------------------------------- +USE conversions_mod, ONLY: rsec_per_day + +IMPLICIT NONE + +REAL, INTENT(IN) :: npp_dr + ! PFT NPP driving RED, normalised per unit PFT area. + ! (kg C m-2 (360d)-1) +REAL, INTENT(IN) :: frac_old, frac_new + ! PFT fraction before/after this call. (-) +REAL, INTENT(IN) :: frac_flux + ! Representative PFT fraction used to normalise the litter + ! flux, taken as the midpoint of frac_old and frac_new. (-) +REAL, INTENT(IN) :: vegCpft_old, vegCpft_new + ! PFT carbon density per PFT area before/after this call. + ! (kg C m-2) +REAL, INTENT(IN) :: dt + ! Vegetation dynamics timestep over which growth was applied + ! and vegCpft changed. (s) +REAL, INTENT(OUT) :: litCpft + ! Implicit litter flux consistent with carbon conservation. + ! (kg C m-2 (360d)-1) + +!End of header + +IF (frac_flux > 0.0) THEN + litCpft = npp_dr - (vegCpft_new * frac_new - vegCpft_old * frac_old) & + / (frac_flux * dt) * rsec_per_day * 360.0 +ELSE + litCpft = 0.0 +END IF + +END SUBROUTINE litCpft_implicit +!----------------------------------------------------------------------------- + END MODULE veg3_field_mod diff --git a/src/io/model_interface/extract_var.inc b/src/io/model_interface/extract_var.inc index 83077bc9..20139acf 100644 --- a/src/io/model_interface/extract_var.inc +++ b/src/io/model_interface/extract_var.inc @@ -44,7 +44,10 @@ USE jules_radiation_mod, ONLY: & wght_alb USE jules_vegetation_mod, ONLY: & - l_vegdrag_surft + l_vegdrag_surft, l_red + +USE veg3_field_mod, ONLY: & + veg_state USE c_z0h_z0m, ONLY: & z0h_z0m @@ -1786,15 +1789,19 @@ CASE ( 'npp_n_gb' ) cube = map_from_land(cube_land) CASE ( 'nbp_gb' ) - cube_land = cube_from_array(trif_vars%npp_n_gb - & - SUM(trif_vars%resp_s_to_atmos_gb(:,:),2) - & - trif_vars%wp_fast_out_gb - & - trif_vars%wp_med_out_gb - & - trif_vars%wp_slow_out_gb - & - trif_vars%harvest_gb - & - trif_vars%veg_c_fire_emission_gb - & - trif_vars%burnt_carbon_dpm - & - trif_vars%burnt_carbon_rpm ) + IF ( l_red ) THEN + cube_land = cube_from_array(veg_state%nbp_gb) + ELSE + cube_land = cube_from_array(trif_vars%npp_n_gb - & + SUM(trif_vars%resp_s_to_atmos_gb(:,:),2) - & + trif_vars%wp_fast_out_gb - & + trif_vars%wp_med_out_gb - & + trif_vars%wp_slow_out_gb - & + trif_vars%harvest_gb - & + trif_vars%veg_c_fire_emission_gb - & + trif_vars%burnt_carbon_dpm - & + trif_vars%burnt_carbon_rpm ) + END IF cube = map_from_land(cube_land) CASE ( 'n_demand' ) diff --git a/src/science/vegetation/veg3_red_dynamic_mod.F90 b/src/science/vegetation/veg3_red_dynamic_mod.F90 index b6a37ef7..ba21f067 100644 --- a/src/science/vegetation/veg3_red_dynamic_mod.F90 +++ b/src/science/vegetation/veg3_red_dynamic_mod.F90 @@ -206,11 +206,11 @@ SUBROUTINE growth_onto_mass_class( & INTEGER :: k REAL :: & -p, & +P, & ! The total PFT carbon assimilate across the gridbox. ! (kgC m-2 s-1) -g, & - ! Total gridbox carbon assimilate devoted to vegetation +G, & + ! Total gridbox carbon assimilate devoted to vegetation ! structural growth. (kgC m-2 s-1) plantNumDensity_g_sum ! Summation of the relative cohort contribution towards the @@ -230,12 +230,12 @@ SUBROUTINE growth_onto_mass_class( & END DO ! Partition the growth into recruitment and structural growth -p= frac * growth -P_s = alpha_recrt * p -g = (1.0 - alpha_recrt) * p +P = frac * growth +P_s = alpha_recrt * P +G = (1.0 - alpha_recrt) * P IF (plantNumDensity_g_sum > 0) THEN - g0 = g / plantNumDensity_g_sum + g0 = G / plantNumDensity_g_sum END IF IF (growth < 0.0) THEN @@ -246,7 +246,7 @@ SUBROUTINE growth_onto_mass_class( & ! Adjust g0 to account for the loss of recruitment IF (plantNumDensity_g_sum > 0) THEN - g0 = p / plantNumDensity_g_sum + g0 = P / plantNumDensity_g_sum END IF END IF @@ -380,7 +380,7 @@ SUBROUTINE update_pft_size_structure( & IF (k == 1) THEN ! Truncate shrinkage at the lowest mass class - individuals cannot ! shrink below the lowest mass class, so we reduce the litterfall flux - ! to account for the negative growth. + ! to account for the negative growth. flux_out(k) = 0.0 mort_litC = mort_litC + plantNumDensity(k) * g_mass(k) @@ -445,8 +445,12 @@ SUBROUTINE update_pft_size_structure( & ! If the resultant vegetation fraction is less than the minimum ! fraction, add trees to the lowest mass class to make up the difference IF (frac_check < frac_min) THEN + ! Take this additional carbon from the litterfall flux + mort_litC = mort_litC - (frac_min - frac_check)/crwn_area_mass(1) * & + mass_mass(1) / dt plantNumDensity(1) = plantNumDensity(1) & +(frac_min - frac_check) / crwn_area_mass(1) + END IF END SUBROUTINE update_pft_size_structure From eff1267659fb084f293703df6f2e23897470ef9f Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Fri, 18 Sep 2026 11:57:41 +0100 Subject: [PATCH 13/14] Fixing veg3_red_dynamic_mod.F90 with umpd3_fixer.py --- src/science/vegetation/veg3_red_dynamic_mod.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/science/vegetation/veg3_red_dynamic_mod.F90 b/src/science/vegetation/veg3_red_dynamic_mod.F90 index ba21f067..94974060 100644 --- a/src/science/vegetation/veg3_red_dynamic_mod.F90 +++ b/src/science/vegetation/veg3_red_dynamic_mod.F90 @@ -206,10 +206,10 @@ SUBROUTINE growth_onto_mass_class( & INTEGER :: k REAL :: & -P, & +P_a, & ! The total PFT carbon assimilate across the gridbox. ! (kgC m-2 s-1) -G, & +G_struct, & ! Total gridbox carbon assimilate devoted to vegetation ! structural growth. (kgC m-2 s-1) plantNumDensity_g_sum @@ -230,12 +230,12 @@ SUBROUTINE growth_onto_mass_class( & END DO ! Partition the growth into recruitment and structural growth -P = frac * growth -P_s = alpha_recrt * P -G = (1.0 - alpha_recrt) * P +P_a = frac * growth +P_s = alpha_recrt * P_a +G_struct = (1.0 - alpha_recrt) * P_a IF (plantNumDensity_g_sum > 0) THEN - g0 = G / plantNumDensity_g_sum + g0 = G_struct / plantNumDensity_g_sum END IF IF (growth < 0.0) THEN @@ -246,7 +246,7 @@ SUBROUTINE growth_onto_mass_class( & ! Adjust g0 to account for the loss of recruitment IF (plantNumDensity_g_sum > 0) THEN - g0 = P / plantNumDensity_g_sum + g0 = G_struct / plantNumDensity_g_sum END IF END IF From 9e4b0f6f792ea6ef82bd58857f0f49ad44b51bda Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Fri, 18 Sep 2026 12:32:09 +0100 Subject: [PATCH 14/14] Undoing change where we were trying to load in veg_state with extract_var.inc. This was causing issues with the rivers only routine. We won't output nbp_gb in this PR/issue. This can be tackled at a later point. --- src/io/model_interface/extract_var.inc | 27 ++++++++++---------------- 1 file changed, 10 insertions(+), 17 deletions(-) diff --git a/src/io/model_interface/extract_var.inc b/src/io/model_interface/extract_var.inc index 20139acf..83077bc9 100644 --- a/src/io/model_interface/extract_var.inc +++ b/src/io/model_interface/extract_var.inc @@ -44,10 +44,7 @@ USE jules_radiation_mod, ONLY: & wght_alb USE jules_vegetation_mod, ONLY: & - l_vegdrag_surft, l_red - -USE veg3_field_mod, ONLY: & - veg_state + l_vegdrag_surft USE c_z0h_z0m, ONLY: & z0h_z0m @@ -1789,19 +1786,15 @@ CASE ( 'npp_n_gb' ) cube = map_from_land(cube_land) CASE ( 'nbp_gb' ) - IF ( l_red ) THEN - cube_land = cube_from_array(veg_state%nbp_gb) - ELSE - cube_land = cube_from_array(trif_vars%npp_n_gb - & - SUM(trif_vars%resp_s_to_atmos_gb(:,:),2) - & - trif_vars%wp_fast_out_gb - & - trif_vars%wp_med_out_gb - & - trif_vars%wp_slow_out_gb - & - trif_vars%harvest_gb - & - trif_vars%veg_c_fire_emission_gb - & - trif_vars%burnt_carbon_dpm - & - trif_vars%burnt_carbon_rpm ) - END IF + cube_land = cube_from_array(trif_vars%npp_n_gb - & + SUM(trif_vars%resp_s_to_atmos_gb(:,:),2) - & + trif_vars%wp_fast_out_gb - & + trif_vars%wp_med_out_gb - & + trif_vars%wp_slow_out_gb - & + trif_vars%harvest_gb - & + trif_vars%veg_c_fire_emission_gb - & + trif_vars%burnt_carbon_dpm - & + trif_vars%burnt_carbon_rpm ) cube = map_from_land(cube_land) CASE ( 'n_demand' )