feat: expose initialization of soil (N, microb, water) pools to R (optional)

Previously, default initial values were derived inside of Fortran code. This has now been moved to R in case the user does not provide initial values for the soil pools.

Bottom line: The default recovers former default initialization.

Author: fabern

NEWS.md

@@ -17,10 +17,14 @@
   * Added optional parameter `init_cohort$init_cohort_age` to specify a corresponding 
   initial age distribution in years (default = 0). Default ensures backwards compatibility.
   Age affects reproduction through `params_species$matural_age`
-  * Added optional parameter to `init_cohort`, to enable specification of initial
+  * Added optional parameter to `init_cohort`, to enable specification of initial plant
   nitrogen pools. Namely: `init_cohort_bl_n14`, `init_cohort_br_n14`,
   `init_cohort_bsw_n14`,`init_cohort_bHW_n14`,`init_cohort_seedC_n14`,
   `init_cohort_nsc_n14``. If not provided, default values ensure backwards compatibility.
+  * Added optional parameter to `init_soil`, to enable specification of initial soil
+  nitrogen pools. Namely: `init_fast_soil_N`, `init_slow_soil_N`, `init_pmicr_C`, 
+  `init_pmicr_d13C`, `init_pmicr_N`, `init_wcl1`, `init_wcl2`, `init_wcl3`, `init_N0_ecosystem`. 
+  If not provided, default values ensure backwards compatibility.
   * Added check that all `species$LMA` >= `params_tile$LMAmin`
 * P-model:
   * no changes

R/data.R

@@ -316,12 +316,12 @@
 #'       \item{init_cohort_bHW}{Initial biomass of heartwood, in kg C per individual.}
 #'       \item{init_cohort_seedC}{Initial biomass of seed, in kg C per individual.}
 #'       \item{init_cohort_nsc}{Initial non-structural biomass, in kg C per individual. (optional, defaults to value derived from br_max and bl_max)}
-#'       \item{init_cohort_bl_n14}{Initial nitorgen in leaf, in kg N per individual.        (optional, defaults to C biomass * CN-ratio)}
-#'       \item{init_cohort_br_n14}{Initial nitorgen in fine root, in Ng C per individual.   (optional, defaults to C biomass * CN-ratio)}
-#'       \item{init_cohort_bsw_n14}{Initial nitorgen in sapwood, in kg N per individual.    (optional, defaults to C biomass * CN-ratio)}
-#'       \item{init_cohort_bHW_n14}{Initial nitorgen in heartwood, in kg N per individual.  (optional, defaults to C biomass * CN-ratio)}
-#'       \item{init_cohort_seedC_n14}{Initial nitorgen in seed, in kg N per individual.     (optional, defaults to C biomass * CN-ratio)}
-#'       \item{init_cohort_nsc_n14}{Initial non-structural biomass, in kg C per individual. (optional, defaults to value derived from br_max and bl_max)}
+#'       \item{init_cohort_bl_n14}{Initial nitrogen in leaf, in kg N per individual.        (optional, defaults to C biomass * CN-ratio)}
+#'       \item{init_cohort_br_n14}{Initial nitrogen in fine root, in kg N per individual.   (optional, defaults to C biomass * CN-ratio)}
+#'       \item{init_cohort_bsw_n14}{Initial nitrogen in sapwood, in kg N per individual.    (optional, defaults to C biomass * CN-ratio)}
+#'       \item{init_cohort_bHW_n14}{Initial nitrogen in heartwood, in kg N per individual.  (optional, defaults to C biomass * CN-ratio)}
+#'       \item{init_cohort_seedC_n14}{Initial nitrogen in seed, in kg N per individual.     (optional, defaults to C biomass * CN-ratio)}
+#'       \item{init_cohort_nsc_n14}{Initial non-structural biomass, in kg N per individual. (optional, defaults to value derived from br_max and bl_max)}
 #'       \item{lu_index}{Land use type this cohorts belongs to (given as index in init_lu aray). (optional)
 #'        Default: 0 (attach to all LU types except thoses which do not accept vegetation -- cf  init_lu.vegetated).}
 #'     }}
@@ -332,6 +332,15 @@
 #'       \item{init_slow_soil_C}{Initial slow soil carbon, in kg C m\eqn{^{-2}}.}
 #'       \item{init_Nmineral}{Mineral nitrogen pool, in kg N m\eqn{^{-2}}.}
 #'       \item{N_input}{Annual nitrogen input to soil N pool, in kg N m\eqn{^{-2}} yr\eqn{^{-1}}.}
+#'       \item{init_fast_soil_N}{Initial fast soil nitrogen, in kg N m\eqn{^{-2}}.                     (optional, defaults to C biomass * CN-ratio)}
+#'       \item{init_slow_soil_N}{Initial slow soil nitrogen, in kg N m\eqn{^{-2}}.                     (optional, defaults to C biomass * CN-ratio)}
+#'       \item{init_pmicr_C}{Initial carbon in soil microbial biomass, in kg C m\eqn{^{-2}}.           (optional, defaults to 0.0)}
+#'       \item{init_pmicr_d13C}{Initial d13C signature in soil microbial biomass, in permil.           (optional, defaults to -9999)}
+#'       \item{init_pmicr_N}{Initial nitrogen in soil microbial biomass, in kg N m\eqn{^{-2}}.         (optional, defaults to 0.0)}
+#'       \item{init_wcl1}{Initial volumetric water content of soil layer 1, in m\eqn{^{3}m\eqn{^{-3}}. (optional, defaults to field capacity)}
+#'       \item{init_wcl2}{Initial volumetric water content of soil layer 2, in m\eqn{^{3}m\eqn{^{-3}}. (optional, defaults to field capacity)}
+#'       \item{init_wcl3}{Initial volumetric water content of soil layer 3, in m\eqn{^{3}m\eqn{^{-3}}. (optional, defaults to field capacity)}
+#'       \item{init_N0_ecosystem}{Initial total amount of nitrogen in ecosystem (only used for nitrogen workaround), in kg N m\eqn{^{-2}}. (optional, defaults to sum of soil and plant pools.) Might not be needed.}
 #'     }}
 #'   \item{init_lu}{A data.frame of initial land unit (LU) specifications, including
 #'     the following data:

R/run_biomee_f_bysite.R

@@ -246,7 +246,7 @@ run_biomee_f_bysite <- function(
   params_tile <- build_params_tile(params_tile)
   params_species <- build_params_species(params_species, params_tile)
   init_cohort    <- build_init_cohort(init_cohort, params_species)
-  init_soil <- build_init_soil(init_soil)
+  init_soil <- build_init_soil(init_soil, init_cohort, params_tile)
   forcing <- build_forcing(forcing)
 
   # Build LULUC parameters
@@ -658,8 +658,41 @@ build_init_cohort <- function(init_cohort, params_species){
   return(init_cohort)
 }
 
-build_init_soil <- function(init_soil){
-  #browser() # currently build_init_soil has no-effect
+build_init_soil <- function(init_soil, init_cohort, params_tile){
+  if ('init_fast_soil_N' %nin% names(init_soil)) {
+    init_soil$init_fast_soil_N = init_soil$init_fast_soil_C / params_tile$CN0metabolicL # former default
+  }
+  if ('init_slow_soil_N' %nin% names(init_soil)) {
+    init_soil$init_slow_soil_N = init_soil$init_slow_soil_C / params_tile$CN0structuralL # former default
+  }
+  if ('init_pmicr_C' %nin% names(init_soil)) {
+    init_soil$init_pmicr_C = 0.0 # former default
+  }
+  if ('init_pmicr_d13C' %nin% names(init_soil)) {
+    init_soil$init_pmicr_d13C = -9999.0 # former default
+  }
+  if ('init_pmicr_N' %nin% names(init_soil)) {
+    init_soil$init_pmicr_N = 0.0 # former default
+  }
+  if ('init_wcl1' %nin% names(init_soil)) {
+    init_soil$init_wcl1 = params_tile$FLDCAP # former default
+  }
+  if ('init_wcl2' %nin% names(init_soil)) {
+    init_soil$init_wcl2 = params_tile$FLDCAP # former default
+  }
+  if ('init_wcl3' %nin% names(init_soil)) {
+    init_soil$init_wcl3 = params_tile$FLDCAP # former default
+  }
+  if ('init_N0_ecosystem' %nin% names(init_soil)) { # this is used for nitrogen workaround
+    Ntot_soil <- init_soil$init_pmicr_N + init_soil$init_fast_soil_N + init_soil$init_slow_soil_N + init_soil$init_Nmineral # kgN/m2
+    N_in_each_cohort <- with(init_cohort,
+         init_cohort_nindivs * # tree/m2
+           (init_cohort_nsc_n14 + init_cohort_bl_n14 + init_cohort_br_n14 +     # kgN per tree
+            init_cohort_bsw_n14 + init_cohort_bHW_n14 + init_cohort_seedC_n14))
+    Ntot_plant <- sum(N_in_each_cohort)
+    Ntot <- Ntot_soil + Ntot_plant
+    init_soil$init_N0_ecosystem = Ntot # former default: sum of the initialized soil and plant pools
+  }
   return(init_soil)
 }
 
@@ -917,7 +950,16 @@ prepare_init_soil <- function(init_soil){
     "init_fast_soil_C",
     "init_slow_soil_C",
     "init_Nmineral",
-    "N_input"
+    "N_input", 
+    "init_fast_soil_N", 
+    "init_slow_soil_N", 
+    "init_pmicr_C", 
+    "init_pmicr_d13C", 
+    "init_pmicr_N", 
+    "init_wcl1", 
+    "init_wcl2", 
+    "init_wcl3",
+    "init_N0_ecosystem"
   )
 }
 

src/interface_in_biosphere_biomee.mod.f90

@@ -11,7 +11,7 @@ module md_interface_in_biomee
   private
   public  inputs, interface_in_biosphere_biomee, params_species_biomee, init_lu_biomee
 
-  !===== Soil water hydrualics
+  !===== Soil water hydraulics
   integer, public, parameter :: MAX_LEVELS = 3  ! Soil layers, for soil water dynamics
   real, public, parameter ::  thksl(MAX_LEVELS) = (/0.05, 0.45, 1.5/)  ! m, thickness of soil layers
 
@@ -22,7 +22,7 @@ module md_interface_in_biomee
   integer, public, parameter :: nvars_params_siml    = 11
   integer, public, parameter :: nvars_site_info      = 4
   integer, public, parameter :: nvars_params_tile    = 20
-  integer, public, parameter :: nvars_init_soil      = 4
+  integer, public, parameter :: nvars_init_soil      = 13
   integer, public, parameter :: nvars_init_cohorts   = 16
   integer, public, parameter :: nvars_params_species = 65
   integer, public, parameter :: nvars_init_lu        = 5
@@ -215,6 +215,15 @@ module md_interface_in_biomee
     real :: init_slow_soil_C
     real :: init_Nmineral
     real :: N_input
+    real :: init_fast_soil_N
+    real :: init_slow_soil_N
+    real :: init_pmicr_C
+    real :: init_pmicr_d13C
+    real :: init_pmicr_N
+    real :: init_wcl1
+    real :: init_wcl2
+    real :: init_wcl3
+    real :: init_N0_ecosystem
 
   contains
 
@@ -405,6 +414,16 @@ subroutine populate_init_soil(self, init_soil)
     self%init_slow_soil_C         = real( init_soil(2) )
     self%init_Nmineral            = real( init_soil(3) )
     self%N_input                  = real( init_soil(4) )
+    self%init_fast_soil_N         = real( init_soil(5) )
+    self%init_slow_soil_N         = real( init_soil(6) )
+    self%init_pmicr_C             = real( init_soil(7) )
+    self%init_pmicr_d13C          = real( init_soil(8) )
+    self%init_pmicr_N             = real( init_soil(9) )
+    self%init_wcl1                = real( init_soil(10))
+    self%init_wcl2                = real( init_soil(11))
+    self%init_wcl3                = real( init_soil(12))
+    self%init_N0_ecosystem        = real( init_soil(13))
+    
   end subroutine populate_init_soil
 
   subroutine populate_params_tile(self, params_tile)

src/vegetation_tile_biomee.mod.f90

@@ -119,7 +119,6 @@ module md_vegetation_tile_biomee
     real    :: N_P2S_yr           ! N turnover (plant to soil) (kg N m-2 yr-1)
 
     !=====  Memory
-    real    :: previousN          ! weighted annual available N
     real    :: initialN0          ! initial available N (kg N m-2)
 
     !=====  Fast fluxes, kg m-2 timestep-1
@@ -225,7 +224,7 @@ pure real function soilwater(self)
     !---------------------------------------------------------------
     class(vegn_tile_type), intent(in) :: self
 
-    soilwater  = SUM(self%wcl(:)*thksl(:)*1000.0)
+    soilwater  = SUM(self%wcl(:)*thksl(:)*1000.0) ! in m*1000 = m3/m2 *1000 = kg/m2
 
   end function soilwater
 
@@ -1079,20 +1078,23 @@ subroutine initialize_vegn_tile( self, lu_index )
     ! Initial Soil pools and environmental conditions
     self%psoil_fs%c12 = inputs%init_soil%init_fast_soil_C  ! fast soil carbon pool, (kg C/m2)
     self%psoil_sl%c12 = inputs%init_soil%init_slow_soil_C  ! slow soil carbon pool, (kg C/m2)
-    self%psoil_fs%n14 = self%psoil_fs%c12 / inputs%params_tile%CN0metabolicL  ! fast soil nitrogen pool, (kg N/m2)
-    self%psoil_sl%n14 = self%psoil_sl%c12 / inputs%params_tile%CN0structuralL ! slow soil nitrogen pool, (kg N/m2)
+    self%psoil_fs%n14 = inputs%init_soil%init_fast_soil_N ! fast soil nitrogen pool, (kg N/m2)
+    self%psoil_sl%n14 = inputs%init_soil%init_slow_soil_N ! slow soil nitrogen pool, (kg N/m2)
     self%inorg%n14    = inputs%init_soil%init_Nmineral     ! Mineral nitrogen pool, (kg N/m2)
-    self%previousN    = self%inorg%n14
 
     ! debug: adding microbial biomass initialisation
-    self%pmicr = orgpool() ! to do: add to: inputs%init_soil%xxxxx
+    self%pmicr = orgpool(inputs%init_soil%init_pmicr_C,   &
+                         inputs%init_soil%init_pmicr_d13C, &
+                         inputs%init_soil%init_pmicr_N)
 
-    ! Initialize soil volumetric water conent with field capacity (maximum soil moisture to start with)
-    self%wcl = inputs%params_tile%FLDCAP
+    ! Initialize soil volumetric water content
+    self%wcl = (/inputs%init_soil%init_wcl1, inputs%init_soil%init_wcl2, inputs%init_soil%init_wcl3/)
 
-    call self%aggregate_cohorts()
+    ! Initialize initialN0, that is used for nitrogen workaround: keep the N in the system constant at this value
+    self%initialN0 = inputs%init_soil%init_N0_ecosystem
 
-    self%initialN0 =  self%totN
+    !call self%aggregate_cohorts()
+    !self%initialN0 =  self%totN
 
   end subroutine initialize_vegn_tile