geco-bern · fabern · Feb 16, 2026 · Feb 16, 2026 · Feb 17, 2026 · Feb 17, 2026
diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml
@@ -35,11 +35,11 @@ jobs:
 
       - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          extra-packages: any::pkgdown, any::devtools, local::.
+          extra-packages: any::pkgdown, any::pak, local::.
           needs: website
 
       - name: Build site
-        run: devtools::install_github("geco-bern/GECO_documentation_template"); pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE)
+        run: pak::pak("geco-bern/GECO_documentation_template"); pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE)
         shell: Rscript {0}
 
       - name: Deploy to GitHub pages 🚀

diff --git a/CITATION.cff b/CITATION.cff
@@ -70,7 +70,7 @@ abstract: '<h1>rsofun</h1>
   <p>To install the latest development release of the package run the following commands
   to install rsofun directly from GitHub:</p>
 
-  <p><code>if(!require(remotes)){install.packages("remotes")}</code><br><code>remotes::install_github("geco-bern/rsofun")</code><br><code>library("rsofun")</code><br><br></p>
+  <p><code>if(!require(pak)){install.packages("pak")}</code><br><code>pak::pak("geco-bern/rsofun")</code><br><code>library("rsofun")</code><br><br></p>
 
   <p>**NOTE:** Installing from GitHub requires compilation of Fortran and C source
   code contained in {rsofun}. To enable compiling source code, install <a href="https://cran.r-project.org/bin/windows/Rtools/">Rtools</a>
@@ -79,12 +79,6 @@ abstract: '<h1>rsofun</h1>
   (see also ''Mandatory tools'' <a href="https://mac.r-project.org/tools/">here</a>).
   On Linux, the gfortran compiler is usually installed already.</p>
 
-  <p>Vignettes are not rendered by default, if you want to include additional documentation
-  please use:</p>
-
-  <p><code>if(!require(remotes)){install.packages("remotes")}</code><br><code>remotes::install_github("geco-bern/rsofun",
-  build_vignettes = TRUE)</code><br><code>library("rsofun")</code></p>
-
   <h3>From source&nbsp;</h3>
 
   <p>Assuming <code>rsofun</code> is the location of the source directory, on can

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -74,7 +74,7 @@ LazyData: true
 LazyDataCompression: xz
 ByteCompile: true
 NeedsCompilation: yes
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Suggests:
     covr,
     constructive,

diff --git a/NEWS.md b/NEWS.md
@@ -1,4 +1,9 @@
 # rsofun (development version)
+* BiomeEP:
+  * Bugfix: annual, cohort-level output had mixed up column names for variables:
+  `NSN`,`seedC`,`leafC`,`rootC`,`sapwoodC`,`heartwoodC`
+  * Cohorts are now less aggressively merged. Merging criteria was simplified from
+  relative to absolute DBH difference. Now merging if difference <= 0.01 m.
 
 # rsofun 5.1.0
 

diff --git a/R/data.R b/R/data.R
@@ -277,7 +277,7 @@
 #'       \item{thetaHT}{Coefficient for allometry (height = alphaHT * DBH_m ** thetaHT), in m m\eqn{^{-thetaHT}}.}
 #'       \item{alphaCA}{Coefficient for allometry (projected crown area = pi * (alphaCA * DBH_m) ** thetaCA), in m\eqn{^{2/thetaCA-1}}.}
 #'       \item{thetaCA}{Coefficient for allometry (projected crown area = pi * (alphaCA * DBH_m) ** thetaCA), unitless. Dybzinski (eq. G1) showed that thetaCA = theatBM - 1.}
-#'       \item{alphaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM), in kg C m\eqn{^{-thetaBM}}.}
+#'       \item{alphaBM}{TODO: unused. alphaBM is overwritten by: (alphaBM = rho_wood * taperfactor * pi/4.*alphaHT). Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM), in kg C m\eqn{^{-thetaBM}}.}
 #'       \item{thetaBM}{Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM), unitless. Dybzinski (eq. G1) showed that thetaCA = theatBM - 1.}
 #'     }}\item{The following columns pertaining to the \bold{reproduction parameters}:}{\describe{
 #'       \item{seedlingsize}{Initial size of seedlings, in kg C per individual.}

