From c297b8d43fc312458ff318c6e01a3044a9000f1d Mon Sep 17 00:00:00 2001
From: Dongchen Zhang <zhangdc@scc-gl2.scc.bu.edu>
Date: Thu, 13 Apr 2023 13:56:11 -0400
Subject: [PATCH 01/32] Adding qsub_parallel module into pecan.

---
 base/workflow/NAMESPACE                       |   2 +
 base/workflow/R/merge_job_files.R             |  42 +++++++
 base/workflow/R/qsub_parallel.R               | 103 ++++++++++++++++++
 base/workflow/man/merge_job_files.Rd          |  26 +++++
 base/workflow/man/qsub_parallel.Rd            |  26 +++++
 .../assim.sequential/R/sda.enkf_MultiSite.R   |   4 +-
 6 files changed, 202 insertions(+), 1 deletion(-)
 create mode 100644 base/workflow/R/merge_job_files.R
 create mode 100644 base/workflow/R/qsub_parallel.R
 create mode 100644 base/workflow/man/merge_job_files.Rd
 create mode 100644 base/workflow/man/qsub_parallel.Rd

diff --git a/base/workflow/NAMESPACE b/base/workflow/NAMESPACE
index e8a4c246a4c..62b40898236 100644
--- a/base/workflow/NAMESPACE
+++ b/base/workflow/NAMESPACE
@@ -2,7 +2,9 @@
 
 export(create_execute_test_xml)
 export(do_conversions)
+export(merge_job_files)
 export(model_specific_tags)
+export(qsub_parallel)
 export(run.write.configs)
 export(runModule.get.trait.data)
 export(runModule.run.write.configs)
diff --git a/base/workflow/R/merge_job_files.R b/base/workflow/R/merge_job_files.R
new file mode 100644
index 00000000000..23c84ebe09f
--- /dev/null
+++ b/base/workflow/R/merge_job_files.R
@@ -0,0 +1,42 @@
+#' merge_job_files
+#' Merge multiple job.sh files into one larger file.
+#' @param settings PEcAn.settings object with host section.
+#' @param jobs_per_file the number of files you want to merge.
+#' @param outdir output directory of merged job files.
+#'
+#' @return
+#' @export
+#' @author Dongchen Zhang
+#'
+#' @examples
+merge_job_files <- function(settings, jobs_per_file = 10, outdir = NULL){
+  # default outdir
+  if(is.null(outdir)){
+    outdir <- file.path(settings$rundir, "merged_jobs")
+  }
+  # create folder or delete previous job files.
+  if(dir.exists(outdir)){
+    unlink(list.files(outdir, recursive = T, full.names = T))
+  }else{
+    dir.create(outdir)
+  }
+  # merge job files.
+  run_list <- readLines(con = file.path(settings$rundir, "runs.txt"))
+  job_file_paths <- list.files(settings$host$rundir, pattern = "*.sh", recursive = T, full.names = T)
+  i <- 0
+  files <- c()
+  while (i < length(job_file_paths)) {
+    jobs <- c()
+    for (j in 1:jobs_per_file) {
+      if((i+j) > length(job_file_paths)){
+        break
+      }
+      jobs <- c(jobs, readLines(job_file_paths[i+j]))
+    }
+    writeLines(jobs, con = file.path(outdir, paste0("job_", i,".sh")))
+    Sys.chmod(file.path(outdir, paste0("job_", i,".sh")))
+    files <- c(files, file.path(outdir, paste0("job_", i,".sh")))
+    i <- i + jobs_per_file
+  }
+  files
+}
diff --git a/base/workflow/R/qsub_parallel.R b/base/workflow/R/qsub_parallel.R
new file mode 100644
index 00000000000..6c33d2b764f
--- /dev/null
+++ b/base/workflow/R/qsub_parallel.R
@@ -0,0 +1,103 @@
+#' qsub_parallel
+#'
+#' @param settings pecan settings object
+#' @param files allow submit jobs based on job.sh file paths.
+#' @param prefix used for detecting if jobs are completed or not.
+#' @export
+#' @examples
+#' \dontrun{
+#'   qsub_parallel(settings)
+#' }
+#' @author Dongchen Zhang
+#'
+qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out") {
+  run_list <- readLines(con = file.path(settings$rundir, "runs.txt"))
+  is_local <- PEcAn.remote::is.localhost(settings$host)
+  # loop through runs and either call start run, or launch job on remote machine
+  #parallel submit jobs
+  cores <- parallel::detectCores()
+  cl <- makeSOCKcluster(cores)
+  registerDoSNOW(cl)
+  #progress bar
+  pb <- txtProgressBar(min=1, max=ifelse(is.null(files), length(run_list), length(files)), style=3)
+  progress <- function(n) setTxtProgressBar(pb, n)
+  opts <- list(progress=progress)
+  PEcAn.logger::logger.info("Submitting jobs!")
+  # if we want to submit jobs separately.
+  if(is.null(files)){
+    jobids <- foreach(run = run_list, .packages="Kendall", .options.snow=opts, settings = rep(settings, length(files))) %dopar% {
+      run_id_string <- format(run, scientific = FALSE)
+      qsub <- settings$host$qsub
+      qsub <- gsub("@NAME@", paste0("PEcAn-", run_id_string), qsub)
+      qsub <- gsub("@STDOUT@", file.path(settings$host$outdir, run_id_string, "stdout.log"), qsub)
+      qsub <- gsub("@STDERR@", file.path(settings$host$outdir, run_id_string, "stderr.log"), qsub)
+      qsub <- strsplit(qsub, " (?=([^\"']*\"[^\"']*\")*[^\"']*$)", perl = TRUE)
+      # start the actual model run
+      cmd <- qsub[[1]]
+      if(PEcAn.remote::is.localhost(settings$host)){
+        out <- system2(cmd, file.path(settings$host$rundir, run_id_string, "job.sh"), stdout = TRUE, stderr = TRUE)
+      }else{
+        out <- remote.execute.cmd(host, cmd, file.path(settings$host$rundir, run_id_string, "job.sh"), stderr = TRUE)
+      }
+      jobid <- PEcAn.remote::qsub_get_jobid(
+        out = out[length(out)],
+        qsub.jobid = settings$host$qsub.jobid,
+        stop.on.error = TRUE)
+      return(jobid)
+    }
+  }else{
+    # if we want to submit merged job files.
+    std_out <- file.path(settings$host$outdir, "merged_stdout")
+    if(!dir.exists(std_out)){
+      dir.create(std_out)
+    }else{
+      unlink(list.files(std_out, recursive = T, full.names = T))
+    }
+    jobids <- foreach(file = files, .packages="Kendall", .options.snow=opts, settings = rep(settings, length(files))) %dopar% {
+      qsub <- settings$host$qsub
+      base_name <- basename(file)
+      num <- gsub("\\D", "", base_name)
+      name <- paste0("SDA", num)
+      qsub <- gsub("@NAME@", name, qsub)
+      qsub <- gsub("@STDOUT@", file.path(std_out, paste0("stdout", num, ".log")), qsub)
+      qsub <- gsub("@STDERR@", file.path(std_out, paste0("stderr", num, ".log")), qsub)
+      qsub <- strsplit(qsub, " (?=([^\"']*\"[^\"']*\")*[^\"']*$)", perl = TRUE)
+      cmd <- qsub[[1]]
+      if(PEcAn.remote::is.localhost(settings$host)){
+        out <- system2(cmd, file, stdout = TRUE, stderr = TRUE)
+      }else{
+        out <- remote.execute.cmd(host, cmd, file, stderr = TRUE)
+      }
+      jobid <- PEcAn.remote::qsub_get_jobid(
+        out = out[length(out)],
+        qsub.jobid = settings$host$qsub.jobid,
+        stop.on.error = TRUE)
+      return(jobid)
+    }
+  }
+  
+  PEcAn.logger::logger.info("Jobs submitted!")
+  #check if jobs are completed
+  PEcAn.logger::logger.info("Checking the qsub jobs status!")
+  ## setup progressbar
+  pb <- utils::txtProgressBar(min = 0, max = length(unlist(run_list)), style = 3)
+  pbi <- 0
+  folders <- file.path(settings$host$outdir, run_list)
+  while (length(folders) > 0) {
+    Sys.sleep(10)
+    completed_folders <- foreach(folder = folders, settings = rep(settings, length(run_list))) %dopar% {
+      if(file.exists(file.path(folder, prefix))){
+        return(folder)
+      }
+    }
+    if(length(unlist(completed_folders)) > 0){
+      Ind <- which(unlist(completed_folders) %in% folders)
+      folders <- folders[-Ind]
+      pbi <- pbi + length(completed_folders)
+      utils::setTxtProgressBar(pb, pbi)
+    }
+  }
+  close(pb)
+  stopCluster(cl)
+  PEcAn.logger::logger.info("Completed!")
+}
\ No newline at end of file
diff --git a/base/workflow/man/merge_job_files.Rd b/base/workflow/man/merge_job_files.Rd
new file mode 100644
index 00000000000..b0cb1b096e6
--- /dev/null
+++ b/base/workflow/man/merge_job_files.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/merge_job_files.R
+\name{merge_job_files}
+\alias{merge_job_files}
+\title{merge_job_files
+Merge multiple job.sh files into one larger file.}
+\usage{
+merge_job_files(settings, jobs_per_file = 10, outdir = NULL)
+}
+\arguments{
+\item{settings}{PEcAn.settings object with host section.}
+
+\item{jobs_per_file}{the number of files you want to merge.}
+
+\item{outdir}{output directory of merged job files.}
+}
+\value{
+
+}
+\description{
+merge_job_files
+Merge multiple job.sh files into one larger file.
+}
+\author{
+Dongchen Zhang
+}
diff --git a/base/workflow/man/qsub_parallel.Rd b/base/workflow/man/qsub_parallel.Rd
new file mode 100644
index 00000000000..cbd2c6614d4
--- /dev/null
+++ b/base/workflow/man/qsub_parallel.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/qsub_parallel.R
+\name{qsub_parallel}
+\alias{qsub_parallel}
+\title{qsub_parallel}
+\usage{
+qsub_parallel(settings, files = NULL, prefix = "sipnet.out")
+}
+\arguments{
+\item{settings}{pecan settings object}
+
+\item{files}{allow submit jobs based on job.sh file paths.}
+
+\item{prefix}{used for detecting if jobs are completed or not.}
+}
+\description{
+qsub_parallel
+}
+\examples{
+\dontrun{
+  qsub_parallel(settings)
+}
+}
+\author{
+Dongchen Zhang
+}
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index d5975c38c7f..07f44aa3604 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -420,7 +420,8 @@ sda.enkf.multisite <- function(settings,
         runs.tmp <- list.dirs(rundir, full.names = F)
         runs.tmp <- runs.tmp[grepl("ENS-*|[0-9]", runs.tmp)] 
         writeLines(runs.tmp[runs.tmp != ''], file.path(rundir, 'runs.txt'))
-        PEcAn.workflow::start_model_runs(settings, write=settings$database$bety$write)
+        # PEcAn.workflow::start_model_runs(settings, write=settings$database$bety$write)
+        PEcAn.workflow::qsub_parallel(settings, files=PEcAn.workflow::merge_job_files(settings, 20))
         
         #------------- Reading - every iteration and for SDA
         
@@ -711,5 +712,6 @@ sda.enkf.multisite <- function(settings,
     {
       unlink(list.files(outdir, "*.nc", recursive = TRUE, full.names = TRUE))
     }
+    PEcAn.utils::sendmail("zhangdc@bu.edu", "zhangdc@bu.edu", "SDA progress report", paste("Time point:", obs.times[t], "has been completed!"))
   } ### end loop over time
 } # sda.enkf
\ No newline at end of file

From ff0a52a3209e14b2d86679de7829243fa82eae05 Mon Sep 17 00:00:00 2001
From: Dongchen Zhang <zhangdc@scc-gb1.scc.bu.edu>
Date: Mon, 24 Apr 2023 12:11:55 -0400
Subject: [PATCH 02/32] Making the function more robust.

---
 .../assim.sequential/R/sda.enkf_MultiSite.R   | 30 +++++++++++--------
 1 file changed, 18 insertions(+), 12 deletions(-)

diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index 07f44aa3604..f558d6cd868 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -341,7 +341,7 @@ sda.enkf.multisite <- function(settings,
         load(file.path(settings$outdir, "samples.Rdata"))
       }
       #reformatting params
