diff --git a/NEWS.md b/NEWS.md
index bec4d5fd..abc96091 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,8 @@
 # rsofun (development version)
 
+* Breaking change: biomee drivers' `init_cohort$init_n_cohorts` column has been phased out and must not be present in  
+drivers to protect against data corruption.
+
 # rsofun v5.0.0
 
 * new BiomeE forcing data matching that of P-model
diff --git a/R/data.R b/R/data.R
index b627d7e6..731e20d1 100644
--- a/R/data.R
+++ b/R/data.R
@@ -305,9 +305,12 @@
 #'       \item{init_cohort_species}{Index of a species described in param_species.}
 #'       \item{init_cohort_nindivs}{Initial individual density, in individuals per
 #'         m\eqn{^{2}}.}
+#'       \item{init_cohort_bl}{Initial biomass of leaf, in kg C per individual.}
+#'       \item{init_cohort_br}{Initial biomass of fine root, in kg C per individual.}
 #'       \item{init_cohort_bsw}{Initial biomass of sapwood, in kg C per individual.}
-#'       \item{init_cohort_bHW}{Initial biomass of heartwood, in kg C per tree.}
-#'       \item{init_cohort_nsc}{Initial non-structural biomass.}
+#'       \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.}
 #'     }}
 #'   \item{init_soil}{A data.frame of initial soil pools, including
 #'   the following data:
diff --git a/R/run_biomee_f_bysite.R b/R/run_biomee_f_bysite.R
index c48faa1d..d229b4f3 100644
--- a/R/run_biomee_f_bysite.R
+++ b/R/run_biomee_f_bysite.R
@@ -343,6 +343,11 @@ run_biomee_f_bysite <- function(
           return(TRUE)
         }
       })
+
+    if ('init_n_cohorts' %in% names(init_cohort)) {
+      warning("Error: column 'init_n_cohorts' under 'init_cohort' has been phased out and must be removed from the drivers.")
+      data_integrity <- append(data_integrity, FALSE)
+    }
     
     # only return true if all checked variables are TRUE
     # suppress warning on coercion of list to single logical