diff --git a/R/run_biomee_f_bysite.R b/R/run_biomee_f_bysite.R
@@ -518,32 +518,32 @@ prepare_init_lu <- function(init_lu){
     init_lu <- init_lu %>% mutate(preset = 'unmanaged')
   }
   if(!'extra_N_input' %in% names(init_lu)) {
-    init_lu <- init_lu %>% mutate('extra_N_input' = case_match(
+    init_lu <- init_lu %>% mutate('extra_N_input' = recode_values(
       'preset',
       "cropland" ~ 0.01,
-      .default = 0.0
+      default = 0.0
     ))
   }
   if(!'extra_turnover_rate' %in% names(init_lu)) {
-    init_lu <- init_lu %>% mutate('extra_turnover_rate' = case_match(
+    init_lu <- init_lu %>% mutate('extra_turnover_rate' = recode_values(
       'preset',
       "cropland" ~ 0.2,
-      .default = 0.0
+      default = 0.0
     ))
   }
   if(!'oxidized_litter_fraction' %in% names(init_lu)) {
-    init_lu <- init_lu %>% mutate('oxidized_litter_fraction' = case_match(
+    init_lu <- init_lu %>% mutate('oxidized_litter_fraction' = recode_values(
       'preset',
       "cropland" ~ 0.9,
       "pasture" ~ 0.4,
-      .default = 0.0
+      default = 0.0
     ))
   }
   init_lu <- init_lu %>% mutate(
-    'vegetated' = case_match(
+    'vegetated' = recode_values(
       'preset',
       "urban" ~ FALSE,
-      .default = TRUE
+      default = TRUE
     )
   )
   init_lu <- init_lu %>% select(
@@ -846,41 +846,41 @@ annual_tile_output <- function(raw_data, aggregated_LU = FALSE){
 
 annual_cohort_output <- function(raw_data){
   annual_values <- c(
-    "cohort",
-    "year",
-    "cID",
-    "PFT",
-    "layer",
-    "density",
-    "flayer",
-    "DBH",
-    "dDBH",
-    "height",
-    "age",
-    "BA",
-    "dBA",
-    "Acrown",
-    "Aleaf",
-    "NSC",
-    "seedC",
-    "leafC",
-    "rootC",
-    "sapwoodC",
-    "heartwoodC",
-    "NSN",
-    "treeG",
-    "fseed",
-    "fleaf",
-    "froot",
-    "fwood",
-    "NPP",
-    "GPP",
-    "Rauto",
-    "N_uptk",
-    "N_fxed",
-    "deathrate",
-    "n_deadtrees",
-    "c_deadtrees"
+    "cohort",      # ANNUAL_COHORTS_ID             =  1
+    "year",        # ANNUAL_COHORTS_YEAR           =  2
+    "cID",         # ANNUAL_COHORTS_CID            =  3
+    "PFT",         # ANNUAL_COHORTS_PFT            =  4
+    "layer",       # ANNUAL_COHORTS_LAYER          =  5
+    "density",     # ANNUAL_COHORTS_DENSITY        =  6
+    "flayer",      # ANNUAL_COHORTS_FLAYER         =  7
+    "DBH",         # ANNUAL_COHORTS_DBH            =  8
+    "dDBH",        # ANNUAL_COHORTS_DDBH           =  9
+    "height",      # ANNUAL_COHORTS_HEIGHT         = 10
+    "age",         # ANNUAL_COHORTS_AGE            = 11
+    "BA",          # ANNUAL_COHORTS_BA             = 12
+    "dBA",         # ANNUAL_COHORTS_DBA            = 13
+    "Acrown",      # ANNUAL_COHORTS_ACROWN         = 14
+    "Aleaf",       # ANNUAL_COHORTS_ALEAF          = 15
+    "NSC",         # ANNUAL_COHORTS_NCS            = 16
+    "NSN",         # ANNUAL_COHORTS_NSN            = 17
+    "seedC",       # ANNUAL_COHORTS_SEED_C         = 18
+    "leafC",       # ANNUAL_COHORTS_LEAF_C         = 19
+    "rootC",       # ANNUAL_COHORTS_ROOT_C         = 20
+    "sapwoodC",    # ANNUAL_COHORTS_SW_C           = 21
+    "heartwoodC",  # ANNUAL_COHORTS_HW_C           = 22
+    "treeG",       # ANNUAL_COHORTS_TREEG          = 23
+    "fseed",       # ANNUAL_COHORTS_FSEED          = 24
+    "fleaf",       # ANNUAL_COHORTS_FLEAF          = 25
+    "froot",       # ANNUAL_COHORTS_FROOT          = 26
+    "fwood",       # ANNUAL_COHORTS_FWOOD          = 27
+    "NPP",         # ANNUAL_COHORTS_NPP            = 28
+    "GPP",         # ANNUAL_COHORTS_GPP            = 29
+    "Rauto",       # ANNUAL_COHORTS_RESP           = 30
+    "N_uptk",      # ANNUAL_COHORTS_N_UPTK         = 31
+    "N_fxed",      # ANNUAL_COHORTS_N_FIX          = 32
+    "deathrate",   # ANNUAL_COHORTS_DEATHRATE      = 33
+    "n_deadtrees", # ANNUAL_COHORTS_N_LOSS         = 34
+    "c_deadtrees"  # ANNUAL_COHORTS_C_LOSS         = 35
   )
 
   dimensions <- dim(raw_data)

diff --git a/README.md b/README.md
@@ -26,21 +26,13 @@ library("rsofun")
 To install the latest development release of the package run the following commands to install rsofun directly from GitHub:
 
 ``` r
-if(!require(remotes)){install.packages("remotes")}
-remotes::install_github("geco-bern/rsofun")
+if(!require(pak)){install.packages("pak")}
+pak::pak("geco-bern/rsofun")
 library("rsofun")
 ```
 
 **NOTE:** Installing from GitHub requires compilation of Fortran and C source code contained in {rsofun}. To enable compiling source code, install [Rtools](https://cran.r-project.org/bin/windows/Rtools/) on Windows, or [Xcode](https://developer.apple.com/xcode/) and the [GNU Fortran compiler on Mac](https://github.com/fxcoudert/gfortran-for-macOS) (see also 'Mandatory tools' [here](https://mac.r-project.org/tools/)). On Linux, the gfortran compiler is usually installed already.
 
-Vignettes are not rendered by default, if you want to include additional documentation please use:
-
-``` r
-if(!require(remotes)){install.packages("remotes")}
-remotes::install_github("geco-bern/rsofun", build_vignettes = TRUE)
-library("rsofun")
-```
-
 ### From source
 
 Assuming `rsofun` is the location of the source directory, on can build the R package (with extension .tar.gz) from the command line using:

diff --git a/data-raw/generate_biomee_drivers.R b/data-raw/generate_biomee_drivers.R
@@ -153,7 +153,7 @@ init_cohort <- tibble(
   init_cohort_bsw     = rep(0.05, 1), # initial biomass of sapwood, kg C/individual
   init_cohort_bHW     = rep(0.0, 1),  # initial biomass of heartwood, kg C/tree
   init_cohort_seedC   = rep(0.0, 1),  # initial biomass of seeds, kg C/individual
-  init_cohort_nsc     = rep(0.05, 1), # initial non-structural biomass
+  init_cohort_nsc     = rep(0.05, 1), # initial non-structural biomass, kg C/individual
   lu_index            = rep(0, 1)     # index land use (LU) containing this cohort. 0 (default) means any vegetated tile will contain a copy.
 )
 

diff --git a/data/biomee_gs_leuning_output.rda b/data/biomee_gs_leuning_output.rda
diff --git a/data/biomee_p_model_luluc_output.rda b/data/biomee_p_model_luluc_output.rda
diff --git a/data/biomee_p_model_output.rda b/data/biomee_p_model_output.rda
diff --git a/src/cohort_biomee.mod.f90 b/src/cohort_biomee.mod.f90
@@ -373,8 +373,7 @@ function can_be_merged_with(self, other) result(res)
 
       sameSizeTree = (spdata(self%species)%lifeform > 0).and.  &
               (spdata(other%species)%lifeform > 0).and.  &
-              ((dbh_diff/(self_dbh + other_dbh) < 0.1 ) .or.  &
-                      (dbh_diff < 0.001))  ! it'll be always true for grasses
+              (dbh_diff <= 0.01)  ! for TREES: if DBH difference <= 1cm    ! it'll be always true for grasses
 
       sameSizeGrass= (spdata(self%species)%lifeform == 0) .and. &
               (spdata(other%species)%lifeform == 0) .and. &

diff --git a/src/interface_in_biosphere_biomee.mod.f90 b/src/interface_in_biosphere_biomee.mod.f90
@@ -445,7 +445,7 @@ subroutine populate_spec_data(self, params_species)
     self%thetaHT            = real( params_species(25)) ! prescribed
     self%alphaCA            = real( params_species(26)) ! prescribed
     self%thetaCA            = real( params_species(27)) ! prescribed
-    self%alphaBM            = real( params_species(28)) ! prescribed
+    self%alphaBM            = real( params_species(28)) ! prescribed ! TODO: not actually not prescribed, it is overwritten by alphaBM
     self%thetaBM            = real( params_species(29)) ! prescribed
     self%seedlingsize       = real( params_species(30))
     self%maturalage         = real( params_species(31))
@@ -487,7 +487,7 @@ subroutine init_pft_data(self)
     ! specific root area
     self%SRA           = 2.0/(self%root_r*self%rho_FR)
 
-    ! calculate alphaBM parameter of allometry. note that rho_wood was re-introduced for this calculation
+    ! calculate alphaBM parameter of allometry. note that rho_wood was re-introduced for this calculation ! TODO: note that this overwrites the parameter alphaBM
     self%alphaBM = self%rho_wood * self%taperfactor * PI/4. * self%alphaHT ! 5200
 
     ! Vmax as a function of LNbase

diff --git a/src/main_biomee.mod.f90 b/src/main_biomee.mod.f90
@@ -80,7 +80,7 @@ subroutine biomee_f(            &
     ! Output arrays (naked) to be passed back to C/R
     real(kind=c_double), dimension(nt_daily,nvars_daily_tile, n_lu), intent(out) :: output_daily_tile
     real(kind=c_double), dimension(nt_annual,nvars_annual_tile, n_lu), intent(out) :: output_annual_tile
-    real(kind=c_double), dimension(out_max_cohorts, nt_annual_trans, nvars_annual_cohorts, n_lu), &
+    real(kind=c_double), dimension(NCohortMax, nt_annual_trans, nvars_annual_cohorts, n_lu), &
             intent(out) :: output_annual_cohorts
     real(kind=c_double), dimension(nt_annual,nvars_aggregated_out), intent(out) :: output_annual_aggregated
 

diff --git a/src/vegetation_processes_biomee.mod.f90 b/src/vegetation_processes_biomee.mod.f90
@@ -735,14 +735,15 @@ subroutine vegn_reproduction( vegn )
     type(cohort_stack_item), pointer :: new ! new cohort
     type(cohort_stack_item), pointer :: it
     integer, dimension(NCohortMax) :: reproPFTs
-    type(orgpool), dimension(NCohortMax) :: seed ! seed pool of productible PFTs
-    integer :: nPFTs ! number of new cohorts to be created
+    type(orgpool), dimension(NCohortMax) :: seed ! seed pool of reproducing PFTs
+    integer :: nPFTs ! number of new cohorts to be created (one new cohort per PFT)
     integer :: pft_idx ! PFT index
     integer :: k ! new cohort indices
 
-    ! Looping through all reproductable cohorts and check if reproduction happens
-    reproPFTs = -999 ! the code of reproductive PFT
+    ! Looping through all cohorts and check if reproduction happens
+    reproPFTs = -999 ! the code of reproducing PFT
     nPFTs = 0
+    ! NOTE: since seed is an orgpool(), it initializes by design to zero
 
     ! We loop through each cohort and add C and N to species-specific seed pools
     ! whenever a cohort can reproduce
@@ -760,7 +761,7 @@ subroutine vegn_reproduction( vegn )
         enddo
 
         if (pft_idx == 0) then ! when it is a new PFT, put it to the next place
-          nPFTs            = nPFTs + 1 ! update the number of reproducible PFTs
+          nPFTs            = nPFTs + 1 ! update the number of reproducing PFTs
           reproPFTs(nPFTs) = cc%species ! PFT number
           pft_idx = nPFTs
         endif
@@ -777,7 +778,7 @@ subroutine vegn_reproduction( vegn )
     ! We build new cohorts for seedlings
     do k = 1, nPFTs
 
-      new => vegn%create_cohort()
+      new => vegn%create_cohort() ! NOTE: this could lead to more cohorts than `NCohortMax`
       cc => new%cohort
 
       ! update child cohort parameters

diff --git a/src/vegetation_tile_biomee.mod.f90 b/src/vegetation_tile_biomee.mod.f90
@@ -18,8 +18,7 @@ module md_vegetation_tile_biomee
   public :: params_species_biomee, cohort_type, vegn_tile_type, dampended_forcing_type
 
   !=============== Parameters ======================================================
-  integer, public, parameter :: NCohortMax = 50 ! maximum number of cohorts
-  integer, public, parameter :: out_max_cohorts      = NCohortMax
+  integer, public, parameter :: NCohortMax = 50 ! maximum number of cohorts that are output (there could be more but they would not be reported)
 
   !===== Model
   integer, public, parameter :: NLAYERS_MAX = 9     ! maximum number of canopy layers to be considered
@@ -165,7 +164,7 @@ module md_vegetation_tile_biomee
     !===== Reporting diagnostics
     real, dimension(nvars_annual_tile) :: out_annual_tile
     real, dimension(ndayyear, nvars_daily_tile) :: out_daily_tile
-    real, dimension(out_max_cohorts, nvars_annual_cohorts) :: out_annual_cohorts
+    real, dimension(NCohortMax, nvars_annual_cohorts) :: out_annual_cohorts
 
   contains