diff --git a/NAMESPACE b/NAMESPACE
index 599002e41..42bc8501b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -15,4 +15,5 @@ import(BayesianTools)
 import(GenSA)
 import(dplyr)
 importFrom(magrittr,"%>%")
+importFrom(stats,setNames)
 useDynLib(rsofun)
diff --git a/NEWS.md b/NEWS.md
index c02240591..530918e1d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,10 @@
 # rsofun (development version)
 
 * New `run_pmodel_onestep_f_bysite()` adds single-step leaf-level simulations of the P-model.
+* `calib_sofun()` now passes parameters as a _named_ vector to cost-functions for 
+  easier processing within cost-functions. Default cost-functions were updated,
+  but currently ignore the names. This is fully backward compatible, but allows
+  to use the names in user-created cost functions.
 
 # rsofun 5.1.0
 
diff --git a/R/calib_sofun.R b/R/calib_sofun.R
index bf066b011..98ea42c0a 100644
--- a/R/calib_sofun.R
+++ b/R/calib_sofun.R
@@ -4,10 +4,13 @@
 #' calibration of SOFUN model parameters. 
 #' 
 #' @param drivers A data frame with driver data. See \code{\link{p_model_drivers}}
-#' for a description of the data structure.
+#' for a description of the data structure. Additional columns can optionally be 
+#' provided to \code{drivers} to control e.g. the processing within a 
+#' personalized cost function.
 #' @param obs A data frame containing observational data used for model
 #'  calibration. See \code{\link{p_model_validation}} for a description of the data
-#'  structure.
+#'  structure. Additional columns can optionally be provided to \code{obs} to 
+#'  control e.g. the processing within a personalized cost function.
 #' @param settings A list containing model calibration settings. 
 #'  See the 'P-model usage' vignette for more information and examples.
 #'  \describe{
@@ -31,9 +34,8 @@
 #'  }
 #' @param optim_out A logical indicating whether the function returns the raw
 #'  output of the optimization functions (defaults to TRUE).
-#' @param ... Optional arguments passed on to the cost function specified as
-#'  \code{settings$metric}. 
-#'  . 
+#' @param ... Optional arguments, simply passed on to the cost function. 
+#' 
 #' @return A named list containing the calibrated parameter vector `par` and
 #' the output object from the optimization `mod`. For more details on this
 #' output and how to evaluate it, see \link[BayesianTools:runMCMC]{runMCMC} (also
@@ -41,6 +43,7 @@
 #' and \link[GenSA]{GenSA}.
 #' @export
 #' @importFrom magrittr %>%
+#' @importFrom stats setNames
 #' @import GenSA BayesianTools
 #' 
 #' @examples
@@ -153,12 +156,16 @@ calib_sofun <- function(
     }
     
     # reformat parameters
-    pars <- as.data.frame(do.call("rbind", settings$par), row.names = FALSE)
+    pars <- as.data.frame(do.call("rbind", settings$par)) 
+         # NOTE: This keeps parameters as row names. 
+         # This does not change anything to previous behavior.
+         # But this agrees better with the example data 
+         # `dput(BayesianTools::VSEMgetDefaults())`
     
     priors  <- BayesianTools::createUniformPrior(
-      unlist(pars$lower),
-      unlist(pars$upper),
-      unlist(pars$init)
+      lower = unlist(pars$lower),
+      upper = unlist(pars$upper),
+      best  = unlist(pars$init)
     )
     
     # setup the bayes run, no message forwarding is provided
@@ -168,14 +175,18 @@ calib_sofun <- function(
         do.call(
           "cost",
           list(
-            par = random_par,
+            par = setNames(random_par, rownames(pars)),
+            # NOTE: if we could make use of setup$names from within the cost
+            # function then we wouldn't need this closure (using `pars`), but it
+            # appears that BayesianTools does not pass the names into the
+            # likelihood.
             obs = obs,
             drivers = drivers
           )
         )
       },
     prior = priors,
-    names = names(settings$par)
+    names = rownames(pars)
     )    
     
     # set bt control parameters
diff --git a/R/cost_likelihood_biomee.R b/R/cost_likelihood_biomee.R
index 40f8d89cb..116fb85ee 100644
--- a/R/cost_likelihood_biomee.R
+++ b/R/cost_likelihood_biomee.R
@@ -4,7 +4,7 @@
 #' computes the log-likelihood for the biomee model fitting several target 
 #' variables for a given set of parameters.
 #' 
-#' @param par A vector containing parameter values for \code{'phiRL',
+#' @param par A named vector containing parameter values for \code{'phiRL',
 #' 'LAI_light', 'tf_base', 'par_mort'} in that order, and for the error terms
 #' corresponding to the target variables, e.g. \code{'err_GPP'} if GPP is a target. 
 #' Make sure that
@@ -34,8 +34,11 @@
 #' # BiomeE model parameter values
 #' # and the example data
 #' cost_likelihood_biomee(
-#'  par = c(3.5, 3.5, 1, 1,    # model params
-#'          0.5),              # err_GPP
+#'  par = c(phiRL = 3.5, 
+#'          LAI_light = 3.5, 
+#'          tf_base = 1, 
+#'          par_mort = 1,    # model params
+#'          err_GPP = 0.5),  # err_GPP
 #'  obs = biomee_validation,
 #'  drivers = biomee_gs_leuning_drivers,
 #'  targets = c("GPP")