diff --git a/data-raw/generate_biomee_drivers.R b/data-raw/generate_biomee_drivers.R
index 9060e406..a1df93c8 100644
--- a/data-raw/generate_biomee_drivers.R
+++ b/data-raw/generate_biomee_drivers.R
@@ -74,81 +74,80 @@ params_tile <- tibble(
 )
 
 params_species <- tibble(
-  lifeform      = rep(1,16),
-  phenotype     = c(0,1,1,rep(1,13)),
-  pt            = rep(0,16),
+  lifeform      = rep(1,4),
+  phenotype     = c(0,1,1,1),
+  pt            = rep(0,4),
   # Root parameters
-  alpha_FR      = rep(1.2,16),
-  rho_FR        = rep(200,16),
-  root_r        = rep(2.9E-4,16), 
-  root_zeta     = rep(0.29,16), 
-  Kw_root       = rep(3.5e-09,16),
-  leaf_size     = rep(0.04,16), 
+  alpha_FR      = rep(1.2,4),
+  rho_FR        = rep(200,4),
+  root_r        = rep(2.9E-4,4),
+  root_zeta     = rep(0.29,4),
+  Kw_root       = rep(3.5e-09,4),
+  leaf_size     = rep(0.04,4),
   # Photosynthesis parameters
-  Vmax          = rep(35.0E-6,16),
-  Vannual       = rep(1.2,16),
-  wet_leaf_dreg = rep(0.3,16),
-  m_cond        = rep(7.0,16), 
-  alpha_phot    = rep(0.06,16), 
-  gamma_L       = rep(0.02,16), 
-  gamma_LN      = rep(70.5 ,16),
-  gamma_SW      = rep(0.08,16),
-  gamma_FR      = rep(12.0,16),
-  tc_crit       = rep(283.16,16),
-  tc_crit_on    = rep(280.16,16),
-  gdd_crit      = rep(280.0,16),
-  betaON        = rep(0,2,16),     
-  betaOFF       = rep(0,1,16), 
+  Vmax          = rep(35.0E-6,4),
+  Vannual       = rep(1.2,4),
+  wet_leaf_dreg = rep(0.3,4),
+  m_cond        = rep(7.0,4),
+  alpha_phot    = rep(0.06,4),
+  gamma_L       = rep(0.02,4),
+  gamma_LN      = rep(70.5,4),
+  gamma_SW      = rep(0.08,4),
+  gamma_FR      = rep(12.0,4),
+  tc_crit       = rep(283.16,4),
+  tc_crit_on    = rep(280.16,4),
+  gdd_crit      = rep(280.0,4),
+  betaON        = rep(0,2,4),
+  betaOFF       = rep(0,1,4),
   # Allometry parameters
-  alphaHT       = rep(36,16),                   
-  thetaHT       = rep(0.5,16),                   
-  alphaCA       = rep(150,16),                   
-  thetaCA       = rep(1.5,16),                   
-  alphaBM       = rep(5200,16),                   
-  thetaBM       = c(2.36,2.30,2.54,rep(2.30,13)), 
+  alphaHT       = rep(36,4),
+  thetaHT       = rep(0.5,4),
+  alphaCA       = rep(150,4),
+  thetaCA       = rep(1.5,4),
+  alphaBM       = rep(5200,4),
+  thetaBM       = c(2.36,2.30,2.54,2.30),
   # Reproduction parameters
-  seedlingsize  = rep(0.05,16),
-  maturalage    = rep(5,16),     
-  v_seed        = rep(0.1,16), 
+  seedlingsize  = rep(0.05,4),
+  maturalage    = rep(5,4),
+  v_seed        = rep(0.1,4),
   # Mortality parameters
-  mortrate_d_c  = rep(0.01,16),                   
-  mortrate_d_u  = rep(0.075,16), 
+  mortrate_d_c  = rep(0.01,4),
+  mortrate_d_u  = rep(0.075,4),
   # Leaf parameters
-  LMA           = c(0.05,0.17,0.11,rep(0.1,13)),  
-  leafLS        = rep(1,16), 
-  LNbase        = rep(0.8E-3,16), 
-  CNleafsupport = rep(80,16),
-  rho_wood      = c(590,370,350,rep(300,13)),   
-  taperfactor   = rep(0.75,16),
-  lAImax        = rep(3.5,16), 
-  tauNSC        = rep(3,16), 
-  fNSNmax       = rep(5,16),                      
-  phiCSA        = rep(0.25E-4,16),    
+  LMA           = c(0.05,0.17,0.11,0.1),
+  leafLS        = rep(1,4),
+  LNbase        = rep(0.8E-3,4),
+  CNleafsupport = rep(80,4),
+  rho_wood      = c(590,370,350,300),
+  taperfactor   = rep(0.75,4),
+  lAImax        = rep(3.5,4),
+  tauNSC        = rep(3,4),
+  fNSNmax       = rep(5,4),
+  phiCSA        = rep(0.25E-4,4),
   # C/N ratios for plant pools
-  CNleaf0      = rep(25,16),  
-  CNsw0        = rep(350,16),  
-  CNwood0      = rep(350,16),  
-  CNroot0      = rep(40,16),  
-  CNseed0      = rep(20,16),  
-  Nfixrate0     = rep(0,16),   
-  NfixCost0     = rep(12,16),
-  internal_gap_frac      = rep(0.1,16),
+  CNleaf0       = rep(25,4),
+  CNsw0         = rep(350,4),
+  CNwood0       = rep(350,4),
+  CNroot0       = rep(40,4),
+  CNseed0       = rep(20,4),
+  Nfixrate0     = rep(0,4),
+  NfixCost0     = rep(12,4),
+  internal_gap_frac      = rep(0.1,4),
   # calibratable params
-  kphio         = rep(0.05,16),
-  phiRL         = rep(3.5,16),
-  LAI_light     = rep(3.5,16)
+  kphio         = rep(0.05,4),
+  phiRL         = rep(3.5,4),
+  LAI_light     = rep(3.5,4)
 )
 
 init_cohort <- tibble(
-  init_n_cohorts = 1,   # number of PFTs
-  init_cohort_species = rep(2, 10),    # indicates sps # 1 - Fagus sylvatica
-  init_cohort_nindivs = rep(0.05,10),  # initial individual density, individual/m2 ! 1 indiv/m2 = 10.000 indiv/ha
-  init_cohort_bl      = rep(0.0,10),   # initial biomass of leaves, kg C/individual
-  init_cohort_br      = rep(0.0, 10),  # initial biomass of fine roots, kg C/individual
-  init_cohort_bsw     = rep(0.05,10),  # initial biomass of sapwood, kg C/individual
-  init_cohort_bHW     = rep(0.0, 10),  # initial biomass of heartwood, kg C/tree
-  init_cohort_seedC   = rep(0.0, 10),  # initial biomass of seeds, kg C/individual
-  init_cohort_nsc     = rep(0.05,10)   # initial non-structural biomass
+  init_cohort_species = 2,    # indicates sps # 1 - Fagus sylvatica
+  init_cohort_nindivs = 0.05, # initial individual density, individual/m2 ! 1 indiv/m2 = 10.000 indiv/ha
+  init_cohort_bl      = 0.0,  # initial biomass of leaves, kg C/individual
+  init_cohort_br      = 0.0,  # initial biomass of fine roots, kg C/individual
+  init_cohort_bsw     = 0.05, # initial biomass of sapwood, kg C/individual
+  init_cohort_bHW     = 0.0,  # initial biomass of heartwood, kg C/tree
+  init_cohort_seedC   = 0.0,  # initial biomass of seeds, kg C/individual
+  init_cohort_nsc     = 0.05  # initial non-structural biomass
 )
 
 init_soil <- tibble( #list
diff --git a/data/biomee_gs_leuning_drivers.rda b/data/biomee_gs_leuning_drivers.rda
index 5df8b5be..b729ddfd 100644
Binary files a/data/biomee_gs_leuning_drivers.rda and b/data/biomee_gs_leuning_drivers.rda differ
diff --git a/data/biomee_p_model_drivers.rda b/data/biomee_p_model_drivers.rda
index 18d2f22e..4a2372ba 100644
Binary files a/data/biomee_p_model_drivers.rda and b/data/biomee_p_model_drivers.rda differ
diff --git a/man/biomee_gs_leuning_drivers.Rd b/man/biomee_gs_leuning_drivers.Rd
index c43ee693..596702d5 100644
--- a/man/biomee_gs_leuning_drivers.Rd
+++ b/man/biomee_gs_leuning_drivers.Rd
@@ -130,9 +130,12 @@ This data structure can be freely used for documenting the dataset, but must inc
       \item{init_cohort_species}{Index of a species described in param_species.}
       \item{init_cohort_nindivs}{Initial individual density, in individuals per
         m\eqn{^{2}}.}
+      \item{init_cohort_bl}{Initial biomass of leaf, in kg C per individual.}
+      \item{init_cohort_br}{Initial biomass of fine root, in kg C per individual.}
       \item{init_cohort_bsw}{Initial biomass of sapwood, in kg C per individual.}
-      \item{init_cohort_bHW}{Initial biomass of heartwood, in kg C per tree.}
-      \item{init_cohort_nsc}{Initial non-structural biomass.}
+      \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.}
     }}
   \item{init_soil}{A data.frame of initial soil pools, including
   the following data:
diff --git a/src/interface_biosphere_biomee.mod.f90 b/src/interface_biosphere_biomee.mod.f90
index 45ad9a9b..e9a3c48f 100644
--- a/src/interface_biosphere_biomee.mod.f90
+++ b/src/interface_biosphere_biomee.mod.f90
@@ -119,8 +119,7 @@ module md_interface_biomee
 
   end type
 
-  type inittype_cohort 
-    integer :: init_n_cohorts
+  type inittype_cohort
     integer :: init_cohort_species
     real :: init_cohort_nindivs
     real :: init_cohort_bl
diff --git a/src/sofun_r.f90 b/src/sofun_r.f90
index 17ff7d39..c9da107f 100644
--- a/src/sofun_r.f90
+++ b/src/sofun_r.f90
@@ -558,15 +558,14 @@ subroutine biomee_f(            &
 
     ! Initial cohort sizes
     allocate(myinterface%init_cohort(n_init_cohort))
-    myinterface%init_cohort(:)%init_n_cohorts      = int(init_cohort(:,1))
-    myinterface%init_cohort(:)%init_cohort_species = int(init_cohort(:,2))
-    myinterface%init_cohort(:)%init_cohort_nindivs = real(init_cohort(:,3))
-    myinterface%init_cohort(:)%init_cohort_bl      = real(init_cohort(:,4))
-    myinterface%init_cohort(:)%init_cohort_br      = real(init_cohort(:,5))
-    myinterface%init_cohort(:)%init_cohort_bsw     = real(init_cohort(:,6))
-    myinterface%init_cohort(:)%init_cohort_bHW     = real(init_cohort(:,7))
-    myinterface%init_cohort(:)%init_cohort_seedC   = real(init_cohort(:,8))
-    myinterface%init_cohort(:)%init_cohort_nsc     = real(init_cohort(:,9))
+    myinterface%init_cohort(:)%init_cohort_species = int(init_cohort(:,1))
+    myinterface%init_cohort(:)%init_cohort_nindivs = real(init_cohort(:,2))
+    myinterface%init_cohort(:)%init_cohort_bl      = real(init_cohort(:,3))
+    myinterface%init_cohort(:)%init_cohort_br      = real(init_cohort(:,4))
+    myinterface%init_cohort(:)%init_cohort_bsw     = real(init_cohort(:,5))
+    myinterface%init_cohort(:)%init_cohort_bHW     = real(init_cohort(:,6))
+    myinterface%init_cohort(:)%init_cohort_seedC   = real(init_cohort(:,7))
+    myinterface%init_cohort(:)%init_cohort_nsc     = real(init_cohort(:,8))
 
     ! Initial soil pools
     myinterface%init_soil%init_fast_soil_C = real( init_fast_soil_C )
diff --git a/src/vegetation_biomee.mod.f90 b/src/vegetation_biomee.mod.f90
index 745ff74c..baea6aab 100755
--- a/src/vegetation_biomee.mod.f90
+++ b/src/vegetation_biomee.mod.f90
@@ -2213,7 +2213,7 @@ subroutine initialize_vegn_tile( vegn )
 
     ! xxx seems new from d-ben - missing if?
     ! Initialize plant cohorts
-    init_n_cohorts = myinterface%init_cohort(1)%init_n_cohorts
+    init_n_cohorts = size(myinterface%init_cohort)
     allocate(cc(1:init_n_cohorts), STAT = istat)
     vegn%cohorts => cc
     vegn%n_cohorts = init_n_cohorts