diff --git a/R/run_biomee_f_bysite.R b/R/run_biomee_f_bysite.R
index 3d08cfe2..f476b27f 100644
--- a/R/run_biomee_f_bysite.R
+++ b/R/run_biomee_f_bysite.R
@@ -561,10 +561,10 @@ run_biomee_f_bysite <- function(
       init_Nmineral    = as.numeric(init_soil$init_Nmineral),
       N_input          = as.numeric(init_soil$N_input),
       n                = as.integer(nrow(forcing)), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
-      n_daily          = as.integer(n_daily), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
-      n_annual         = as.integer(runyears), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
-      #n_annual_cohorts = as.integer(params_siml$nyeartrend), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
-      n_annual_cohorts = as.integer(runyears), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
+      n_daily          = as.integer(n_daily), 
+      n_annual         = as.integer(runyears), 
+      n_annual_cohorts = as.integer(params_siml$nyeartrend), # to get cohort outputs after spinup year
+      #n_annual_cohorts = as.integer(runyears), # to get cohort outputs from year 1
       forcing          = as.matrix(forcing)
       )
     
diff --git a/src/sofun_r.f90 b/src/sofun_r.f90
index a59dc381..0d0ce793 100644
--- a/src/sofun_r.f90
+++ b/src/sofun_r.f90
@@ -768,10 +768,18 @@ subroutine biomee_f(            &
       !----------------------------------------------------------------
       ! Output output_annual_cohorts (without subroutine)
       !----------------------------------------------------------------
-      if (.not. myinterface%steering%spinup) then  ! To get outputs only after spinupyears
+
+      ! To get outputs only after spinupyears make if below and 
+      ! also in run_biomee_f_bysite.R make n_annual_cohorts = as.integer(params_siml$nyeartrend)
+
+      if (.not. myinterface%steering%spinup) then  
 
         idx =  yr - myinterface%params_siml%spinupyears
 
+      ! To get outputs for all runyears idx=yr and also in run_biomee_f_bysite.R make n_annual_cohorts = as.integer(runyears)
+       
+       !idx =  yr
+
         output_annual_cohorts_year(idx, :)        = dble(out_biosphere%annual_cohorts(:)%year)
         output_annual_cohorts_cID(idx, :)         = dble(out_biosphere%annual_cohorts(:)%cID)
         output_annual_cohorts_PFT(idx, :)         = dble(out_biosphere%annual_cohorts(:)%PFT)