diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 9d443117..40d3fce6 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -8,3 +8,4 @@ | maggiehendry | Maggie Hendry | Met Office | 2026-08-21 | | DanCopsey | Dan Copsey | Met Office | 2026-08-24 | | Pierre-siddall | Pierre Siddall | Met Office | 2026-08-26 | +| mo-arthurargles | Arthur Argles | Met Office | 2026-09-04 | \ No newline at end of file diff --git a/src/control/shared/next_gen_biogeochem_mod.F90 b/src/control/shared/next_gen_biogeochem_mod.F90 index b5619d3a..34096a36 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 !----------------------------------------------------------------------------- @@ -188,32 +194,44 @@ 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), & ! growth (kg C/m2/s). -mort_add(land_pts,nnpft,nmasst), & +mort_add(land_pts,nnpft,nmasst) ! mortality above baseline (/m2) -demographic_lit(land_pts,nnpft) - ! demographic litter aggregated to PFT !End of headers !Initialise Arrays -npp_dr(:,:) = 0.0 -mort_add(:,:,:) = 0.0 +npp_dr(:,:) = 0.0 +g_leaf_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_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_dr, & !IN state veg_state, & ! OUT Fields @@ -221,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 @@ -229,13 +250,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 +275,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 +304,179 @@ 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 Diagnostics + g_leaf_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_dr(land_pts,nnpft) + ! Mean phenology-driven leaf turnover rate for driving + ! litterfall and 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. + +INTEGER :: l,n,k + ! Loop counters. + +!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 +! 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.0 + lai_bal_dummy(:) = 0.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) + 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, & + 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 + 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, 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. + 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_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_dr(l,n) * veg3_ctrl%dt_red > 1.0) THEN + g_leaf_dr(l,n) = 1.0 / veg3_ctrl%dt_red + + END IF + ELSE + g_leaf_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 76d4d279..d3f0d439 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..bc44a065 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,67 @@ 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) + 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) + 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 +119,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_deallocate, & + 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 +152,49 @@ 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%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) ) 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%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 +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 +208,153 @@ 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_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%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) +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) + +! 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. + +! 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 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 + + ! 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 + ! 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%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 + 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 @@ -160,39 +362,71 @@ 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 +INTEGER :: l,n !End of header !----------------------------------------------------------------------------- IF (l_red .AND. l_triffid) THEN + veg_state%phen(:,:) = 1.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 - - !----------------------------------------------------------------------------- - ! AJW The following code section temporarily sets up the veg3 model - ! In further revisions these fields will imported from the dump - !----------------------------------------------------------------------------- + DO l = 1, land_pts + DO n = 1, nnpft + + ! 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, 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) - 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 + ! 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 - 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 - ! 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 < 0.0) THEN + vegCpft = vegCpft / frac + lai_bal = lai_bal / frac + canht = canht / frac +ELSE + + ! 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 pft_mean_from_mass_class +!----------------------------------------------------------------------------- + END MODULE veg3_field_mod diff --git a/src/initialisation/standalone/init_mod.F90 b/src/initialisation/standalone/init_mod.F90 index b7850e9f..d4dd9879 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,trifctl_data,trif_vars_data) CALL chemvars_assoc(chemvars,chemvars_data) CALL water_resources_assoc(water_resources,water_resources_data) CALL wtrac_jls_assoc(wtrac_jls,wtrac_jls_data) diff --git a/src/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..9f821f6d 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(KIND=real_jlslsm), 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..bb7b285f 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.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) END DO diff --git a/src/science/vegetation/veg3_red_dynamic_mod.F90 b/src/science/vegetation/veg3_red_dynamic_mod.F90 index 84b43607..6896138d 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 ) @@ -61,190 +62,395 @@ 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) -!----------------------------------------------------------------------------- -! 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 m-2 s-1) 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 s-1) +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 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 - - END IF + ! 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), & + !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) & + ) + + !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 - END DO -END DO -! Loop updates the demographic distribution. -DO k = 1, nmasst - DO n = 1, nnpft - DO l = 1, land_pts - - g_mass(l,n,k) = g0(l,n) * red_state%g_mass_scale(n,k) - - 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 - - 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 - - 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 + ! 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), & + 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 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) + ELSE + veg_state%mort_litC(l,n) = 0.0 + END IF + END DO END DO END SUBROUTINE veg3_red_dynamic !----------------------------------------------------------------------------- +!----------------------------------------------------------------------------- +SUBROUTINE growth_onto_mass_class( & + !IN sizing + mclass, & + !IN PFT parameters + alpha_recrt, & + !IN fields + frac,growth, & + !IN mass-cohort properties + plantNumDensity,g_mass_scale, & + !OUT fields + P_s, g0 & + ) + +IMPLICIT NONE + +!----------------------------------------------------------------------------- +! Integers with INTENT IN +!----------------------------------------------------------------------------- +INTEGER, INTENT(IN) :: mclass + ! Number of mass classes for this PFT. + +!----------------------------------------------------------------------------- +! 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 m-2 s-1) +plantNumDensity(mclass), & + ! Population density within each mass cohort. (m-2) +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 m-2 s-1) +g0 + ! Boundary growth for an individual member of the smallest mass + ! cohort. (kgC s-1) + +!----------------------------------------------------------------------------- +!Local Vars +!----------------------------------------------------------------------------- +INTEGER :: k + +REAL :: & +p, & + ! The total PFT carbon assimilate across the gridbox. + ! (kgC m-2 s-1) +g, & + ! Total gridbox carbon assimilate devoted to vegetation + ! structural growth. (kgC m-2 s-1) +plantNumDensity_g_sum + ! Summation of the relative cohort contribution towards the + ! total PFT assimilate (m-2) + +!End of headers + +! Initialise vars +P_s = 0.0 +g0 = 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 + +! 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 + ! 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 + + ! Adjust g0 to account for the loss of recruitment + IF (plantNumDensity_g_sum > 0) THEN + g0 = p / plantNumDensity_g_sum + 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,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-1) +frac_min, & + ! Minimum vegetation fraction for this PFT. (-) +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 s-1) +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(IN OUT) :: & +plantNumDensity(mclass), & + ! Population density within each mass cohort. (m-2) +mort(mclass) + ! 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 gridbox area, veg3_red_dynamics renormalises to + ! PFT canopy area. (kgC m-2 s-1) + +!----------------------------------------------------------------------------- +!Local Vars +!----------------------------------------------------------------------------- +INTEGER :: k + +REAL :: & +g_mass(mclass), & + ! Individual growth across the mass cohorts. (kgC s-1) +dplantNumDensity_dt(mclass), & + ! Net rate of change of population density within each mass + ! cohort. (m-2 s-1) +flux_in(mclass), & + ! 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. + ! (m-2 s-1) +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 + + ! 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 + ! nowhere lower to shrink into, so its share of the loss is instead matched + ! by an equivalent mortality rate. + mort(k) = mort_base + mort_add(k) + + IF (k == 1) THEN + mort(k) = mort(k) - g_mass(k) / mass_mass(k) + END IF + + IF (k < mclass) THEN + ! Flux shrinking down into this class from the class above + flux_in(k) = - plantNumDensity(k+1) * g_mass(k+1) & + / (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) * g_mass(k) & + / (mass_mass(k) - mass_mass(k-1)) + ELSE + flux_out(k) = 0.0 + END IF + + ELSE + 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