From 7228c5202014a972b2f61bb55a220f6394726a80 Mon Sep 17 00:00:00 2001 From: marcadella Date: Fri, 10 Jan 2025 09:27:54 +0100 Subject: [PATCH] Refactoring: removed static variables in gpp --- src/biosphere_biomee.mod.f90 | 4 +- src/datatypes_biomee.mod.f90 | 17 ++++++++- src/gpp_biomee.mod.f90 | 72 ++++++++++++++++------------------- src/sofunutils.mod.f90 | 2 +- src/vegetation_biomee.mod.f90 | 5 +-- 5 files changed, 52 insertions(+), 48 deletions(-) diff --git a/src/biosphere_biomee.mod.f90 b/src/biosphere_biomee.mod.f90 index e2245fb8..79a2d24e 100644 --- a/src/biosphere_biomee.mod.f90 +++ b/src/biosphere_biomee.mod.f90 @@ -50,7 +50,6 @@ subroutine biosphere_annual( & integer :: dayloop_idx, fastloop_idx, simu_steps real, dimension(ndayyear) :: daily_temp ! Daily temperatures (average) real :: tsoil - logical :: first_simu_step !---------------------------------------------------------------- ! INITIALISATIONS @@ -105,7 +104,6 @@ subroutine biosphere_annual( & fastloop: do fastloop_idx = 1,myinterface%steps_per_day simu_steps = simu_steps + 1 - first_simu_step = state%init .and. (simu_steps == 1) vegn%thetaS = (vegn%wcl(2) - myinterface%params_tile%WILTPT) & / (myinterface%params_tile%FLDCAP - myinterface%params_tile%WILTPT) @@ -113,7 +111,7 @@ subroutine biosphere_annual( & !---------------------------------------------------------------- ! Sub-daily time step at resolution given by forcing (can be 1 = daily) !---------------------------------------------------------------- - call vegn_CNW_budget( vegn, myinterface%climate(simu_steps), first_simu_step, tsoil ) + call vegn_CNW_budget( vegn, myinterface%climate(simu_steps), tsoil ) call hourly_diagnostics( vegn, myinterface%climate(simu_steps) ) diff --git a/src/datatypes_biomee.mod.f90 b/src/datatypes_biomee.mod.f90 index 7a27517e..5d849c0d 100644 --- a/src/datatypes_biomee.mod.f90 +++ b/src/datatypes_biomee.mod.f90 @@ -12,7 +12,7 @@ module datatypes_biomee ! define data types and constants implicit none !=============== Public types =========================================================== - public :: spec_data_type, cohort_type, vegn_tile_type + public :: spec_data_type, cohort_type, vegn_tile_type, dampended_forcing_type !=============== Public subroutines ===================================================== public :: Zero_diagnostics, hourly_diagnostics, daily_diagnostics, & @@ -40,6 +40,9 @@ module datatypes_biomee integer, public, parameter :: nvars_params_species = 55 integer, public, parameter :: nvars_init_lu = 1 + !===== Model + integer, public, parameter :: NLAYERS_MAX = 9 ! maximum number of canopy layers to be considered + !===== Photosynthesis real, public, parameter :: extinct = 0.75 ! light extinction coefficient in the canopy for photosynthesis @@ -60,6 +63,15 @@ module datatypes_biomee !===== Leaf life span real, parameter :: c_LLS = 28.57143 ! yr/ (kg C m-2), c_LLS=1/LMAs, where LMAs = 0.035 + type :: dampended_forcing_type + logical :: initialized = .False. + real :: co2 = 0.0 + real :: vpd = 0.0 + real :: temp = 0.0 + real :: patm = 0.0 + real, dimension(NLAYERS_MAX) :: par = dummy ! Initialization to a dummy value. Important! Keep it. + end type dampended_forcing_type + !=============== Tile level data type ============================================================ type :: vegn_tile_type @@ -172,6 +184,9 @@ module datatypes_biomee real :: totNewCC = 0.0 ! New cohort C (all compartments but seed), kg C yr-1 m-2 real :: totNewCN = 0.0 ! New cohort N (all compartments but seed), kg N yr-1 m-2 + ! Scrap variable used by gpp() + type(dampended_forcing_type) :: dampended_forcing + end type vegn_tile_type contains diff --git a/src/gpp_biomee.mod.f90 b/src/gpp_biomee.mod.f90 index 38dab4cb..f69fd2e2 100644 --- a/src/gpp_biomee.mod.f90 +++ b/src/gpp_biomee.mod.f90 @@ -37,7 +37,7 @@ module md_gpp_biomee contains - subroutine gpp( forcing, vegn, first_simu_step ) + subroutine gpp( forcing, vegn ) !////////////////////////////////////////////////////////////////////// ! GPP ! Calculates light availability and photosynthesis for each cohort @@ -58,11 +58,9 @@ subroutine gpp( forcing, vegn, first_simu_step ) type(climate_type), intent(in):: forcing type(vegn_tile_type), intent(inout) :: vegn - logical, intent(in) :: first_simu_step ! is true on the very first simulation step ! local variables used for BiomeE-Allocation part type(cohort_type), pointer :: cc - integer, parameter :: nlayers_max = 9 ! maximum number of canopy layers to be considered real :: rad_top ! downward radiation at the top of the canopy, W/m2 real :: rad_net ! net radiation absorbed by the canopy, W/m2 real :: Tair, TairK ! air temperature, degC and degK @@ -82,15 +80,10 @@ subroutine gpp( forcing, vegn, first_simu_step ) real, dimension(vegn%n_cohorts) :: fapar_tree ! tree-level fAPAR based on LAI within the crown real, dimension(nlayers_max-1) :: fapar_layer - integer:: i, layer=0 + integer:: i ! local variables used for P-model part real :: tk, kphio_temp - real, save :: co2_memory - real, save :: vpd_memory - real, save :: temp_memory - real, save :: patm_memory - real, dimension(nlayers_max), save :: par_memory type(outtype_pmodel) :: out_pmodel ! list of P-model output variables !----------------------------------------------------------- @@ -109,10 +102,9 @@ subroutine gpp( forcing, vegn, first_simu_step ) fapar_layer(:) = 0.0 do i = 1, vegn%n_cohorts cc => vegn%cohorts(i) - layer = Max(1, Min(cc%layer, nlayers_max)) - LAIlayer(layer) = LAIlayer(layer) + cc%leafarea * cc%nindivs / (1.0 - f_gap) + LAIlayer(cc%layer) = LAIlayer(cc%layer) + cc%leafarea * cc%nindivs / (1.0 - f_gap) fapar_tree(i) = 1.0 - exp(-kappa * cc%leafarea / cc%crownarea) ! at individual-level: cc%leafarea represents leaf area index within the crown - fapar_layer(layer) = fapar_layer(layer) + fapar_tree(i) * cc%crownarea * cc%nindivs + fapar_layer(cc%layer) = fapar_layer(cc%layer) + fapar_tree(i) * cc%crownarea * cc%nindivs enddo ! Get light received at each crown layer as a fraction of top-of-canopy -> f_light(layer) @@ -143,13 +135,12 @@ subroutine gpp( forcing, vegn, first_simu_step ) if (cc%status == LEAF_ON) then !.and. cc%lai > 0.1 ! Convert forcing data - layer = Max (1, Min(cc%layer, nlayers_max)) - rad_top = f_light(layer) * forcing%radiation ! downward radiation at the top of the canopy, W/m2 + rad_top = f_light(cc%layer) * forcing%radiation ! downward radiation at the top of the canopy, W/m2 !=============================== ! ORIGINAL !=============================== - rad_net = f_light(layer) * forcing%radiation * 0.9 ! net radiation absorbed by the canopy, W/m2 + rad_net = f_light(cc%layer) * forcing%radiation * 0.9 ! net radiation absorbed by the canopy, W/m2 p_surf = forcing%P_air ! Pa TairK = forcing%Tair ! K Tair = forcing%Tair - 273.16 ! degC @@ -194,19 +185,22 @@ subroutine gpp( forcing, vegn, first_simu_step ) ! Calculate environmental conditions with memory, time scale ! relevant for Rubisco turnover !---------------------------------------------------------------- - if (first_simu_step) then - co2_memory = forcing%CO2 * 1.0e6 - temp_memory = (forcing%Tair - kTkelvin) - vpd_memory = forcing%vpd - patm_memory = forcing%P_air - par_memory = -1.0 ! We initialize par_memory to a dummy value used to detect the need for initialization, - ! as the initialization process using init flags is error prone given the dunmaic the adjuction of layers. - end if - - co2_memory = dampen_variability( forcing%CO2 * 1.0e6, params_gpp%tau_acclim, co2_memory ) - temp_memory = dampen_variability( (forcing%Tair - kTkelvin), params_gpp%tau_acclim, temp_memory) - vpd_memory = dampen_variability( forcing%vpd, params_gpp%tau_acclim, vpd_memory ) - patm_memory = dampen_variability( forcing%P_air, params_gpp%tau_acclim, patm_memory ) + if (.not. vegn%dampended_forcing%initialized) then + vegn%dampended_forcing%co2 = forcing%CO2 * 1.0e6 + vegn%dampended_forcing%temp = (forcing%Tair - kTkelvin) + vegn%dampended_forcing%vpd = forcing%vpd + vegn%dampended_forcing%patm = forcing%P_air + vegn%dampended_forcing%initialized = .True. + else + vegn%dampended_forcing%co2 = dampen_variability( forcing%CO2 * 1.0e6, params_gpp%tau_acclim, & + vegn%dampended_forcing%co2 ) + vegn%dampended_forcing%temp = dampen_variability( (forcing%Tair - kTkelvin), params_gpp%tau_acclim, & + vegn%dampended_forcing%temp) + vegn%dampended_forcing%vpd = dampen_variability( forcing%vpd, & + params_gpp%tau_acclim, vegn%dampended_forcing%vpd ) + vegn%dampended_forcing%patm = dampen_variability( forcing%P_air, & + params_gpp%tau_acclim, vegn%dampended_forcing%patm ) + end if tk = forcing%Tair + kTkelvin @@ -218,9 +212,6 @@ subroutine gpp( forcing, vegn, first_simu_step ) cc => vegn%cohorts(i) associate ( sp => myinterface%params_species(cc%species) ) - ! get canopy layer of this cohort - layer = max(1, min(cc%layer, nlayers_max)) - !---------------------------------------------------------------- ! Instantaneous temperature effect on quantum yield efficiency !---------------------------------------------------------------- @@ -231,12 +222,13 @@ subroutine gpp( forcing, vegn, first_simu_step ) params_gpp%kphio_par_b ) ! photosynthetically active radiation level at this layer - par = f_light(layer) * forcing%radiation * kfFEC * 1.0e-6 + par = f_light(cc%layer) * forcing%radiation * kfFEC * 1.0e-6 ! slowly varying light conditions per layer, relevant for acclimation (P-model quantities) - if (par_memory(layer) <= -1.0) then - par_memory(layer) = par + if (vegn%dampended_forcing%par(cc%layer) == dummy) then + vegn%dampended_forcing%par(cc%layer) = par else - par_memory(layer) = dampen_variability(par, params_gpp%tau_acclim, par_memory(layer)) + vegn%dampended_forcing%par(cc%layer) = dampen_variability(par, params_gpp%tau_acclim, & + vegn%dampended_forcing%par(cc%layer)) end if if (cc%status == LEAF_ON) then @@ -249,11 +241,11 @@ subroutine gpp( forcing, vegn, first_simu_step ) kphio = kphio_temp, & beta = params_gpp%beta, & kc_jmax = params_gpp%kc_jmax, & - ppfd = par_memory(layer), & - co2 = co2_memory, & - tc = temp_memory, & - vpd = vpd_memory, & - patm = patm_memory, & + ppfd = vegn%dampended_forcing%par(cc%layer), & + co2 = vegn%dampended_forcing%co2, & + tc = vegn%dampended_forcing%temp, & + vpd = vegn%dampended_forcing%vpd, & + patm = vegn%dampended_forcing%patm, & c4 = .false., & method_optci = "prentice14", & method_jmaxlim = "wang17" & diff --git a/src/sofunutils.mod.f90 b/src/sofunutils.mod.f90 index 9fb73044..c1ef470b 100644 --- a/src/sofunutils.mod.f90 +++ b/src/sofunutils.mod.f90 @@ -13,7 +13,7 @@ function dampen_variability( var, tau, var_memory ) result( out_memory ) ! Calculates the updated variable accounting for a memory time scale tau. ! Following Eq. 5 in Makela et al. (2004) Tree Physiology 24, 369–376 ! - ! d(var_memory) / dt = (1 / tau) * var - var_memory + ! d(var_memory) / dt = (1 / tau) * (var - var_memory) ! !------------------------------------------------------------------------- ! arguments diff --git a/src/vegetation_biomee.mod.f90 b/src/vegetation_biomee.mod.f90 index b0739bcf..f6f788c1 100755 --- a/src/vegetation_biomee.mod.f90 +++ b/src/vegetation_biomee.mod.f90 @@ -25,7 +25,7 @@ module md_vegetation_biomee !============= Carbon, nitrogen and water budget ===================== !======================================================================== - subroutine vegn_CNW_budget( vegn, forcing, first_simu_step, tsoil ) + subroutine vegn_CNW_budget( vegn, forcing, tsoil ) !//////////////////////////////////////////////////////////////// ! Fast loop carbon, nitrogen, and water dynamics, Weng 2016-11-25 ! include Nitrogen uptake and carbon budget @@ -37,7 +37,6 @@ subroutine vegn_CNW_budget( vegn, forcing, first_simu_step, tsoil ) type(vegn_tile_type), intent(inout) :: vegn type(climate_type), intent(in) :: forcing ! is true on the very first simulation step (first subroutine call of each gridcell) - logical, intent(in) :: first_simu_step real, intent(in) :: tsoil ! Soil temperature in K ! local variables @@ -45,7 +44,7 @@ subroutine vegn_CNW_budget( vegn, forcing, first_simu_step, tsoil ) integer:: i ! Photosynsthesis - call gpp( forcing, vegn, first_simu_step ) + call gpp( forcing, vegn ) ! Update soil water call SoilWaterDynamicsLayer( forcing, vegn )