From c193294a65579cbce3113030ceebefe3440daa76 Mon Sep 17 00:00:00 2001 From: mo-arthurargles Date: Fri, 21 Aug 2026 16:46:09 +0100 Subject: [PATCH 01/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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