@@ -52,6 +55,7 @@ cost_likelihood_biomee <- function(
   # predefine variables for CRAN check compliance
   GPP <- LAI <- Density12 <- plantC <- NULL
   
+  par <- unname(par) # reproduces previous behavior, when par was unnamed
   # Add changed model parameters to drivers, overwriting where necessary.
   drivers$params_species[[1]]$phiRL[]  <- par[1]
   drivers$params_species[[1]]$LAI_light[]  <- par[2]
diff --git a/R/cost_likelihood_pmodel.R b/R/cost_likelihood_pmodel.R
index 74b16a52a..a82a9b42c 100644
--- a/R/cost_likelihood_pmodel.R
+++ b/R/cost_likelihood_pmodel.R
@@ -6,7 +6,7 @@
 #' observed values and with standard deviation given as an input parameter 
 #' (calibratable).
 #' 
-#' @param par A vector of values for the parameters to be calibrated, including
+#' @param par A named vector of values for the parameters to be calibrated, including
 #' a subset of model parameters (described in \code{\link{runread_pmodel_f}}),
 #' in order, and error terms
 #' for each target variable (for example \code{'gpp_err'}), in the same order as
@@ -56,8 +56,10 @@
 #' # temperature dependence of kphio 
 #' # and example data
 #' cost_likelihood_pmodel(
-#'  par = c(0.05, -0.01, 1,     # model parameters
-#'          2),                # err_gpp
+#' par = c(kphio       = 0.05, 
+#'         kphio_par_a = -0.01, 
+#'         kphio_par_b = 1,     # model parameters
+#'         err_gpp     = 2),    # err_gpp
 #'  obs = p_model_validation,
 #'  drivers = p_model_drivers,
 #'  targets = c('gpp'),
@@ -82,6 +84,7 @@ cost_likelihood_pmodel <- function(
 ){
   # predefine variables for CRAN check compliance
   sitename <- data <- gpp_mod <- NULL
+  par <- unname(par) # reproduces previous behavior, when par was unnamed
   
   ## check input parameters
   if( (length(par) + length(par_fixed)) != (9 + length(targets)) ){
@@ -95,7 +98,7 @@ cost_likelihood_pmodel <- function(
                          'beta_unitcostratio', 'rd_to_vcmax', 
                          'tau_acclim', 'kc_jmax')
   
-  if(!is.null(par_fixed)){
+  if(!is.null(par_fixed) && length(par)>0){
     params_modl <- list()
     # complete with calibrated values
     i <- 1 # start counter
@@ -107,9 +110,12 @@ cost_likelihood_pmodel <- function(
         params_modl[[par_name]] <- par_fixed[[par_name]]  # use fixed par value
       }
     }
+  }else if(length(par)==0){                # no parameters calibrated
+    params_modl <- as.list(par_fixed[calib_param_names])
+    par <- par_fixed[grepl("err_",names(par_fixed))]
   }else{
     params_modl <- as.list(par[1:9])       # all parameters calibrated
-    names(params_modl) <- calib_param_names
+    names(params_modl) <- calib_param_names# TODO: problematic, since it assumes they are in the right order
   }
   
   ## run the model
diff --git a/analysis/01-sensitivity-analysis.R b/analysis/01-sensitivity-analysis.R
index 836bf1c4f..9f954dab0 100644
--- a/analysis/01-sensitivity-analysis.R
+++ b/analysis/01-sensitivity-analysis.R
@@ -13,7 +13,7 @@ set.seed(432)
 
 # Define log-likelihood function
 ll_pmodel <- function(
-    par_v                 # a vector of all calibratable parameters including errors
+    par_v                 # a named vector of all calibratable parameters including errors
 ){
   rsofun::cost_likelihood_pmodel(        # likelihood cost function from package
     par_v,
diff --git a/man/calib_sofun.Rd b/man/calib_sofun.Rd
index 4ee95d41a..4081d76fb 100644
--- a/man/calib_sofun.Rd
+++ b/man/calib_sofun.Rd
@@ -8,11 +8,14 @@ calib_sofun(drivers, obs, settings, optim_out = TRUE, ...)
 }
 \arguments{
 \item{drivers}{A data frame with driver data. See \code{\link{p_model_drivers}}
-for a description of the data structure.}
+for a description of the data structure. Additional columns can optionally be 
+provided to \code{drivers} to control e.g. the processing within a 
+personalized cost function.}
 
 \item{obs}{A data frame containing observational data used for model
 calibration. See \code{\link{p_model_validation}} for a description of the data
-structure.}
+structure. Additional columns can optionally be provided to \code{obs} to 
+control e.g. the processing within a personalized cost function.}
 
 \item{settings}{A list containing model calibration settings. 
 See the 'P-model usage' vignette for more information and examples.
@@ -39,9 +42,7 @@ See the 'P-model usage' vignette for more information and examples.
 \item{optim_out}{A logical indicating whether the function returns the raw
 output of the optimization functions (defaults to TRUE).}
 
-\item{...}{Optional arguments passed on to the cost function specified as
-\code{settings$metric}. 
-.}
+\item{...}{Optional arguments, simply passed on to the cost function.}
 }
 \value{
 A named list containing the calibrated parameter vector `par` and
diff --git a/man/cost_likelihood_biomee.Rd b/man/cost_likelihood_biomee.Rd
index 2ba354ba2..e010a7b72 100644
--- a/man/cost_likelihood_biomee.Rd
+++ b/man/cost_likelihood_biomee.Rd
@@ -7,7 +7,7 @@
 cost_likelihood_biomee(par, obs, drivers, targets)
 }
 \arguments{
-\item{par}{A vector containing parameter values for \code{'phiRL',
+\item{par}{A named vector containing parameter values for \code{'phiRL',
 'LAI_light', 'tf_base', 'par_mort'} in that order, and for the error terms
 corresponding to the target variables, e.g. \code{'err_GPP'} if GPP is a target. 
 Make sure that
@@ -45,8 +45,11 @@ should be run using \code{BayesianTools}, so the likelihood is maximized.
 # BiomeE model parameter values
 # and the example data
 cost_likelihood_biomee(
- par = c(3.5, 3.5, 1, 1,    # model params
-         0.5),              # err_GPP
+ par = c(phiRL = 3.5, 
+         LAI_light = 3.5, 
+         tf_base = 1, 
+         par_mort = 1,    # model params
+         err_GPP = 0.5),  # err_GPP
  obs = biomee_validation,
  drivers = biomee_gs_leuning_drivers,
  targets = c("GPP")
diff --git a/man/cost_likelihood_pmodel.Rd b/man/cost_likelihood_pmodel.Rd
index d206d8126..302b66530 100644
--- a/man/cost_likelihood_pmodel.Rd
+++ b/man/cost_likelihood_pmodel.Rd
@@ -16,7 +16,7 @@ cost_likelihood_pmodel(
 )
 }
 \arguments{
-\item{par}{A vector of values for the parameters to be calibrated, including
+\item{par}{A named vector of values for the parameters to be calibrated, including
 a subset of model parameters (described in \code{\link{runread_pmodel_f}}),
 in order, and error terms
 for each target variable (for example \code{'gpp_err'}), in the same order as
@@ -78,8 +78,10 @@ trait value predicted on that date.
 # temperature dependence of kphio 
 # and example data
 cost_likelihood_pmodel(
- par = c(0.05, -0.01, 1,     # model parameters
-         2),                # err_gpp
+par = c(kphio       = 0.05, 
+        kphio_par_a = -0.01, 
+        kphio_par_b = 1,     # model parameters
+        err_gpp     = 2),    # err_gpp
  obs = p_model_validation,
  drivers = p_model_drivers,
  targets = c('gpp'),
diff --git a/tests/testthat/test-likelihoods.R b/tests/testthat/test-likelihoods.R
new file mode 100644
index 000000000..ceb924ea4
--- /dev/null
+++ b/tests/testthat/test-likelihoods.R
@@ -0,0 +1,257 @@
+context("test P-model and BiomeE likelihood frameworks")
+set.seed(10)
+
+test_that("test likelihood/RMSE calculations with pmodel", {
+  
+  test_params_pmodel <- data.frame( # test_params_pmodel was generated with dput(test_params_pmodel)
+    kphio              = c(0.04998, 0.101645649550483, 0.142234791386873,0.136563287638128, 0.0529878485854715), 
+    kphio_par_a        = c(0.01   , -0.00398199869785458, -0.00144068546174094, -0.00302413401333615, -0.00134308318863623), 
+    kphio_par_b        = c(1.0    , 15.59719867073, 14.5120448060334, 20.2669072570279, 18.2723819604144), 
+    soilm_thetastar    = c(0.6*240, 21.4685604907572, 129.50339589268, 36.7124842107296, 209.40362457186), 
+    soilm_betao        = c(0.01   , 0.508695727679878, 0.69720912235789, 0.306312796194106, 0.546225873986259), 
+    beta_unitcostratio = c(146.0  , 109.196580422577, 113.911265798379, 146.538958174642, 114.256114442833), 
+    rd_to_vcmax        = c(0.014  , 0.0468210919597186, 0.0725045789754949, 0.0594035412720405, 0.0491376937157475), 
+    tau_acclim         = c(30.0   , 53.9236626827624, 26.6468327937182, 34.3328710102942, 52.5165054232348), 
+    kc_jmax            = c(0.41   , 0.624640251602978, 0.492042974522337, 0.23459618492052, 0.502756303641945), 
+    err_gpp            = c(0.5    , 3.13616614692146, 2.82713630301412, 0.282233132501133, 3.00473784686066))
+  # test_params_pmodel was created with: par_cal_best <- c(
+  # test_params_pmodel was created with:   kphio              = 0.09423773,
+  # test_params_pmodel was created with:   kphio_par_a        = -0.0025,
+  # test_params_pmodel was created with:   kphio_par_b        = 20,
+  # test_params_pmodel was created with:   soilm_thetastar    = 0.6*240,
+  # test_params_pmodel was created with:   soilm_betao        = 0.2,
+  # test_params_pmodel was created with:   beta_unitcostratio = 146.0,
+  # test_params_pmodel was created with:   rd_to_vcmax        = 0.014,
+  # test_params_pmodel was created with:   tau_acclim         = 30.0,
+  # test_params_pmodel was created with:   kc_jmax            = 0.41,
+  # test_params_pmodel was created with:   err_gpp          = 1
+  # test_params_pmodel was created with: )
+  # test_params_pmodel was created with: par_cal_min <- c(
+  # test_params_pmodel was created with:   kphio              = 0.03,
+  # test_params_pmodel was created with:   kphio_par_a        = -0.004,
+  # test_params_pmodel was created with:   kphio_par_b        = 10,
+  # test_params_pmodel was created with:   soilm_thetastar    = 0,
+  # test_params_pmodel was created with:   soilm_betao        = 0,
+  # test_params_pmodel was created with:   beta_unitcostratio = 50.0,
+  # test_params_pmodel was created with:   rd_to_vcmax        = 0.01,
+  # test_params_pmodel was created with:   tau_acclim         = 7.0,
+  # test_params_pmodel was created with:   kc_jmax            = 0.2,
+  # test_params_pmodel was created with:   err_gpp          = 0.01
+  # test_params_pmodel was created with: )
+  # test_params_pmodel was created with: par_cal_max <- c(
+  # test_params_pmodel was created with:   kphio              = 0.15,
+  # test_params_pmodel was created with:   kphio_par_a        = -0.001,
+  # test_params_pmodel was created with:   kphio_par_b        = 30,
+  # test_params_pmodel was created with:   soilm_thetastar    = 240,
+  # test_params_pmodel was created with:   soilm_betao        = 1,
+  # test_params_pmodel was created with:   beta_unitcostratio = 200.0,
+  # test_params_pmodel was created with:   rd_to_vcmax        = 0.1,
+  # test_params_pmodel was created with:   tau_acclim         = 60.0,
+  # test_params_pmodel was created with:   kc_jmax            = 0.8,
+  # test_params_pmodel was created with:   err_gpp          = 4
+  # test_params_pmodel was created with: )
+  # test_params_pmodel was created with: library(BayesianTools) # for prior sampling
+  # test_params_pmodel was created with: prior <- createUniformPrior(lower = par_cal_min, upper = par_cal_max, best = par_cal_best)
+  # test_params_pmodel was created with: test_params_pmodel <- prior$sampler(4) |> as.data.frame() |>
+  # test_params_pmodel was created with:   stats::setNames(nm = names(par_cal_min))
+  # test_params_pmodel was created with: test_params_pmodel <- bind_rows(par_cal_best, test_params_pmodel)
+  
+  # also add error model for vcmax25
+  test_params_pmodel <- dplyr::mutate(test_params_pmodel, err_vcmax25 = 0.5) 
+  
+  
+  # Test rsofun::cost_likelihood_pmodel()
+  ll_values <- apply(test_params_pmodel |> dplyr::select(-err_vcmax25), 1, function(par_v) { # par_v is a named vector
+    # TODO: when rewriting cost_likelihood_pmodel: activate more difficult check to ignore unneded error parameter 'err_vcmax25'
+    #       ll_values <- apply(test_params_pmodel, 1, function(par_v) {...})
+    rsofun::cost_likelihood_pmodel(     # likelihood cost function from package
+      par = par_v,                      # par: should be a named vector
+      obs = rsofun::p_model_validation, # obs: example data from package
+      drivers = rsofun::p_model_drivers,# drivers: example data from package
+      targets = c('gpp'),
+      par_fixed = NULL)
+  })
+  testthat::expect_equal(
+    tolerance = 1e-4,
+    object = ll_values, 
+    # expected was generated with dput(ll_values)
+    expected = c(
+      -13706.5058738304,
+      -4109.97422397591,
+      -11148.2269071033,
+      -2255167.68663634,
+      -3846.02937529864)
+  )
+  
+  # Test rsofun::cost_rmse_pmodel()
+  rmse_values <- apply(dplyr::select(test_params_pmodel,-err_gpp, -err_vcmax25), 1, function(par_v) { # par_v is a named vector
+    rsofun::cost_rmse_pmodel(
+      par = par_v,                      # par: should be a named vector
+      obs = rsofun::p_model_validation, # obs: example data from package
+      drivers = rsofun::p_model_drivers,
+      targets = c('gpp'),
+      par_fixed = NULL
+    )
+  })
+  testthat::expect_equal(
+    tolerance = 1e-4,
+    object = rmse_values, 
+    # expected was generated with dput(rmse_values)
+    expected = c(
+      1.91661972744907, 
+      2.02647969696678,
+      8.19483248539096,
+      14.0907280919599,
+      1.38184787783589
+    )
+  )
+  
+  # Also test vcmax25 target and multi-target loglikelihoods:
+  ll_values2 <- apply(dplyr::select(test_params_pmodel, -err_gpp), 1, function(par_v) { # par_v is a named vector
+    # TODO: when rewriting cost_likelihood_pmodel: activate more difficult check to ignore unneded error parameter 'err_gpp'
+    #       ll_values2 <- apply(test_params_pmodel, 1, function(par_v) {...})
+    rsofun::cost_likelihood_pmodel(         # likelihood cost function from package
+      par     = par_v,                      # par: should be a named vector
+      obs     = p_model_validation_vcmax25, # obs: example data from package
+      drivers = p_model_drivers_vcmax25,    # drivers: example data from package
+      targets = c('vcmax25'))
+  })
+  ll_values3 <- apply(test_params_pmodel, 1, function(par_v) { # par_v is a named vector
+    rsofun::cost_likelihood_pmodel(                                    # likelihood cost function from package
+      par     = par_v,                                                 # par: should be a named vector
+      obs     = rbind(p_model_validation, p_model_validation_vcmax25), # obs: example data from package 
+      drivers = rbind(p_model_drivers, p_model_drivers_vcmax25),       # drivers: example data from package
+      targets = c('gpp', 'vcmax25'))
+  })
+  
+  testthat::expect_equal(ll_values3, ll_values + ll_values2) # loglikelihood of multiple targets is additive
+  testthat::expect_equal(
+    tolerance = 0.5, #tolerance = 1e-4,
+    object = ll_values3, 
+    # expected was generated with dput(ll_values3)
+    expected = c(
+      -13707.4063840446,
+      -4110.8773900703,
+      -11149.1301175005,
+      -2255168.59027969,
+      -3846.93254413504
+    )
+  )
+  testthat::expect_equal(
+    tolerance = 1e-4,
+    object = ll_values2, 
+    # expected was generated with dput(ll_values2)
+    expected = c(
+      -0.903165435731823,
+      -0.903165445893757,
+      -0.903165590656233,
+      -0.903165645611412, 
+      -0.90316542572855
+    )
+  )
+  
+  # test p-model likelihood with only fixed parameters
+  ll_pmodel_fixed <- rsofun::cost_likelihood_pmodel(
+    par     = c(),                                                   # par: should be a named vector
+    obs     = rbind(p_model_validation, p_model_validation_vcmax25), # obs: example data from package 
+    drivers = rbind(p_model_drivers, p_model_drivers_vcmax25),       # drivers: example data from package
+    
+    # additional arguments for the cost function
+    par_fixed = c(         # fix parameter value from previous calibration
+      kc_jmax            = 0.8,
+      kphio              = 0.041,
+      kphio_par_a        = 0.0,
+      kphio_par_b        = 16,
+      soilm_thetastar    = 0.6 * 240,  # to recover paper setup with soil moisture stress
+      soilm_betao        = 0.0,
+      beta_unitcostratio = 146.0,
+      rd_to_vcmax        = 0.014,      # value from Atkin et al. 2015 for C3 herbaceous
+      tau_acclim         = 30.0,
+      err_gpp = 0.2, err_vcmax25 = 3.0
+    ), 
+    targets = c('gpp', 'vcmax25')
+  )
+  testthat::expect_equal(tolerance = 0.5, #tolerance = 1e-4,
+                         object = ll_pmodel_fixed,
+                         # expected was generated with dput(ll_pmodel_fixed)
+                         expected = -336583.32327482)
+})
+
+test_that("test likelihood/RMSE calculations with BiomeE", {
+  test_params_BiomeE <- data.frame( # test_params_BiomeE was generated with dput(test_params_BiomeE)
+    phiRL     = c(6.59158648136072, 2.41828079945408, 4.51794087081216, 0.323927985038608), 
+    LAI_light = c(4.83413460890297, 4.89137732107192, 6.25084221335128, 1.65691818702035), 
+    tf_base   = c(0.986252965405583, 1.52580757206306, 0.278885046485811, 0.125027264398523), 
+    par_mort  = c(1.64211843877565, 0.579043845250271, 1.28934027748182, 1.11228716920596), 
+    err_GPP   = c(2.9679689736967, 3.70911861001514, 1.16307689385489, 0.195016647893935)
+  ) # TODO: in BiomeE output is uppercase GPP, but in p-model it is lowercase
+  
+  # test_params_BiomeE was created with:  Test cost_likelihood_biomee()
+  # test_params_BiomeE was created with:  parBiomeE_cal_best <- c(
+  # test_params_BiomeE was created with:    phiRL              = 3.5,
+  # test_params_BiomeE was created with:    LAI_light          = 3.5,
+  # test_params_BiomeE was created with:    tf_base            = 1,
+  # test_params_BiomeE was created with:    par_mort           = 1,
+  # test_params_BiomeE was created with:    err_GPP          = 1
+  # test_params_BiomeE was created with:  )
+  # test_params_BiomeE was created with:  parBiomeE_cal_min <- c(
+  # test_params_BiomeE was created with:    phiRL              = 0.1,
+  # test_params_BiomeE was created with:    LAI_light          = 0.1,
+  # test_params_BiomeE was created with:    tf_base            = 0.1,
+  # test_params_BiomeE was created with:    par_mort           = 0.1,
+  # test_params_BiomeE was created with:    err_GPP          = 0.01
+  # test_params_BiomeE was created with:  )
+  # test_params_BiomeE was created with:  parBiomeE_cal_max <- c(
+  # test_params_BiomeE was created with:    phiRL              = 7.0,
+  # test_params_BiomeE was created with:    LAI_light          = 7.0,
+  # test_params_BiomeE was created with:    tf_base            = 2.0,
+  # test_params_BiomeE was created with:    par_mort           = 2.0,
+  # test_params_BiomeE was created with:    err_GPP          = 4
+  # test_params_BiomeE was created with:  )
+  # test_params_BiomeE was created with:  prior_BiomeE <- createUniformPrior(lower = parBiomeE_cal_min, upper = parBiomeE_cal_max, best = parBiomeE_cal_best)
+  # test_params_BiomeE was created with:  test_params_BiomeE <- prior_BiomeE$sampler(4) |> as.data.frame() |>
+  # test_params_BiomeE was created with:    stats::setNames(nm = names(par_cal_min))
+  
+  # Test rsofun::cost_likelihood_biomee()
+  ll_values_BiomeE <- apply(test_params_BiomeE, 1, function(par_v) { # par_v is a named vector
+    rsofun::cost_likelihood_biomee(    # likelihood cost function from package
+      par = par_v,                     # par: should be a named vector
+      obs = rsofun::biomee_validation, # obs: example data from package
+      drivers = rsofun::biomee_p_model_drivers,
+      targets = c('GPP')) # TODO: in BiomeE output is uppercase GPP, but in p-model it is lowercase
+  })
+  testthat::expect_equal(
+    tolerance = 1e-4, 
+    object = ll_values_BiomeE, 
+    # expected was generated with dput(ll_values_BiomeE)
+    expected = c(
+      -2.02016968026715,
+      -2.23968996093749,
+      -1.24162843355044,
+      -0.383538022003261
+    )
+  )
+  
+  
+  # Test rsofun::cost_rmse_biomee()
+  relError_values_BiomeE <- apply(dplyr::select(test_params_BiomeE, -err_GPP), 1, function(par_v) { # par_v is a named vector
+    rsofun::cost_rmse_biomee(          # cost function for RMSE (actually relative error) from package
+      par = par_v,                     # par: should be a named vector
+      obs = rsofun::biomee_validation, # obs: example data from package
+      drivers = rsofun::biomee_p_model_drivers)
+  })
+  testthat::expect_equal(
+    tolerance = 1e-4,
+    object = relError_values_BiomeE,
+    # expected was generated with dput(relError_values_BiomeE)
+    # NOTE: these errors are relative Errors, not RMSE:
+    expected = c(
+      0.435809001728683, 
+      0.174141992595677, 
+      0.152546040596932, 
+      0.316161901910987
+    )
+  )
+})
+
+
diff --git a/vignettes/new_cost_function.Rmd b/vignettes/new_cost_function.Rmd
index 137bd123d..70fd336fa 100644
--- a/vignettes/new_cost_function.Rmd
+++ b/vignettes/new_cost_function.Rmd
@@ -19,22 +19,21 @@ library(dplyr)
 library(ggplot2)
 ```
 
-The `rsofun` package allows to calibrate parameters of the `pmodel` and `biomee` models via the `calib_sofun()` function. The implementation of the calibration is fairly flexible and can be adapted to a specific use-case via a cost function (used as metrics for the optimization routines in `calib_sofun()`). The package provides a set of standard cost functions named `cost_*`, which can be used for a variety of calibrations (different sets of model parameters, using various target variables, etc.). Alternatively, it's possible to write a more specific new cost function to be used together with `calib_sofun()`.
+The `rsofun` package allows to calibrate parameters of the `pmodel` and `biomee` models via the `calib_sofun()` function. The implementation of the calibration is fairly flexible and can be adapted to a specific use-case via a tailor-made cost function (used as metrics for the optimization routines in `calib_sofun()`). The package provides a set of standard cost functions named `cost_*`, which can be used for a variety of calibrations (different sets of model parameters, using various target variables, etc.). Alternatively, it's possible to write a more specific new cost function to be used together with `calib_sofun()`.
 
-In this vignette, we go over some examples on how to use the `rsofun` cost functions for parameter calibration and how to write your own custom one from scratch.
+In this vignette, we go over some examples on how to use the `rsofun` cost functions for parameter calibration with `calib_sofun()` and how to write your own custom one from scratch.
 
 ### Calibration to GPP using RMSE and GenSA optimizer
 
 A simple approach to parameter calibration is to find the parameter values that lead to the best prediction performance, in terms of the RMSE (root mean squared error). The function `cost_rmse_pmodel()` runs the P-model internally to calculate the RMSE between predicted target values (in this case GPP) and the corresponding observations.
 
-The implementation of `cost_rmse_pmodel()` allows flexibility in various ways. We can simultaneously calibrate a subset of model parameters and also replicate the different calibration setups in Stocker et al., 2020 GMD. For example, following the `ORG` setup, only parameter `kphio` is calibrated. Furthermore, the standard cost functions allow to calibrate to several targets (fluxes and leaf traits predicted by the P-model) simultaneously and to parallelize the simulations. Since the P-model is run internally to make predictions, we must always specify which values the model parameters should take, i.e. the parameters that aren't calibrated (via argument `par_fixed`). 
+The implementations of `cost_rmse_pmodel()` and `calib_sofun()` allows flexibility in various ways. We can simultaneously calibrate a subset of model parameters and also replicate the different calibration setups in Stocker et al., 2020 GMD, simply by providing the appropriate inputs as `settings`, `par_fixed` and `targets` to `calib_sofun()`. Since the P-model is run internally to make predictions, we must always specify the values of the model parameters that aren't calibrated (via argument `par_fixed`). For example, following the `ORG` setup, only parameter `kphio` is calibrated with below code:
 
-The syntax to run the calibration routine is as follows:
-```{r eval = FALSE}
+```{r run GenSA calibration, eval = TRUE}
 # Define calibration settings and parameter ranges from previous work
 settings_rmse <- list(
   method = 'GenSA',                   # minimizes the RMSE
-  metric = cost_rmse_pmodel,          # our cost function
+  metric = cost_rmse_pmodel,          # our cost function returning the RMSE
   control = list(                     # control parameters for optimizer GenSA
     maxit = 100),                     
   par = list(                         # bounds for the parameter space
@@ -63,16 +62,20 @@ pars_calib_rmse <- calib_sofun(
   ),
   targets = "gpp"           # define target variable GPP
 )
+
+pars_calib_rmse
 ```
 
-The output of `calib_sofun()` is a list containing the calibrated parameter values and the raw optimization output from the optimizer (here from `GenSA` or, as we see next, from `BayesianTools::runMCMC`).
+The output of `calib_sofun()` is a list containing the calibrated parameter values (element `par`) and the raw optimization output from the optimizer (element `mod`; here from `GenSA` or, as we see next, from `BayesianTools::runMCMC`).
+
+Note that the standard cost functions allow to calibrate to several targets (fluxes and leaf traits predicted by the P-model) simultaneously and to parallelize the simulations. 
 
 ### Calibration to GPP using a simple likelihood function and BayesianTools
 
-Let's calibrate the parameters involved in the temperature dependency of the quantum yield efficiency, `kphio`, `kphio_par_a` and `kphio_par_b`, taking a Bayesian calibration approach. We assume that the target variable (`'gpp'`) follows a normal distribution centered at the observations and with its standard deviation being a new calibratable parameter (`'err_gpp'`). We also assume a uniform prior distribution for all calibratable parameters.
-By maximizing the normal log-likelihood, the MAP (maximum a posteriori) estimators for all 4 parameters are computed. With the function `cost_likelihood_pmodel()`, we can easily perform this calibration, as follows:
+Let's calibrate the parameters involved in the temperature dependency of the quantum yield efficiency, `kphio`, `kphio_par_a` and `kphio_par_b`. Taking a Bayesian calibration approach, we need to define a likelihood as cost function (we'll use `cost_likelihood_pmodel()`). We assume that the target variable (`'gpp'`) follows a normal distribution centered at the observations and with its _unknown_ standard deviation (`'err_gpp'`), that we need to add to the calibratable parameters. We also assume a uniform prior distribution for all calibratable parameters.
+By maximizing the log-likelihood, the MAP (maximum a posteriori) estimators for all 4 parameters are computed. With the functions `cost_likelihood_pmodel()` and `calib_sofun()`, we can easily perform this type of calibration:
 
-```{r eval = FALSE}
+```{r run Bayesian calibration, eval = TRUE}
 # Define calibration settings
 settings_likelihood <- list(
   method = 'BayesianTools',
@@ -109,15 +112,17 @@ pars_calib_likelihood <- calib_sofun(
   ),
   targets = "gpp"
 )
+
+pars_calib_likelihood
 ```
 
 Furthermore, there are equivalent cost functions available for the BiomeE model. Check out the reference pages for more details on how to use `cost_likelihood_biomee()` and `cost_rmse_biomee()`.
 
 ### Calibration to GPP and Vcmax25 using the joint log-likelihood and BayesianTools
 
-You may be interested in calibrating the model to different target variables simultaneously, like flux and leaf trait measurements. Here we present an example, where we use `cost_likelihood_pmodel()` to compute the joint normal likelihood of all the targets specified (that is, by summing the log-likelihoods of GPP and Vcmax25) and ultimately calibrate the `kc_jmax` parameter. It would be possible to follow this workflow for several-target calibration also with RMSE as the optimization metric, using `cost_rmse_pmodel()` and `GenSA` optimization.
+You may be interested in calibrating the model to different target variables simultaneously, like flux and leaf trait measurements. Here we present an example, where we use `cost_likelihood_pmodel()` to compute the joint likelihood of all the targets specified (that is, by summing the log-likelihoods of GPP and Vcmax25) and ultimately calibrate the `kc_jmax` parameter. It would be possible to follow this workflow for several-target calibration also with RMSE as the optimization metric, using `cost_rmse_pmodel()` and `GenSA` optimization.
 
-```{r eval = FALSE}
+```{r run concatenated calibration, eval = FALSE}
 # Define calibration settings for two targets
 settings_joint_likelihood <- list(
   method = "BayesianTools",
@@ -153,42 +158,98 @@ par_calib_join <- calib_sofun(
   ),    
   targets = c('gpp', 'vcmax25')
 )
+par_calib_join
 ```
 
-Note that GPP predictions are directly compared to GPP observations on that day, but Vcmax25 predicted by the P-model (being a leaf trait) is averaged over the growing season and compared to a single Vcmax25 observation taken per site. The cost functions provided in the package tell apart fluxes and leaf traits by the presence of a `"date"` column in the nested validation data frames `p_model_validation` and `p_model_validation_vcmax25`.
+Note that GPP predictions are directly compared to GPP observations on that day, but Vcmax25 predicted by the P-model (being a leaf trait) is averaged over the growing season and compared to a single Vcmax25 observation taken per site. 
+
+The cost functions provided in the package tell apart fluxes and leaf traits by the presence of a `"date"` column in the nested validation data frames `p_model_validation` and `p_model_validation_vcmax25`.
 
 ### Write your custom cost function 
 
-If the RMSE or normal log-likelihood (for one or several targets) cost functions that we provide do not fit your use case, you can easily write a custom one. In this section, we drive you through the main ideas with an example.
+If the RMSE or log-likelihood (for one or several targets) cost functions that we provide do not fit your use case, you can easily write a custom one. In this section, we drive you through the main ideas with an example.
+To run the calibration, you can still use `calib_sofun()` in combination with 
+your custom cost function. 
+
+The routine `calib_sofun()` requires `drivers`, `obs` 
+and `settings` as mandatory arguments. These provide data.frames with driver and
+observational data, as well as settings for the calibration.
+The optional argument `optim_out` defines if the raw optimization output should 
+be returned.
+All other (optional) arguments to `calib_sofun()` are passed through to the cost 
+function (e.g. `par_fixed` in above example). 
+They can be used freely inside of your custom cost function, e.g. to control 
+the simulation setup or the processing. On top of these optional arguments, it 
+is also possible to extend the `drivers` and `obs` data.frames with additional 
+columns that can be used freely for fine-grained control within your custom cost 
+function.
 
 All cost functions must take at least three arguments:
 
-* `par`: A vector of calibratable model parameters. In each iteration of the optimization, a new set of values of `par` is used to run the model and compute the cost.
+* `par`: A named vector of calibratable model parameters. In each iteration of the optimization, a new set of values of `par` is used to run the model and compute the cost.
 
 * `obs`: A data frame of observations, against which to compare the simulation results.
 
 * `drivers`: A data frame of driver data, used to run the simulations.
 
-* Additional optional arguments can be used, like for example the model parameter values that should be fixed across simulations, etc.
+* Additional optional arguments can be used, An example would be model parameter values that should be fixed across simulations, etc.
 
-Since we are calibrating the parameters based on model outputs, the cost function runs the P-model and compare its output to observed validation data.
-```{r, eval = FALSE}
-function(par, obs, drivers){
+Below we'll walk you through the definition of a custom cost function.
+In this example, we'll calibrate the soil moisture stress parameters and use the 
+mean absolute error (MAE) as custom cost function.
+
+Since we are calibrating the parameters based on model outputs, the cost 
+function will eventually need to run the P-model and compare its output to 
+observed validation data. 
+
+To get started we suggest to write a dummy cost function and use it together 
+with `calib_sofun()` as shown below. Note that one way of developing the cost
+function would be to use a `browser()` statement during the development. It 
+allows you to explore the variables that you have access to from within the cost
+function.
+```{r eval = FALSE}
+# Define the custom cost function to be used
+cost_mae <- function(par, obs, drivers, my_own_message){
   # Your code
+  browser() # can facilitate the development, remove afterwards
 }
+
+# Define calibration settings and parameter ranges
+settings_mae <- list(
+  method = 'GenSA',
+  metric = cost_mae, # directly uses the custom cost function
+  control = list(
+    maxit = 100
+  ),
+  par = list(
+    soilm_thetastar = list(lower=0.0, upper=3000, init=0.6*240),
+    soilm_betao = list(lower=0, upper=1, init=0.2)
+  )
+)
+
+# Calibrate the model and optimize the free parameters
+pars_calib_mae <- calib_sofun(
+  drivers = p_model_drivers,
+  obs = p_model_validation,
+  settings = settings_mae,
+  # optional arguments if needed in the cost function
+  my_own_message = "Hi from inside the cost_mae function."
+)
+
+pars_calib_mae
 ```
 
-In the optimization procedure, the cost function only takes as argument the parameters `par` that are fed to `calib_sofun()` via `settings$par` (see previous sections). Nevertheless, within the cost function we call `runread_pmodel_f()` and this function needs a full set of model parameters. Therefore, the parameters that aren't being calibrated must be hard coded inside the cost function (or passed as an argument like the `rsofun` cost functions). In this example, we only want to calibrate the soil moisture stress parameters.
+During the optimization procedure, the cost function receives as argument a suggestion of the parameters `par`. This might be just a subset of all needed parameters (defined via `settings$par`). Thus to call `runread_pmodel_f()` within the cost function, a full set of model parameters is needed. Here, we'll hardcode the parameters that aren't being calibrated inside the cost function. (Note, that in above examples they were passed as an additional argument `par_fixed`).
 ```{r, eval = FALSE}
-function(par, obs, drivers){
+cost_mae <- function(par, obs, drivers, my_own_message){
   
   # Set values for the list of calibrated and non-calibrated model parameters
   params_modl <- list(
     kphio              = 0.09423773,
     kphio_par_a        = 0.0,
     kphio_par_b        = 25,
-    soilm_thetastar    = par[1],
-    soilm_betao        = par[2],
+    soilm_thetastar    = par[["soilm_thetastar"]],
+    soilm_betao        = par[["soilm_betao"]],
     beta_unitcostratio = 146.0,
     rd_to_vcmax        = 0.014,
     tau_acclim         = 30.0,
@@ -204,11 +265,13 @@ function(par, obs, drivers){
   )
   
   # Your code to compute the cost
+  print(my_own_message) # useless, but showcases how to use additional arguments
+  browser() # can facilitate the development, remove afterwards
 }
 ``` 
 
-The following chunk defines the final function. We clean the observations and model output and align the data according to site and date, to compute the mean absolute error (MAE) on GPP. Finally, the function should return a scalar value, in this case the MAE, which we want to minimize. Keep in mind that the GenSA optimization will minimize the cost, but with the BayesianTools method the cost is always maximized.
-```{r}
+The following chunk defines the final function. We clean the observations and model output and align the data according to site and date, to compute the mean absolute error (MAE) on GPP. Finally, the function should return a scalar value, in this case the MAE, which we want to minimize. Keep in mind that the GenSA optimization will minimize the cost, but with the BayesianTools method the cost (i.e. the likelihood) is always maximized.
+```{r define custom cost function, eval = TRUE}
 cost_mae <- function(par, obs, drivers){
 
   # Set values for the list of calibrated and non-calibrated model parameters
@@ -216,8 +279,8 @@ cost_mae <- function(par, obs, drivers){
     kphio              = 0.09423773,
     kphio_par_a        = 0.0,
     kphio_par_b        = 25,
-    soilm_thetastar    = par[1],
-    soilm_betao        = par[2],
+    soilm_thetastar    = par[["soilm_thetastar"]],
+    soilm_betao        = par[["soilm_betao"]],
     beta_unitcostratio = 146.0,
     rd_to_vcmax        = 0.014,
     tau_acclim         = 30.0,
@@ -229,7 +292,7 @@ cost_mae <- function(par, obs, drivers){
   df <- runread_pmodel_f(
     drivers = drivers,
     par = params_modl,
-    makecheck = TRUE,
+    makecheck = FALSE,
     parallel = FALSE
   )
   
@@ -252,11 +315,12 @@ cost_mae <- function(par, obs, drivers){
   
   # Return the computed cost
   return(cost)
+  # browser() # can facilitate the development, remove afterwards
 }
 ``` 
 
 As a last step, let's verify that the calibration procedure runs using this cost function.
-```{r eval = FALSE}
+```{r run custom calibration, eval=TRUE}
 # Define calibration settings and parameter ranges
 settings_mae <- list(
   method = 'GenSA',
@@ -275,4 +339,7 @@ pars_calib_mae <- calib_sofun(
   obs = p_model_validation,
   settings = settings_mae
 )
-```
\ No newline at end of file
+
+pars_calib_mae
+```
+
diff --git a/vignettes/sensitivity_analysis.Rmd b/vignettes/sensitivity_analysis.Rmd
index 4bf0793b0..a19ffdbd7 100644
--- a/vignettes/sensitivity_analysis.Rmd
+++ b/vignettes/sensitivity_analysis.Rmd
@@ -56,7 +56,7 @@ We want to see how sensitive this function is to changes in the parameter values
 ```{r}
 # Define log-likelihood function
 ll_pmodel <- function(
-    par_v                 # a vector of all calibratable parameters including errors
+    par_v                 # a named vector of all calibratable parameters including errors
 ){
   rsofun::cost_likelihood_pmodel(        # reuse likelihood cost function
     par_v,