-      new.params <- sda_matchparam(settings, ensemble.samples, site.ids, nens)
+      new.params <- PEcAnAssimSequential:::sda_matchparam(settings, ensemble.samples, site.ids, nens)
     }
       #sample met ensemble members
       inputs <- conf.settings %>% purrr::map(function(setting) {
@@ -421,21 +421,25 @@ sda.enkf.multisite <- function(settings,
         runs.tmp <- runs.tmp[grepl("ENS-*|[0-9]", runs.tmp)] 
         writeLines(runs.tmp[runs.tmp != ''], file.path(rundir, 'runs.txt'))
         # PEcAn.workflow::start_model_runs(settings, write=settings$database$bety$write)
-        PEcAn.workflow::qsub_parallel(settings, files=PEcAn.workflow::merge_job_files(settings, 20))
+        PEcAn.workflow::qsub_parallel(settings, files=PEcAn.workflow::merge_job_files(settings, 10), prefix = paste0(obs.year, ".nc"))
         
         #------------- Reading - every iteration and for SDA
         
         #put building of X into a function that gets called
-        reads <- build_X(out.configs = out.configs, 
-                         settings = settings, 
-                         new.params = new.params, 
-                         nens = nens, 
-                         read_restart_times = read_restart_times, 
-                         outdir = outdir, 
-                         t = t, 
-                         var.names = var.names, 
-                         my.read_restart = my.read_restart,
-                         restart_flag = restart_flag)
+        while (is.character(try(reads <- PEcAnAssimSequential:::build_X(out.configs = out.configs, 
+                                                                        settings = settings, 
+                                                                        new.params = new.params, 
+                                                                        nens = nens, 
+                                                                        read_restart_times = read_restart_times, 
+                                                                        outdir = outdir, 
+                                                                        t = t, 
+                                                                        var.names = var.names, 
+                                                                        my.read_restart = my.read_restart,
+                                                                        restart_flag = restart_flag), silent = T))) {
+          Sys.sleep(5)
+          PEcAn.logger::logger.info("Empty folder, try again!")
+        }
+        
         
         if (control$debug) browser()
         #let's read the parameters of each site/ens
@@ -500,6 +504,8 @@ sda.enkf.multisite <- function(settings,
         }
         if(!is.null(pre_enkf_params)){
           Pf <- pre_enkf_params[[t]]$Pf
+        }else{
+          Pf <- NULL
         }
         
         if(!exists('Cmcmc_tobit2space')) {

From b8f97d3f5be6bdb153cf1bc12e5529fbf1d4fe4c Mon Sep 17 00:00:00 2001
From: Dongchen Zhang <zhangdc@scc-gb1.scc.bu.edu>
Date: Mon, 24 Apr 2023 12:12:15 -0400
Subject: [PATCH 03/32] Adding the gamma parameters' names.

---
 modules/assim.sequential/R/Nimble_codes.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/modules/assim.sequential/R/Nimble_codes.R b/modules/assim.sequential/R/Nimble_codes.R
index 38929c17e25..798e96c6cbd 100644
--- a/modules/assim.sequential/R/Nimble_codes.R
+++ b/modules/assim.sequential/R/Nimble_codes.R
@@ -178,7 +178,7 @@ tobit.model <-  nimbleCode({
 GEF.MultiSite.Nimble <-  nimbleCode({
   if (q.type == 1) {
     # Sorting out qs
-    qq ~ dgamma(aq, bq) ## aq and bq are estimated over time
+    qq ~ dgamma(shape = aq, rate = bq) ## aq and bq are estimated over time
     q[1:YN, 1:YN] <- qq * diag(YN)
   } else if (q.type == 2) {
     # Sorting out qs

From 4333f11958508d1f4feb821f701e6a0a42503b77 Mon Sep 17 00:00:00 2001
From: Dongchen Zhang <zhangdc@scc-gb1.scc.bu.edu>
Date: Mon, 24 Apr 2023 12:12:41 -0400
Subject: [PATCH 04/32] Correct the way of updating Q.

---
 modules/assim.sequential/R/Analysis_sda_multiSite.R | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/modules/assim.sequential/R/Analysis_sda_multiSite.R b/modules/assim.sequential/R/Analysis_sda_multiSite.R
index d832a7b41d8..c27a0a4d670 100644
--- a/modules/assim.sequential/R/Analysis_sda_multiSite.R
+++ b/modules/assim.sequential/R/Analysis_sda_multiSite.R
@@ -496,8 +496,9 @@ GEF.MultiSite<-function(settings, Forecast, Observed, H, extraArg,...){
   # Setting up the prior for the next step from the posterior of this step
   if (t < nt){
     if (q.type == single.q){ #if it's a gamma case
-      aqq[1, 1, t + 1] <- mean(mq)
-      bqq[t + 1] <- stats::var(mq  %>%  as.numeric())
+      qq <- dat[,  grep("qq", colnames(dat))]
+      aqq[1, 1, t + 1] <- (mean(qq))^2/stats::var(qq)
+      bqq[t + 1] <- mean(qq)/stats::var(qq)
     } else { # if it's a wish case
       col <- matrix(1:length(elements.W.Data) ^ 2,
                     length(elements.W.Data),

From 4f4630b7575db2241f72a0abc469724b5874ed57 Mon Sep 17 00:00:00 2001
From: Dongchen Zhang <zhangdc@scc-gb1.scc.bu.edu>
Date: Mon, 24 Apr 2023 12:12:56 -0400
Subject: [PATCH 05/32] Avoid extremely high SoilC.

---
 models/sipnet/R/write.configs.SIPNET.R | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 91c5224c9a0..92cefc7cebe 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -441,7 +441,11 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     plant_wood_vars <- c("AbvGrndWood", "abvGrndWoodFrac", "coarseRootFrac", "fineRootFrac")
     if (all(plant_wood_vars %in% ic.names)) {
       # reconstruct total wood C
-      wood_total_C <- IC$AbvGrndWood / IC$abvGrndWoodFrac
+      if(IC$abvGrndWoodFrac < 0.1){
+        wood_total_C <- 0
+      }else{
+        wood_total_C <- IC$AbvGrndWood / IC$abvGrndWoodFrac
+      }
       #Sanity check
       if (is.infinite(wood_total_C) | is.nan(wood_total_C) | wood_total_C < 0) {
         wood_total_C <- 0

From b88db3329641a841783f0c967a2bffea75085b8d Mon Sep 17 00:00:00 2001
From: Dongchen Zhang <zhangdc@scc-gb1.scc.bu.edu>
Date: Mon, 24 Apr 2023 12:13:10 -0400
Subject: [PATCH 06/32] Update.

---
 base/workflow/R/qsub_parallel.R | 1 +
 1 file changed, 1 insertion(+)

diff --git a/base/workflow/R/qsub_parallel.R b/base/workflow/R/qsub_parallel.R
index 6c33d2b764f..0318f479951 100644
--- a/base/workflow/R/qsub_parallel.R
+++ b/base/workflow/R/qsub_parallel.R
@@ -79,6 +79,7 @@ qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out") {
   PEcAn.logger::logger.info("Jobs submitted!")
   #check if jobs are completed
   PEcAn.logger::logger.info("Checking the qsub jobs status!")
+  PEcAn.logger::logger.info(paste0("Checking the file ", prefix))
   ## setup progressbar
   pb <- utils::txtProgressBar(min = 0, max = length(unlist(run_list)), style = 3)
   pbi <- 0

From 657339f7a3f646f05da45a6aa3da372e93aa6405 Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Sat, 24 Jun 2023 01:09:24 -0400
Subject: [PATCH 07/32] Use Alexis's script for reading adjusted soil
 respiration parameter

---
 modules/uncertainty/R/get.parameter.samples.R | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/modules/uncertainty/R/get.parameter.samples.R b/modules/uncertainty/R/get.parameter.samples.R
index 26c136faef0..dec77cdd5d8 100644
--- a/modules/uncertainty/R/get.parameter.samples.R
+++ b/modules/uncertainty/R/get.parameter.samples.R
@@ -58,9 +58,14 @@ get.parameter.samples <- function(settings,
     if (!is.na(posterior.files[i])) {
       # Load specified file
       load(posterior.files[i])
-      if (!exists("prior.distns") & exists("post.distns")) {
+      fname <- file.path(outdirs[i], "post.distns.Rdata")
+      if (file.exists(fname)) {
+        load(fname)
         prior.distns <- post.distns
       }
+      if (exists("trait.mcmc") & exists("post.distns")) {
+        post.distns <- post.distns[!(row.names(post.distns) %in% names(trait.mcmc)),]
+      }
     } else {
       # Default to most recent posterior in the workflow, or the prior if there is none
       fname <- file.path(outdirs[i], "post.distns.Rdata")
@@ -104,7 +109,13 @@ get.parameter.samples <- function(settings,
     
     ### When no ma for a trait, sample from prior
     ### Trim all chains to shortest mcmc chain, else 20000 samples
-    priors <- rownames(prior.distns)
+    if(exists("post.distns") & exists("trait.mcmc")){
+      priors <- c(names(trait.mcmc), rownames(post.distns))
+      prior.distns <- post.distns
+    }else{
+      priors <- rownames(prior.distns)
+    }
+    #priors <- names(trait.mcmc)
     
     if (exists("trait.mcmc")) {
       param.names[[i]] <- names(trait.mcmc)

From f267550c4e3951a49966393fff957e979e7ff206 Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Mon, 17 Jul 2023 00:23:39 -0400
Subject: [PATCH 08/32] Enable monthly SDA

---
 base/workflow/man/merge_job_files.Rd          |  3 --
 models/sipnet/R/model2netcdf.SIPNET.R         |  8 +++--
 models/sipnet/R/write_restart.SIPNET.R        |  2 +-
 models/sipnet/man/model2netcdf.SIPNET.Rd      |  2 +-
 .../assim.sequential/R/sda.enkf_MultiSite.R   | 32 +++++++------------
 .../man/sda.enkf.multisite.Rd                 | 15 ++++-----
 6 files changed, 27 insertions(+), 35 deletions(-)

diff --git a/base/workflow/man/merge_job_files.Rd b/base/workflow/man/merge_job_files.Rd
index b0cb1b096e6..6f6e2f4fd89 100644
--- a/base/workflow/man/merge_job_files.Rd
+++ b/base/workflow/man/merge_job_files.Rd
@@ -13,9 +13,6 @@ merge_job_files(settings, jobs_per_file = 10, outdir = NULL)
 \item{jobs_per_file}{the number of files you want to merge.}
 
 \item{outdir}{output directory of merged job files.}
-}
-\value{
-
 }
 \description{
 merge_job_files
diff --git a/models/sipnet/R/model2netcdf.SIPNET.R b/models/sipnet/R/model2netcdf.SIPNET.R
index 33ca6480dc7..0cc58b8f1eb 100644
--- a/models/sipnet/R/model2netcdf.SIPNET.R
+++ b/models/sipnet/R/model2netcdf.SIPNET.R
@@ -103,7 +103,7 @@ sipnet2datetime <- function(sipnet_tval, base_year, base_month = 1,
 ##' @export
 ##' @author Shawn Serbin, Michael Dietze
 model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, delete.raw = FALSE, revision, prefix = "sipnet.out",
-                                overwrite = FALSE, conflict = FALSE) {
+                                overwrite = FALSE, conflict = TRUE) {
 
   ### Read in model output in SIPNET format
   sipnet_out_file <- file.path(outdir, prefix)
@@ -303,6 +303,10 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
         files <- file.path(outdir, "current.nc")
       }
       mergeNC(files = files, outfile = file.path(outdir, paste(y, "nc", sep = ".")))
+      nc<- ncdf4::nc_open(file.path(outdir, paste(y, "nc", sep = ".")),write=TRUE)
+      nc<-ncdf4::ncvar_rename(nc,"time_bnds","time_bounds")
+      ncdf4::ncatt_put(nc, "time", "bounds","time_bounds", prec=NA)
+      ncdf4::nc_close(nc)
       unlink(files, recursive = T)
     }else{
       nc      <- ncdf4::nc_create(file.path(outdir, paste(y, "nc", sep = ".")), nc_var)
@@ -323,4 +327,4 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
   }
 } # model2netcdf.SIPNET
 #--------------------------------------------------------------------------------------------------#
-### EOF
\ No newline at end of file
+### EOF
diff --git a/models/sipnet/R/write_restart.SIPNET.R b/models/sipnet/R/write_restart.SIPNET.R
index 18387a9f893..fadbd27353c 100755
--- a/models/sipnet/R/write_restart.SIPNET.R
+++ b/models/sipnet/R/write_restart.SIPNET.R
@@ -109,7 +109,7 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
   if ("LAI" %in% variables) {
     analysis.save[[length(analysis.save) + 1]] <- new.state$LAI  
     if (new.state$LAI < 0) analysis.save[[length(analysis.save)]] <- 0
-    names(analysis.save[[length(analysis.save)]]) <- c("lai")
+   names(analysis.save[[length(analysis.save)]]) <- c("lai")
   }
   
   if (!is.null(analysis.save) && length(analysis.save)>0){
diff --git a/models/sipnet/man/model2netcdf.SIPNET.Rd b/models/sipnet/man/model2netcdf.SIPNET.Rd
index 3853c307e6a..697ec32d1bd 100644
--- a/models/sipnet/man/model2netcdf.SIPNET.Rd
+++ b/models/sipnet/man/model2netcdf.SIPNET.Rd
@@ -14,7 +14,7 @@ model2netcdf.SIPNET(
   revision,
   prefix = "sipnet.out",
   overwrite = FALSE,
-  conflict = FALSE
+  conflict = TRUE
 )
 }
 \arguments{
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index 802f8ffc385..0b21c8a60ee 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -38,9 +38,9 @@ sda.enkf.multisite <- function(settings,
                                run_parallel = TRUE,
                                ensemble.samples = NULL,
                                parallel_qsub = TRUE,
-                               free_run = FALSE,
-                               send_email = NULL,
                                control=list(trace = TRUE,
+					    free_run = FALSE,
+					    send_email = NULL,
                                             FF = FALSE,
                                             interactivePlot = FALSE,
                                             TimeseriesPlot = FALSE,
@@ -131,9 +131,9 @@ sda.enkf.multisite <- function(settings,
   #Here I'm trying to make a temp config list name and put it into map to iterate
   if(multi.site.flag){
     conf.settings<-settings
-    site.ids <- conf.settings$run %>% purrr::map('site') %>% purrr::map('id') %>% base::unlist() %>% base::as.character()
+    site.ids <- conf.settings %>% map(~.x[['run']] ) %>% purrr::map('site') %>% purrr::map('id') %>% base::unlist() %>% base::as.character()
     # a matrix ready to be sent to spDistsN1 in sp package - first col is the long second is the lat and row names are the site ids
-    site.locs <- conf.settings$run %>% purrr::map('site') %>% purrr::map_dfr(~c(.x[['lon']],.x[['lat']]) %>% as.numeric)%>% 
+    site.locs <- conf.settings %>% map(~.x[['run']] ) %>% purrr::map('site') %>% purrr::map_dfr(~c(.x[['lon']],.x[['lat']]) %>% as.numeric)%>% 
       t %>%
       `colnames<-`(c("Lon","Lat")) %>%
       `rownames<-`(site.ids)
@@ -441,9 +441,10 @@ sda.enkf.multisite <- function(settings,
         #------------- Reading - every iteration and for SDA
         
         #put building of X into a function that gets called
-        max_t <- 0
-        while("try-error" %in% class(
-          try(reads <- build_X(out.configs = out.configs, 
+        #max_t <- 0
+        #while("try-error" %in% class(
+         # try(
+	reads <-  build_X(out.configs = out.configs, 
                                settings = settings, 
                                new.params = new.params, 
                                nens = nens, 
@@ -452,16 +453,7 @@ sda.enkf.multisite <- function(settings,
                                t = t, 
                                var.names = var.names, 
                                my.read_restart = my.read_restart,
-                               restart_flag = restart_flag), silent = T))
-        ){
-          Sys.sleep(10)
-          max_t <- max_t + 1
-          if(max_t > 20){
-            PEcAn.logger::logger.info("Can't find outputed NC file! Please rerun the code!")
-            break
-          }
-          PEcAn.logger::logger.info("Empty folder, try again!")
-        }
+                               restart_flag = restart_flag)
         
         if (control$debug) browser()
         #let's read the parameters of each site/ens
@@ -721,9 +713,9 @@ sda.enkf.multisite <- function(settings,
       unlink(mailfile)
     }
     # useful for debugging to keep .nc files for assimilated years. T = 2, because this loops removes the files that were run when starting the next loop
-#    if (keepNC && t == 1){
-#      unlink(list.files(outdir, "*.nc", recursive = TRUE, full.names = TRUE))
-#    }
+   if (keepNC && t == 1){
+    unlink(list.files(outdir, "*.nc", recursive = TRUE, full.names = TRUE))
+   }
       ## MCD: I commented the above "if" out because if you are restarting from a previous forecast, this might delete the files in that earlier folder
   } ### end loop over time
 } # sda.enkf
diff --git a/modules/assim.sequential/man/sda.enkf.multisite.Rd b/modules/assim.sequential/man/sda.enkf.multisite.Rd
index 3f216df4803..81e6477faa8 100644
--- a/modules/assim.sequential/man/sda.enkf.multisite.Rd
+++ b/modules/assim.sequential/man/sda.enkf.multisite.Rd
@@ -16,11 +16,10 @@ sda.enkf.multisite(
   run_parallel = TRUE,
   ensemble.samples = NULL,
   parallel_qsub = TRUE,
-  free_run = FALSE,
-  send_email = NULL,
-  control = list(trace = TRUE, FF = FALSE, interactivePlot = FALSE, TimeseriesPlot =
-    FALSE, BiasPlot = FALSE, plot.title = NULL, facet.plots = FALSE, debug = FALSE, pause
-    = FALSE, Profiling = FALSE, OutlierDetection = FALSE),
+  control = list(trace = TRUE, free_run = FALSE, send_email = NULL, FF = FALSE,
+    interactivePlot = FALSE, TimeseriesPlot = FALSE, BiasPlot = FALSE, plot.title = NULL,
+    facet.plots = FALSE, debug = FALSE, pause = FALSE, Profiling = FALSE,
+    OutlierDetection = FALSE),
   ...
 )
 }
@@ -47,12 +46,12 @@ sda.enkf.multisite(
 
 \item{parallel_qsub}{Bool variable decide if you want to submit the bash jobs under the parallel mode, the default value is TRUE.}
 
+\item{control}{List of flags controlling the behaviour of the SDA. trace for reporting back the SDA outcomes, interactivePlot for plotting the outcomes after each step, 
+TimeseriesPlot for post analysis examination, BiasPlot for plotting the correlation between state variables, plot.title is the title of post analysis plots and debug mode allows for pausing the code and examining the variables inside the function.}
+
 \item{free_run}{Bool variable decide if the run is a free run with no analysis been used, the default value is FALSE.}
 
 \item{send_email}{List object containing the "from", "to", and "body", the default value is NULL.}
-
-\item{control}{List of flags controlling the behaviour of the SDA. trace for reporting back the SDA outcomes, interactivePlot for plotting the outcomes after each step, 
-TimeseriesPlot for post analysis examination, BiasPlot for plotting the correlation between state variables, plot.title is the title of post analysis plots and debug mode allows for pausing the code and examining the variables inside the function.}
 }
 \value{
 NONE

From 5c9a9872372efbe9cdc91157afeeac97be349909 Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Tue, 22 Aug 2023 19:57:02 -0400
Subject: [PATCH 09/32] Add update_phenology function

---
 base/remote/R/qsub_parallel.R                 |  16 +--
 models/sipnet/R/write.configs.SIPNET.R        |  26 +++-
 models/sipnet/R/write_restart.SIPNET.R        |   7 +-
 models/sipnet/man/write.config.SIPNET.Rd      |   5 +-
 models/sipnet/man/write_restart.SIPNET.Rd     |   5 +-
 .../assim.sequential/R/Analysis_sda_block.R   | 123 +++++++++++++-----
 modules/assim.sequential/R/Nimble_codes.R     |   4 +
 .../assim.sequential/R/SDA_OBS_Assembler.R    |   3 +-
 modules/assim.sequential/R/matrix_operation.R |   2 +-
 .../assim.sequential/R/sda.enkf_MultiSite.R   |  58 +++++----
 .../man/analysis_sda_block.Rd                 |   3 +-
 .../man/sda.enkf.multisite.Rd                 |   3 +-
 modules/assim.sequential/man/update_q.Rd      |  10 +-
 modules/uncertainty/R/ensemble.R              |  12 +-
 .../uncertainty/man/write.ensemble.configs.Rd |   5 +-
 15 files changed, 204 insertions(+), 78 deletions(-)

diff --git a/base/remote/R/qsub_parallel.R b/base/remote/R/qsub_parallel.R
index aca657d1ba8..b1f05b16354 100644
--- a/base/remote/R/qsub_parallel.R
+++ b/base/remote/R/qsub_parallel.R
@@ -91,19 +91,15 @@ qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out") {
   pb <- utils::txtProgressBar(min = 0, max = length(unlist(run_list)), style = 3)
   pbi <- 0
   folders <- file.path(settings$host$outdir, run_list)
-  while (length(folders) > 0) {
-    Sys.sleep(10)
-    completed_folders <- foreach::foreach(folder = folders, settings = rep(settings, length(run_list))) %dopar% {
+  completed_folders <- c()
+  while (length(completed_folders) < length(folders)) {
+    completed_folders <- foreach::foreach(folder = folders) %dopar% {
       if(file.exists(file.path(folder, prefix))){
         return(folder)
       }
-    }
-    if(length(unlist(completed_folders)) > 0){
-      Ind <- which(unlist(completed_folders) %in% folders)
-      folders <- folders[-Ind]
-      pbi <- pbi + length(completed_folders)
-      utils::setTxtProgressBar(pb, pbi)
-    }
+    } %>% unlist
+    pbi <- length(completed_folders)
+    utils::setTxtProgressBar(pb, pbi)
   }
   close(pb)
   parallel::stopCluster(cl)
diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 92cefc7cebe..119afd0db8c 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -14,7 +14,7 @@
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
-                                restart = NULL, spinup = NULL) {
+                                restart = NULL, spinup = NULL,obs.mean=NULL,time_t=NULL,update_phenology=NULL) {
   ### WRITE sipnet.in
   template.in <- system.file("sipnet.in", package = "PEcAn.SIPNET")
   config.text <- readLines(con = template.in, n = -1)
@@ -417,6 +417,27 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     if ("leafGrowth" %in% pft.names) {
       param[which(param[, 1] == "leafGrowth"), 2] <- pft.traits[which(pft.names == "leafGrowth")]
     }
+    if (update_phenology){
+    #update LeafOnday
+      PEcAn.logger::logger.warn("time_t")
+      print(time_t)
+      obs.times<-names(obs.mean)
+      if (time_t>1){
+	  tracker <- 1
+          print((obs.mean[[time_t]][[settings$run$site$id]][["LAI"]]-obs.mean[[time_t-1]][[settings$run$site$id]][["LAI"]])/obs.mean[[time_t-1]][[settings$run$site$id]][["LAI"]])
+      if (((obs.mean[[time_t]][[settings$run$site$id]][["LAI"]]-obs.mean[[time_t-1]][[settings$run$site$id]])[["LAI"]]/obs.mean[[time_t-1]][[settings$run$site$id]][["LAI"]])>1){
+          leafOnDay<-lubridate::yday(obs.times[[time_t]])
+          param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
+	  tracker <- 0
+      }
+      if (tracker){
+          param_leafonday<-utils::read.table(file.path(settings$rundir, run.id, "sipnet.param")) #read from previous sipnet.param
+          param[which(param[, 1] == "leafOnDay"), 2] <- param_leafonday[which(param_leafonday[, 1] == "leafOnDay"), 2]
+       }
+      }
+      PEcAn.logger::logger.warn("leafonday after")
+      print(param[which(param[, 1] == "leafOnDay"), 2])
+    }
   }  ## end loop over PFTS
   ####### end parameter update
   #working on reading soil file (only working for 1 soil file)
@@ -545,6 +566,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       if (!is.na(leafOnDay) && is.numeric(leafOnDay)) {
         param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
       }
+
       ## leafOffDay
       leafOffDay <- try(ncdf4::ncvar_get(IC.nc,"date_of_senescence"),silent = TRUE)
       if (!is.na(leafOffDay) && is.numeric(leafOffDay)) {
@@ -609,4 +631,4 @@ remove.config.SIPNET <- function(main.outdir, settings) {
   } else {
     print("*** WARNING: Removal of files on remote host not yet implemented ***")
   }
-} # remove.config.SIPNET 
\ No newline at end of file
+} # remove.config.SIPNET 
diff --git a/models/sipnet/R/write_restart.SIPNET.R b/models/sipnet/R/write_restart.SIPNET.R
index fadbd27353c..c0d1caa9a8b 100755
--- a/models/sipnet/R/write_restart.SIPNET.R
+++ b/models/sipnet/R/write_restart.SIPNET.R
@@ -26,7 +26,7 @@
 ##' @return NONE
 ##' @export
 write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings, new.state,
-                                 RENAME = TRUE, new.params = FALSE, inputs) {
+                                 RENAME = TRUE, new.params = FALSE, inputs,obs.mean=NULL,time=NULL,update_phenology=NULL) {
   
   rundir <- settings$host$rundir
   variables <- colnames(new.state)
@@ -127,6 +127,9 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
                                            settings = settings,
                                            run.id = runid,
                                            inputs = inputs,
-                                           IC = analysis.save.mat))
+                                           IC = analysis.save.mat,
+					   obs.mean = obs.mean,
+                                           time_t=time,
+					   update_phenology=update_phenology))
   print(runid)
 } # write_restart.SIPNET
diff --git a/models/sipnet/man/write.config.SIPNET.Rd b/models/sipnet/man/write.config.SIPNET.Rd
index 889b9237cda..86019f657c7 100644
--- a/models/sipnet/man/write.config.SIPNET.Rd
+++ b/models/sipnet/man/write.config.SIPNET.Rd
@@ -12,7 +12,10 @@ write.config.SIPNET(
   inputs = NULL,
   IC = NULL,
   restart = NULL,
-  spinup = NULL
+  spinup = NULL,
+  obs.mean = NULL,
+  time_t = NULL,
+  update_phenology = NULL
 )
 }
 \description{
diff --git a/models/sipnet/man/write_restart.SIPNET.Rd b/models/sipnet/man/write_restart.SIPNET.Rd
index d080faabd93..895b1857d91 100644
--- a/models/sipnet/man/write_restart.SIPNET.Rd
+++ b/models/sipnet/man/write_restart.SIPNET.Rd
@@ -13,7 +13,10 @@ write_restart.SIPNET(
   new.state,
   RENAME = TRUE,
   new.params = FALSE,
-  inputs
+  inputs,
+  obs.mean = NULL,
+  time = NULL,
+  update_phenology = NULL
 )
 }
 \arguments{
diff --git a/modules/assim.sequential/R/Analysis_sda_block.R b/modules/assim.sequential/R/Analysis_sda_block.R
index e573280108d..7b90dc08974 100644
--- a/modules/assim.sequential/R/Analysis_sda_block.R
+++ b/modules/assim.sequential/R/Analysis_sda_block.R
@@ -16,7 +16,7 @@
 ##' 
 ##' @return It returns the `build.block.xy` object and the analysis results.
 ##' @importFrom dplyr %>%
-analysis_sda_block <- function (settings, block.list.all, X, obs.mean, obs.cov, t, nt, MCMC.args) {
+analysis_sda_block <- function (settings, block.list.all, X, obs.mean, obs.cov, t, nt, MCMC.args, block.list.all.pre = NULL) {
   #convert from vector values to block lists.
   if ("try-error" %in% class(try(block.results <- build.block.xy(settings = settings, 
                                                                  block.list.all = block.list.all, 
@@ -33,7 +33,10 @@ analysis_sda_block <- function (settings, block.list.all, X, obs.mean, obs.cov,
   R <- block.results[[4]]
   
   #update q.
-  if ("try-error" %in% class(try(block.list.all <- update_q(block.list.all, t, nt)))) {
+  if ("try-error" %in% class(try(block.list.all <- update_q(block.list.all, t, nt, aqq.Init = as.numeric(settings$state.data.assimilation$aqq.Init),
+                                                            bqq.Init = as.numeric(settings$state.data.assimilation$bqq.Init),
+                                                            MCMC_dat = NULL,
+                                                            block.list.all.pre)))) {
     PEcAn.logger::logger.error("Something wrong within the update_q function.")
   }
   
@@ -51,8 +54,26 @@ analysis_sda_block <- function (settings, block.list.all, X, obs.mean, obs.cov,
   
   #parallel for loop over each block.
   PEcAn.logger::logger.info(paste0("Running MCMC ", "for ", length(block.list.all[[t]]), " blocks"))
-  if ("try-error" %in% class(try(block.list.all[[t]] <- furrr::future_map(block.list.all[[t]], MCMC_block_function, .progress = T)))) {
-    PEcAn.logger::logger.error("Something wrong within the MCMC_block_function function.")
+ # if ("try-error" %in% class(try(block.list.all[[t]] <- furrr::future_map(block.list.all[[t]], MCMC_block_function, .progress = T)))) {
+ #  PEcAn.logger::logger.error("Something wrong within the MCMC_block_function function.")
+ # }
+
+  if  (length(block.list.all[[t]]) >1){
+  cores <- parallel::detectCores()
+  cl <- parallel::makeCluster(cores)
+  doSNOW::registerDoSNOW(cl)
+  #progress bar
+  pb <- utils::txtProgressBar(min=1, max=length(block.list.all[[t]]), style=3)
+  progress <- function(n) utils::setTxtProgressBar(pb, n)
+  opts <- list(progress=progress)
+  PEcAn.logger::logger.info("Submitting jobs!")
+   block.list.all[[t]]<-foreach::foreach(block = block.list.all[[t]], .packages="Kendall", .options.snow=opts) %dopar% {
+    library(nimble)
+    library(purrr)
+    MCMC_block_function(block)
+}
+}else{
+    block.list.all[[t]][[1]]<-MCMC_block_function(block=block.list.all[[t]][[1]])
   }
   PEcAn.logger::logger.info("Completed!")
   
@@ -177,14 +198,25 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
       
       #fill in constants.
       block.h <- Construct.H.multisite(site.ids[i], var.names, obs.mean[[t]])
+      block.list[[i]]$H <- block.h
       block.list[[i]]$constant$H <- which(apply(block.h, 2, sum) == 1)
       block.list[[i]]$constant$N <- length(f.start:f.end)
       block.list[[i]]$constant$YN <- length(y.start:y.end)
       block.list[[i]]$constant$q.type <- q.type
     }
+    names(block.list) <- site.ids
   } else {
     #find networks given TRUE/FALSE matrix representing sites' interactions.
     block.vec <- matrix_network(dis.matrix <= as.numeric(settings$state.data.assimilation$scalef))
+    #check if the matrix_network function is working correctly.
+    #check if the blocks are calculated correctly.
+    if (block.vec %>% 
+        purrr::map(function(l){length(l)}) %>%
+        unlist %>%
+        sum() != length(site.ids)) {
+      PEcAn.logger::logger.error("Block calculation failed, please check the matrix_network function!")
+      return(0)
+    }
     block.list <- vector("list", length(block.vec))
     #loop over sites
     for (i in seq_along(block.vec)) {#i is site index
@@ -213,6 +245,7 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
       
       #fill in constants
       block.h <- Construct.H.multisite(site.ids[ids], var.names, obs.mean[[t]])
+      block.list[[i]]$H <- block.h
       block.list[[i]]$constant$H <- which(apply(block.h, 2, sum) == 1)
       block.list[[i]]$constant$N <- length(f.ind)
       block.list[[i]]$constant$YN <- length(y.ind)
@@ -237,7 +270,7 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
 MCMC_Init <- function (block.list, X) {
   var.names <- unique(attributes(X)$dimnames[[2]])
   #sample mu.f from X.
-  sample.mu.f <- X[sample(seq_along(1:nrow(X)), 1),]
+  sample.mu.f <- colMeans(X)
   for (i in seq_along(block.list)) {
     #number of observations.
     num.obs <- length(block.list[[i]]$data$y.censored)
@@ -248,7 +281,7 @@ MCMC_Init <- function (block.list, X) {
       end <- (block.list[[i]]$sites.per.block[j]) * length(var.names)
       block.list[[i]]$Inits$X.mod <- c(block.list[[i]]$Inits$X.mod, sample.mu.f[start:end])
       #initialize X
-      block.list[[i]]$Inits$X <- block.list[[i]]$Inits$X.mod[block.list[[i]]$constant$H]
+      block.list[[i]]$Inits$X <- block.list[[i]]$data$y.censored
       #initialize Xs
       block.list[[i]]$Inits$Xs <- block.list[[i]]$Inits$X.mod[block.list[[i]]$constant$H]
     }
@@ -293,6 +326,7 @@ MCMC_block_function <- function(block) {
   #hear we change the RW_block sampler to the ess sampler 
   #because it has a better performance of MVN sampling
   samplerLists <- conf$getSamplers()
+  samplerNumberOffset <- length(samplerLists)
   if (block$constant$q.type == 1) {
     #if we have vector q
     #only X.mod should be sampled with ess sampler.
@@ -305,11 +339,21 @@ MCMC_block_function <- function(block) {
   }
   conf$setSamplers(samplerLists)
   
+  #add toggle Y sampler.
+  for (i in 1:block$constant$YN) {
+    conf$addSampler(paste0("y.censored[", i, "]"), 'toggle', control=list(type='RW'))
+  }
+  conf$printSamplers()
   #compile MCMC
   Rmcmc <- nimble::buildMCMC(conf)
   Cmodel <- nimble::compileNimble(model_pred)
   Cmcmc <- nimble::compileNimble(Rmcmc, project = model_pred, showCompilerOutput = FALSE)
   
+  #add toggle Y sampler.
+  for(i in 1:block$constant$YN) {
+    valueInCompiledNimbleFunction(Cmcmc$samplerFunctions[[samplerNumberOffset+i]], 'toggle', 0)
+  }
+  
   #run MCMC
   dat <- runMCMC(Cmcmc, niter = block$MCMC$niter, nburnin = block$MCMC$nburnin, thin = block$MCMC$nthin, nchains = block$MCMC$nchain)
   
@@ -355,27 +399,35 @@ MCMC_block_function <- function(block) {
     aq <- solve(q.bar) * bq
     block$aqq[,,block$t+1] <- GrabFillMatrix(block$aqq[,,block$t], block$constant$H, aq)
     block$bqq[block$t+1] <- bq
-    
-    # #if it's a wishart case
-    # bq <- block$data$bq
-    # aq <- block$data$aq
-    # for (i in 1:dim(aq)[1]) {
-    #   CHAR <- paste0("[", i, "]")
-    #   for (j in 1:dim(aq)[2]) {
-    #     aq[i, j] <- M[paste0("q[", i, ", ", j, "]")]/bq
-    #   }
-    # }
-    # #update aqq and bqq
-    # block$aqq[,,block$t+1] <- GrabFillMatrix(block$aqq[,,block$t], block$constant$H, aq)
-    # block$bqq[block$t+1] <- block$bqq[block$t]
   }
   #update mua and pa; mufa, and pfa
   iX <- grep("X[", colnames(dat), fixed = TRUE)
   iX.mod <- grep("X.mod[", colnames(dat), fixed = TRUE)
-  mua <- colMeans(dat[, iX])
-  pa <- stats::cov(dat[, iX])
-  mufa <- colMeans(dat[, iX.mod])
-  pfa <- stats::cov(dat[, iX.mod])
+  if (length(iX) == 1) {
+    mua <- mean(dat[, iX])
+    pa <- var(dat[, iX])
+  } else {
+    mua <- colMeans(dat[, iX])
+    pa <- stats::cov(dat[, iX])
+  }
+  
+  if (length(iX.mod) == 1) {
+    mufa <- mean(dat[, iX.mod])
+    pfa <- var(dat[, iX.mod])
+  } else {
+    mufa <- colMeans(dat[, iX.mod])
+    pfa <- stats::cov(dat[, iX.mod])
+  }
+  
+  #
+  # if (length(iX) == 1) {
+  #   K <- pfa %*% t(block$H) %*% solve(1/block$data$r + block$H %*% pfa %*% t(block$H))
+  # } else {
+  #   K <- pfa %*% t(block$H) %*% solve((diag(1/diag(block$data$r)) + block$H %*% pfa %*% t(block$H)))
+  # }
+  # 
+  # pfa <- (diag(block$constant$N) - K %*% block$H) %*% pfa
+  # pa <- GrabFillMatrix(kf.pa, block$constant$H)
   
   #return values.
   block$update <- list(aq = aq, bq = bq, mua = mua, pa = pa, mufa = mufa, pfa = pfa)
@@ -392,7 +444,7 @@ MCMC_block_function <- function(block) {
 ##' @param MCMC_dat data frame of MCMC samples, the default it NULL.
 ##' 
 ##' @return It returns the `block.list.all` object with initialized/updated Q filled in.
-update_q <- function (block.list.all, t, nt, MCMC_dat = NULL) {
+update_q <- function (block.list.all, t, nt, aqq.Init = NULL, bqq.Init = NULL, MCMC_dat = NULL, block.list.all.pre = NULL) {
   block.list <- block.list.all[[t]]
   #if it's an update.
   if (is.null(MCMC_dat)) {
@@ -403,24 +455,33 @@ update_q <- function (block.list.all, t, nt, MCMC_dat = NULL) {
         nobs <- length(block.list[[i]]$data$y.censored)
         if (block.list[[i]]$constant$q.type == 1) {
           #initialize aqq and bqq for nt
-          block.list[[i]]$aqq <- array(1, dim = c(nvar, nt))
-          block.list[[i]]$bqq <- array(1, dim = c(nvar, nt))
+          if (!is.null(aqq.Init) && !is.null(bqq.Init)) {
+            block.list[[i]]$aqq <- array(aqq.Init, dim = c(nvar, nt + 1))
+            block.list[[i]]$bqq <- array(bqq.Init, dim = c(nvar, nt + 1))
+          } else {
+            block.list[[i]]$aqq <- array(1, dim = c(nvar, nt + 1))
+            block.list[[i]]$bqq <- array(1, dim = c(nvar, nt + 1))
+          }
           #update aq and bq based on aqq and bqq
           block.list[[i]]$data$aq <- block.list[[i]]$aqq[block.list[[i]]$constant$H, t]
           block.list[[i]]$data$bq <- block.list[[i]]$bqq[block.list[[i]]$constant$H, t]
         } else if (block.list[[i]]$constant$q.type == 2) {
           #initialize aqq and bqq for nt
-          block.list[[i]]$aqq <- array(1, dim = c(nvar, nvar, nt))
+          block.list[[i]]$aqq <- array(1, dim = c(nvar, nvar, nt + 1))
           block.list[[i]]$aqq[,,t] <- stats::toeplitz((nvar:1)/nvar)
-          block.list[[i]]$bqq <- rep(nobs, nt)
+          block.list[[i]]$bqq <- rep(nobs, nt + 1)
           #update aq and bq based on aqq and bqq
           block.list[[i]]$data$aq <- GrabFillMatrix(block.list[[i]]$aqq[,,t], block.list[[i]]$constant$H)
           block.list[[i]]$data$bq <- block.list[[i]]$bqq[t]
         }
       }
     } else if (t > 1) {
-      #if we want to update q from previous SDA runs.
-      block.list.pre <- block.list.all[[t - 1]]
+      if (!is.null(block.list.all.pre)) {
+        block.list.pre <- block.list.all.pre[[t - 1]]
+      } else {
+        #if we want to update q from previous SDA runs.
+        block.list.pre <- block.list.all[[t - 1]]
+      }
       for (i in seq_along(block.list)) {
         nvar <- length(block.list[[i]]$data$muf)
         nobs <- length(block.list[[i]]$data$y.censored)
@@ -484,4 +545,4 @@ block.2.vector <- function (block.list, X, H) {
               Pf = Pf,
               mu.a = mu.a,
               Pa = Pa))
-}
\ No newline at end of file
+}
diff --git a/modules/assim.sequential/R/Nimble_codes.R b/modules/assim.sequential/R/Nimble_codes.R
index 1a97239b340..236f0351428 100644
--- a/modules/assim.sequential/R/Nimble_codes.R
+++ b/modules/assim.sequential/R/Nimble_codes.R
@@ -238,8 +238,12 @@ GEF.Block.Nimble <-  nimbleCode({
     for (i in 1:YN) {
       #sample Q.
       q[i] ~ dgamma(shape = aq[i], rate = bq[i])
+      if (length(H) == 1){
+        X[i] ~ dnorm(X.mod[H], sd = 1/sqrt(q[i]))
+      } else{
       #sample latent variable X.
       X[i]  ~ dnorm(X.mod[H[i]], sd = 1/sqrt(q[i]))
+      }
       #likelihood
       y.censored[i] ~ dnorm(X[i], sd = 1/sqrt(r[i, i]))
     }
diff --git a/modules/assim.sequential/R/SDA_OBS_Assembler.R b/modules/assim.sequential/R/SDA_OBS_Assembler.R
index 9510a2ad89a..233f62353a2 100644
--- a/modules/assim.sequential/R/SDA_OBS_Assembler.R
+++ b/modules/assim.sequential/R/SDA_OBS_Assembler.R
@@ -173,12 +173,13 @@ SDA_OBS_Assembler <- function(settings){
     for (j in seq_along(obs.mean[[i]])) {
       if (sum(is.na(obs.mean[[i]][[j]]))){
         na_ind <- which(is.na(obs.mean[[i]][[j]]))
-        obs.mean[[i]][[j]] <- obs.mean[[i]][[j]][-na_ind]
+        #obs.mean[[i]][[j]] <- obs.mean[[i]][[j]][-na_ind]
         if(length(obs.mean[[i]][[j]]) == 1){
           obs.cov[[i]][[j]] <- obs.cov[[i]][[j]][-na_ind]
         }else{
           obs.cov[[i]][[j]] <- obs.cov[[i]][[j]][-na_ind, -na_ind]
         }
+	obs.mean[[i]][[j]] <- obs.mean[[i]][[j]][-na_ind]
       }
       SoilC_ind <- which(names(obs.mean[[i]][[j]]) == "TotSoilCarb")
       if (length(SoilC_ind) > 0){
diff --git a/modules/assim.sequential/R/matrix_operation.R b/modules/assim.sequential/R/matrix_operation.R
index 5567e153aac..c5177c21ece 100644
--- a/modules/assim.sequential/R/matrix_operation.R
+++ b/modules/assim.sequential/R/matrix_operation.R
@@ -59,7 +59,7 @@ matrix_network <- function (mat) {
         Inits <- c(Inits, which(mat[init,]))
       }
       Inits <- Inits[which(!Inits %in% vec)]
-      vec <- c(vec, Inits)
+      vec <- sort(unique(c(vec, Inits)))
       #if we don't have any new site that belongs to this network.
       if (length(Inits) == 0) {
         #then stop.
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index 20b3339bdf7..3fa4e58d4dd 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -29,8 +29,6 @@ sda.enkf.multisite <- function(settings,
                                pre_enkf_params = NULL,
                                ensemble.samples = NULL,
                                control=list(trace = TRUE,
-					    free_run = FALSE,
-					    send_email = NULL,
                                             FF = FALSE,
                                             interactivePlot = FALSE,
                                             TimeseriesPlot = FALSE,
@@ -46,7 +44,8 @@ sda.enkf.multisite <- function(settings,
                                             send_email = NULL,
                                             keepNC = TRUE,
                                             forceRun = TRUE,
-                                            run_parallel = TRUE),
+                                            run_parallel = TRUE,
+					    update_phenology = FALSE),
                                ...) {
   #add if/else for when restart points to folder instead if T/F set restart as T
   if(is.list(restart)){
@@ -416,7 +415,10 @@ sda.enkf.multisite <- function(settings,
               model = settings$model$type,
               write.to.db = settings$database$bety$write,
               restart = restart.arg,
-              rename = TRUE
+              rename = TRUE,
+	      obs.mean = obs.mean,
+	      time=t,
+              update_phenology=control$update_phenology
             )
           }) %>%
           stats::setNames(site.ids)
@@ -428,16 +430,15 @@ sda.enkf.multisite <- function(settings,
         paste(file.path(rundir, 'runs.txt'))  ## testing
         Sys.sleep(0.01)                       ## testing
         if(control$parallel_qsub){
-          PEcAn.remote::qsub_parallel(settings, files=PEcAn.remote::merge_job_files(settings, 10), prefix = paste0(obs.year, ".nc"))
+          PEcAn.remote::qsub_parallel(settings, files=PEcAn.remote::merge_job_files(settings, control$jobs.per.file), prefix = paste0(obs.year, ".nc"))
         }else{
           PEcAn.workflow::start_model_runs(settings, write=settings$database$bety$write)
         }
         #------------- Reading - every iteration and for SDA
         #put building of X into a function that gets called
-        #max_t <- 0
-        #while("try-error" %in% class(
-         # try(
-	reads <-  build_X(out.configs = out.configs, 
+        max_t <- 0
+        while("try-error" %in% class(
+          try(reads <- PEcAnAssimSequential:::build_X(out.configs = out.configs, 
                                settings = settings, 
                                new.params = new.params, 
                                nens = nens, 
@@ -446,7 +447,16 @@ sda.enkf.multisite <- function(settings,
                                t = t, 
                                var.names = var.names, 
                                my.read_restart = my.read_restart,
-                               restart_flag = restart_flag)
+                               restart_flag = restart_flag), silent = T))
+        ){
+          Sys.sleep(10)
+          max_t <- max_t + 1
+          if(max_t > 1000){
+            PEcAn.logger::logger.info("Can't find outputed NC file! Please rerun the code!")
+            break
+          }
+          PEcAn.logger::logger.info("Empty folder, try again!")
+        }
         
         if (control$debug) browser()
         #let's read the parameters of each site/ens
@@ -497,7 +507,7 @@ sda.enkf.multisite <- function(settings,
         #decide if we want the block analysis function or multi-site analysis function.
         if (processvar == TRUE && settings$state.data.assimilation$q.type %in% c("vector", "wishart")) {
           #initialize block.list.all.
-          if (t == 1) {
+          if (t == 1 | !exists("block.list.all")) {
             block.list.all <- obs.mean %>% purrr::map(function(l){NULL})
           }
           #initialize MCMC arguments.
@@ -505,10 +515,13 @@ sda.enkf.multisite <- function(settings,
             MCMC.args <- list(niter = 1e5,
                               nthin = 10,
                               nchain = 1,
-                              nburnin = 1e4)
+                              nburnin = 5e4)
           }
           #running analysis function.
-          enkf.params[[obs.t]] <- analysis_sda_block(settings, block.list.all, X, obs.mean, obs.cov, t, nt, MCMC.args)
+	  save(X,file = file.path(settings$outdir,"X.Rdata"))
+          PEcAn.logger::logger.info("The total number of time is:")
+	  print(nt)
+          enkf.params[[obs.t]] <- analysis_sda_block(settings, block.list.all, X, obs.mean, obs.cov, t, nt, MCMC.args, pre_enkf_params)
           enkf.params[[obs.t]] <- c(enkf.params[[obs.t]], RestartList = list(restart.list %>% stats::setNames(site.ids)))
           block.list.all <- enkf.params[[obs.t]]$block.list.all
           #Forecast
@@ -546,7 +559,7 @@ sda.enkf.multisite <- function(settings,
           if(!is.null(pre_enkf_params)){
             Pf <- pre_enkf_params[[t]]$Pf
           }
-          recompileTobit = !exists('Cmcmc_tobit2space')
+          recompileTobit = !exists('Cmcmc_tobit2space')     
           recompileGEF   = !exists('Cmcmc')
           #weight list
           # This reads ensemble weights generated by `get_ensemble_weights` function from assim.sequential package
@@ -657,11 +670,12 @@ sda.enkf.multisite <- function(settings,
           #will throw an error when q.bar and Pf are different sizes i.e. when you are running with no obs and do not variance for all state variables
           #Pa <- Pf + solve(q.bar)
           #hack have Pa = Pf for now
-          if(!is.null(pre_enkf_params)){
-            Pf <- pre_enkf_params[[t]]$Pf
-          }else{
-            Pf <- stats::cov(X) # Cov Forecast - This is used as an initial condition
-          }
+          # if(!is.null(pre_enkf_params)){
+          #   Pf <- pre_enkf_params[[t]]$Pf
+          # }else{
+          #   Pf <- stats::cov(X) # Cov Forecast - This is used as an initial condition
+          # }
+          Pf <- stats::cov(X)
           Pa <- Pf
         }
         enkf.params[[obs.t]] <- list(mu.f = mu.f, Pf = Pf, mu.a = mu.a, Pa = Pa)
@@ -724,9 +738,9 @@ sda.enkf.multisite <- function(settings,
       unlink(mailfile)
     }
     # useful for debugging to keep .nc files for assimilated years. T = 2, because this loops removes the files that were run when starting the next loop
-   if (keepNC && t == 1){
-    unlink(list.files(outdir, "*.nc", recursive = TRUE, full.names = TRUE))
-   }
+#    if (keepNC && t == 1){
+#      unlink(list.files(outdir, "*.nc", recursive = TRUE, full.names = TRUE))
+#    }
       ## MCD: I commented the above "if" out because if you are restarting from a previous forecast, this might delete the files in that earlier folder
   } ### end loop over time
 } # sda.enkf
diff --git a/modules/assim.sequential/man/analysis_sda_block.Rd b/modules/assim.sequential/man/analysis_sda_block.Rd
index b5f836aa932..cc8346fb0c4 100644
--- a/modules/assim.sequential/man/analysis_sda_block.Rd
+++ b/modules/assim.sequential/man/analysis_sda_block.Rd
@@ -12,7 +12,8 @@ analysis_sda_block(
   obs.cov,
   t,
   nt,
-  MCMC.args
+  MCMC.args,
+  block.list.all.pre = NULL
 )
 }
 \arguments{
diff --git a/modules/assim.sequential/man/sda.enkf.multisite.Rd b/modules/assim.sequential/man/sda.enkf.multisite.Rd
index 5017d63df8f..b3523f109fa 100644
--- a/modules/assim.sequential/man/sda.enkf.multisite.Rd
+++ b/modules/assim.sequential/man/sda.enkf.multisite.Rd
@@ -15,7 +15,8 @@ sda.enkf.multisite(
   control = list(trace = TRUE, FF = FALSE, interactivePlot = FALSE, TimeseriesPlot =
     FALSE, BiasPlot = FALSE, plot.title = NULL, facet.plots = FALSE, debug = FALSE, pause
     = FALSE, Profiling = FALSE, OutlierDetection = FALSE, parallel_qsub = TRUE, free_run
-    = FALSE, send_email = NULL, keepNC = TRUE, forceRun = TRUE, run_parallel = TRUE),
+    = FALSE, send_email = NULL, keepNC = TRUE, forceRun = TRUE, run_parallel = TRUE,
+    update_phenology = FALSE),
   ...
 )
 }
diff --git a/modules/assim.sequential/man/update_q.Rd b/modules/assim.sequential/man/update_q.Rd
index bdf6acef08f..9d44af37253 100644
--- a/modules/assim.sequential/man/update_q.Rd
+++ b/modules/assim.sequential/man/update_q.Rd
@@ -4,7 +4,15 @@
 \alias{update_q}
 \title{update_q}
 \usage{
-update_q(block.list.all, t, nt, MCMC_dat = NULL)
+update_q(
+  block.list.all,
+  t,
+  nt,
+  aqq.Init = NULL,
+  bqq.Init = NULL,
+  MCMC_dat = NULL,
+  block.list.all.pre = NULL
+)
 }
 \arguments{
 \item{block.list.all}{each block within the `block.list` lists.}
diff --git a/modules/uncertainty/R/ensemble.R b/modules/uncertainty/R/ensemble.R
index 6bde5171414..62c058aa057 100644
--- a/modules/uncertainty/R/ensemble.R
+++ b/modules/uncertainty/R/ensemble.R
@@ -210,7 +210,7 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 ##' @export
 ##' @author David LeBauer, Carl Davidson, Hamze Dokoohaki
 write.ensemble.configs <- function(defaults, ensemble.samples, settings, model, 
-                                   clean = FALSE, write.to.db = TRUE, restart=NULL, rename = FALSE) {
+                                   clean = FALSE, write.to.db = TRUE, restart=NULL, rename = FALSE,obs.mean=NULL,time=NULL,update_phenology=NULL) {
   
   con <- NULL
   my.write.config <- paste("write.config.", model, sep = "")
@@ -408,7 +408,10 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
       do.call(my.write.config, args = list( defaults = defaults, 
                                             trait.values = lapply(samples$parameters$samples, function(x, n) { x[n, , drop=FALSE] }, n=i), # this is the params
                                             settings = settings, 
-                                            run.id = run.id
+                                            run.id = run.id,
+					    obs.mean = obs.mean,
+					    time_t=time,
+					    update_phenology=update_phenology
       )
       )
       cat(format(run.id, scientific = FALSE), file = file.path(settings$rundir, "runs.txt"), sep = "\n", append = TRUE)
@@ -458,7 +461,10 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
                            new.state = new.state[i, ], 
                            new.params = new.params[[i]], #new.params$`646`[[i]] for debugging
                            inputs =list(met=list(path=inputs$samples[[i]])), 
-                           RENAME = rename)#for restart from previous model runs, not sharing the same outdir
+                           RENAME = rename,
+			   obs.mean=obs.mean,
+			   time=time,
+			   update_phenology=update_phenology)#for restart from previous model runs, not sharing the same outdir
       )
     }
     params<-new.params
diff --git a/modules/uncertainty/man/write.ensemble.configs.Rd b/modules/uncertainty/man/write.ensemble.configs.Rd
index 45d67cdcb7e..4872a8202e8 100644
--- a/modules/uncertainty/man/write.ensemble.configs.Rd
+++ b/modules/uncertainty/man/write.ensemble.configs.Rd
@@ -12,7 +12,10 @@ write.ensemble.configs(
   clean = FALSE,
   write.to.db = TRUE,
   restart = NULL,
-  rename = FALSE
+  rename = FALSE,
+  obs.mean = NULL,
+  time = NULL,
+  update_phenology = NULL
 )
 }
 \arguments{

From b80deff378167b17f09d9e05cec42d5aedd08828 Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Fri, 1 Sep 2023 00:11:43 -0400
Subject: [PATCH 10/32] Adjust update_phenology function based on MODIS
 phenology data

---
 models/sipnet/R/write.configs.SIPNET.R        | 28 ++++++++-----------
 models/sipnet/R/write_restart.SIPNET.R        |  7 ++---
 models/sipnet/man/write.config.SIPNET.Rd      |  3 +-
 models/sipnet/man/write_restart.SIPNET.Rd     |  3 +-
 .../assim.sequential/R/sda.enkf_MultiSite.R   |  3 +-
 .../man/sda.enkf.multisite.Rd                 |  3 +-
 modules/uncertainty/R/ensemble.R              | 12 ++++----
 modules/uncertainty/R/get.parameter.samples.R |  7 ++++-
 .../uncertainty/man/write.ensemble.configs.Rd |  1 -
 9 files changed, 30 insertions(+), 37 deletions(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 77935ac27e6..6847f2cba12 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -14,7 +14,7 @@
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
-                                restart = NULL, spinup = NULL,obs.mean=NULL,time_t=NULL,update_phenology=NULL) {
+                                restart = NULL, spinup = NULL,obs_time=NULL,update_phenology=NULL) {
   ### WRITE sipnet.in
   template.in <- system.file("sipnet.in", package = "PEcAn.SIPNET")
   config.text <- readLines(con = template.in, n = -1)
@@ -417,24 +417,18 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     if ("leafGrowth" %in% pft.names) {
       param[which(param[, 1] == "leafGrowth"), 2] <- pft.traits[which(pft.names == "leafGrowth")]
     }
+
+    #update LeafOnday and LeafOffDay
     if (update_phenology){
-    #update LeafOnday
-      PEcAn.logger::logger.warn("time_t")
-      print(time_t)
-      obs.times<-names(obs.mean)
-      if (time_t>1){
-	  tracker <- 1
-          print((obs.mean[[time_t]][[settings$run$site$id]][["LAI"]]-obs.mean[[time_t-1]][[settings$run$site$id]][["LAI"]])/obs.mean[[time_t-1]][[settings$run$site$id]][["LAI"]])
-      if (((obs.mean[[time_t]][[settings$run$site$id]][["LAI"]]-obs.mean[[time_t-1]][[settings$run$site$id]])[["LAI"]]/obs.mean[[time_t-1]][[settings$run$site$id]][["LAI"]])>1){
-          leafOnDay<-lubridate::yday(obs.times[[time_t]])
-          param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
-	  tracker <- 0
-      }
-      if (tracker){
-          param_leafonday<-utils::read.table(file.path(settings$rundir, run.id, "sipnet.param")) #read from previous sipnet.param
-          param[which(param[, 1] == "leafOnDay"), 2] <- param_leafonday[which(param_leafonday[, 1] == "leafOnDay"), 2]
+       leafphdata <- utils::read.csv("/data2/Projects/NASA_CMS/neon_soc_sda/leaf_phenology.csv")
+       leafon <- leafphdata$leafonday[leafphdata$year==obs_time & leafphdata$site_id==settings$run$site$id]
+       leafoff<- leafphdata$leafoffday[leafphdata$year==obs_time & leafphdata$site_id==settings$run$site$id]
+       if (!is.na(leafon)){
+	 param[which(param[, 1] == "leafOnDay"), 2] <- as.numeric(leafon)
+       }
+       if (!is.na(leafoff)){
+        param[which(param[, 1] == "leafOffDay"), 2] <- as.numeric(leafoff)
        }
-      }
       PEcAn.logger::logger.warn("leafonday after")
       print(param[which(param[, 1] == "leafOnDay"), 2])
     }
diff --git a/models/sipnet/R/write_restart.SIPNET.R b/models/sipnet/R/write_restart.SIPNET.R
index c0d1caa9a8b..bbf8b3ef33d 100755
--- a/models/sipnet/R/write_restart.SIPNET.R
+++ b/models/sipnet/R/write_restart.SIPNET.R
@@ -26,7 +26,7 @@
 ##' @return NONE
 ##' @export
 write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings, new.state,
-                                 RENAME = TRUE, new.params = FALSE, inputs,obs.mean=NULL,time=NULL,update_phenology=NULL) {
+                                 RENAME = TRUE, new.params = FALSE, inputs, obs_time=NULL, update_phenology=NULL) {
   
   rundir <- settings$host$rundir
   variables <- colnames(new.state)
@@ -128,8 +128,7 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
                                            run.id = runid,
                                            inputs = inputs,
                                            IC = analysis.save.mat,
-					   obs.mean = obs.mean,
-                                           time_t=time,
-					   update_phenology=update_phenology))
+                                           obs_time=obs_time,
+					                                 update_phenology=update_phenology))
   print(runid)
 } # write_restart.SIPNET
diff --git a/models/sipnet/man/write.config.SIPNET.Rd b/models/sipnet/man/write.config.SIPNET.Rd
index 86019f657c7..5748aa60866 100644
--- a/models/sipnet/man/write.config.SIPNET.Rd
+++ b/models/sipnet/man/write.config.SIPNET.Rd
@@ -13,8 +13,7 @@ write.config.SIPNET(
   IC = NULL,
   restart = NULL,
   spinup = NULL,
-  obs.mean = NULL,
-  time_t = NULL,
+  obs_time = NULL,
   update_phenology = NULL
 )
 }
diff --git a/models/sipnet/man/write_restart.SIPNET.Rd b/models/sipnet/man/write_restart.SIPNET.Rd
index 895b1857d91..62d5019d279 100644
--- a/models/sipnet/man/write_restart.SIPNET.Rd
+++ b/models/sipnet/man/write_restart.SIPNET.Rd
@@ -14,8 +14,7 @@ write_restart.SIPNET(
   RENAME = TRUE,
   new.params = FALSE,
   inputs,
-  obs.mean = NULL,
-  time = NULL,
+  obs_time = NULL,
   update_phenology = NULL
 )
 }
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index 76ad05e323e..c83c58dbd7c 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -415,8 +415,7 @@ sda.enkf.multisite <- function(settings,
               write.to.db = settings$database$bety$write,
               restart = restart.arg,
               rename = TRUE,
-	      obs.mean = obs.mean,
-	      time=t,
+              time = obs.year,
               update_phenology=control$update_phenology
             )
           }) %>%
diff --git a/modules/assim.sequential/man/sda.enkf.multisite.Rd b/modules/assim.sequential/man/sda.enkf.multisite.Rd
index 62a06a1821d..f3a87ce6a2b 100644
--- a/modules/assim.sequential/man/sda.enkf.multisite.Rd
+++ b/modules/assim.sequential/man/sda.enkf.multisite.Rd
@@ -15,7 +15,8 @@ sda.enkf.multisite(
   control = list(trace = TRUE, FF = FALSE, interactivePlot = FALSE, TimeseriesPlot =
     FALSE, BiasPlot = FALSE, plot.title = NULL, facet.plots = FALSE, debug = FALSE, pause
     = FALSE, Profiling = FALSE, OutlierDetection = FALSE, parallel_qsub = TRUE,
-    send_email = NULL, keepNC = TRUE, forceRun = TRUE, run_parallel = TRUE),
+    send_email = NULL, keepNC = TRUE, forceRun = TRUE, run_parallel = TRUE,
+    update_phenology = FALSE),
   ...
 )
 }
diff --git a/modules/uncertainty/R/ensemble.R b/modules/uncertainty/R/ensemble.R
index 62c058aa057..24fb76a9e38 100644
--- a/modules/uncertainty/R/ensemble.R
+++ b/modules/uncertainty/R/ensemble.R
@@ -210,7 +210,7 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 ##' @export
 ##' @author David LeBauer, Carl Davidson, Hamze Dokoohaki
 write.ensemble.configs <- function(defaults, ensemble.samples, settings, model, 
-                                   clean = FALSE, write.to.db = TRUE, restart=NULL, rename = FALSE,obs.mean=NULL,time=NULL,update_phenology=NULL) {
+                                   clean = FALSE, write.to.db = TRUE, restart=NULL, rename = FALSE,time=NULL,update_phenology=NULL) {
   
   con <- NULL
   my.write.config <- paste("write.config.", model, sep = "")
@@ -409,9 +409,8 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
                                             trait.values = lapply(samples$parameters$samples, function(x, n) { x[n, , drop=FALSE] }, n=i), # this is the params
                                             settings = settings, 
                                             run.id = run.id,
-					    obs.mean = obs.mean,
-					    time_t=time,
-					    update_phenology=update_phenology
+					                                  obs_time = time,
+					                                  update_phenology=update_phenology
       )
       )
       cat(format(run.id, scientific = FALSE), file = file.path(settings$rundir, "runs.txt"), sep = "\n", append = TRUE)
@@ -462,9 +461,8 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
                            new.params = new.params[[i]], #new.params$`646`[[i]] for debugging
                            inputs =list(met=list(path=inputs$samples[[i]])), 
                            RENAME = rename,
-			   obs.mean=obs.mean,
-			   time=time,
-			   update_phenology=update_phenology)#for restart from previous model runs, not sharing the same outdir
+			                     obs_time = time,
+			                     update_phenology = update_phenology)#for restart from previous model runs, not sharing the same outdir
       )
     }
     params<-new.params
diff --git a/modules/uncertainty/R/get.parameter.samples.R b/modules/uncertainty/R/get.parameter.samples.R
index cc323b69fb2..7853530f292 100644
--- a/modules/uncertainty/R/get.parameter.samples.R
+++ b/modules/uncertainty/R/get.parameter.samples.R
@@ -163,8 +163,13 @@ get.parameter.samples <- function(settings,
     }
     for (prior in priors) {
       if (prior %in% param.names[[i]]) {
+        if (prior=="som_respiration_rate") {
+         samples <- trait.mcmc[[prior]] %>% purrr::map(~ .x[,'beta.o']) %>% unlist() %>% as.matrix()
+         samples <- samples[(samples>0)&(samples<0.007)]
+    } else {
         samples <- trait.mcmc[[prior]] %>% purrr::map(~ .x[,'beta.o']) %>% unlist() %>% as.matrix()
-      } else {
+      } }
+      else {
         samples <- PEcAn.priors::get.sample(prior.distns[prior, ], samples.num, q_samples[ , priors==prior])
       }
       trait.samples[[pft.name]][[prior]] <- samples
diff --git a/modules/uncertainty/man/write.ensemble.configs.Rd b/modules/uncertainty/man/write.ensemble.configs.Rd
index 4872a8202e8..d888588019a 100644
--- a/modules/uncertainty/man/write.ensemble.configs.Rd
+++ b/modules/uncertainty/man/write.ensemble.configs.Rd
@@ -13,7 +13,6 @@ write.ensemble.configs(
   write.to.db = TRUE,
   restart = NULL,
   rename = FALSE,
-  obs.mean = NULL,
   time = NULL,
   update_phenology = NULL
 )

From 8ff4f4da0d8145f7b1c00fc837b6b33ef2c2492d Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Fri, 1 Sep 2023 19:41:49 -0400
Subject: [PATCH 11/32] Fix some bugs

---
 models/sipnet/R/write.configs.SIPNET.R          | 4 ++--
 modules/assim.sequential/R/Analysis_sda_block.R | 6 ++++--
 modules/assim.sequential/R/sda.enkf_MultiSite.R | 2 +-
 modules/uncertainty/R/get.parameter.samples.R   | 2 +-
 4 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 6847f2cba12..57470d4d4bd 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -424,10 +424,10 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
        leafon <- leafphdata$leafonday[leafphdata$year==obs_time & leafphdata$site_id==settings$run$site$id]
        leafoff<- leafphdata$leafoffday[leafphdata$year==obs_time & leafphdata$site_id==settings$run$site$id]
        if (!is.na(leafon)){
-	 param[which(param[, 1] == "leafOnDay"), 2] <- as.numeric(leafon)
+	      param[which(param[, 1] == "leafOnDay"), 2] <- leafon
        }
        if (!is.na(leafoff)){
-        param[which(param[, 1] == "leafOffDay"), 2] <- as.numeric(leafoff)
+        param[which(param[, 1] == "leafOffDay"), 2] <- leafoff
        }
       PEcAn.logger::logger.warn("leafonday after")
       print(param[which(param[, 1] == "leafOnDay"), 2])
diff --git a/modules/assim.sequential/R/Analysis_sda_block.R b/modules/assim.sequential/R/Analysis_sda_block.R
index 2608f3df376..e603f1a24eb 100644
--- a/modules/assim.sequential/R/Analysis_sda_block.R
+++ b/modules/assim.sequential/R/Analysis_sda_block.R
@@ -213,8 +213,10 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
         block.list[[i]]$data$r <- diag(1, length(var.names))
         block.h <- matrix(1, 1, length(var.names))
       } else {
-        y.start <- obs_per_site[i] * (i - 1) + 1
-        y.end <- obs_per_site[i] * i
+        #y.start <- obs_per_site[i] * (i - 1) + 1
+        #y.end <- obs_per_site[i] * i
+        y.start <- sum(obs_per_site[1:i-1]) + 1  #obs_per_site[i] * (i - 1) + 1
+        y.end <- y.start + obs_per_site[i] - 1  #obs_per_site[i] * i
         block.list[[i]]$data$y.censored <- y.censored[y.start:y.end]
         block.list[[i]]$data$r <- solve(R[y.start:y.end, y.start:y.end])
         block.h <- Construct.H.multisite(site.ids[i], var.names, obs.mean[[t]])
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index c83c58dbd7c..08a05c49379 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -492,7 +492,7 @@ sda.enkf.multisite <- function(settings,
       ###-------------------------------------------------------------------###---- 
       #To trigger the analysis function with free run, you need to first specify the control$forceRun as TRUE,
       #Then specify the settings$state.data.assimilation$by.site as TRUE, and finally
-      if (!is.null(obs.mean[[t]][[1]]) & as.logical(settings$state.data.assimilation$free.run) | control$forceRun) {
+      if (!is.null(obs.mean[[t]][[1]]) | as.logical(settings$state.data.assimilation$free.run) & control$forceRun) {
         # TODO: as currently configured, Analysis runs even if all obs are NA, 
         #  which clearly should be triggering the `else` of this if, but the
         #  `else` has not been invoked in a while an may need updating
diff --git a/modules/uncertainty/R/get.parameter.samples.R b/modules/uncertainty/R/get.parameter.samples.R
index 7853530f292..3aadb3d9ff4 100644
--- a/modules/uncertainty/R/get.parameter.samples.R
+++ b/modules/uncertainty/R/get.parameter.samples.R
@@ -165,7 +165,7 @@ get.parameter.samples <- function(settings,
       if (prior %in% param.names[[i]]) {
         if (prior=="som_respiration_rate") {
          samples <- trait.mcmc[[prior]] %>% purrr::map(~ .x[,'beta.o']) %>% unlist() %>% as.matrix()
-         samples <- samples[(samples>0)&(samples<0.007)]
+         samples <- samples[(samples>0)&(samples<0.003)]
     } else {
         samples <- trait.mcmc[[prior]] %>% purrr::map(~ .x[,'beta.o']) %>% unlist() %>% as.matrix()
       } }

From 8d5ba2fdc1f5cf17dd419a86cb81d896b15f970d Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Thu, 18 Jan 2024 17:29:48 -0500
Subject: [PATCH 12/32] Fix some bugs of block-based SDA run and update
 phenology_MODIS_extract code

---
 models/sipnet/R/model2netcdf.SIPNET.R         |   5 +-
 models/sipnet/man/model2netcdf.SIPNET.Rd      |   3 +-
 .../assim.sequential/R/Analysis_sda_block.R   |  47 +--
 .../assim.sequential/R/sda.enkf_MultiSite.R   |   3 +-
 .../qianyu(cherry)/pecan_HARV_site.xml        | 331 ++++++++++++++++++
 modules/data.remote/NAMESPACE                 |   1 +
 .../data.remote/R/phenology_MODIS_extract.R   | 120 ++++---
 .../man/phenology_MODIS_extract.Rd            |  37 ++
 8 files changed, 461 insertions(+), 86 deletions(-)
 create mode 100644 modules/assim.sequential/inst/sda_backup/qianyu(cherry)/pecan_HARV_site.xml
 create mode 100644 modules/data.remote/man/phenology_MODIS_extract.Rd

diff --git a/models/sipnet/R/model2netcdf.SIPNET.R b/models/sipnet/R/model2netcdf.SIPNET.R
index 0cc58b8f1eb..da755a61bd2 100644
--- a/models/sipnet/R/model2netcdf.SIPNET.R
+++ b/models/sipnet/R/model2netcdf.SIPNET.R
@@ -97,6 +97,7 @@ sipnet2datetime <- function(sipnet_tval, base_year, base_month = 1,
 ##' @param revision model revision
 ##' @param overwrite Flag for overwriting nc files or not
 ##' @param conflict Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next.
+##' conflict is set to TRUE to enable the assimilation of monthly data.
 ##' @param prefix prefix to read the output files
 ##' @param delete.raw Flag to remove sipnet.out files, FALSE = do not remove files TRUE = remove files
 ##'
@@ -296,13 +297,15 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
       close(varfile)
       ncdf4::nc_close(nc)
       
-      #merge nc files
+      #merge nc files of the same year together to enable the assimilation of subannual data
       if(file.exists(file.path(outdir, "previous.nc"))){
         files <- c(file.path(outdir, "previous.nc"), file.path(outdir, "current.nc"))
       }else{
         files <- file.path(outdir, "current.nc")
       }
       mergeNC(files = files, outfile = file.path(outdir, paste(y, "nc", sep = ".")))
+      #The command "cdo" in mergeNC will automatically rename "time_bounds" to "time_bnds". However, "time_bounds" is used 
+      #in read_restart codes later. So we need to read the new NetCDF file and convert the variable name back. 
       nc<- ncdf4::nc_open(file.path(outdir, paste(y, "nc", sep = ".")),write=TRUE)
       nc<-ncdf4::ncvar_rename(nc,"time_bnds","time_bounds")
       ncdf4::ncatt_put(nc, "time", "bounds","time_bounds", prec=NA)
diff --git a/models/sipnet/man/model2netcdf.SIPNET.Rd b/models/sipnet/man/model2netcdf.SIPNET.Rd
index 697ec32d1bd..67a2b89466a 100644
--- a/models/sipnet/man/model2netcdf.SIPNET.Rd
+++ b/models/sipnet/man/model2netcdf.SIPNET.Rd
@@ -36,7 +36,8 @@ model2netcdf.SIPNET(
 
 \item{overwrite}{Flag for overwriting nc files or not}
 
-\item{conflict}{Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next.}
+\item{conflict}{Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next.
+conflict is set to TRUE to enable the assimilation of monthly data.}
 }
 \description{
 Convert SIPNET output to netCDF
diff --git a/modules/assim.sequential/R/Analysis_sda_block.R b/modules/assim.sequential/R/Analysis_sda_block.R
index b99df223235..6cb603f5538 100644
--- a/modules/assim.sequential/R/Analysis_sda_block.R
+++ b/modules/assim.sequential/R/Analysis_sda_block.R
@@ -58,27 +58,9 @@ analysis_sda_block <- function (settings, block.list.all, X, obs.mean, obs.cov,
   
   #parallel for loop over each block.
   PEcAn.logger::logger.info(paste0("Running MCMC ", "for ", length(block.list.all[[t]]), " blocks"))
- # if ("try-error" %in% class(try(block.list.all[[t]] <- furrr::future_map(block.list.all[[t]], MCMC_block_function, .progress = T)))) {
- #  PEcAn.logger::logger.error("Something wrong within the MCMC_block_function function.")
-  # return(0)
- # }
-
-  if  (length(block.list.all[[t]]) >1){
-  cores <- parallel::detectCores()
-  cl <- parallel::makeCluster(cores)
-  doSNOW::registerDoSNOW(cl)
-  #progress bar
-  pb <- utils::txtProgressBar(min=1, max=length(block.list.all[[t]]), style=3)
-  progress <- function(n) utils::setTxtProgressBar(pb, n)
-  opts <- list(progress=progress)
-  PEcAn.logger::logger.info("Submitting jobs!")
-   block.list.all[[t]]<-foreach::foreach(block = block.list.all[[t]], .packages="Kendall", .options.snow=opts) %dopar% {
-    library(nimble)
-    library(purrr)
-    MCMC_block_function(block)
-}
-}else{
-    block.list.all[[t]][[1]]<-MCMC_block_function(block=block.list.all[[t]][[1]])
+  if ("try-error" %in% class(try(block.list.all[[t]] <- furrr::future_map(block.list.all[[t]], MCMC_block_function, .progress = T)))) {
+   PEcAn.logger::logger.severe("Something wrong within the MCMC_block_function function.")
+   return(0)
   }
   PEcAn.logger::logger.info("Completed!")
   
@@ -132,7 +114,7 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
   }
   
   #distance calculations and localization
-  site.locs <- settings$run %>%
+  site.locs <- settings  %>% map(~.x[['run']] ) %>%
     purrr::map('site') %>%
     purrr::map_dfr(~c(.x[['lon']],.x[['lat']]) %>% as.numeric)%>%
     t %>%
@@ -318,11 +300,11 @@ MCMC_Init <- function (block.list, X) {
     }
     #initialize q.
     #if we want the vector q.
-    if (block.list[[i]]$constant$q.type == 1) {
+    if (block.list[[i]]$constant$q.type == 3) {
       for (j in seq_along(block.list[[i]]$data$y.censored)) {
         block.list[[i]]$Inits$q <- c(block.list[[i]]$Inits$q, stats::rgamma(1, shape = block.list[[i]]$data$aq[j], rate = block.list[[i]]$data$bq[j]))
       }
-    } else if (block.list[[i]]$constant$q.type == 2) {
+    } else if (block.list[[i]]$constant$q.type == 4) {
       #if we want the wishart Q.
       if ("try-error" %in% class(try(block.list[[i]]$Inits$q <- 
                                      stats::rWishart(1, df = block.list[[i]]$data$bq, Sigma = block.list[[i]]$data$aq)[,,1], silent = T))) {
@@ -359,12 +341,12 @@ MCMC_block_function <- function(block) {
   #because it has a better performance of MVN sampling
   samplerLists <- conf$getSamplers()
   samplerNumberOffset <- length(samplerLists)
-  if (block$constant$q.type == 1) {
+  if (block$constant$q.type == 3) {
     #if we have vector q
     #only X.mod should be sampled with ess sampler.
     X.mod.ind <- which(grepl("X.mod", samplerLists %>% purrr::map(~ .x$target) %>% unlist()))
     samplerLists[[X.mod.ind]]$setName("ess")
-  } else if (block$constant$q.type == 2) {
+  } else if (block$constant$q.type == 4) {
     #if we have wishart q
     #everything should be sampled with ess sampler.
     samplerLists %>% purrr::map(function(l){l$setName("ess")})
@@ -391,11 +373,10 @@ MCMC_block_function <- function(block) {
   
   #run MCMC
   dat <- runMCMC(Cmcmc, niter = block$MCMC$niter, nburnin = block$MCMC$nburnin, thin = block$MCMC$nthin, nchains = block$MCMC$nchain)
-  
   #update aq, bq, mua, and pa
   M <- colMeans(dat)
   block$update$aq <- block$Inits$q
-  if (block$constant$q.type == 1) {
+  if (block$constant$q.type == 3) {
     #if it's a vector q case
     aq <- bq <- rep(NA, length(block$data$y.censored))
     for (i in seq_along(aq)) {
@@ -408,7 +389,7 @@ MCMC_block_function <- function(block) {
     block$aqq[block$constant$H, block$t+1] <- aq
     block$bqq[,block$t+1] <- block$bqq[, block$t]
     block$bqq[block$constant$H, block$t+1] <- bq
-  } else if (block$constant$q.type == 2) {
+  } else if (block$constant$q.type == 4) {
     #previous updates
     mq <- dat[,  grep("q", colnames(dat))]  # Omega, Precision
     q.bar <- matrix(apply(mq, 2, mean),
@@ -481,7 +462,7 @@ update_q <- function (block.list.all, t, nt, aqq.Init = NULL, bqq.Init = NULL, M
       for (i in seq_along(block.list)) {
         nvar <- length(block.list[[i]]$data$muf)
         nobs <- length(block.list[[i]]$data$y.censored)
-        if (block.list[[i]]$constant$q.type == 1) {
+        if (block.list[[i]]$constant$q.type == 3) {
           #initialize aqq and bqq for nt
           if (!is.null(aqq.Init) && !is.null(bqq.Init)) {
             block.list[[i]]$aqq <- array(aqq.Init, dim = c(nvar, nt + 1))
@@ -493,7 +474,7 @@ update_q <- function (block.list.all, t, nt, aqq.Init = NULL, bqq.Init = NULL, M
           #update aq and bq based on aqq and bqq
           block.list[[i]]$data$aq <- block.list[[i]]$aqq[block.list[[i]]$constant$H, t]
           block.list[[i]]$data$bq <- block.list[[i]]$bqq[block.list[[i]]$constant$H, t]
-        } else if (block.list[[i]]$constant$q.type == 2) {
+        } else if (block.list[[i]]$constant$q.type == 4) {
           #initialize aqq and bqq for nt
           block.list[[i]]$aqq <- array(1, dim = c(nvar, nvar, nt + 1))
           block.list[[i]]$aqq[,,t] <- stats::toeplitz((nvar:1)/nvar)
@@ -513,14 +494,14 @@ update_q <- function (block.list.all, t, nt, aqq.Init = NULL, bqq.Init = NULL, M
       for (i in seq_along(block.list)) {
         nvar <- length(block.list[[i]]$data$muf)
         nobs <- length(block.list[[i]]$data$y.censored)
-        if (block.list[[i]]$constant$q.type == 1) {
+        if (block.list[[i]]$constant$q.type == 3) {
           #copy previous aqq and bqq to the current t
           block.list[[i]]$aqq <- block.list.pre[[i]]$aqq
           block.list[[i]]$bqq <- block.list.pre[[i]]$bqq
           #update aq and bq
           block.list[[i]]$data$aq <- block.list[[i]]$aqq[block.list[[i]]$constant$H, t]
           block.list[[i]]$data$bq <- block.list[[i]]$bqq[block.list[[i]]$constant$H, t]
-        } else if (block.list[[i]]$constant$q.type == 2) {
+        } else if (block.list[[i]]$constant$q.type == 4) {
           #initialize aqq and bqq for nt
           block.list[[i]]$aqq <- block.list.pre[[i]]$aqq
           block.list[[i]]$bqq <- block.list.pre[[i]]$bqq
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index b044026fa1a..79c56e38b26 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -523,13 +523,12 @@ sda.enkf.multisite <- function(settings,
           if (is.null(control$MCMC.args)) {
             MCMC.args <- list(niter = 1e5,
                               nthin = 10,
-                              nchain = 3,
+                              nchain = 1,
                               nburnin = 5e4)
           } else {
             MCMC.args <- control$MCMC.args
           }
           #running analysis function.
-          save(X,file = file.path(settings$outdir, "X.data"))
           enkf.params[[obs.t]] <- analysis_sda_block(settings, block.list.all, X, obs.mean, obs.cov, t, nt, MCMC.args, pre_enkf_params)
           enkf.params[[obs.t]] <- c(enkf.params[[obs.t]], RestartList = list(restart.list %>% stats::setNames(site.ids)))
           block.list.all <- enkf.params[[obs.t]]$block.list.all
diff --git a/modules/assim.sequential/inst/sda_backup/qianyu(cherry)/pecan_HARV_site.xml b/modules/assim.sequential/inst/sda_backup/qianyu(cherry)/pecan_HARV_site.xml
new file mode 100644
index 00000000000..18fb501126e
--- /dev/null
+++ b/modules/assim.sequential/inst/sda_backup/qianyu(cherry)/pecan_HARV_site.xml
@@ -0,0 +1,331 @@
+<?xml version="1.0"?>
+<pecan.multi>
+ <state.data.assimilation>
+  <process.variance>TRUE</process.variance>
+  <adjustment>FALSE</adjustment>
+  <censored.data>FALSE</censored.data>
+  <FullYearNC>TRUE</FullYearNC>
+  <NC.Overwrite>FALSE</NC.Overwrite>
+  <NC.Prefix>sipnet.out</NC.Prefix>
+  <q.type>vector</q.type>
+  <by.site>TRUE</by.site>
+  <free.run>FALSE</free.run>
+  <aqq.Init>1</aqq.Init>
+  <bqq.Init>1</bqq.Init>
+  <scalef>1</scalef>
+  <chains>5</chains>
+  <data>
+   <format_id>1000000040</format_id>
+   <input.id>1000013298</input.id>
+  </data>
+  <state.variables>
+   <variable>
+    <variable.name>AbvGrndWood</variable.name>
+    <unit>MgC/ha</unit>
+    <min_value>0</min_value>
+    <max_value>9999</max_value>
+   </variable>
+   <variable>
+    <variable.name>LAI</variable.name>
+    <unit></unit>
+    <min_value>0</min_value>
+    <max_value>9999</max_value>
+   </variable>
+   <variable>
+    <variable.name>SoilMoistFrac</variable.name>
+    <unit></unit>
+    <min_value>0</min_value>
+    <max_value>100</max_value>
+   </variable>
+   <variable>
+    <variable.name>TotSoilCarb</variable.name>
+    <unit>kg/m^2</unit>
+    <min_value>0</min_value>
+    <max_value>9999</max_value>
+   </variable>
+  </state.variables>
+  <forecast.time.step>month</forecast.time.step>
+  <start.date>2017/01/01</start.date>
+  <end.date>2021/12/31</end.date>
+  <Obs_Prep>
+   <Landtrendr_AGB>
+    <AGB_indir>/data2/RS_GIS_Data/LandTrendr/LandTrendr_AGB_data/v1</AGB_indir>
+    <timestep>
+     <unit>year</unit>
+     <num>1</num>
+    </timestep>
+    <allow_download>TRUE</allow_download>
+    <export_csv>TRUE</export_csv>
+   </Landtrendr_AGB>
+   <MODIS_LAI>
+    <search_window>30</search_window>
+    <timestep>
+     <unit>year</unit>
+     <num>1</num>
+    </timestep>
+    <export_csv>TRUE</export_csv>
+    <run_parallel>TRUE</run_parallel>
+   </MODIS_LAI>
+   <SMAP_SMP>
+    <search_window>30</search_window>
+    <timestep>
+     <unit>year</unit>
+     <num>1</num>
+    </timestep>
+    <export_csv>TRUE</export_csv>
+    <update_csv>FALSE</update_csv>
+   </SMAP_SMP>
+   <Soilgrids_SoilC>
+    <timestep>
+     <unit>year</unit>
+     <num>1</num>
+    </timestep>
+    <export_csv>TRUE</export_csv>
+   </Soilgrids_SoilC>
+   <start.date>2010-07-15</start.date>
+   <end.date>2021-07-15</end.date>
+   <outdir>/data2/Projects/NASA_CMS/neon_soc_sda/OBS_4var</outdir>
+   <timestep>
+    <unit>year</unit>
+    <num>1</num>
+   </timestep>
+  </Obs_Prep>
+ </state.data.assimilation>
+ <info>
+  <notes></notes>
+  <userid>-1</userid>
+  <username></username>
+  <date>2017/12/06 21:19:33 +0000</date>
+ </info>
+ <outdir>/data2/Projects/NASA_CMS/neon_soc_sda/SDA_single_site_LAI_SWC_monthly_HARV_PR</outdir>
+ <rundir>/data2/Projects/NASA_CMS/neon_soc_sda/SDA_single_site_LAI_SWC_monthly_HARV_PR/run</rundir>
+ <modeloutdir>/data2/Projects/NASA_CMS/neon_soc_sda/SDA_single_site_LAI_SWC_monthly_HARV_PR/out</modeloutdir>
+ <database>
+  <bety>
+   <user>bety</user>
+   <password>bety</password>
+   <host>localhost</host>
+   <port>5432</port>
+   <dbname>betydb</dbname>
+   <driver>PostgreSQL</driver>
+   <write>FALSE</write>
+  </bety>
+  <dbfiles>/data/pecan_dbfiles/</dbfiles>
+ </database>
+ <pfts>
+  <pft>
+   <name>temperate.deciduous.HPDA</name>
+   <constants>
+    <num>1</num>
+   </constants>
+   <outdir>/data2/Projects/NASA_CMS/neon_soc_sda/setup_run/PFT_Files_new/temperate</outdir>
+  </pft>
+  <pft>
+   <name>boreal.coniferous</name>
+   <constants>
+    <num>1</num>
+   </constants>
+   <outdir>/data2/Projects/NASA_CMS/neon_soc_sda/setup_run/PFT_Files_new/conifer</outdir>
+  </pft>
+  <pft>
+   <name>semiarid.grassland_HPDA</name>
+   <constants>
+    <num>1</num>
+   </constants>
+   <outdir>/data2/Projects/NASA_CMS/neon_soc_sda/setup_run/PFT_Files_new/grass</outdir>
+  </pft>
+ </pfts>
+ <meta.analysis>
+  <iter>3000</iter>
+  <random.effects>FALSE</random.effects>
+ </meta.analysis>
+ <ensemble>
+  <size>100</size>
+  <variable>NPP</variable>
+  <samplingspace>
+   <parameters>
+    <method>lhc</method>
+   </parameters>
+   <met>
+    <method>sampling</method>
+   </met>
+  </samplingspace>
+ </ensemble>
+ <model>
+  <id>1000000014</id>
+  <type>SIPNET</type>
+  <default.param>/data2/Projects/NASA_CMS/neon_soc_sda/Bartlett.param</default.param>
+  <revision>r136</revision>
+  <delete.raw>FALSE</delete.raw>
+  <binary>/data2/qli1/sipnet_prescribed_LAI.r136</binary>
+  <leaf_phenology>
+   <outdir>/data2/Projects/NASA_CMS/neon_soc_sda</outdir>
+  </leaf_phenology>  
+ </model>
+ <host>
+  <name>localhost</name>
+  <scratchdir>/scratch</scratchdir>
+  <prerun>module load geos/3.10.2-gcc850 geotiff/1.7.1-gcc850 libtiff/4.3.0 hdf4/4.2.15-gcc850 hdf5/1.12.1-gcc850 jags/4.3.0-gcc850 libxml2/2.9.13 libxslt/1.1.35 netcdf/4.8.1-gcc850 proj/7.2.1-gcc850 sqlite/3.38.1 udunits/2.2.28-gcc850 gdal/3.4.2-gcc850</prerun>
+  <cdosetup>module load cdo/2.1.0-gcc850</cdosetup>
+  <qsub>sbatch -t 36:00:00 -v -J @NAME@ -o @STDOUT@ -e @STDERR@</qsub>
+  <qsub.jobid>Submitted batch job ([[:digit:]]+)</qsub.jobid>
+  <qstat>squeue --job @JOBID@ || echo DONE</qstat>
+  <modellauncher>
+   <binary>/data/software/modellauncher/modellauncher</binary>
+   <qsub.extra>--ntasks 45</qsub.extra>
+  </modellauncher>
+  <outdir>/data2/Projects/NASA_CMS/neon_soc_sda/SDA_single_site_LAI_SWC_monthly_HARV_PR/out</outdir>
+  <rundir>/data2/Projects/NASA_CMS/neon_soc_sda/SDA_single_site_LAI_SWC_monthly_HARV_PR/run</rundir>
+ </host>
+ <settings.info>
+  <deprecated.settings.fixed>TRUE</deprecated.settings.fixed>
+  <settings.updated>TRUE</settings.updated>
+  <checked>TRUE</checked>
+ </settings.info>
+ <run>
+  <settings.1000004945>
+   <site>
+    <id>1000004945</id>
+    <met.start>2017/01/01</met.start>
+    <met.end>2021/12/31</met.end>
+    <lat>42.5369</lat>
+    <lon>-72.17266</lon>
+    <name>Harvard Forest (NEON-D01-HARV)</name>
+   </site>
+   <start.date>2017/01/01</start.date>
+   <end.date>2021/12/31</end.date>
+   <inputs>
+    <met>
+     <source>ERA5</source>
+     <output>SIPNET</output>
+     <id></id>
+     <path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.1.2010-01-01.2021-12-31.clim</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.10.2010-01-01.2021-12-31.clim</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.2.2010-01-01.2021-12-31.clim</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.3.2010-01-01.2021-12-31.clim</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.4.2010-01-01.2021-12-31.clim</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.5.2010-01-01.2021-12-31.clim</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.6.2010-01-01.2021-12-31.clim</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.7.2010-01-01.2021-12-31.clim</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.8.2010-01-01.2021-12-31.clim</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/ERA5_forcing/1000004945/ERA5.9.2010-01-01.2021-12-31.clim</path>
+     </path>
+    </met>
+    <pft.site>
+     <path>/data2/Projects/NASA_CMS/neon_soc_sda/site_pft.csv</path>
+    </pft.site>
+    <poolinitcond>
+     <source>NEON_veg</source>
+     <output>poolinitcond</output>
+     <ensemble>100</ensemble>
+     <path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_100.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_99.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_98.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_97.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_96.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_95.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_94.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_93.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_92.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_91.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_90.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_89.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_88.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_87.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_86.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_85.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_84.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_83.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_82.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_81.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_80.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_79.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_78.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_77.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_76.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_75.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_74.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_73.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_72.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_71.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_70.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_69.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_68.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_67.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_66.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_65.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_64.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_63.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_62.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_61.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_60.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_59.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_58.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_57.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_56.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_55.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_54.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_53.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_52.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_51.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_50.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_49.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_48.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_47.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_46.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_45.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_44.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_43.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_42.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_41.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_40.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_39.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_38.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_37.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_36.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_35.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_34.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_33.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_32.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_31.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_30.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_29.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_28.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_27.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_26.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_25.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_24.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_23.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_22.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_21.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_20.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_19.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_18.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_17.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_16.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_15.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_14.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_13.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_12.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_11.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_10.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_9.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_8.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_7.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_6.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_5.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_4.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_3.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_2.nc</path>
+      <path>/data2/Projects/NASA_CMS/neon_soc_sda/neon_datasets/output_ic_data_0330/1000004945/IC_site_1-4945_1.nc</path>
+     </path>
+    </poolinitcond>
+   </inputs>
+    </settings.1000004945>
+ </run>
+ <multisettings>
+  <multisettings>run</multisettings>
+ </multisettings>
+</pecan.multi>
diff --git a/modules/data.remote/NAMESPACE b/modules/data.remote/NAMESPACE
index 2e0c8d51c0b..a45db9f3968 100644
--- a/modules/data.remote/NAMESPACE
+++ b/modules/data.remote/NAMESPACE
@@ -9,6 +9,7 @@ export(download.NLCD)
 export(download.thredds.AGB)
 export(extract.LandTrendr.AGB)
 export(extract_NLCD)
+export(phenology_MODIS_extract)
 export(remote_process)
 importFrom(foreach,"%do%")
 importFrom(foreach,"%dopar%")
diff --git a/modules/data.remote/R/phenology_MODIS_extract.R b/modules/data.remote/R/phenology_MODIS_extract.R
index a46876fe47f..93d10aef988 100644
--- a/modules/data.remote/R/phenology_MODIS_extract.R
+++ b/modules/data.remote/R/phenology_MODIS_extract.R
@@ -1,5 +1,30 @@
-phenology_MODIS_extract<- function(settings,NCore = NULL, run_parallel = TRUE,update_csv = TRUE){
-  outdir<-settings$model$leaf_phenology$outdir
+##' Get MODIS phenology data by date and location  
+##' 
+##' @name phenology_MODIS_extract
+##' @title phenology_MODIS_extract
+##' @export
+##' @param settings PEcAn settings
+##' @param run_parallel optional method to download data parallely. Only works if more than 1 
+##' site is needed and there are >1 CPUs available.
+##' @param ncores number of cpus to use if run_parallel is set to TRUE. If you do not know the 
+##' number of CPU's available, enter NULL. 
+##' 
+##' The output file will be added to the folder defined in settings.
+##' 
+##' @examples
+##' \dontrun{
+##' settings <-  PEcAn.settings::read.settings("pecan.xml")
+##' test_modis_phenology <- PEcAn.data.remote::phenology_MODIS_extract(
+##'   settings = settings,
+##'   run_parallel = TRUE,
+##'   ncores = NULL)
+##' }
+##' @return NONE
+##' @author Qianyu Li
+##'
+
+phenology_MODIS_extract<- function(settings,run_parallel = TRUE,ncores = NULL){ 
+  #Get site id, lat, lon and site name for each site
   site_info <- settings %>% 
     purrr::map(~.x[['run']] ) %>% 
     purrr::map('site')%>% 
@@ -8,86 +33,83 @@ phenology_MODIS_extract<- function(settings,NCore = NULL, run_parallel = TRUE,up
       site.list$lat <- as.numeric(site.list$lat)
       site.list$lon <- as.numeric(site.list$lon)
       list(site_id=site.list$id, lat=site.list$lat, lon=site.list$lon, site_name=site.list$name)
-    })%>% 
+      })%>% 
     dplyr::bind_rows() %>% 
     as.list()
-  
-  start_date <- as.Date(settings$state.data.assimilation$start.date, tz="UTC")
-  end_date <- as.Date(settings$state.data.assimilation$end.date, tz="UTC")
-  
-  if(file.exists(file.path(outdir, "leaf_phenology_pecan.csv"))){
-    Previous_CSV <- utils::read.csv(file.path(outdir, "leaf_phenology_pecan.csv"))
-  }
-
-  start_YEARDOY <- paste0(lubridate::year(start_date),sprintf("%03d", 1))#using 1 and 365 DOY to avoid any possible missing data.
-  end_YEARDOY <- paste0(lubridate::year(end_date),sprintf("%03d", 365))
-  
-  leafon_data <- PEcAn.data.remote::call_MODIS(outdir = NULL, 
+  #Set up the outdir to store leaf phenology data file           
+  outdir<-settings$model$leaf_phenology$outdir
+  if (!is.null(outdir)) {
+    #Set up the start and end date of the extraction
+    start_date <- as.Date(settings$state.data.assimilation$start.date, tz="UTC")
+    end_date <- as.Date(settings$state.data.assimilation$end.date, tz="UTC")
+    start_YEARDOY <- paste0(lubridate::year(start_date),sprintf("%03d", 1))
+    end_YEARDOY <- paste0(lubridate::year(end_date),sprintf("%03d", 365))
+    #Extracting leaf-on, leaf-off, and Quality Assurance (QA) from MODIS MCD12Q2 product
+    leafon_data <- PEcAn.data.remote::call_MODIS(outdir = NULL, 
                                                var = "Greenup", 
                                                site_info = site_info, 
                                                product_dates = c(start_YEARDOY, end_YEARDOY),
                                                run_parallel = as.logical(run_parallel), 
-                                               ncores = NCore, 
+                                               ncores = ncores, 
                                                product = "MCD12Q2", 
                                                band = "MidGreenup.Num_Modes_01",
                                                package_method = "MODISTools", 
                                                QC_filter = FALSE, 
                                                progress = TRUE)
-  leafoff_data <- PEcAn.data.remote::call_MODIS(outdir = NULL, 
+
+    leafoff_data <- PEcAn.data.remote::call_MODIS(outdir = NULL, 
                                                 var = "Greendown", 
                                                 site_info = site_info, 
                                                 product_dates = c(start_YEARDOY, end_YEARDOY),
                                                 run_parallel = as.logical(run_parallel), 
-                                                ncores = NCore, 
+                                                ncores = ncores, 
                                                 product = "MCD12Q2", 
                                                 band = "MidGreendown.Num_Modes_01",
                                                 package_method = "MODISTools", 
                                                 QC_filter = FALSE, 
                                                 progress = TRUE)
-  
-  
-  qc_data <- PEcAn.data.remote::call_MODIS(outdir = NULL, 
+
+    qc_data <- PEcAn.data.remote::call_MODIS(outdir = NULL, 
                                            var = "QA", 
                                            site_info = site_info, 
                                            product_dates = c(start_YEARDOY, end_YEARDOY),
                                            run_parallel = as.logical(run_parallel), 
-                                           ncores = NCore, 
+                                           ncores = ncores, 
                                            product = "MCD12Q2", 
                                            band = "QA_Detailed.Num_Modes_01",
                                            package_method = "MODISTools", 
                                            QC_filter = FALSE, 
                                            progress = TRUE)
-  UnpackDetailedQA <- function(x){
+
+    #Example function for unpacking QA for each date phenometric from MCD12Q2_Collection6_UserGuide
+    #in the order: Greenup, MidGreenup, Maturity, Peak, Senescence, MidGreendown, Dormancy                                        
+    UnpackDetailedQA <- function(x){
     bits <- as.integer(intToBits(x))
     quals <- sapply(seq(1, 16, by=2), function(i) sum(bits[i:(i+1)] * 2^c(0, 1)))[1:7]
     return(quals)
-  }
-  
-  qa_ind<-matrix(NA,nrow=length(qc_data[,1]),ncol=2)
-  for (i in seq_along(qc_data$data)){
-    qa_ind[i,]<-cbind(UnpackDetailedQA(qc_data$data[i])[2],UnpackDetailedQA(qc_data$data[i])[6])
-  }
+    }
+    #Extract QA for leaf-on and leaf-off dates. Values are in the range 0–3 corresponding to “best”, “good”, “fair”, and “poor”.
+    qa_ind<-matrix(NA,nrow=length(qc_data[,1]),ncol=2)
+    for (i in seq_along(qc_data$data)){
+      qa_ind[i,]<-cbind(UnpackDetailedQA(qc_data$data[i])[2],UnpackDetailedQA(qc_data$data[i])[6])
+    }
   
-  leafphdata <- cbind(leafon_data %>% as.data.frame %>% dplyr::select("calendar_date", "site_id", "lat", "lon", "data"), leafoff_data$data,qa_ind)%>% 
+    leafphdata <- cbind(leafon_data %>% as.data.frame %>% dplyr::select("calendar_date", "site_id", "lat", "lon", "data"), leafoff_data$data,qa_ind)%>% 
     `colnames<-`(c("year", "site_id", "lat", "lon", "leafonday","leafoffday","leafon_qa","leafoff_qa"))
-  leafphdata$leafonday[leafphdata$leafonday==32767]<-NA
-  leafphdata$leafoffday[leafphdata$leafoffday==32767]<-NA
-  leafphdata$leafonday[leafphdata$leafon_qa==3]<-NA
-  leafphdata$leafoffday[leafphdata$leafoff_qa==3]<-NA
-  leafphdata$leafonday[leafphdata$leafonday>leafphdata$leafoffday]<-NA
-  leafphdata$leafonday<-lubridate::yday(as.Date(leafphdata$leafonday))
-  leafphdata$leafoffday<-lubridate::yday(as.Date(leafphdata$leafoffday))
-  leafphdata$year<-lubridate::year(leafphdata$year)
-  leafphdata$site_id<-as.character(leafphdata$site_id)
-  
-  if(as.logical(update_csv)){
-    if(exists("Previous_CSV")){#we already read the csv file previously.
-      Current_CSV <- rbind(Previous_CSV, leafphdata)
-      Current_CSV <- Current_CSV[!duplicated(paste0(Current_CSV$site_id, Current_CSV$year)),]#using site_id and date to remove duplicated records.
-      utils::write.csv(Current_CSV, file = file.path(outdir, "leaf_phenology_pecan.csv"), row.names = FALSE)
-    }else{
-      Current_CSV <- leafphdata
-      utils::write.csv(Current_CSV, file = file.path(outdir, "leaf_phenology_pecan.csv"), row.names = FALSE)
-    }
+    leafphdata$leafonday[leafphdata$leafonday==32767]<-NA #exclude the data with fill values
+    leafphdata$leafoffday[leafphdata$leafoffday==32767]<-NA
+    leafphdata$leafonday[leafphdata$leafon_qa==3]<-NA #exclude the data when QA is poor
+    leafphdata$leafoffday[leafphdata$leafoff_qa==3]<-NA
+    leafphdata$leafonday[leafphdata$leafonday>leafphdata$leafoffday]<-NA #exclude the data when leaf-on date is larger than leaf-off date
+    leafphdata$leafonday<-lubridate::yday(as.Date(leafphdata$leafonday)) #convert the dates to Day-of-Year format
+    leafphdata$leafoffday<-lubridate::yday(as.Date(leafphdata$leafoffday))
+    leafphdata$year<-lubridate::year(leafphdata$year)
+    leafphdata$site_id<-as.character(leafphdata$site_id)
+    
+    PEcAn.logger::logger.info(paste0("Storing results in: ",file.path(outdir,"leaf_phenology.csv")))
+    utils::write.csv(leafphdata, file = file.path(outdir, "leaf_phenology.csv"), row.names = FALSE)
+  }
+  else {
+    PEcAn.logger::logger.error("No output directory found. Please set it up in pecan.xml under the model tag.")
   }
 }
\ No newline at end of file
diff --git a/modules/data.remote/man/phenology_MODIS_extract.Rd b/modules/data.remote/man/phenology_MODIS_extract.Rd
new file mode 100644
index 00000000000..0beaf172f72
--- /dev/null
+++ b/modules/data.remote/man/phenology_MODIS_extract.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/phenology_MODIS_extract.R
+\name{phenology_MODIS_extract}
+\alias{phenology_MODIS_extract}
+\title{phenology_MODIS_extract}
+\usage{
+phenology_MODIS_extract(settings, run_parallel = TRUE, ncores = NULL)
+}
+\arguments{
+\item{settings}{PEcAn settings}
+
+\item{run_parallel}{optional method to download data parallely. Only works if more than 1 
+site is needed and there are >1 CPUs available.}
+
+\item{ncores}{number of cpus to use if run_parallel is set to TRUE. If you do not know the 
+number of CPU's available, enter NULL. 
+
+The output file will be added to the folder defined in settings.}
+}
+\value{
+NONE
+}
+\description{
+Get MODIS phenology data by date and location
+}
+\examples{
+\dontrun{
+settings <-  PEcAn.settings::read.settings("pecan.xml")
+test_modis_phenology <- PEcAn.data.remote::phenology_MODIS_extract(
+  settings = settings,
+  run_parallel = TRUE,
+  ncores = NULL)
+}
+}
+\author{
+Qianyu Li
+}

From 84c9d6be01b6e059139165f822ef6f1b8e1c871d Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Fri, 19 Jan 2024 18:30:57 -0500
Subject: [PATCH 13/32] Add documention and adjust some spaces

---
 models/sipnet/R/write.configs.SIPNET.R           |  4 ++--
 models/sipnet/R/write_restart.SIPNET.R           |  4 +++-
 models/sipnet/man/write_restart.SIPNET.Rd        |  4 ++++
 modules/assim.sequential/R/Analysis_sda_block.R  | 16 ++++++++--------
 modules/assim.sequential/R/Nimble_codes.R        |  3 +++
 modules/assim.sequential/R/sda.enkf_MultiSite.R  |  6 +++---
 modules/uncertainty/R/ensemble.R                 |  2 ++
 modules/uncertainty/R/get.parameter.samples.R    |  5 ++---
 .../uncertainty/man/write.ensemble.configs.Rd    |  4 ++++
 9 files changed, 31 insertions(+), 17 deletions(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index f57d5046ed1..eecd9696eff 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -465,6 +465,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       }else{
         wood_total_C <- IC$AbvGrndWood / IC$abvGrndWoodFrac
       }
+
       #Sanity check
       if (is.infinite(wood_total_C) | is.nan(wood_total_C) | wood_total_C < 0) {
         wood_total_C <- 0
@@ -564,7 +565,6 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       if (!is.na(leafOnDay) && is.numeric(leafOnDay)) {
         param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
       }
-
       ## leafOffDay
       leafOffDay <- try(ncdf4::ncvar_get(IC.nc,"date_of_senescence"),silent = TRUE)
       if (!is.na(leafOffDay) && is.numeric(leafOffDay)) {
@@ -629,4 +629,4 @@ remove.config.SIPNET <- function(main.outdir, settings) {
   } else {
     print("*** WARNING: Removal of files on remote host not yet implemented ***")
   }
-} # remove.config.SIPNET 
+} # remove.config.SIPNET
\ No newline at end of file
diff --git a/models/sipnet/R/write_restart.SIPNET.R b/models/sipnet/R/write_restart.SIPNET.R
index bbf8b3ef33d..1b06c23c63f 100755
--- a/models/sipnet/R/write_restart.SIPNET.R
+++ b/models/sipnet/R/write_restart.SIPNET.R
@@ -20,6 +20,8 @@
 ##' @param RENAME flag to either rename output file or not
 ##' @param new.params list of parameters to convert between different states 
 ##' @param inputs list of model inputs to use in write.configs.SIPNET
+##' @param obs_time obervation timepoints
+##' @param update_phenology TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA
 ##'
 ##' @description Write restart files for SIPNET. WARNING: Some variables produce illegal values < 0 and have been hardcoded to correct these values!!
 ##' 
@@ -109,7 +111,7 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
   if ("LAI" %in% variables) {
     analysis.save[[length(analysis.save) + 1]] <- new.state$LAI  
     if (new.state$LAI < 0) analysis.save[[length(analysis.save)]] <- 0
-   names(analysis.save[[length(analysis.save)]]) <- c("lai")
+    names(analysis.save[[length(analysis.save)]]) <- c("lai")
   }
   
   if (!is.null(analysis.save) && length(analysis.save)>0){
diff --git a/models/sipnet/man/write_restart.SIPNET.Rd b/models/sipnet/man/write_restart.SIPNET.Rd
index 62d5019d279..8ed94f5d2bd 100644
--- a/models/sipnet/man/write_restart.SIPNET.Rd
+++ b/models/sipnet/man/write_restart.SIPNET.Rd
@@ -36,6 +36,10 @@ write_restart.SIPNET(
 \item{new.params}{list of parameters to convert between different states}
 
 \item{inputs}{list of model inputs to use in write.configs.SIPNET}
+
+\item{obs_time}{obervation timepoints}
+
+\item{update_phenology}{TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA}
 }
 \value{
 NONE
diff --git a/modules/assim.sequential/R/Analysis_sda_block.R b/modules/assim.sequential/R/Analysis_sda_block.R
index 6cb603f5538..953d8cd6977 100644
--- a/modules/assim.sequential/R/Analysis_sda_block.R
+++ b/modules/assim.sequential/R/Analysis_sda_block.R
@@ -4,9 +4,9 @@
 ##' 
 ##' @param settings  pecan standard multi-site settings list.  
 ##' @param block.list.all Lists of forecast and analysis outputs for each time point of each block. If t=1, we initialize those outputs of each block with NULL from the `sda.enkf.multisite` function.
-##' @param X A matrix contains ensemble forecasts with the dimensions of `[ensemble number, site number * number of state variables]`. The columns are matched with the site.ids and state variable names of the inside the `FORECAST` object in the `sda.enkf.multisite` script.
-##' @param obs.mean Lists of date times named by time points, which contains lists of sites named by site ids, which contains observation means for each state variables of each site for each time point.
-##' @param obs.cov   Lists of date times named by time points, which contains lists of sites named by site ids, which contains observation covariances for all state variables of each site for each time point.
+##' @param X A matrix contains ensemble forecasts with the dimensions of `[ensemble number, site number * number of state variables]`. The columns are matched with the site.ids and state variable names of the inside the `FORECAST` object in the `sda.enkf.multisite` script. 
+##' @param obs.mean Lists of date times named by time points, which contains lists of sites named by site ids, which contains observation means for each state variables of each site for each time point. 
+##' @param obs.cov   Lists of date times named by time points, which contains lists of sites named by site ids, which contains observation covariances for all state variables of each site for each time point. 
 ##' @param t time point in format of YYYY-MM-DD.
 ##' @param nt total length of time steps, corresponding to the `nt` variable in the `sda.enkf.multisite` function.
 ##' @param MCMC.args arguments for the MCMC sampling, details can be found in the roxygen strucutre for control list in the `sda.enkf.multisite` function.
@@ -59,8 +59,8 @@ analysis_sda_block <- function (settings, block.list.all, X, obs.mean, obs.cov,
   #parallel for loop over each block.
   PEcAn.logger::logger.info(paste0("Running MCMC ", "for ", length(block.list.all[[t]]), " blocks"))
   if ("try-error" %in% class(try(block.list.all[[t]] <- furrr::future_map(block.list.all[[t]], MCMC_block_function, .progress = T)))) {
-   PEcAn.logger::logger.severe("Something wrong within the MCMC_block_function function.")
-   return(0)
+    PEcAn.logger::logger.severe("Something wrong within the MCMC_block_function function.")
+    return(0)
   }
   PEcAn.logger::logger.info("Completed!")
   
@@ -151,11 +151,11 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
       }
     }
     #create matrix the describes the support for each observed state variable at time t
-    min_max <- settings$state.data.assimilation$state.variables %>%
+    min_max <- settings$state.data.assimilation$state.variables %>% 
       purrr::map(function(state.variable){
         c(as.numeric(state.variable$min_value),
           as.numeric(state.variable$max_value))
-      }) %>% unlist() %>% as.vector() %>%
+      }) %>% unlist() %>% as.vector() %>% 
       matrix(length(settings$state.data.assimilation$state.variables), 2, byrow = T) %>%
       `rownames<-`(var.names)
     #Create y.censored and y.ind
@@ -554,4 +554,4 @@ block.2.vector <- function (block.list, X, H) {
               Pf = Pf,
               mu.a = mu.a,
               Pa = Pa))
-}
+}
\ No newline at end of file
diff --git a/modules/assim.sequential/R/Nimble_codes.R b/modules/assim.sequential/R/Nimble_codes.R
index f52da620db4..9f0de482a44 100644
--- a/modules/assim.sequential/R/Nimble_codes.R
+++ b/modules/assim.sequential/R/Nimble_codes.R
@@ -187,14 +187,17 @@ GEF.MultiSite.Nimble <-  nimbleCode({
       # Sorting out qs
       q[1:YN, 1:YN] ~ dwish(R = aq[1:YN, 1:YN], df = bq) ## aq and bq are estimated over time
     }
+
     for (i in 1:nH) {
       tmpX[i]  <- X.mod[H[i]]
       Xs[i] <- tmpX[i]
     }
     ## add process error to x model but just for the state variables that we have data and H knows who
     X[1:YN]  ~ dmnorm(Xs[1:YN], prec = q[1:YN, 1:YN])
+
     ## Likelihood
     y.censored[1:YN] ~ dmnorm(X[1:YN], prec = r[1:YN, 1:YN])
+    
     # #puting the ones that they don't have q in Xall - They come from X.model
     # # If I don't have data on then then their q is not identifiable, so we use the same Xs as Xmodel
     if(nNotH > 0){
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index 79c56e38b26..35a6994b0e0 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -350,7 +350,7 @@ sda.enkf.multisite <- function(settings,
         load(file.path(settings$outdir, "samples.Rdata"))
       }
       #reformatting params
-      new.params <- PEcAnAssimSequential:::sda_matchparam(settings, ensemble.samples, site.ids, nens)
+      new.params <- sda_matchparam(settings, ensemble.samples, site.ids, nens)
     }
       #sample met ensemble members
       #TODO: incorporate Phyllis's restart work
@@ -444,7 +444,7 @@ sda.enkf.multisite <- function(settings,
         #put building of X into a function that gets called
         max_t <- 0
         while("try-error" %in% class(
-          try(reads <- PEcAnAssimSequential:::build_X(out.configs = out.configs, 
+          try(reads <- build_X(out.configs = out.configs, 
                                settings = settings, 
                                new.params = new.params, 
                                nens = nens, 
@@ -755,4 +755,4 @@ sda.enkf.multisite <- function(settings,
 #    }
       ## MCD: I commented the above "if" out because if you are restarting from a previous forecast, this might delete the files in that earlier folder
   } ### end loop over time
-} # sda.enkf
+} # sda.enkf
\ No newline at end of file
diff --git a/modules/uncertainty/R/ensemble.R b/modules/uncertainty/R/ensemble.R
index 24fb76a9e38..8c582c01814 100644
--- a/modules/uncertainty/R/ensemble.R
+++ b/modules/uncertainty/R/ensemble.R
@@ -197,6 +197,8 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 ##' @param write.to.db logical: Record this run in BETY?
 ##' @param restart In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.
 ##' @param rename Decide if we want to rename previous output files, for example convert from sipnet.out to sipnet.2020-07-16.out.
+##' @param time obervation timepoints
+##' @param update_phenology TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA
 ##'
 ##' @return list, containing $runs = data frame of runids, $ensemble.id = the ensemble ID for these runs and $samples with ids and samples used for each tag.  Also writes sensitivity analysis configuration files as a side effect
 ##' @details The restart functionality is developed using model specific functions by calling write_restart.modelname function. First, you need to make sure that this function is already exist for your desired model.See here \url{https://pecanproject.github.io/pecan-documentation/master/pecan-models.html}
diff --git a/modules/uncertainty/R/get.parameter.samples.R b/modules/uncertainty/R/get.parameter.samples.R
index 454e1fb7d2d..cc323b69fb2 100644
--- a/modules/uncertainty/R/get.parameter.samples.R
+++ b/modules/uncertainty/R/get.parameter.samples.R
@@ -163,9 +163,8 @@ get.parameter.samples <- function(settings,
     }
     for (prior in priors) {
       if (prior %in% param.names[[i]]) {
-         samples <- trait.mcmc[[prior]] %>% purrr::map(~ .x[,'beta.o']) %>% unlist() %>% as.matrix()
-      }
-      else {
+        samples <- trait.mcmc[[prior]] %>% purrr::map(~ .x[,'beta.o']) %>% unlist() %>% as.matrix()
+      } else {
         samples <- PEcAn.priors::get.sample(prior.distns[prior, ], samples.num, q_samples[ , priors==prior])
       }
       trait.samples[[pft.name]][[prior]] <- samples
diff --git a/modules/uncertainty/man/write.ensemble.configs.Rd b/modules/uncertainty/man/write.ensemble.configs.Rd
index d888588019a..4986fc99d48 100644
--- a/modules/uncertainty/man/write.ensemble.configs.Rd
+++ b/modules/uncertainty/man/write.ensemble.configs.Rd
@@ -33,6 +33,10 @@ write.ensemble.configs(
 \item{restart}{In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.}
 
 \item{rename}{Decide if we want to rename previous output files, for example convert from sipnet.out to sipnet.2020-07-16.out.}
+
+\item{time}{obervation timepoints}
+
+\item{update_phenology}{TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA}
 }
 \value{
 list, containing $runs = data frame of runids, $ensemble.id = the ensemble ID for these runs and $samples with ids and samples used for each tag.  Also writes sensitivity analysis configuration files as a side effect

From 59947a3ddc26e3f6a73bd417a66c4c6bfa344836 Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Fri, 19 Jan 2024 20:17:54 -0500
Subject: [PATCH 14/32] Update Create_Multi_settings.R to include the path for
 setting up leaf phenology

---
 .../inst/MultiSite-Exs/SDA/Create_Multi_settings.R             | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R b/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R
index eda487b4cf1..cd86b71df03 100644
--- a/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R
+++ b/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R
@@ -198,7 +198,8 @@ template <- PEcAn.settings::Settings(list(
                          revision = "ssr",
                          delete.raw = FALSE,
                          binary = "/usr2/postdoc/istfer/SIPNET/trunk//sipnet_if",
-                         jobtemplate = "~/sipnet_geo.job"
+                         jobtemplate = "~/sipnet_geo.job",
+                         leaf_phenology=structure(list(outdir="..."))
   )),
   
   ###########################################################################

From 9136c439356d33a092888fa25a4126abf6626c1e Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-gl3.scc.bu.edu>
Date: Fri, 8 Mar 2024 17:58:11 -0500
Subject: [PATCH 15/32] Revert to multisite mapping approach

---
 modules/assim.sequential/R/Analysis_sda_block.R | 2 +-
 modules/assim.sequential/R/sda.enkf_MultiSite.R | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/modules/assim.sequential/R/Analysis_sda_block.R b/modules/assim.sequential/R/Analysis_sda_block.R
index fdff11b660d..e38684eec31 100644
--- a/modules/assim.sequential/R/Analysis_sda_block.R
+++ b/modules/assim.sequential/R/Analysis_sda_block.R
@@ -112,7 +112,7 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
     PEcAn.logger::logger.warn("The zero variances in Pf is being replaced by one fifth of the minimum variance in those matrices respectively.")
   }
   #distance calculations and localization
-  site.locs <- settings  %>% map(~.x[['run']] ) %>%
+  site.locs <- settings$run %>%
     purrr::map('site') %>%
     purrr::map_dfr(~c(.x[['lon']],.x[['lat']]) %>% as.numeric)%>%
     t %>%
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index 22d55c95db8..c82a308df3f 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -132,9 +132,9 @@ sda.enkf.multisite <- function(settings,
   #Here I'm trying to make a temp config list name and put it into map to iterate
   if(multi.site.flag){
     conf.settings<-settings
-    site.ids <- conf.settings %>% map(~.x[['run']] ) %>% purrr::map('site') %>% purrr::map('id') %>% base::unlist() %>% base::as.character()
+    site.ids <- conf.settings$run %>% purrr::map('site') %>% purrr::map('id') %>% base::unlist() %>% base::as.character()
     # a matrix ready to be sent to spDistsN1 in sp package - first col is the long second is the lat and row names are the site ids
-    site.locs <- conf.settings %>% map(~.x[['run']] ) %>% purrr::map('site') %>% purrr::map_dfr(~c(.x[['lon']],.x[['lat']]) %>% as.numeric)%>% 
+    site.locs <- conf.settings$run %>% purrr::map('site') %>% purrr::map_dfr(~c(.x[['lon']],.x[['lat']]) %>% as.numeric)%>% 
       t %>%
       `colnames<-`(c("Lon","Lat")) %>%
       `rownames<-`(site.ids)

From bb33db15fdd1bef489fa45d2d74fb6379a724252 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-ta4.scc.bu.edu>
Date: Fri, 22 Mar 2024 13:38:34 -0400
Subject: [PATCH 16/32] Fix bugs in preparing initial LAI

---
 models/sipnet/R/write.configs.SIPNET.R |  2 +-
 modules/data.land/R/prepare_pools.R    | 13 ++++++++++---
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 9d740b2a366..d5c1d084572 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -541,7 +541,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     ICs_num <- length(settings$run$inputs$poolinitcond$path)
     IC.path <- settings$run$inputs$poolinitcond$path[[sample(1:ICs_num, 1)]]
 
-    IC.pools <- PEcAn.data.land::prepare_pools(IC.path, constants = list(sla = SLA))
+    IC.pools <- PEcAn.data.land::prepare_pools(IC.path, settings, constants = list(sla = SLA))
     
     if(!is.null(IC.pools)){
       IC.nc <- ncdf4::nc_open(IC.path) #for additional variables specific to SIPNET
diff --git a/modules/data.land/R/prepare_pools.R b/modules/data.land/R/prepare_pools.R
index 258a9dc9949..6d9924ea176 100644
--- a/modules/data.land/R/prepare_pools.R
+++ b/modules/data.land/R/prepare_pools.R
@@ -7,7 +7,7 @@
 ##' @param constants list of constants; must include SLA in m2 / kg C if providing LAI for leaf carbon
 ##' @return list of pool values in kg C / m2 with generic names
 ##' @author Anne Thomas
-prepare_pools <- function(nc.path, constants = NULL){
+prepare_pools <- function(nc.path, settings, constants = NULL){
   #function to check that var was loaded (numeric) and has a valid value (not NA or negative)
   is.valid <- function(var){
     return(all(is.numeric(var) && !is.na(var) &&  var >= 0)) 
@@ -77,11 +77,18 @@ prepare_pools <- function(nc.path, constants = NULL){
           PEcAn.logger::logger.error("TotLivBiom is less than sum of AbvGrndWood and roots; will use default for leaf biomass")
         }
       }
+      #Initial LAI at the beginning of year is set as 0 for deciduous forests and grasslands
+      site_pft <- utils::read.csv(settings$run$inputs$pft.site$path)
+      site.pft.name <- site_pft$pft[site_pft$site == settings$run$site$id]
+      if (site.pft.name!="boreal.coniferous") { 
+        leaf <- 0
+        IC.params[["leaf"]] <- leaf        
+      }
+      
       #Calculate LAI given leaf and sla
       sla <- constants$sla
       if (!is.null(sla) && is.valid(leaf)) {
-        LAI <- PEcAn.utils::ud_convert(leaf * sla, "g", "kg") #convert from g to kg
-        
+        LAI <- leaf * sla
         IC.params[["LAI"]] <- LAI
       }
       

From 1e1fde65621698396dc3d74ad89a65e0ac0b7c06 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-gc2.scc.bu.edu>
Date: Mon, 25 Mar 2024 15:08:32 -0400
Subject: [PATCH 17/32] Delete files that were placed by mistake

---
 base/workflow/NAMESPACE                |   2 -
 base/workflow/R/merge_job_files.R      |  42 ----------
 base/workflow/R/qsub_parallel.R        | 104 -------------------------
 base/workflow/man/merge_job_files.Rd   |  23 ------
 base/workflow/man/qsub_parallel.Rd     |  26 -------
 modules/data.land/man/prepare_pools.Rd |   2 +-
 6 files changed, 1 insertion(+), 198 deletions(-)
 delete mode 100644 base/workflow/R/merge_job_files.R
 delete mode 100644 base/workflow/R/qsub_parallel.R
 delete mode 100644 base/workflow/man/merge_job_files.Rd
 delete mode 100644 base/workflow/man/qsub_parallel.Rd

diff --git a/base/workflow/NAMESPACE b/base/workflow/NAMESPACE
index 62b40898236..e8a4c246a4c 100644
--- a/base/workflow/NAMESPACE
+++ b/base/workflow/NAMESPACE
@@ -2,9 +2,7 @@
 
 export(create_execute_test_xml)
 export(do_conversions)
-export(merge_job_files)
 export(model_specific_tags)
-export(qsub_parallel)
 export(run.write.configs)
 export(runModule.get.trait.data)
 export(runModule.run.write.configs)
diff --git a/base/workflow/R/merge_job_files.R b/base/workflow/R/merge_job_files.R
deleted file mode 100644
index 23c84ebe09f..00000000000
--- a/base/workflow/R/merge_job_files.R
+++ /dev/null
@@ -1,42 +0,0 @@
-#' merge_job_files
-#' Merge multiple job.sh files into one larger file.
-#' @param settings PEcAn.settings object with host section.
-#' @param jobs_per_file the number of files you want to merge.
-#' @param outdir output directory of merged job files.
-#'
-#' @return
-#' @export
-#' @author Dongchen Zhang
-#'
-#' @examples
-merge_job_files <- function(settings, jobs_per_file = 10, outdir = NULL){
-  # default outdir
-  if(is.null(outdir)){
-    outdir <- file.path(settings$rundir, "merged_jobs")
-  }
-  # create folder or delete previous job files.
-  if(dir.exists(outdir)){
-    unlink(list.files(outdir, recursive = T, full.names = T))
-  }else{
-    dir.create(outdir)
-  }
-  # merge job files.
-  run_list <- readLines(con = file.path(settings$rundir, "runs.txt"))
-  job_file_paths <- list.files(settings$host$rundir, pattern = "*.sh", recursive = T, full.names = T)
-  i <- 0
-  files <- c()
-  while (i < length(job_file_paths)) {
-    jobs <- c()
-    for (j in 1:jobs_per_file) {
-      if((i+j) > length(job_file_paths)){
-        break
-      }
-      jobs <- c(jobs, readLines(job_file_paths[i+j]))
-    }
-    writeLines(jobs, con = file.path(outdir, paste0("job_", i,".sh")))
-    Sys.chmod(file.path(outdir, paste0("job_", i,".sh")))
-    files <- c(files, file.path(outdir, paste0("job_", i,".sh")))
-    i <- i + jobs_per_file
-  }
-  files
-}
diff --git a/base/workflow/R/qsub_parallel.R b/base/workflow/R/qsub_parallel.R
deleted file mode 100644
index 0318f479951..00000000000
--- a/base/workflow/R/qsub_parallel.R
+++ /dev/null
@@ -1,104 +0,0 @@
-#' qsub_parallel
-#'
-#' @param settings pecan settings object
-#' @param files allow submit jobs based on job.sh file paths.
-#' @param prefix used for detecting if jobs are completed or not.
-#' @export
-#' @examples
-#' \dontrun{
-#'   qsub_parallel(settings)
-#' }
-#' @author Dongchen Zhang
-#'
-qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out") {
-  run_list <- readLines(con = file.path(settings$rundir, "runs.txt"))
-  is_local <- PEcAn.remote::is.localhost(settings$host)
-  # loop through runs and either call start run, or launch job on remote machine
-  #parallel submit jobs
-  cores <- parallel::detectCores()
-  cl <- makeSOCKcluster(cores)
-  registerDoSNOW(cl)
-  #progress bar
-  pb <- txtProgressBar(min=1, max=ifelse(is.null(files), length(run_list), length(files)), style=3)
-  progress <- function(n) setTxtProgressBar(pb, n)
-  opts <- list(progress=progress)
-  PEcAn.logger::logger.info("Submitting jobs!")
-  # if we want to submit jobs separately.
-  if(is.null(files)){
-    jobids <- foreach(run = run_list, .packages="Kendall", .options.snow=opts, settings = rep(settings, length(files))) %dopar% {
-      run_id_string <- format(run, scientific = FALSE)
-      qsub <- settings$host$qsub
-      qsub <- gsub("@NAME@", paste0("PEcAn-", run_id_string), qsub)
-      qsub <- gsub("@STDOUT@", file.path(settings$host$outdir, run_id_string, "stdout.log"), qsub)
-      qsub <- gsub("@STDERR@", file.path(settings$host$outdir, run_id_string, "stderr.log"), qsub)
-      qsub <- strsplit(qsub, " (?=([^\"']*\"[^\"']*\")*[^\"']*$)", perl = TRUE)
-      # start the actual model run
-      cmd <- qsub[[1]]
-      if(PEcAn.remote::is.localhost(settings$host)){
-        out <- system2(cmd, file.path(settings$host$rundir, run_id_string, "job.sh"), stdout = TRUE, stderr = TRUE)
-      }else{
-        out <- remote.execute.cmd(host, cmd, file.path(settings$host$rundir, run_id_string, "job.sh"), stderr = TRUE)
-      }
-      jobid <- PEcAn.remote::qsub_get_jobid(
-        out = out[length(out)],
-        qsub.jobid = settings$host$qsub.jobid,
-        stop.on.error = TRUE)
-      return(jobid)
-    }
-  }else{
-    # if we want to submit merged job files.
-    std_out <- file.path(settings$host$outdir, "merged_stdout")
-    if(!dir.exists(std_out)){
-      dir.create(std_out)
-    }else{
-      unlink(list.files(std_out, recursive = T, full.names = T))
-    }
-    jobids <- foreach(file = files, .packages="Kendall", .options.snow=opts, settings = rep(settings, length(files))) %dopar% {
-      qsub <- settings$host$qsub
-      base_name <- basename(file)
-      num <- gsub("\\D", "", base_name)
-      name <- paste0("SDA", num)
-      qsub <- gsub("@NAME@", name, qsub)
-      qsub <- gsub("@STDOUT@", file.path(std_out, paste0("stdout", num, ".log")), qsub)
-      qsub <- gsub("@STDERR@", file.path(std_out, paste0("stderr", num, ".log")), qsub)
-      qsub <- strsplit(qsub, " (?=([^\"']*\"[^\"']*\")*[^\"']*$)", perl = TRUE)
-      cmd <- qsub[[1]]
-      if(PEcAn.remote::is.localhost(settings$host)){
-        out <- system2(cmd, file, stdout = TRUE, stderr = TRUE)
-      }else{
-        out <- remote.execute.cmd(host, cmd, file, stderr = TRUE)
-      }
-      jobid <- PEcAn.remote::qsub_get_jobid(
-        out = out[length(out)],
-        qsub.jobid = settings$host$qsub.jobid,
-        stop.on.error = TRUE)
-      return(jobid)
-    }
-  }
-  
-  PEcAn.logger::logger.info("Jobs submitted!")
-  #check if jobs are completed
-  PEcAn.logger::logger.info("Checking the qsub jobs status!")
-  PEcAn.logger::logger.info(paste0("Checking the file ", prefix))
-  ## setup progressbar
-  pb <- utils::txtProgressBar(min = 0, max = length(unlist(run_list)), style = 3)
-  pbi <- 0
-  folders <- file.path(settings$host$outdir, run_list)
-  while (length(folders) > 0) {
-    Sys.sleep(10)
-    completed_folders <- foreach(folder = folders, settings = rep(settings, length(run_list))) %dopar% {
-      if(file.exists(file.path(folder, prefix))){
-        return(folder)
-      }
-    }
-    if(length(unlist(completed_folders)) > 0){
-      Ind <- which(unlist(completed_folders) %in% folders)
-      folders <- folders[-Ind]
-      pbi <- pbi + length(completed_folders)
-      utils::setTxtProgressBar(pb, pbi)
-    }
-  }
-  close(pb)
-  stopCluster(cl)
-  PEcAn.logger::logger.info("Completed!")
-}
\ No newline at end of file
diff --git a/base/workflow/man/merge_job_files.Rd b/base/workflow/man/merge_job_files.Rd
deleted file mode 100644
index 6f6e2f4fd89..00000000000
--- a/base/workflow/man/merge_job_files.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/merge_job_files.R
-\name{merge_job_files}
-\alias{merge_job_files}
-\title{merge_job_files
-Merge multiple job.sh files into one larger file.}
-\usage{
-merge_job_files(settings, jobs_per_file = 10, outdir = NULL)
-}
-\arguments{
-\item{settings}{PEcAn.settings object with host section.}
-
-\item{jobs_per_file}{the number of files you want to merge.}
-
-\item{outdir}{output directory of merged job files.}
-}
-\description{
-merge_job_files
-Merge multiple job.sh files into one larger file.
-}
-\author{
-Dongchen Zhang
-}
diff --git a/base/workflow/man/qsub_parallel.Rd b/base/workflow/man/qsub_parallel.Rd
deleted file mode 100644
index cbd2c6614d4..00000000000
--- a/base/workflow/man/qsub_parallel.Rd
+++ /dev/null
@@ -1,26 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/qsub_parallel.R
-\name{qsub_parallel}
-\alias{qsub_parallel}
-\title{qsub_parallel}
-\usage{
-qsub_parallel(settings, files = NULL, prefix = "sipnet.out")
-}
-\arguments{
-\item{settings}{pecan settings object}
-
-\item{files}{allow submit jobs based on job.sh file paths.}
-
-\item{prefix}{used for detecting if jobs are completed or not.}
-}
-\description{
-qsub_parallel
-}
-\examples{
-\dontrun{
-  qsub_parallel(settings)
-}
-}
-\author{
-Dongchen Zhang
-}
diff --git a/modules/data.land/man/prepare_pools.Rd b/modules/data.land/man/prepare_pools.Rd
index 3fbdf7c0ac1..8a1290404f0 100644
--- a/modules/data.land/man/prepare_pools.Rd
+++ b/modules/data.land/man/prepare_pools.Rd
@@ -4,7 +4,7 @@
 \alias{prepare_pools}
 \title{prepare_pools}
 \usage{
-prepare_pools(nc.path, constants = NULL)
+prepare_pools(nc.path, settings, constants = NULL)
 }
 \arguments{
 \item{nc.path}{path to netcdf file containing standard dimensions and variables; currently supports these variables: TotLivBiom, leaf_carbon_content, LAI, AbvGrndWood, root_carbon_content, fine_root_carbon_content, coarse_root_carbon_content, litter_carbon_content, soil_organic_carbon_content, soil_carbon_content, wood_debris_carbon_content}

From 2c77fd9242e9adc5f5f97a92545f2ce59bb98ac5 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-gc2.scc.bu.edu>
Date: Mon, 25 Mar 2024 20:26:27 -0400
Subject: [PATCH 18/32] Update per Mike and Chris' comments

---
 models/sipnet/R/write.configs.SIPNET.R        | 10 ++++++++++
 models/sipnet/man/write.config.SIPNET.Rd      | 19 +++++++++++++++++++
 modules/data.land/R/prepare_pools.R           |  1 +
 modules/data.land/man/prepare_pools.Rd        |  2 ++
 .../data.remote/R/phenology_MODIS_extract.R   | 15 +++++++--------
 .../man/phenology_MODIS_extract.Rd            |  7 ++++---
 6 files changed, 43 insertions(+), 11 deletions(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index d5c1d084572..70e828d9fdc 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -11,6 +11,16 @@
 ##' Writes a configuration files for your model
 ##' @name write.config.SIPNET
 ##' @title Writes a configuration files for SIPNET model
+##' @param defaults pft
+##' @param trait.values vector of samples for a given trait
+##' @param settings PEcAn settings object
+##' @param run.id run ID
+##' @param inputs list of model inputs
+##' @param IC initial condition
+##' @param restart In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.
+##' @param spinup
+##' @param obs_time obervation timepoints
+##' @param update_phenology TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
diff --git a/models/sipnet/man/write.config.SIPNET.Rd b/models/sipnet/man/write.config.SIPNET.Rd
index 5748aa60866..f06e358e8f0 100644
--- a/models/sipnet/man/write.config.SIPNET.Rd
+++ b/models/sipnet/man/write.config.SIPNET.Rd
@@ -17,6 +17,25 @@ write.config.SIPNET(
   update_phenology = NULL
 )
 }
+\arguments{
+\item{defaults}{pft}
+
+\item{trait.values}{vector of samples for a given trait}
+
+\item{settings}{PEcAn settings object}
+
+\item{run.id}{run ID}
+
+\item{inputs}{list of model inputs}
+
+\item{IC}{initial condition}
+
+\item{restart}{In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.}
+
+\item{obs_time}{obervation timepoints}
+
+\item{update_phenology}{TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA}
+}
 \description{
 Writes a configuration files for your model
 }
diff --git a/modules/data.land/R/prepare_pools.R b/modules/data.land/R/prepare_pools.R
index 6d9924ea176..5d82bdd1db2 100644
--- a/modules/data.land/R/prepare_pools.R
+++ b/modules/data.land/R/prepare_pools.R
@@ -4,6 +4,7 @@
 ##' @export
 ##'
 ##' @param nc.path path to netcdf file containing standard dimensions and variables; currently supports these variables: TotLivBiom, leaf_carbon_content, LAI, AbvGrndWood, root_carbon_content, fine_root_carbon_content, coarse_root_carbon_content, litter_carbon_content, soil_organic_carbon_content, soil_carbon_content, wood_debris_carbon_content
+##' @param settings PEcAn settings object
 ##' @param constants list of constants; must include SLA in m2 / kg C if providing LAI for leaf carbon
 ##' @return list of pool values in kg C / m2 with generic names
 ##' @author Anne Thomas
diff --git a/modules/data.land/man/prepare_pools.Rd b/modules/data.land/man/prepare_pools.Rd
index 8a1290404f0..d2eba2f6829 100644
--- a/modules/data.land/man/prepare_pools.Rd
+++ b/modules/data.land/man/prepare_pools.Rd
@@ -9,6 +9,8 @@ prepare_pools(nc.path, settings, constants = NULL)
 \arguments{
 \item{nc.path}{path to netcdf file containing standard dimensions and variables; currently supports these variables: TotLivBiom, leaf_carbon_content, LAI, AbvGrndWood, root_carbon_content, fine_root_carbon_content, coarse_root_carbon_content, litter_carbon_content, soil_organic_carbon_content, soil_carbon_content, wood_debris_carbon_content}
 
+\item{settings}{PEcAn settings object}
+
 \item{constants}{list of constants; must include SLA in m2 / kg C if providing LAI for leaf carbon}
 }
 \value{
diff --git a/modules/data.remote/R/phenology_MODIS_extract.R b/modules/data.remote/R/phenology_MODIS_extract.R
index 93d10aef988..ab1af859057 100644
--- a/modules/data.remote/R/phenology_MODIS_extract.R
+++ b/modules/data.remote/R/phenology_MODIS_extract.R
@@ -1,7 +1,5 @@
 ##' Get MODIS phenology data by date and location  
 ##' 
-##' @name phenology_MODIS_extract
-##' @title phenology_MODIS_extract
 ##' @export
 ##' @param settings PEcAn settings
 ##' @param run_parallel optional method to download data parallely. Only works if more than 1 
@@ -9,7 +7,7 @@
 ##' @param ncores number of cpus to use if run_parallel is set to TRUE. If you do not know the 
 ##' number of CPU's available, enter NULL. 
 ##' 
-##' The output file will be added to the folder defined in settings.
+##' The output file will be saved to the folder defined in settings.
 ##' 
 ##' @examples
 ##' \dontrun{
@@ -19,7 +17,8 @@
 ##'   run_parallel = TRUE,
 ##'   ncores = NULL)
 ##' }
-##' @return NONE
+##' @return a dataframe containing phenology dates and corresponding QA. 
+##' Output column names are "year", "site_id", "lat", "lon", "leafonday","leafoffday","leafon_qa","leafoff_qa"
 ##' @author Qianyu Li
 ##'
 
@@ -38,7 +37,9 @@ phenology_MODIS_extract<- function(settings,run_parallel = TRUE,ncores = NULL){
     as.list()
   #Set up the outdir to store leaf phenology data file           
   outdir<-settings$model$leaf_phenology$outdir
-  if (!is.null(outdir)) {
+  if (is.null(outdir)) {
+    PEcAn.logger::logger.error("No output directory found. Please set it up in pecan.xml under the model tag.")
+  } else {
     #Set up the start and end date of the extraction
     start_date <- as.Date(settings$state.data.assimilation$start.date, tz="UTC")
     end_date <- as.Date(settings$state.data.assimilation$end.date, tz="UTC")
@@ -108,8 +109,6 @@ phenology_MODIS_extract<- function(settings,run_parallel = TRUE,ncores = NULL){
     
     PEcAn.logger::logger.info(paste0("Storing results in: ",file.path(outdir,"leaf_phenology.csv")))
     utils::write.csv(leafphdata, file = file.path(outdir, "leaf_phenology.csv"), row.names = FALSE)
-  }
-  else {
-    PEcAn.logger::logger.error("No output directory found. Please set it up in pecan.xml under the model tag.")
+    return(leafphdata)
   }
 }
\ No newline at end of file
diff --git a/modules/data.remote/man/phenology_MODIS_extract.Rd b/modules/data.remote/man/phenology_MODIS_extract.Rd
index 0beaf172f72..486f9d5d5cd 100644
--- a/modules/data.remote/man/phenology_MODIS_extract.Rd
+++ b/modules/data.remote/man/phenology_MODIS_extract.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/phenology_MODIS_extract.R
 \name{phenology_MODIS_extract}
 \alias{phenology_MODIS_extract}
-\title{phenology_MODIS_extract}
+\title{Get MODIS phenology data by date and location}
 \usage{
 phenology_MODIS_extract(settings, run_parallel = TRUE, ncores = NULL)
 }
@@ -15,10 +15,11 @@ site is needed and there are >1 CPUs available.}
 \item{ncores}{number of cpus to use if run_parallel is set to TRUE. If you do not know the 
 number of CPU's available, enter NULL. 
 
-The output file will be added to the folder defined in settings.}
+The output file will be saved to the folder defined in settings.}
 }
 \value{
-NONE
+a dataframe containing phenology dates and corresponding QA. 
+Output column names are "year", "site_id", "lat", "lon", "leafonday","leafoffday","leafon_qa","leafoff_qa"
 }
 \description{
 Get MODIS phenology data by date and location

From e6c683acedb1b6ae8fb9b42279788869f48e778f Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-gc2.scc.bu.edu>
Date: Tue, 26 Mar 2024 02:01:06 -0400
Subject: [PATCH 19/32] load the files into a specific envir

---
 modules/uncertainty/R/get.parameter.samples.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/modules/uncertainty/R/get.parameter.samples.R b/modules/uncertainty/R/get.parameter.samples.R
index e4821652ff2..d61cd951f09 100644
--- a/modules/uncertainty/R/get.parameter.samples.R
+++ b/modules/uncertainty/R/get.parameter.samples.R
@@ -58,10 +58,10 @@ get.parameter.samples <- function(settings,
     ## Load posteriors
     if (!is.na(posterior.files[i])) {
       # Load specified file
-      load(posterior.files[i])
+      load(posterior.files[i], envir = distns)
       fname <- file.path(outdirs[i], "post.distns.Rdata")
       if (file.exists(fname)) {
-        load(fname)
+        load(fname, envir = distns)
         prior.distns <- post.distns
       }
       if (exists("trait.mcmc") & exists("post.distns")) {

From f56eed244f88913da29da12a479e97325d6924b8 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-zp3.scc.bu.edu>
Date: Wed, 27 Mar 2024 17:33:52 -0400
Subject: [PATCH 20/32] Update default arguments in functions

---
 models/sipnet/R/write.configs.SIPNET.R |  2 +-
 models/sipnet/R/write_restart.SIPNET.R |  4 ++--
 modules/uncertainty/R/ensemble.R       | 12 ++++++------
 3 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 70e828d9fdc..6cc1f8a959d 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -24,7 +24,7 @@
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
-                                restart = NULL, spinup = NULL,obs_time=NULL,update_phenology=NULL) {
+                                restart = NULL, spinup = NULL,obs_time=NULL,update_phenology=FALSE) {
   ### WRITE sipnet.in
   template.in <- system.file("sipnet.in", package = "PEcAn.SIPNET")
   config.text <- readLines(con = template.in, n = -1)
diff --git a/models/sipnet/R/write_restart.SIPNET.R b/models/sipnet/R/write_restart.SIPNET.R
index aadb44602aa..38662852037 100755
--- a/models/sipnet/R/write_restart.SIPNET.R
+++ b/models/sipnet/R/write_restart.SIPNET.R
@@ -29,7 +29,7 @@
 ##' @return NONE
 ##' @export
 write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings, new.state,
-                                 RENAME = TRUE, new.params = FALSE, inputs, obs_time=NULL, update_phenology=NULL, verbose = FALSE) {
+                                 RENAME = TRUE, new.params = FALSE, inputs, obs_time=NULL, update_phenology=FALSE, verbose = FALSE) {
 
   rundir <- settings$host$rundir
   variables <- colnames(new.state)
@@ -134,6 +134,6 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
                                            inputs = inputs,
                                            IC = analysis.save.mat,
                                            obs_time=obs_time,
-					                                 update_phenology=update_phenology))
+                                           update_phenology=update_phenology))
   print(runid)
 } # write_restart.SIPNET
\ No newline at end of file
diff --git a/modules/uncertainty/R/ensemble.R b/modules/uncertainty/R/ensemble.R
index caa73e70fcc..1f8b836bc71 100644
--- a/modules/uncertainty/R/ensemble.R
+++ b/modules/uncertainty/R/ensemble.R
@@ -213,7 +213,7 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 ##' @export
 ##' @author David LeBauer, Carl Davidson, Hamze Dokoohaki
 write.ensemble.configs <- function(defaults, ensemble.samples, settings, model, 
-                                   clean = FALSE, write.to.db = TRUE, restart=NULL, rename = FALSE,time=NULL,update_phenology=NULL) {
+                                   clean = FALSE, write.to.db = TRUE, restart = NULL, rename = FALSE, time = NULL, update_phenology = FALSE) {
   
   con <- NULL
   my.write.config <- paste("write.config.", model, sep = "")
@@ -410,10 +410,10 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
       
       do.call(my.write.config, args = list( defaults = defaults, 
                                             trait.values = lapply(samples$parameters$samples, function(x, n) { x[n, , drop=FALSE] }, n=i), # this is the params
-                                            settings = settings, 
+                                            settings = settings,
                                             run.id = run.id,
-					                                  obs_time = time,
-					                                  update_phenology=update_phenology
+                                            obs_time = time,
+                                            update_phenology=update_phenology
       )
       )
       cat(format(run.id, scientific = FALSE), file = file.path(settings$rundir, "runs.txt"), sep = "\n", append = TRUE)
@@ -464,8 +464,8 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
                            new.params = new.params[[i]], #new.params$`646`[[i]] for debugging
                            inputs =list(met=list(path=inputs$samples[[i]])), 
                            RENAME = rename,
-			                     obs_time = time,
-			                     update_phenology = update_phenology)#for restart from previous model runs, not sharing the same outdir
+                           obs_time = time,
+                           update_phenology = update_phenology)#for restart from previous model runs, not sharing the same outdir
       )
     }
     params<-new.params

From 45007d27e0975b469556bceb0a3b006822b96e29 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-wo3.scc.bu.edu>
Date: Thu, 30 May 2024 19:40:06 -0400
Subject: [PATCH 21/32] remove timepoint arguement in write_config

---
 models/sipnet/R/write.configs.SIPNET.R            | 6 +++---
 models/sipnet/R/write_restart.SIPNET.R            | 3 +--
 models/sipnet/man/write.config.SIPNET.Rd          | 7 +++----
 models/sipnet/man/write_restart.SIPNET.Rd         | 7 +++----
 modules/assim.sequential/R/sda.enkf_MultiSite.R   | 2 +-
 modules/uncertainty/R/ensemble.R                  | 4 +---
 modules/uncertainty/R/get.parameter.samples.R     | 4 ++--
 modules/uncertainty/man/write.ensemble.configs.Rd | 7 +++----
 8 files changed, 17 insertions(+), 23 deletions(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 6cc1f8a959d..77666c736ae 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -24,7 +24,7 @@
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
-                                restart = NULL, spinup = NULL,obs_time=NULL,update_phenology=FALSE) {
+                                restart = NULL, spinup = NULL,update_phenology=FALSE) {
   ### WRITE sipnet.in
   template.in <- system.file("sipnet.in", package = "PEcAn.SIPNET")
   config.text <- readLines(con = template.in, n = -1)
@@ -454,8 +454,8 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
      if (!is.null(leaf_pheno_outdir)) {
     ##read data
        leafphdata <- utils::read.csv(paste0(leaf_pheno_outdir,"/leaf_phenology.csv"))
-       leafOnDay <- leafphdata$leafonday[leafphdata$year==obs_time & leafphdata$site_id==settings$run$site$id]
-       leafOffDay<- leafphdata$leafoffday[leafphdata$year==obs_time & leafphdata$site_id==settings$run$site$id]
+       leafOnDay <- leafphdata$leafonday[leafphdata$year==obs.year & leafphdata$site_id==settings$run$site$id]
+       leafOffDay<- leafphdata$leafoffday[leafphdata$year==obs.year & leafphdata$site_id==settings$run$site$id]
        if (!is.na(leafOnDay)){
 	      param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
        }
diff --git a/models/sipnet/R/write_restart.SIPNET.R b/models/sipnet/R/write_restart.SIPNET.R
index 38662852037..04aec87a8a0 100755
--- a/models/sipnet/R/write_restart.SIPNET.R
+++ b/models/sipnet/R/write_restart.SIPNET.R
@@ -29,7 +29,7 @@
 ##' @return NONE
 ##' @export
 write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings, new.state,
-                                 RENAME = TRUE, new.params = FALSE, inputs, obs_time=NULL, update_phenology=FALSE, verbose = FALSE) {
+                                 RENAME = TRUE, new.params = FALSE, inputs, update_phenology=FALSE, verbose = FALSE) {
 
   rundir <- settings$host$rundir
   variables <- colnames(new.state)
@@ -133,7 +133,6 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
                                            run.id = runid,
                                            inputs = inputs,
                                            IC = analysis.save.mat,
-                                           obs_time=obs_time,
                                            update_phenology=update_phenology))
   print(runid)
 } # write_restart.SIPNET
\ No newline at end of file
diff --git a/models/sipnet/man/write.config.SIPNET.Rd b/models/sipnet/man/write.config.SIPNET.Rd
index f06e358e8f0..fa366757088 100644
--- a/models/sipnet/man/write.config.SIPNET.Rd
+++ b/models/sipnet/man/write.config.SIPNET.Rd
@@ -13,8 +13,7 @@ write.config.SIPNET(
   IC = NULL,
   restart = NULL,
   spinup = NULL,
-  obs_time = NULL,
-  update_phenology = NULL
+  update_phenology = FALSE
 )
 }
 \arguments{
@@ -32,9 +31,9 @@ write.config.SIPNET(
 
 \item{restart}{In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.}
 
-\item{obs_time}{obervation timepoints}
-
 \item{update_phenology}{TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA}
+
+\item{obs_time}{obervation timepoints}
 }
 \description{
 Writes a configuration files for your model
diff --git a/models/sipnet/man/write_restart.SIPNET.Rd b/models/sipnet/man/write_restart.SIPNET.Rd
index 7896deefb1d..127ec027d96 100644
--- a/models/sipnet/man/write_restart.SIPNET.Rd
+++ b/models/sipnet/man/write_restart.SIPNET.Rd
@@ -14,8 +14,7 @@ write_restart.SIPNET(
   RENAME = TRUE,
   new.params = FALSE,
   inputs,
-  obs_time = NULL,
-  update_phenology = NULL,
+  update_phenology = FALSE,
   verbose = FALSE
 )
 }
@@ -38,11 +37,11 @@ write_restart.SIPNET(
 
 \item{inputs}{list of model inputs to use in write.configs.SIPNET}
 
-\item{obs_time}{obervation timepoints}
-
 \item{update_phenology}{TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA}
 
 \item{verbose}{decide if we want to print the outputs.}
+
+\item{obs_time}{obervation timepoints}
 }
 \value{
 NONE
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index c82a308df3f..9b4ce7b6845 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -403,6 +403,7 @@ sda.enkf.multisite <- function(settings,
       } else { ## t == 1
         restart.list <- vector("list", length(conf.settings))
       }
+      save(restart.list,file = file.path(settings$outdir, "restart_list.Rdata"))
       #add flag for restart t=1 to skip model runs
       if(restart_flag & t == 1){
         #for restart when t=1 do not need to do model runs and X should already exist in environment by this point
@@ -423,7 +424,6 @@ sda.enkf.multisite <- function(settings,
               write.to.db = settings$database$bety$write,
               restart = restart.arg,
               rename = TRUE,
-              time = obs.year,
               update_phenology=control$update_phenology
             )
           }) %>%
diff --git a/modules/uncertainty/R/ensemble.R b/modules/uncertainty/R/ensemble.R
index 1f8b836bc71..334d9373dd7 100644
--- a/modules/uncertainty/R/ensemble.R
+++ b/modules/uncertainty/R/ensemble.R
@@ -213,7 +213,7 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 ##' @export
 ##' @author David LeBauer, Carl Davidson, Hamze Dokoohaki
 write.ensemble.configs <- function(defaults, ensemble.samples, settings, model, 
-                                   clean = FALSE, write.to.db = TRUE, restart = NULL, rename = FALSE, time = NULL, update_phenology = FALSE) {
+                                   clean = FALSE, write.to.db = TRUE, restart = NULL, rename = FALSE, update_phenology = FALSE) {
   
   con <- NULL
   my.write.config <- paste("write.config.", model, sep = "")
@@ -412,7 +412,6 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
                                             trait.values = lapply(samples$parameters$samples, function(x, n) { x[n, , drop=FALSE] }, n=i), # this is the params
                                             settings = settings,
                                             run.id = run.id,
-                                            obs_time = time,
                                             update_phenology=update_phenology
       )
       )
@@ -464,7 +463,6 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
                            new.params = new.params[[i]], #new.params$`646`[[i]] for debugging
                            inputs =list(met=list(path=inputs$samples[[i]])), 
                            RENAME = rename,
-                           obs_time = time,
                            update_phenology = update_phenology)#for restart from previous model runs, not sharing the same outdir
       )
     }
diff --git a/modules/uncertainty/R/get.parameter.samples.R b/modules/uncertainty/R/get.parameter.samples.R
index d61cd951f09..e4821652ff2 100644
--- a/modules/uncertainty/R/get.parameter.samples.R
+++ b/modules/uncertainty/R/get.parameter.samples.R
@@ -58,10 +58,10 @@ get.parameter.samples <- function(settings,
     ## Load posteriors
     if (!is.na(posterior.files[i])) {
       # Load specified file
-      load(posterior.files[i], envir = distns)
+      load(posterior.files[i])
       fname <- file.path(outdirs[i], "post.distns.Rdata")
       if (file.exists(fname)) {
-        load(fname, envir = distns)
+        load(fname)
         prior.distns <- post.distns
       }
       if (exists("trait.mcmc") & exists("post.distns")) {
diff --git a/modules/uncertainty/man/write.ensemble.configs.Rd b/modules/uncertainty/man/write.ensemble.configs.Rd
index 4986fc99d48..ec4ca464e09 100644
--- a/modules/uncertainty/man/write.ensemble.configs.Rd
+++ b/modules/uncertainty/man/write.ensemble.configs.Rd
@@ -13,8 +13,7 @@ write.ensemble.configs(
   write.to.db = TRUE,
   restart = NULL,
   rename = FALSE,
-  time = NULL,
-  update_phenology = NULL
+  update_phenology = FALSE
 )
 }
 \arguments{
@@ -34,9 +33,9 @@ write.ensemble.configs(
 
 \item{rename}{Decide if we want to rename previous output files, for example convert from sipnet.out to sipnet.2020-07-16.out.}
 
-\item{time}{obervation timepoints}
-
 \item{update_phenology}{TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA}
+
+\item{time}{obervation timepoints}
 }
 \value{
 list, containing $runs = data frame of runids, $ensemble.id = the ensemble ID for these runs and $samples with ids and samples used for each tag.  Also writes sensitivity analysis configuration files as a side effect

From 68e44045e078e16e06e84f84e4ca9ea31e15a531 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-wo3.scc.bu.edu>
Date: Mon, 3 Jun 2024 00:34:41 -0400
Subject: [PATCH 22/32] Update per Mike's comments

---
 base/remote/R/qsub_parallel.R                 |  2 +-
 base/remote/man/qsub_parallel.Rd              |  2 +-
 base/workflow/R/do_conversions.R              |  7 +++
 .../03_topical_pages/04_R_workflow.Rmd        |  8 +++
 models/sipnet/R/model2netcdf.SIPNET.R         |  7 +--
 models/sipnet/R/write.configs.SIPNET.R        | 30 ++++++----
 models/sipnet/R/write_restart.SIPNET.R        |  7 +--
 models/sipnet/man/model2netcdf.SIPNET.Rd      |  5 +-
 models/sipnet/man/write.config.SIPNET.Rd      |  7 +--
 models/sipnet/man/write_restart.SIPNET.Rd     |  5 --
 .../assim.sequential/R/sda.enkf_MultiSite.R   | 10 +---
 .../MultiSite-Exs/SDA/Create_Multi_settings.R |  8 ++-
 .../man/sda.enkf.multisite.Rd                 |  3 +-
 modules/data.land/R/prepare_pools.R           | 11 +---
 modules/data.land/man/prepare_pools.Rd        |  4 +-
 modules/data.remote/NAMESPACE                 |  2 +-
 ...IS_extract.R => extract_phenology_MODIS.R} | 59 ++++++-------------
 .../man/extract_phenology_MODIS.Rd            | 42 +++++++++++++
 .../man/phenology_MODIS_extract.Rd            | 38 ------------
 modules/uncertainty/R/ensemble.R              | 10 +---
 modules/uncertainty/R/get.parameter.samples.R | 28 ++++-----
 .../uncertainty/man/write.ensemble.configs.Rd |  7 +--
 22 files changed, 131 insertions(+), 171 deletions(-)
 rename modules/data.remote/R/{phenology_MODIS_extract.R => extract_phenology_MODIS.R} (75%)
 create mode 100644 modules/data.remote/man/extract_phenology_MODIS.Rd
 delete mode 100644 modules/data.remote/man/phenology_MODIS_extract.Rd

diff --git a/base/remote/R/qsub_parallel.R b/base/remote/R/qsub_parallel.R
index 064672f11a1..2c855588c84 100644
--- a/base/remote/R/qsub_parallel.R
+++ b/base/remote/R/qsub_parallel.R
@@ -14,7 +14,7 @@
 #' 
 #' @importFrom foreach %dopar%
 #' @importFrom dplyr %>%
-qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out", sleep = 10, hybrid = TRUE) {
+qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out", sleep = 10, hybrid = FALSE) {
   if("try-error" %in% class(try(find.package("doSNOW"), silent = T))){
     PEcAn.logger::logger.info("Package doSNOW is not installed! Please install it and rerun the function!")
     return(0)
diff --git a/base/remote/man/qsub_parallel.Rd b/base/remote/man/qsub_parallel.Rd
index 274104b8139..3e5176379ef 100644
--- a/base/remote/man/qsub_parallel.Rd
+++ b/base/remote/man/qsub_parallel.Rd
@@ -9,7 +9,7 @@ qsub_parallel(
   files = NULL,
   prefix = "sipnet.out",
   sleep = 10,
-  hybrid = TRUE
+  hybrid = FALSE
 )
 }
 \arguments{
diff --git a/base/workflow/R/do_conversions.R b/base/workflow/R/do_conversions.R
index 2acf93a4b06..d6fef8d771d 100644
--- a/base/workflow/R/do_conversions.R
+++ b/base/workflow/R/do_conversions.R
@@ -73,6 +73,13 @@ do_conversions <- function(settings, overwrite.met = FALSE, overwrite.fia = FALS
       ## which is done locally in rundir and then rsync'ed to remote
       ## rather than having a model-format soils file that is processed remotely
     }
+
+    # Phenology data extraction
+    if(input.tag == "leaf_phenology" && is.null(input$path)){
+      settings$run$inputs[[i]]$path <- PEcAn.data.remote::extract_phenology_MODIS(site_info,start_date,end_date,outdir,run_parallel = TRUE,ncores = NULL)
+      needsave <- TRUE
+    }
+
     # met conversion
     
     if (input.tag == "met") {
diff --git a/book_source/03_topical_pages/04_R_workflow.Rmd b/book_source/03_topical_pages/04_R_workflow.Rmd
index dea190f6b6d..2c5c0706be8 100644
--- a/book_source/03_topical_pages/04_R_workflow.Rmd
+++ b/book_source/03_topical_pages/04_R_workflow.Rmd
@@ -121,6 +121,13 @@ The main script that handles Met Processing, is [`met.process`](https://github.c
 
 ### Converting from PEcAn standard to model-specific format {#workflow-met-model}
 
+## Input phenological Data {#workflow-input-phenology}
+
+To enable the use of MODIS phenology data (MODIS Land Cover Dynamics (MCD12Q2)) to update the phenological parameters (leaf-on date) and (leaf-off date) at each restart timepoint: 
+ 1. Generate the phenological parameter CSV file by running `PEcAn.data.remote::extract_phenology_MODIS`.
+ 2. Turn the tag under `settings$model$leaf_phenology` in xml file to be `TRUE`.
+ 3. Provide the generated phenological parameter CSV file to `settings$run$inputs$leaf_phenology$path`.
+
 ## Traits {#workflow-traits}
 
 (TODO: Under construction)
@@ -130,6 +137,7 @@ The main script that handles Met Processing, is [`met.process`](https://github.c
 (TODO: Under construction)
 
 ## Model Configuration {#workflow-modelconfig}
+To enable the state data assimilation with sub-annual data, default `conflict` in `model2netcdf` should be `TRUE`. 
 
 (TODO: Under construction)
 
diff --git a/models/sipnet/R/model2netcdf.SIPNET.R b/models/sipnet/R/model2netcdf.SIPNET.R
index da755a61bd2..ca01fff0814 100644
--- a/models/sipnet/R/model2netcdf.SIPNET.R
+++ b/models/sipnet/R/model2netcdf.SIPNET.R
@@ -97,14 +97,13 @@ sipnet2datetime <- function(sipnet_tval, base_year, base_month = 1,
 ##' @param revision model revision
 ##' @param overwrite Flag for overwriting nc files or not
 ##' @param conflict Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next.
-##' conflict is set to TRUE to enable the assimilation of monthly data.
 ##' @param prefix prefix to read the output files
 ##' @param delete.raw Flag to remove sipnet.out files, FALSE = do not remove files TRUE = remove files
 ##'
 ##' @export
 ##' @author Shawn Serbin, Michael Dietze
 model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, delete.raw = FALSE, revision, prefix = "sipnet.out",
-                                overwrite = FALSE, conflict = TRUE) {
+                                overwrite = FALSE, conflict = FALSE) {
 
   ### Read in model output in SIPNET format
   sipnet_out_file <- file.path(outdir, prefix)
@@ -141,7 +140,7 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
   for (y in year_seq) {
     #initialize the conflicted as FALSE
     conflicted <- FALSE
-    
+    conflict <- TRUE    #conflict is set to TRUE to enable the rename of yearly nc file for merging SDA results with sub-annual data
     #if we have conflicts on this file.
     if (file.exists(file.path(outdir, paste(y, "nc", sep = "."))) & overwrite == FALSE & conflict == FALSE) {
       next
@@ -297,7 +296,7 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
       close(varfile)
       ncdf4::nc_close(nc)
       
-      #merge nc files of the same year together to enable the assimilation of subannual data
+      #merge nc files of the same year together to enable the assimilation of sub-annual data
       if(file.exists(file.path(outdir, "previous.nc"))){
         files <- c(file.path(outdir, "previous.nc"), file.path(outdir, "current.nc"))
       }else{
diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index f206c715c9e..702991ca076 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -19,12 +19,10 @@
 ##' @param IC initial condition
 ##' @param restart In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.
 ##' @param spinup
-##' @param obs_time obervation timepoints
-##' @param update_phenology TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
-                                restart = NULL, spinup = NULL,update_phenology=FALSE) {
+                                restart = NULL, spinup = NULL) {
   ### WRITE sipnet.in
   template.in <- system.file("sipnet.in", package = "PEcAn.SIPNET")
   config.text <- readLines(con = template.in, n = -1)
@@ -450,13 +448,14 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     }
 
     #update LeafOnday and LeafOffDay
-    if (update_phenology){
-     leaf_pheno_outdir <- settings$model$leaf_phenology$outdir  ## read from settings
-     if (!is.null(leaf_pheno_outdir)) {
+    if (settings$model$leaf_phenology){
+     obs_year <- lubridate::year(settings$run$start.date)
+     leaf_pheno_data <- settings$run$inputs$leaf_phenology$path  ## read from settings
+     if (!is.null(leaf_pheno_data)) {
     ##read data
-       leafphdata <- utils::read.csv(paste0(leaf_pheno_outdir,"/leaf_phenology.csv"))
-       leafOnDay <- leafphdata$leafonday[leafphdata$year==obs.year & leafphdata$site_id==settings$run$site$id]
-       leafOffDay<- leafphdata$leafoffday[leafphdata$year==obs.year & leafphdata$site_id==settings$run$site$id]
+       leafphdata <- utils::read.csv(leaf_pheno_data)
+       leafOnDay <- leafphdata$leafonday[leafphdata$year == obs_year & leafphdata$site_id==settings$run$site$id]
+       leafOffDay<- leafphdata$leafoffday[leafphdata$year== obs_year & leafphdata$site_id==settings$run$site$id]
        if (!is.na(leafOnDay)){
 	      param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
        }
@@ -464,7 +463,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
         param[which(param[, 1] == "leafOffDay"), 2] <- leafOffDay
        }
       }else {
-      PEcAn.logger::logger.info("No phenology data were found. Please consider running modules/data.remote/R/phenology_MODIS_extract.R to get the parameter file.")
+      PEcAn.logger::logger.info("No phenology data were found. Please consider running `PEcAn.data.remote::extract_phenology_MODIS` to get the parameter file.")
       }
     }
   } ## end loop over PFTS
@@ -552,7 +551,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     ICs_num <- length(settings$run$inputs$poolinitcond$path)
     IC.path <- settings$run$inputs$poolinitcond$path[[sample(1:ICs_num, 1)]]
 
-    IC.pools <- PEcAn.data.land::prepare_pools(IC.path, settings, constants = list(sla = SLA))
+    IC.pools <- PEcAn.data.land::prepare_pools(IC.path, constants = list(sla = SLA))
     
     if(!is.null(IC.pools)){
       IC.nc <- ncdf4::nc_open(IC.path) #for additional variables specific to SIPNET
@@ -565,6 +564,15 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       if (!is.na(lai) && is.numeric(lai)) {
         param[which(param[, 1] == "laiInit"), 2] <- lai
       }
+
+      #Initial LAI is set as 0 for deciduous forests and grasslands for non-growing seasons
+      if (!(lubridate::month(settings$run$start.date) %in% seq(5,9))){ #Growing seasons are coarsely defined as months from May to September for non-conifers in the US
+         site_pft <- utils::read.csv(settings$run$inputs$pft.site$path)
+         site.pft.name <- site_pft$pft[site_pft$site == settings$run$site$id]
+         if (site.pft.name!="boreal.coniferous") {   #Currently only excluding boreal conifers. Other evergreen PFTs could be added here later.
+              param[which(param[, 1] == "laiInit"), 2] <- 0       
+          }
+      }
       ## neeInit gC/m2
       nee <- try(ncdf4::ncvar_get(IC.nc,"nee"),silent = TRUE)
       if (!is.na(nee) && is.numeric(nee)) {
diff --git a/models/sipnet/R/write_restart.SIPNET.R b/models/sipnet/R/write_restart.SIPNET.R
index 04aec87a8a0..455266f971e 100755
--- a/models/sipnet/R/write_restart.SIPNET.R
+++ b/models/sipnet/R/write_restart.SIPNET.R
@@ -20,8 +20,6 @@
 ##' @param RENAME flag to either rename output file or not
 ##' @param new.params list of parameters to convert between different states 
 ##' @param inputs list of model inputs to use in write.configs.SIPNET
-##' @param obs_time obervation timepoints
-##' @param update_phenology TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA
 ##' @param verbose decide if we want to print the outputs.
 ##'
 ##' @description Write restart files for SIPNET. WARNING: Some variables produce illegal values < 0 and have been hardcoded to correct these values!!
@@ -29,7 +27,7 @@
 ##' @return NONE
 ##' @export
 write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings, new.state,
-                                 RENAME = TRUE, new.params = FALSE, inputs, update_phenology=FALSE, verbose = FALSE) {
+                                 RENAME = TRUE, new.params = FALSE, inputs, verbose = FALSE) {
 
   rundir <- settings$host$rundir
   variables <- colnames(new.state)
@@ -132,7 +130,6 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
                                            settings = settings,
                                            run.id = runid,
                                            inputs = inputs,
-                                           IC = analysis.save.mat,
-                                           update_phenology=update_phenology))
+                                           IC = analysis.save.mat))
   print(runid)
 } # write_restart.SIPNET
\ No newline at end of file
diff --git a/models/sipnet/man/model2netcdf.SIPNET.Rd b/models/sipnet/man/model2netcdf.SIPNET.Rd
index 67a2b89466a..3853c307e6a 100644
--- a/models/sipnet/man/model2netcdf.SIPNET.Rd
+++ b/models/sipnet/man/model2netcdf.SIPNET.Rd
@@ -14,7 +14,7 @@ model2netcdf.SIPNET(
   revision,
   prefix = "sipnet.out",
   overwrite = FALSE,
-  conflict = TRUE
+  conflict = FALSE
 )
 }
 \arguments{
@@ -36,8 +36,7 @@ model2netcdf.SIPNET(
 
 \item{overwrite}{Flag for overwriting nc files or not}
 
-\item{conflict}{Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next.
-conflict is set to TRUE to enable the assimilation of monthly data.}
+\item{conflict}{Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next.}
 }
 \description{
 Convert SIPNET output to netCDF
diff --git a/models/sipnet/man/write.config.SIPNET.Rd b/models/sipnet/man/write.config.SIPNET.Rd
index fa366757088..8b829282c8e 100644
--- a/models/sipnet/man/write.config.SIPNET.Rd
+++ b/models/sipnet/man/write.config.SIPNET.Rd
@@ -12,8 +12,7 @@ write.config.SIPNET(
   inputs = NULL,
   IC = NULL,
   restart = NULL,
-  spinup = NULL,
-  update_phenology = FALSE
+  spinup = NULL
 )
 }
 \arguments{
@@ -30,10 +29,6 @@ write.config.SIPNET(
 \item{IC}{initial condition}
 
 \item{restart}{In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.}
-
-\item{update_phenology}{TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA}
-
-\item{obs_time}{obervation timepoints}
 }
 \description{
 Writes a configuration files for your model
diff --git a/models/sipnet/man/write_restart.SIPNET.Rd b/models/sipnet/man/write_restart.SIPNET.Rd
index 127ec027d96..22d7f885d7b 100644
--- a/models/sipnet/man/write_restart.SIPNET.Rd
+++ b/models/sipnet/man/write_restart.SIPNET.Rd
@@ -14,7 +14,6 @@ write_restart.SIPNET(
   RENAME = TRUE,
   new.params = FALSE,
   inputs,
-  update_phenology = FALSE,
   verbose = FALSE
 )
 }
@@ -37,11 +36,7 @@ write_restart.SIPNET(
 
 \item{inputs}{list of model inputs to use in write.configs.SIPNET}
 
-\item{update_phenology}{TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA}
-
 \item{verbose}{decide if we want to print the outputs.}
-
-\item{obs_time}{obervation timepoints}
 }
 \value{
 NONE
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index bc371c92765..99811e9787b 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -50,8 +50,7 @@ sda.enkf.multisite <- function(settings,
                                             keepNC = TRUE,
                                             forceRun = TRUE,
                                             run_parallel = TRUE,
-                                            MCMC.args = NULL,
-                                            update_phenology = TRUE),
+                                            MCMC.args = NULL),
                                ...) {
   #add if/else for when restart points to folder instead if T/F set restart as T
   if(is.list(restart)){
@@ -297,8 +296,7 @@ sda.enkf.multisite <- function(settings,
             settings = settings,
             model = settings$model$type,
             write.to.db = settings$database$bety$write,
-            restart = restart.arg,
-            update_phenology=control$update_phenology
+            restart = restart.arg
           )
         }) %>%
         stats::setNames(site.ids)
@@ -403,7 +401,6 @@ sda.enkf.multisite <- function(settings,
       } else { ## t == 1
         restart.list <- vector("list", length(conf.settings))
       }
-      save(restart.list,file = file.path(settings$outdir, "restart_list.Rdata"))
       #add flag for restart t=1 to skip model runs
       if(restart_flag & t == 1){
         #for restart when t=1 do not need to do model runs and X should already exist in environment by this point
@@ -423,8 +420,7 @@ sda.enkf.multisite <- function(settings,
               model = settings$model$type,
               write.to.db = settings$database$bety$write,
               restart = restart.arg,
-              rename = TRUE,
-              update_phenology=control$update_phenology
+              rename = TRUE
             )
           }) %>%
           stats::setNames(site.ids)
diff --git a/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R b/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R
index a788edb188f..5c92803dc20 100644
--- a/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R
+++ b/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R
@@ -17,6 +17,7 @@ ERA5_dir <- "/projectnb/dietzelab/dongchen/anchorSites/ERA5_2012_2021/"
 XML_out_dir <- "/projectnb/dietzelab/dongchen/anchorSites/SDA/pecan.xml"
 
 pft_csv_dir <- "/projectnb/dietzelab/dongchen/anchorSites/site_pft.csv"
+modis_phenology_dir <- "/projectnb/dietzelab/Cherry/pft_files/leaf_phenology.csv"
 
 #Obs_prep part
 #AGB
@@ -239,9 +240,9 @@ template <- PEcAn.settings::Settings(list(
                          type = "SIPNET",
                          revision = "ssr",
                          delete.raw = FALSE,
-                         binary = "/usr2/postdoc/istfer/SIPNET/trunk//sipnet_if",
+                         binary = model_binary,
                          jobtemplate = "~/sipnet_geo.job",
-                         leaf_phenology=structure(list(outdir="USER_DEFINED"))
+                         leaf_phenology= TRUE
   )),
   
   ###########################################################################
@@ -291,7 +292,8 @@ template <- PEcAn.settings::Settings(list(
       #                               )),
       # soilinitcond = structure(list(path = "/projectnb/dietzelab/ahelgeso/EFI_Forecast_Challenge/"
       #                               )),
-      pft.site = structure(list(path = pft_csv_dir))
+      pft.site = structure(list(path = pft_csv_dir)),
+      leaf_phenology = structure(list(path = modis_phenology_dir))
     ))
   ))
 ))
diff --git a/modules/assim.sequential/man/sda.enkf.multisite.Rd b/modules/assim.sequential/man/sda.enkf.multisite.Rd
index 96ade951efc..9c3560b3e89 100644
--- a/modules/assim.sequential/man/sda.enkf.multisite.Rd
+++ b/modules/assim.sequential/man/sda.enkf.multisite.Rd
@@ -14,8 +14,7 @@ sda.enkf.multisite(
   ensemble.samples = NULL,
   control = list(trace = TRUE, TimeseriesPlot = FALSE, debug = FALSE, pause = FALSE,
     Profiling = FALSE, OutlierDetection = FALSE, parallel_qsub = TRUE, send_email = NULL,
-    keepNC = TRUE, forceRun = TRUE, run_parallel = TRUE, MCMC.args = NULL,
-    update_phenology = TRUE),
+    keepNC = TRUE, forceRun = TRUE, run_parallel = TRUE, MCMC.args = NULL),
   ...
 )
 }
diff --git a/modules/data.land/R/prepare_pools.R b/modules/data.land/R/prepare_pools.R
index 5d82bdd1db2..ecb8c2f3b44 100644
--- a/modules/data.land/R/prepare_pools.R
+++ b/modules/data.land/R/prepare_pools.R
@@ -4,11 +4,10 @@
 ##' @export
 ##'
 ##' @param nc.path path to netcdf file containing standard dimensions and variables; currently supports these variables: TotLivBiom, leaf_carbon_content, LAI, AbvGrndWood, root_carbon_content, fine_root_carbon_content, coarse_root_carbon_content, litter_carbon_content, soil_organic_carbon_content, soil_carbon_content, wood_debris_carbon_content
-##' @param settings PEcAn settings object
 ##' @param constants list of constants; must include SLA in m2 / kg C if providing LAI for leaf carbon
 ##' @return list of pool values in kg C / m2 with generic names
 ##' @author Anne Thomas
-prepare_pools <- function(nc.path, settings, constants = NULL){
+prepare_pools <- function(nc.path, constants = NULL){
   #function to check that var was loaded (numeric) and has a valid value (not NA or negative)
   is.valid <- function(var){
     return(all(is.numeric(var) && !is.na(var) &&  var >= 0)) 
@@ -78,14 +77,6 @@ prepare_pools <- function(nc.path, settings, constants = NULL){
           PEcAn.logger::logger.error("TotLivBiom is less than sum of AbvGrndWood and roots; will use default for leaf biomass")
         }
       }
-      #Initial LAI at the beginning of year is set as 0 for deciduous forests and grasslands
-      site_pft <- utils::read.csv(settings$run$inputs$pft.site$path)
-      site.pft.name <- site_pft$pft[site_pft$site == settings$run$site$id]
-      if (site.pft.name!="boreal.coniferous") { 
-        leaf <- 0
-        IC.params[["leaf"]] <- leaf        
-      }
-      
       #Calculate LAI given leaf and sla
       sla <- constants$sla
       if (!is.null(sla) && is.valid(leaf)) {
diff --git a/modules/data.land/man/prepare_pools.Rd b/modules/data.land/man/prepare_pools.Rd
index d2eba2f6829..3fbdf7c0ac1 100644
--- a/modules/data.land/man/prepare_pools.Rd
+++ b/modules/data.land/man/prepare_pools.Rd
@@ -4,13 +4,11 @@
 \alias{prepare_pools}
 \title{prepare_pools}
 \usage{
-prepare_pools(nc.path, settings, constants = NULL)
+prepare_pools(nc.path, constants = NULL)
 }
 \arguments{
 \item{nc.path}{path to netcdf file containing standard dimensions and variables; currently supports these variables: TotLivBiom, leaf_carbon_content, LAI, AbvGrndWood, root_carbon_content, fine_root_carbon_content, coarse_root_carbon_content, litter_carbon_content, soil_organic_carbon_content, soil_carbon_content, wood_debris_carbon_content}
 
-\item{settings}{PEcAn settings object}
-
 \item{constants}{list of constants; must include SLA in m2 / kg C if providing LAI for leaf carbon}
 }
 \value{
diff --git a/modules/data.remote/NAMESPACE b/modules/data.remote/NAMESPACE
index 1ed0362a310..e54e6c5d778 100644
--- a/modules/data.remote/NAMESPACE
+++ b/modules/data.remote/NAMESPACE
@@ -13,7 +13,7 @@ export(download.NLCD)
 export(download.thredds.AGB)
 export(extract.LandTrendr.AGB)
 export(extract_NLCD)
-export(phenology_MODIS_extract)
+export(extract_phenology_MODIS)
 export(remote_process)
 importFrom(foreach,"%do%")
 importFrom(foreach,"%dopar%")
diff --git a/modules/data.remote/R/phenology_MODIS_extract.R b/modules/data.remote/R/extract_phenology_MODIS.R
similarity index 75%
rename from modules/data.remote/R/phenology_MODIS_extract.R
rename to modules/data.remote/R/extract_phenology_MODIS.R
index ab1af859057..3853103693e 100644
--- a/modules/data.remote/R/phenology_MODIS_extract.R
+++ b/modules/data.remote/R/extract_phenology_MODIS.R
@@ -1,48 +1,27 @@
 ##' Get MODIS phenology data by date and location  
 ##' 
 ##' @export
-##' @param settings PEcAn settings
+##' @param site_info A dataframe of site info containing the BETYdb site ID, 
+##' site name, latitude, and longitude, e.g. 
+##' @param start_date Start date to download data
+##' @param end_date End date to download data
+##' @param outdir Path to store the outputs
 ##' @param run_parallel optional method to download data parallely. Only works if more than 1 
 ##' site is needed and there are >1 CPUs available.
 ##' @param ncores number of cpus to use if run_parallel is set to TRUE. If you do not know the 
 ##' number of CPU's available, enter NULL. 
-##' 
-##' The output file will be saved to the folder defined in settings.
-##' 
-##' @examples
-##' \dontrun{
-##' settings <-  PEcAn.settings::read.settings("pecan.xml")
-##' test_modis_phenology <- PEcAn.data.remote::phenology_MODIS_extract(
-##'   settings = settings,
-##'   run_parallel = TRUE,
-##'   ncores = NULL)
-##' }
-##' @return a dataframe containing phenology dates and corresponding QA. 
+##' @return a dataframe containing phenology dates and corresponding QA.  
+##' The output file will be saved as a CSV file to the outdir.
 ##' Output column names are "year", "site_id", "lat", "lon", "leafonday","leafoffday","leafon_qa","leafoff_qa"
 ##' @author Qianyu Li
-##'
 
-phenology_MODIS_extract<- function(settings,run_parallel = TRUE,ncores = NULL){ 
-  #Get site id, lat, lon and site name for each site
-  site_info <- settings %>% 
-    purrr::map(~.x[['run']] ) %>% 
-    purrr::map('site')%>% 
-    purrr::map(function(site.list){
-      #conversion from string to number
-      site.list$lat <- as.numeric(site.list$lat)
-      site.list$lon <- as.numeric(site.list$lon)
-      list(site_id=site.list$id, lat=site.list$lat, lon=site.list$lon, site_name=site.list$name)
-      })%>% 
-    dplyr::bind_rows() %>% 
-    as.list()
-  #Set up the outdir to store leaf phenology data file           
-  outdir<-settings$model$leaf_phenology$outdir
+
+extract_phenology_MODIS<- function(site_info,start_date,end_date,outdir,run_parallel = TRUE,ncores = NULL){ 
+ 
   if (is.null(outdir)) {
-    PEcAn.logger::logger.error("No output directory found. Please set it up in pecan.xml under the model tag.")
+    PEcAn.logger::logger.error("No output directory found. Please provide it.")
   } else {
     #Set up the start and end date of the extraction
-    start_date <- as.Date(settings$state.data.assimilation$start.date, tz="UTC")
-    end_date <- as.Date(settings$state.data.assimilation$end.date, tz="UTC")
     start_YEARDOY <- paste0(lubridate::year(start_date),sprintf("%03d", 1))
     end_YEARDOY <- paste0(lubridate::year(end_date),sprintf("%03d", 365))
     #Extracting leaf-on, leaf-off, and Quality Assurance (QA) from MODIS MCD12Q2 product
@@ -70,7 +49,7 @@ phenology_MODIS_extract<- function(settings,run_parallel = TRUE,ncores = NULL){
                                                 QC_filter = FALSE, 
                                                 progress = TRUE)
 
-    qc_data <- PEcAn.data.remote::call_MODIS(outdir = NULL, 
+    qa_data <- PEcAn.data.remote::call_MODIS(outdir = NULL, 
                                            var = "QA", 
                                            site_info = site_info, 
                                            product_dates = c(start_YEARDOY, end_YEARDOY),
@@ -85,14 +64,14 @@ phenology_MODIS_extract<- function(settings,run_parallel = TRUE,ncores = NULL){
     #Example function for unpacking QA for each date phenometric from MCD12Q2_Collection6_UserGuide
     #in the order: Greenup, MidGreenup, Maturity, Peak, Senescence, MidGreendown, Dormancy                                        
     UnpackDetailedQA <- function(x){
-    bits <- as.integer(intToBits(x))
-    quals <- sapply(seq(1, 16, by=2), function(i) sum(bits[i:(i+1)] * 2^c(0, 1)))[1:7]
-    return(quals)
-    }
+       bits <- as.integer(intToBits(x))
+       quals <- sapply(seq(1, 16, by=2), function(i) sum(bits[i:(i+1)] * 2^c(0, 1)))[1:7]
+       return(quals)
+      }
     #Extract QA for leaf-on and leaf-off dates. Values are in the range 0–3 corresponding to “best”, “good”, “fair”, and “poor”.
-    qa_ind<-matrix(NA,nrow=length(qc_data[,1]),ncol=2)
-    for (i in seq_along(qc_data$data)){
-      qa_ind[i,]<-cbind(UnpackDetailedQA(qc_data$data[i])[2],UnpackDetailedQA(qc_data$data[i])[6])
+    qa_ind<-matrix(NA,nrow=length(qa_data[,1]),ncol=2)
+    for (i in seq_along(qa_data$data)){
+      qa_ind[i,]<-cbind(UnpackDetailedQA(qa_data$data[i])[2],UnpackDetailedQA(qa_data$data[i])[6])
     }
   
     leafphdata <- cbind(leafon_data %>% as.data.frame %>% dplyr::select("calendar_date", "site_id", "lat", "lon", "data"), leafoff_data$data,qa_ind)%>% 
diff --git a/modules/data.remote/man/extract_phenology_MODIS.Rd b/modules/data.remote/man/extract_phenology_MODIS.Rd
new file mode 100644
index 00000000000..cf73e693e7e
--- /dev/null
+++ b/modules/data.remote/man/extract_phenology_MODIS.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/extract_phenology_MODIS.R
+\name{extract_phenology_MODIS}
+\alias{extract_phenology_MODIS}
+\title{Get MODIS phenology data by date and location}
+\usage{
+extract_phenology_MODIS(
+  site_info,
+  start_date,
+  end_date,
+  outdir,
+  run_parallel = TRUE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{site_info}{A dataframe of site info containing the BETYdb site ID, 
+site name, latitude, and longitude, e.g.}
+
+\item{start_date}{Start date to download data}
+
+\item{end_date}{End date to download data}
+
+\item{outdir}{Path to store the outputs}
+
+\item{run_parallel}{optional method to download data parallely. Only works if more than 1 
+site is needed and there are >1 CPUs available.}
+
+\item{ncores}{number of cpus to use if run_parallel is set to TRUE. If you do not know the 
+number of CPU's available, enter NULL.}
+}
+\value{
+a dataframe containing phenology dates and corresponding QA.  
+The output file will be saved as a CSV file to the outdir.
+Output column names are "year", "site_id", "lat", "lon", "leafonday","leafoffday","leafon_qa","leafoff_qa"
+}
+\description{
+Get MODIS phenology data by date and location
+}
+\author{
+Qianyu Li
+}
diff --git a/modules/data.remote/man/phenology_MODIS_extract.Rd b/modules/data.remote/man/phenology_MODIS_extract.Rd
deleted file mode 100644
index 486f9d5d5cd..00000000000
--- a/modules/data.remote/man/phenology_MODIS_extract.Rd
+++ /dev/null
@@ -1,38 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/phenology_MODIS_extract.R
-\name{phenology_MODIS_extract}
-\alias{phenology_MODIS_extract}
-\title{Get MODIS phenology data by date and location}
-\usage{
-phenology_MODIS_extract(settings, run_parallel = TRUE, ncores = NULL)
-}
-\arguments{
-\item{settings}{PEcAn settings}
-
-\item{run_parallel}{optional method to download data parallely. Only works if more than 1 
-site is needed and there are >1 CPUs available.}
-
-\item{ncores}{number of cpus to use if run_parallel is set to TRUE. If you do not know the 
-number of CPU's available, enter NULL. 
-
-The output file will be saved to the folder defined in settings.}
-}
-\value{
-a dataframe containing phenology dates and corresponding QA. 
-Output column names are "year", "site_id", "lat", "lon", "leafonday","leafoffday","leafon_qa","leafoff_qa"
-}
-\description{
-Get MODIS phenology data by date and location
-}
-\examples{
-\dontrun{
-settings <-  PEcAn.settings::read.settings("pecan.xml")
-test_modis_phenology <- PEcAn.data.remote::phenology_MODIS_extract(
-  settings = settings,
-  run_parallel = TRUE,
-  ncores = NULL)
-}
-}
-\author{
-Qianyu Li
-}
diff --git a/modules/uncertainty/R/ensemble.R b/modules/uncertainty/R/ensemble.R
index 8cd0af38ede..a91b7385381 100644
--- a/modules/uncertainty/R/ensemble.R
+++ b/modules/uncertainty/R/ensemble.R
@@ -198,8 +198,6 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 ##' @param write.to.db logical: Record this run in BETY?
 ##' @param restart In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.
 ##' @param rename Decide if we want to rename previous output files, for example convert from sipnet.out to sipnet.2020-07-16.out.
-##' @param time obervation timepoints
-##' @param update_phenology TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA
 ##'
 ##' @return list, containing $runs = data frame of runids, $ensemble.id = the ensemble ID for these runs and $samples with ids and samples used for each tag.  Also writes sensitivity analysis configuration files as a side effect
 ##' @details The restart functionality is developed using model specific functions by calling write_restart.modelname function. First, you need to make sure that this function is already exist for your desired model.See here \url{https://pecanproject.github.io/pecan-documentation/master/pecan-models.html}
@@ -213,7 +211,7 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 ##' @export
 ##' @author David LeBauer, Carl Davidson, Hamze Dokoohaki
 write.ensemble.configs <- function(defaults, ensemble.samples, settings, model, 
-                                   clean = FALSE, write.to.db = TRUE, restart = NULL, rename = FALSE, update_phenology = FALSE) {
+                                   clean = FALSE, write.to.db = TRUE, restart = NULL, rename = FALSE) {
   
   con <- NULL
   my.write.config <- paste("write.config.", model, sep = "")
@@ -411,8 +409,7 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
       do.call(my.write.config, args = list( defaults = defaults, 
                                             trait.values = lapply(samples$parameters$samples, function(x, n) { x[n, , drop=FALSE] }, n=i), # this is the params
                                             settings = settings,
-                                            run.id = run.id,
-                                            update_phenology=update_phenology
+                                            run.id = run.id
       )
       )
       cat(format(run.id, scientific = FALSE), file = file.path(settings$rundir, "runs.txt"), sep = "\n", append = TRUE)
@@ -462,8 +459,7 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model,
                            new.state = new.state[i, ], 
                            new.params = new.params[[i]], #new.params$`646`[[i]] for debugging
                            inputs =list(met=list(path=inputs$samples[[i]])), 
-                           RENAME = rename,
-                           update_phenology = update_phenology)#for restart from previous model runs, not sharing the same outdir
+                           RENAME = rename)#for restart from previous model runs, not sharing the same outdir
       )
     }
     params<-new.params
diff --git a/modules/uncertainty/R/get.parameter.samples.R b/modules/uncertainty/R/get.parameter.samples.R
index e4821652ff2..0dd16d09461 100644
--- a/modules/uncertainty/R/get.parameter.samples.R
+++ b/modules/uncertainty/R/get.parameter.samples.R
@@ -58,14 +58,9 @@ get.parameter.samples <- function(settings,
     ## Load posteriors
     if (!is.na(posterior.files[i])) {
       # Load specified file
-      load(posterior.files[i])
-      fname <- file.path(outdirs[i], "post.distns.Rdata")
-      if (file.exists(fname)) {
-        load(fname)
-        prior.distns <- post.distns
-      }
-      if (exists("trait.mcmc") & exists("post.distns")) {
-        post.distns <- post.distns[!(row.names(post.distns) %in% names(trait.mcmc)),]
+      load(posterior.files[i], envir = distns)
+      if (is.null(distns$prior.distns) & !is.null(distns$post.distns)) {
+        distns$prior.distns <- distns$post.distns
       }
     } else {
       # Default to most recent posterior in the workflow, or the prior if there is none
@@ -111,16 +106,13 @@ get.parameter.samples <- function(settings,
     
     ### When no ma for a trait, sample from prior
     ### Trim all chains to shortest mcmc chain, else 20000 samples
-    if(exists("post.distns") & exists("trait.mcmc")){
-      priors <- c(names(trait.mcmc), rownames(post.distns))
-      prior.distns <- post.distns
-    }else{
-      priors <- rownames(prior.distns)
-    }
-    #priors <- names(trait.mcmc)
-  
-    if (exists("trait.mcmc")) {
-      param.names[[i]] <- names(trait.mcmc)
+    if(!is.null(distns$prior.distns)){
+      priors <- rownames(distns$prior.distns)
+    } else {
+      priors <- NULL
+    }  
+    if (!is.null(distns$trait.mcmc)) {
+      param.names[[i]] <- names(distns$trait.mcmc)
       names(param.names)[i] <- pft.name
       
       samples.num <- min(sapply(distns$trait.mcmc, function(x) nrow(as.matrix(x))))
diff --git a/modules/uncertainty/man/write.ensemble.configs.Rd b/modules/uncertainty/man/write.ensemble.configs.Rd
index ec4ca464e09..45d67cdcb7e 100644
--- a/modules/uncertainty/man/write.ensemble.configs.Rd
+++ b/modules/uncertainty/man/write.ensemble.configs.Rd
@@ -12,8 +12,7 @@ write.ensemble.configs(
   clean = FALSE,
   write.to.db = TRUE,
   restart = NULL,
-  rename = FALSE,
-  update_phenology = FALSE
+  rename = FALSE
 )
 }
 \arguments{
@@ -32,10 +31,6 @@ write.ensemble.configs(
 \item{restart}{In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.}
 
 \item{rename}{Decide if we want to rename previous output files, for example convert from sipnet.out to sipnet.2020-07-16.out.}
-
-\item{update_phenology}{TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA}
-
-\item{time}{obervation timepoints}
 }
 \value{
 list, containing $runs = data frame of runids, $ensemble.id = the ensemble ID for these runs and $samples with ids and samples used for each tag.  Also writes sensitivity analysis configuration files as a side effect

From 7b4839d6e015eb87ac0aedaf05f9ae6b05ac98f0 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-wo3.scc.bu.edu>
Date: Mon, 3 Jun 2024 00:52:17 -0400
Subject: [PATCH 23/32] Small tweak

---
 base/remote/R/qsub_parallel.R    | 2 +-
 base/remote/man/qsub_parallel.Rd | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/base/remote/R/qsub_parallel.R b/base/remote/R/qsub_parallel.R
index 2c855588c84..064672f11a1 100644
--- a/base/remote/R/qsub_parallel.R
+++ b/base/remote/R/qsub_parallel.R
@@ -14,7 +14,7 @@
 #' 
 #' @importFrom foreach %dopar%
 #' @importFrom dplyr %>%
-qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out", sleep = 10, hybrid = FALSE) {
+qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out", sleep = 10, hybrid = TRUE) {
   if("try-error" %in% class(try(find.package("doSNOW"), silent = T))){
     PEcAn.logger::logger.info("Package doSNOW is not installed! Please install it and rerun the function!")
     return(0)
diff --git a/base/remote/man/qsub_parallel.Rd b/base/remote/man/qsub_parallel.Rd
index 3e5176379ef..274104b8139 100644
--- a/base/remote/man/qsub_parallel.Rd
+++ b/base/remote/man/qsub_parallel.Rd
@@ -9,7 +9,7 @@ qsub_parallel(
   files = NULL,
   prefix = "sipnet.out",
   sleep = 10,
-  hybrid = FALSE
+  hybrid = TRUE
 )
 }
 \arguments{

From c3abe88b4486a6dcd5ae6bc69d7a51890c0753a7 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-un3.scc.bu.edu>
Date: Tue, 4 Jun 2024 13:50:47 -0400
Subject: [PATCH 24/32] Apply suggestions

---
 .../03_topical_pages/04_R_workflow.Rmd        |  3 +--
 models/sipnet/R/write.configs.SIPNET.R        | 20 +++++++++++--------
 .../assim.sequential/R/sda.enkf_MultiSite.R   |  2 +-
 .../MultiSite-Exs/SDA/Create_Multi_settings.R |  3 +--
 .../data.remote/R/extract_phenology_MODIS.R   |  9 +++++----
 .../man/extract_phenology_MODIS.Rd            |  2 +-
 6 files changed, 21 insertions(+), 18 deletions(-)

diff --git a/book_source/03_topical_pages/04_R_workflow.Rmd b/book_source/03_topical_pages/04_R_workflow.Rmd
index 2c5c0706be8..c1d8fed0145 100644
--- a/book_source/03_topical_pages/04_R_workflow.Rmd
+++ b/book_source/03_topical_pages/04_R_workflow.Rmd
@@ -125,8 +125,7 @@ The main script that handles Met Processing, is [`met.process`](https://github.c
 
 To enable the use of MODIS phenology data (MODIS Land Cover Dynamics (MCD12Q2)) to update the phenological parameters (leaf-on date) and (leaf-off date) at each restart timepoint: 
  1. Generate the phenological parameter CSV file by running `PEcAn.data.remote::extract_phenology_MODIS`.
- 2. Turn the tag under `settings$model$leaf_phenology` in xml file to be `TRUE`.
- 3. Provide the generated phenological parameter CSV file to `settings$run$inputs$leaf_phenology$path`.
+ 2. Provide the generated phenological parameter CSV file to `settings$run$inputs$leaf_phenology$path`.
 
 ## Traits {#workflow-traits}
 
diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 702991ca076..e1107f12b57 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -448,21 +448,25 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     }
 
     #update LeafOnday and LeafOffDay
-    if (settings$model$leaf_phenology){
-     obs_year <- lubridate::year(settings$run$start.date)
-     leaf_pheno_data <- settings$run$inputs$leaf_phenology$path  ## read from settings
-     if (!is.null(leaf_pheno_data)) {
+     if (!is.null(settings$run$inputs$leaf_phenology)){
+     obs_year_start <- lubridate::year(settings$run$start.date)
+     obs_year_end <- lubridate::year(settings$run$end.date)
+     if (obs_year_start != obs_year_end) {
+      PEcAn.logger::logger.info("Start.date and end.date are not in the same year. Currently start.date is used for refering phenological data")
+     }
+     leaf_pheno_path <- settings$run$inputs$leaf_phenology$path  ## read from settings
+      if (!is.null(leaf_pheno_path)){
     ##read data
-       leafphdata <- utils::read.csv(leaf_pheno_data)
-       leafOnDay <- leafphdata$leafonday[leafphdata$year == obs_year & leafphdata$site_id==settings$run$site$id]
-       leafOffDay<- leafphdata$leafoffday[leafphdata$year== obs_year & leafphdata$site_id==settings$run$site$id]
+       leafphdata <- utils::read.csv(leaf_pheno_path)
+       leafOnDay <-  dplyr::filter(leafphdata, year== obs_year_start & site_id==settings$run$site$id) %>%  dplyr::select (leafonday)
+       leafOffDay<-  dplyr::filter(leafphdata, year== obs_year_start & site_id==settings$run$site$id) %>%  dplyr::select (leafoffday)
        if (!is.na(leafOnDay)){
 	      param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
        }
        if (!is.na(leafOffDay)){
         param[which(param[, 1] == "leafOffDay"), 2] <- leafOffDay
        }
-      }else {
+      } else {
       PEcAn.logger::logger.info("No phenology data were found. Please consider running `PEcAn.data.remote::extract_phenology_MODIS` to get the parameter file.")
       }
     }
diff --git a/modules/assim.sequential/R/sda.enkf_MultiSite.R b/modules/assim.sequential/R/sda.enkf_MultiSite.R
index 99811e9787b..b46d954ab95 100644
--- a/modules/assim.sequential/R/sda.enkf_MultiSite.R
+++ b/modules/assim.sequential/R/sda.enkf_MultiSite.R
@@ -530,7 +530,7 @@ sda.enkf.multisite <- function(settings,
           if (is.null(control$MCMC.args)) {
             MCMC.args <- list(niter = 1e5,
                               nthin = 10,
-                              nchain = 1,
+                              nchain = 3,
                               nburnin = 5e4)
           } else {
             MCMC.args <- control$MCMC.args
diff --git a/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R b/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R
index 5c92803dc20..0202d56cede 100644
--- a/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R
+++ b/modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R
@@ -241,8 +241,7 @@ template <- PEcAn.settings::Settings(list(
                          revision = "ssr",
                          delete.raw = FALSE,
                          binary = model_binary,
-                         jobtemplate = "~/sipnet_geo.job",
-                         leaf_phenology= TRUE
+                         jobtemplate = "~/sipnet_geo.job"
   )),
   
   ###########################################################################
diff --git a/modules/data.remote/R/extract_phenology_MODIS.R b/modules/data.remote/R/extract_phenology_MODIS.R
index 3853103693e..d99620950b5 100644
--- a/modules/data.remote/R/extract_phenology_MODIS.R
+++ b/modules/data.remote/R/extract_phenology_MODIS.R
@@ -10,7 +10,7 @@
 ##' site is needed and there are >1 CPUs available.
 ##' @param ncores number of cpus to use if run_parallel is set to TRUE. If you do not know the 
 ##' number of CPU's available, enter NULL. 
-##' @return a dataframe containing phenology dates and corresponding QA.  
+##' @return the path for output file 
 ##' The output file will be saved as a CSV file to the outdir.
 ##' Output column names are "year", "site_id", "lat", "lon", "leafonday","leafoffday","leafon_qa","leafoff_qa"
 ##' @author Qianyu Li
@@ -86,8 +86,9 @@ extract_phenology_MODIS<- function(site_info,start_date,end_date,outdir,run_para
     leafphdata$year<-lubridate::year(leafphdata$year)
     leafphdata$site_id<-as.character(leafphdata$site_id)
     
-    PEcAn.logger::logger.info(paste0("Storing results in: ",file.path(outdir,"leaf_phenology.csv")))
-    utils::write.csv(leafphdata, file = file.path(outdir, "leaf_phenology.csv"), row.names = FALSE)
-    return(leafphdata)
+    file_path<-file.path(outdir,"leaf_phenology.csv")
+    PEcAn.logger::logger.info(paste0("Storing results in: ",file_path))
+    utils::write.csv(leafphdata, file = file_path, row.names = FALSE)
+    return(file_path)
   }
 }
\ No newline at end of file
diff --git a/modules/data.remote/man/extract_phenology_MODIS.Rd b/modules/data.remote/man/extract_phenology_MODIS.Rd
index cf73e693e7e..0a134d32069 100644
--- a/modules/data.remote/man/extract_phenology_MODIS.Rd
+++ b/modules/data.remote/man/extract_phenology_MODIS.Rd
@@ -30,7 +30,7 @@ site is needed and there are >1 CPUs available.}
 number of CPU's available, enter NULL.}
 }
 \value{
-a dataframe containing phenology dates and corresponding QA.  
+the path for output file 
 The output file will be saved as a CSV file to the outdir.
 Output column names are "year", "site_id", "lat", "lon", "leafonday","leafoffday","leafon_qa","leafoff_qa"
 }

From bb5d40c7bfd8465dab5ffc26bf578d031eeea336 Mon Sep 17 00:00:00 2001
From: Michael Dietze <dietze@bu.edu>
Date: Tue, 4 Jun 2024 14:07:00 -0400
Subject: [PATCH 25/32] Update base/workflow/R/do_conversions.R

---
 base/workflow/R/do_conversions.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/base/workflow/R/do_conversions.R b/base/workflow/R/do_conversions.R
index d6fef8d771d..2d886ba4cbf 100644
--- a/base/workflow/R/do_conversions.R
+++ b/base/workflow/R/do_conversions.R
@@ -76,7 +76,7 @@ do_conversions <- function(settings, overwrite.met = FALSE, overwrite.fia = FALS
 
     # Phenology data extraction
     if(input.tag == "leaf_phenology" && is.null(input$path)){
-      settings$run$inputs[[i]]$path <- PEcAn.data.remote::extract_phenology_MODIS(site_info,start_date,end_date,outdir,run_parallel = TRUE,ncores = NULL)
+      #settings$run$inputs[[i]]$path <- PEcAn.data.remote::extract_phenology_MODIS(site_info,start_date,end_date,outdir,run_parallel = TRUE,ncores = NULL)
       needsave <- TRUE
     }
 

From 5ce785cc990c7cf445074be7a0097990613adb6e Mon Sep 17 00:00:00 2001
From: Michael Dietze <dietze@bu.edu>
Date: Tue, 4 Jun 2024 14:07:12 -0400
Subject: [PATCH 26/32] Update models/sipnet/R/write.configs.SIPNET.R

---
 models/sipnet/R/write.configs.SIPNET.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index e1107f12b57..89ee492bf57 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -18,7 +18,7 @@
 ##' @param inputs list of model inputs
 ##' @param IC initial condition
 ##' @param restart In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.
-##' @param spinup
+##' @param spinup currently unused, included for compatibility with other models
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,

From 82a1713b157b72c51b2d0e0d2b9a53b42fb3f52c Mon Sep 17 00:00:00 2001
From: Michael Dietze <dietze@bu.edu>
Date: Tue, 4 Jun 2024 14:26:52 -0400
Subject: [PATCH 27/32] Update models/sipnet/man/write.config.SIPNET.Rd

---
 models/sipnet/man/write.config.SIPNET.Rd | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/models/sipnet/man/write.config.SIPNET.Rd b/models/sipnet/man/write.config.SIPNET.Rd
index 8b829282c8e..47b54d67af9 100644
--- a/models/sipnet/man/write.config.SIPNET.Rd
+++ b/models/sipnet/man/write.config.SIPNET.Rd
@@ -30,6 +30,8 @@ write.config.SIPNET(
 
 \item{restart}{In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.}
 }
+
+\item{spinup}{currently unused, included for compatibility with other models}
 \description{
 Writes a configuration files for your model
 }

From 47ef5c47db28091e0fb5757a4528b0c09ce76106 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-un3.scc.bu.edu>
Date: Tue, 4 Jun 2024 15:03:00 -0400
Subject: [PATCH 28/32] Update document

---
 models/sipnet/man/write.config.SIPNET.Rd | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/sipnet/man/write.config.SIPNET.Rd b/models/sipnet/man/write.config.SIPNET.Rd
index 47b54d67af9..32c3e334238 100644
--- a/models/sipnet/man/write.config.SIPNET.Rd
+++ b/models/sipnet/man/write.config.SIPNET.Rd
@@ -29,9 +29,9 @@ write.config.SIPNET(
 \item{IC}{initial condition}
 
 \item{restart}{In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.}
-}
 
 \item{spinup}{currently unused, included for compatibility with other models}
+}
 \description{
 Writes a configuration files for your model
 }

From 212cdfa810691ca75b5d6e5f7042d8526a1e8acd Mon Sep 17 00:00:00 2001
From: Michael Dietze <dietze@bu.edu>
Date: Tue, 4 Jun 2024 15:45:50 -0400
Subject: [PATCH 29/32] Apply suggestions from code review

---
 models/sipnet/R/write.configs.SIPNET.R | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 89ee492bf57..41f568c8cfa 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -19,6 +19,8 @@
 ##' @param IC initial condition
 ##' @param restart In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.
 ##' @param spinup currently unused, included for compatibility with other models
+##' @importFrom rlang .data
+##' @importFrom dplyr %>%
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
@@ -458,8 +460,8 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       if (!is.null(leaf_pheno_path)){
     ##read data
        leafphdata <- utils::read.csv(leaf_pheno_path)
-       leafOnDay <-  dplyr::filter(leafphdata, year== obs_year_start & site_id==settings$run$site$id) %>%  dplyr::select (leafonday)
-       leafOffDay<-  dplyr::filter(leafphdata, year== obs_year_start & site_id==settings$run$site$id) %>%  dplyr::select (leafoffday)
+       leafOnDay <-  dplyr::filter(leafphdata, .data$yea r== obs_year_start & .data$site_id==settings$run$site$id) %>%  dplyr::select (.data$leafonday)
+       leafOffDay<-  dplyr::filter(leafphdata, .data$year == obs_year_start & .data$site_id==settings$run$site$id) %>%  dplyr::select (.data$leafoffday)
        if (!is.na(leafOnDay)){
 	      param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
        }

From d6a4a3ddc4386aef5ca54c54dbcf5bbdf67d4479 Mon Sep 17 00:00:00 2001
From: Michael Dietze <dietze@bu.edu>
Date: Tue, 4 Jun 2024 16:23:13 -0400
Subject: [PATCH 30/32] Update models/sipnet/R/write.configs.SIPNET.R

Co-authored-by: Chris Black <chris@ckblack.org>
---
 models/sipnet/R/write.configs.SIPNET.R | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 41f568c8cfa..3bd113806cb 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -460,8 +460,10 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       if (!is.null(leaf_pheno_path)){
     ##read data
        leafphdata <- utils::read.csv(leaf_pheno_path)
-       leafOnDay <-  dplyr::filter(leafphdata, .data$yea r== obs_year_start & .data$site_id==settings$run$site$id) %>%  dplyr::select (.data$leafonday)
-       leafOffDay<-  dplyr::filter(leafphdata, .data$year == obs_year_start & .data$site_id==settings$run$site$id) %>%  dplyr::select (.data$leafoffday)
+       leafOnDay <-  dplyr::filter(leafphdata, .data$year == obs_year_start & .data$site_id==settings$run$site$id) %>%
+         dplyr::select (.data$leafonday)
+       leafOffDay<-  dplyr::filter(leafphdata, .data$year == obs_year_start & .data$site_id==settings$run$site$id) %>%
+         dplyr::select (.data$leafoffday)
        if (!is.na(leafOnDay)){
 	      param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
        }

From a6f38347bba109538903b5d6592cfe58e0cdb556 Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-un3.scc.bu.edu>
Date: Tue, 4 Jun 2024 17:55:12 -0400
Subject: [PATCH 31/32] Update document

---
 models/sipnet/NAMESPACE | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/models/sipnet/NAMESPACE b/models/sipnet/NAMESPACE
index cdef11b5367..e2cdb0e98d0 100644
--- a/models/sipnet/NAMESPACE
+++ b/models/sipnet/NAMESPACE
@@ -11,3 +11,5 @@ export(split_inputs.SIPNET)
 export(veg2model.SIPNET)
 export(write.config.SIPNET)
 export(write_restart.SIPNET)
+importFrom(dplyr,"%>%")
+importFrom(rlang,.data)

From 5280bfd4ed55940776dc12322bc1b6603cd76dae Mon Sep 17 00:00:00 2001
From: Qianyu Li <ql1@scc-un3.scc.bu.edu>
Date: Wed, 5 Jun 2024 17:34:43 -0400
Subject: [PATCH 32/32] Revert back to filtering

---
 models/sipnet/NAMESPACE                | 2 --
 models/sipnet/R/write.configs.SIPNET.R | 8 ++------
 2 files changed, 2 insertions(+), 8 deletions(-)

diff --git a/models/sipnet/NAMESPACE b/models/sipnet/NAMESPACE
index e2cdb0e98d0..cdef11b5367 100644
--- a/models/sipnet/NAMESPACE
+++ b/models/sipnet/NAMESPACE
@@ -11,5 +11,3 @@ export(split_inputs.SIPNET)
 export(veg2model.SIPNET)
 export(write.config.SIPNET)
 export(write_restart.SIPNET)
-importFrom(dplyr,"%>%")
-importFrom(rlang,.data)
diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 3bd113806cb..ffa65a5ee62 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -19,8 +19,6 @@
 ##' @param IC initial condition
 ##' @param restart In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details.
 ##' @param spinup currently unused, included for compatibility with other models
-##' @importFrom rlang .data
-##' @importFrom dplyr %>%
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
@@ -460,10 +458,8 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       if (!is.null(leaf_pheno_path)){
     ##read data
        leafphdata <- utils::read.csv(leaf_pheno_path)
-       leafOnDay <-  dplyr::filter(leafphdata, .data$year == obs_year_start & .data$site_id==settings$run$site$id) %>%
-         dplyr::select (.data$leafonday)
-       leafOffDay<-  dplyr::filter(leafphdata, .data$year == obs_year_start & .data$site_id==settings$run$site$id) %>%
-         dplyr::select (.data$leafoffday)
+       leafOnDay <-  leafphdata$leafonday[leafphdata$year == obs_year_start & leafphdata$site_id==settings$run$site$id]
+       leafOffDay<-  leafphdata$leafoffday[leafphdata$year== obs_year_start & leafphdata$site_id==settings$run$site$id]
        if (!is.na(leafOnDay)){
 	      param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
        }