diff --git a/models/lpjguess/NAMESPACE b/models/lpjguess/NAMESPACE
index 56bfda9e366..e814d6a2aec 100644
--- a/models/lpjguess/NAMESPACE
+++ b/models/lpjguess/NAMESPACE
@@ -6,10 +6,14 @@ export(met2model.LPJGUESS)
 export(model2netcdf.LPJGUESS)
 export(pecan2lpjguess)
 export(readStateBinary)
+export(read_binary_LPJGUESS)
+export(read_restart.LPJGUESS)
 export(split_inputs.LPJGUESS)
 export(update_state_LPJGUESS)
 export(write.config.LPJGUESS)
 export(write.insfile.LPJGUESS)
+export(write_binary_LPJGUESS)
+export(write_restart.LPJGUESS)
 importFrom(PEcAn.utils,days_in_year)
 importFrom(Rcpp,sourceCpp)
 importFrom(ncdf4,nc_close)
diff --git a/models/lpjguess/R/read_restart.LPJGUESS.R b/models/lpjguess/R/read_restart.LPJGUESS.R
index d828137c8f1..a8bd857f302 100644
--- a/models/lpjguess/R/read_restart.LPJGUESS.R
+++ b/models/lpjguess/R/read_restart.LPJGUESS.R
@@ -1,13 +1,25 @@
-
-# developing
-# outdir = "/fs/data2/output//PEcAn_1000010473/out"
-# runid = 1002656839
-# stop.time = "1960-12-31 23:59:59 UTC"
-# load("/fs/data2/output/PEcAn_1000010473/SDAsettings_develop.Rdata")
-# var.names = c("AGB.pft", "TotSoilCarb")
-# load("/fs/data2/output/PEcAn_1000010473/SDAparams_develop.Rdata")
-
-
+#' Read Restart for LPJGUESS
+#'
+#' @param outdir      output directory
+#' @param runid       run ID
+#' @param stop.time   year that is being read
+#' @param settings    PEcAn settings object
+#' @param var.names   var.names to be extracted
+#' @param params      passed on to return value
+#'
+#' @return X_tmp      vector of forecasts
+#' @export
+#' @examples
+#' \dontrun{
+#'   rx <- read_restart.LPJGUESS(
+#'            outdir   = "/projectnb/…/LPJ_output",
+#'            runid    = "123456",
+#'            stop.time = as.POSIXct("2001-12-31 23:59:59", tz = "UTC"),
+#'            settings = settings,
+#'            var.names = c("AGB.pft"),
+#'            params = params)
+#' }
+#' @author Istem Fer, Yinghao Sun
 read_restart.LPJGUESS <- function(outdir, runid, stop.time, settings, var.names, params){
   
   # which LPJ-GUESS version, the structure of state file depends a lot on version
@@ -24,7 +36,6 @@ read_restart.LPJGUESS <- function(outdir, runid, stop.time, settings, var.names,
   # read binary state file, takes a couple of minutes
   Gridcell_container <- read_binary_LPJGUESS(outdir  = file.path(outdir, runid), 
                                              version = lpjguess_ver)
-  
   forecast <- list()
   
   # additional varnames for LPJ-GUESS?
@@ -33,8 +44,8 @@ read_restart.LPJGUESS <- function(outdir, runid, stop.time, settings, var.names,
     
     if (var_name == "AGB.pft") {
       
-      cmass_sap_perpft   <- calculateGridcellVariablePerPFT(model.state = Gridcell_container, variable = "cmass_sap")
-      cmass_heart_perpft <- calculateGridcellVariablePerPFT(model.state = Gridcell_container, variable = "cmass_heart")
+      cmass_sap_perpft   <- calculateGridcellVariablePerPFT(model.state = Gridcell_container$state, variable = "cmass_sap")
+      cmass_heart_perpft <- calculateGridcellVariablePerPFT(model.state = Gridcell_container$state, variable = "cmass_heart")
       
       cmass_wood <- cmass_sap_perpft + cmass_heart_perpft
       cmass_wood <- PEcAn.utils::ud_convert(cmass_wood, "kg/m^2", "Mg/ha")
@@ -45,11 +56,12 @@ read_restart.LPJGUESS <- function(outdir, runid, stop.time, settings, var.names,
       cmass_abvg_wood <- cmass_wood - cmass_blwg_wood
       
       forecast[[length(forecast) + 1]]    <- cmass_abvg_wood
-      names(forecast[[length(forecast)]]) <- paste0("AGB.pft.", unlist(Gridcell_container$meta_data$pft))
+      names(forecast[[length(forecast)]]) <- paste0("AGB.pft.", unlist(Gridcell_container$state$meta_data$pft))
       
     }
   }
   
+  # params$LPJGUESS_state include state, pos_list, siz_list
   params$LPJGUESS_state <- Gridcell_container
   
   PEcAn.logger::logger.info("Finished --", runid)
diff --git a/models/lpjguess/R/read_state.R b/models/lpjguess/R/read_state.R
index c1f0ee149bc..7d54c4a8525 100644
--- a/models/lpjguess/R/read_state.R
+++ b/models/lpjguess/R/read_state.R
@@ -1,6 +1,3 @@
-
-######################## Helper functions ########################
-
 #' Find Stream Variable
 #'
 #' A helper function that lists streamed variables. It returns the names of streamed variables.
@@ -8,6 +5,7 @@
 #' @param file_in A character vector representing the file content to search through.
 #' @param line_nos A numeric vector of length 2, specifying the start and end lines to search for streamed variables.
 #' @return A character vector of streamed variable names.
+#' @keywords internal
 # helper function that lists streamed variables, it just returns the names, types are checked by other fucntion
 find_stream_var <- function(file_in, line_nos){
   
@@ -85,6 +83,7 @@ find_stream_var <- function(file_in, line_nos){
 #' @param pattern A character string pattern to look for in the file.
 #' @return A numeric vector of length 2, giving the start and end line numbers.
 #' @importFrom stringr str_match
+#' @keywords internal
 # helper function that scans LPJ-GUESS that returns the beginning and the ending lines of serialized object
 serialize_starts_ends <- function(file_in, pattern = "void Gridcell::serialize"){
   # find the starting line from the given pattern
@@ -159,6 +158,7 @@ find_closing <- function(find = "}", line_no, file_in, if_else_check = FALSE){
 #' @return A numeric value representing the size (number of streamed variables).
 #' @importFrom stringr str_match
 #' @importFrom utils glob2rx
+#' @keywords internal
 # helper function that determines the stream size to read
 find_stream_size <- function(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS){
   
@@ -322,6 +322,7 @@ read_state <- function(file_path) {
 #' @param LPJ_GUESS_TYPES A character vector of recognized LPJ-GUESS types.
 #' @param guessh_in A character vector of LPJ-GUESS header file content.
 #' @return A character string indicating the stream type.
+#' @keywords internal
 # helper function to decide the type of the stream
 # this function relies on the architecture of LPJ-GUESS and has bunch of harcoded checks, see model documentation
 find_stream_type <- function(class = NULL, current_stream_var, LPJ_GUESS_CLASSES, LPJ_GUESS_TYPES, guessh_in){
@@ -429,13 +430,8 @@ find_stream_type <- function(class = NULL, current_stream_var, LPJ_GUESS_CLASSES
   return(list(type = gsub(" ", "", stream_type), name = stream_name, substring = sub_string))
 } # find_stream_type
 
-
-###################################### READ STATE
-
-
 # this fcn is for potential natural vegetation only
 # when there is landcover, there will be more stand types
-
 # also for cohort mode only
 
 # Gridcell: Top-level object containing all dynamic and static data for a particular gridcell
@@ -450,26 +446,31 @@ find_stream_type <- function(class = NULL, current_stream_var, LPJ_GUESS_CLASSES
 # Soil : Stores state variables for soils and the snow pack. One object of class Soil is defined for each patch.
 # Fluxes : The Fluxes class stores accumulated monthly and annual fluxes. One object of type Fluxes is defined for each patch.
 # Individual : Stores state variables for an average individual plant. In cohort mode, it is the average individual of a cohort of plants approximately the same age and from the same patch.
-
-
-# test path
-#outdir <- "/fs/data2/output/PEcAn_1000010473/out/1002656304"
-
-# outdir, at least model version, maybe also settings
+# 
 #' Read Binary File for LPJ-GUESS
 #'
 #' Reads a binary file formatted for LPJ-GUESS and extracts relevant data.
 #'
-#' @param outdir A character string specifying the output directory containing the binary state files.
+#' @param outdir  The output directory where ".state" and "meta.bin" will be written
 #' @param version A character string specifying the LPJ-GUESS version (default is "PalEON").
 #' @importFrom stringr str_match
 #' @importFrom utils glob2rx
 #' @return A matrix or list containing the extracted data.
+#' @export
+#' @author Istem Fer, Yinghao Sun
 read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
   
+  # ## FOR TEST
+  # outdir <- "/projectnb/dietzelab/yinghao/try/write_test/out"
+  # rundir <- "/projectnb/dietzelab/yinghao/try/write_test/run"
+  
   # find rundir too, params.ins is in there and we need to get some values from there
   rundir <- file.path(dirname(dirname(outdir)), "run", basename(outdir))
   
+  # create lists to store byte offset and byte size for each variable
+  pos_list <- list()
+  siz_list <- list()
+  
   # guess.cpp has the info of what is being written
   guesscpp_name <- paste0("guess.", version, ".cpp")  # these are gonna be in the package guess.VERSION.cpp
   guesscpp_in   <- readLines(con = system.file(guesscpp_name, package = "PEcAn.LPJGUESS"), n = -1)
@@ -629,6 +630,10 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
   Gridcell <- list()
   level <- "Gridcell"
   for(g_i in seq_along(streamed_vars_gridcell)){ # Gridcell-loop starts
+    
+    # # Debug for empty nstands
+    # if(g_i == 7) browser()
+    
     current_stream <- streamed_vars_gridcell[g_i]
     # weird, it doesn't go into Gridcell st
     if(current_stream == "st[i]")   next #current_stream <- "Gridcellst" 
@@ -643,7 +648,12 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
       # note that this is streamed under Gridcell, not Stand in guess.cpp, 
       # but I think this info needs to go together with the Stand sublist
       # so prepend landcovertype to the streamed_vars_stand EDIT: I'll actually just read it here
-      Gridcell[["Stand"]][["landcovertype"]] <- readBin(zz, what = integer(), n = 1, size = 4)
+      
+      ## Past version
+      #Gridcell[["Stand"]][["landcovertype"]] <- readBin(zz, what = integer(), n = 1, size = 4)
+      
+      # # Landcover will be read again under stand. So "landcovertype" here is meaningless but we need to read/write.  
+      Gridcell[["landcovertype"]] <- readBin(zz, what = integer(), n = 1, size = 4)
       num_stnd <- as.numeric(Gridcell$nstands)
       Gridcell[["Stand"]] <- vector("list", num_stnd) 
       
@@ -652,7 +662,6 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
     
     # "(*this)[*]" points to different things under different levels, here it is stand
     if(grepl(utils::glob2rx("(*this)[*]"), current_stream)){ # note that first else-part will be evaluated considering the order in guess.cpp
-      
       # STAND
       level <- "Stand"
       current_stream <- "Stand"
@@ -669,15 +678,16 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
       
       for(stnd_i in seq_len(num_stnd)){ #looping over the stands
         for(svs_i in seq_along(streamed_vars_stand)){ # looping over the streamed stand vars
-          
           current_stream <- streamed_vars_stand[svs_i]
           if(grepl(utils::glob2rx("pft[*]"), current_stream)) current_stream <- paste0(level, "pft") # i counter might change, using wildcard
           
           if(current_stream == "nobj" & level == "Stand"){
-            # nobj points to different things under different levels, here it is the number of patches
+            # nobj: Number of Patches
             # number of patches is set through insfiles, read by write.configs and passed to this fcn
             # but it's also written to the state file, need to move bytes
+            pos <- seek(zz)
             nofpatch <- readBin(zz, integer(), 1, size = 4)  
+            # browser()
             if(npatches == nofpatch){ # also not a bad place to check if everything is going fine so far
               Gridcell[["Stand"]][[stnd_i]]$npatches <- npatches
               #Gridcell[["Stand"]] <- vector("list", npatches) 
@@ -687,8 +697,8 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
             next
           }
           
-          # "(*this)[*]" points to different things under different levels, here it is patch
-          if(grepl(utils::glob2rx("(*this)[*]"), current_stream)){ 
+          ##### "(*this)[*]" points to different things under different levels, here it is PATCH ####
+          if(grepl(utils::glob2rx("(*this)[*]"), current_stream)){
             # PATCH
             level <- "Patch"
             current_stream <- "Patch"
@@ -704,6 +714,7 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
             
             for(ptch_i in seq_len(npatches)){ #looping over the patches
               for(svp_i in seq_along(streamed_vars_patch)){ #looping over the streamed patch vars
+                # if(svp_i == 17) browser()
                 current_stream <- streamed_vars_patch[svp_i]
                 if(grepl(utils::glob2rx("pft[*]"), current_stream)) current_stream <- paste0(level, "pft") # i counter might change, using wildcard
                 
@@ -734,16 +745,20 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                     streamed_vars_veg <- find_stream_var(file_in = guesscpp_in, line_nos = beg_end)
                     
                     # NOTE : Unlike other parts, this bit is a lot less generalized!!!
-                    # I'm gonna asumme Vegetation class won't change much in the future
+                    # I'm gonna assume Vegetation class won't change much in the future
                     # indiv.pft.id and indiv needs to be looped over nobj times
                     if(!setequal(streamed_vars_veg, c("nobj", "indiv.pft.id", "indiv"))){
                       PEcAn.logger::logger.severe("Vegetation class object changed in this model version, you need to fix read.state")
                     }
                     
                     # nobj points to different things under different levels, here it is the number of individuals
+                    pos <- seek(zz)
                     number_of_individuals <- readBin(zz, integer(), 1, size = 4) 
                     Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Vegetation"]] <- list()
                     Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Vegetation"]][["number_of_individuals"]] <- number_of_individuals
+                    key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Vegetation", "number_of_individuals", fsep = "/")
+                    pos_list[[key]] <- pos
+                    siz_list[[key]] <- 4
                     
                     # few checks for sensible vals
                     if(number_of_individuals < 0 | number_of_individuals > 10000){ # should there be an upper limit here too?
@@ -764,7 +779,12 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                     for(indv_i in seq_len(number_of_individuals)){
                       Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Vegetation"]][["Individuals"]][[indv_i]] <- list()
                       # which PFT is this?
+                      pos <- seek(zz)
                       Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Vegetation"]][["Individuals"]][[indv_i]][["indiv.pft.id"]] <- readBin(zz, integer(), 1, size = 4)
+                      key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Vegetation", "Individuals", indv_i, "indiv.pft.id", fsep = "/")
+                      pos_list[[key]] <- pos
+                      siz_list[[key]] <- 4
+                      
                       # read all the individual class
                       for(svi_i in seq_along(streamed_vars_indv)){ # 
                         
@@ -797,11 +817,14 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                             
                             current_stream_specs <- find_stream_size(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS)
                             # and read!
-                            
+                            pos <- seek(zz)
                             Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Vegetation"]][["Individuals"]][[indv_i]][["PhotosynthesisResult"]][[current_stream_type$name]] <- readBin(con  = zz, 
                                                                                                                                                                                            what = current_stream_specs$what, 
                                                                                                                                                                                            n    = current_stream_specs$n, 
                                                                                                                                                                                            size = current_stream_specs$size)
+                            key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Vegetation", "Individuals", indv_i, "PhotosynthesisResult", current_stream_type$name, fsep = "/")
+                            pos_list[[key]] <- pos
+                            siz_list[[key]] <- current_stream_specs$size
                             
                           }# streamed_vars_photo-loop ends
                           
@@ -812,16 +835,24 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                           current_stream_specs <- find_stream_size(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS)
                           
                           if(current_stream_specs$single){
+                            pos <- seek(zz)
                             Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Vegetation"]][["Individuals"]][[indv_i]][[current_stream_type$name]] <- readBin(con  = zz, 
                                                                                                                                                                  what = current_stream_specs$what, 
                                                                                                                                                                  n    = current_stream_specs$n, 
                                                                                                                                                                  size = current_stream_specs$size)
+                            key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Vegetation", "Individuals", indv_i, current_stream_type$name, fsep = "/")
+                            pos_list[[key]] <- pos
+                            siz_list[[key]] <- current_stream_specs$size 
                           }else{
                             for(css.i in seq_along(current_stream_specs$what)){
+                              pos <- seek(zz)
                               Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Vegetation"]][["Individuals"]][[indv_i]][[current_stream_specs$names[css.i]]]<- readBin(con  = zz, 
                                                                                                                                                                            what = current_stream_specs$what[css.i], 
                                                                                                                                                                            n    = current_stream_specs$n[css.i], 
                                                                                                                                                                            size = current_stream_specs$size[css.i])
+                              key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Vegetation", "Individuals", indv_i, current_stream_specs$names[css.i], fsep = "/")
+                              pos_list[[key]] <- pos
+                              siz_list[[key]] <- current_stream_specs$size[css.i] 
                             }
                           }
                         }
@@ -845,18 +876,33 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                     # parse from guess.h
                     PerPFTFluxType <- c("NPP", "GPP", "RA", "ISO", "MON")
                     Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["annual_fluxes_per_pft"]] <- list()
-                    key1 <- readBin(zz, "integer", 1, 8)
+                    # The number of PFTS
+                    pos <- seek(zz)
+                    key1 <- readBin(zz, "integer", 1, 8) 
                     Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["annual_fluxes_per_pft"]][["n_pft"]] <- key1
-                    for(fpft_i in seq_len(key1)){ # key1 11 PFTs
-                      Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["annual_fluxes_per_pft"]][[fpft_i]] <- list()
+                    key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Fluxes", "annual_fluxes_per_pft", "n_pft", fsep = "/")
+                    pos_list[[key]] <- pos
+                    siz_list[[key]] <- 8
+                    
+                    for(fpft_i in seq_len(key1)){ # key1 12 PFTs
+                      Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["annual_fluxes_per_pft"]][[paste0("pft", fpft_i)]] <- list()
+                      pos <- seek(zz)
                       key2 <- readBin(zz, "integer", 1, 8)
                       if(key2 > 10000){ #make sure you dind't read a weird number, this is supposed to be number of fluxes per pft, can't have too many
                         PEcAn.logger::logger.severe("Number of fluxes per pft read from the state file is too high. Check read.state function")
                       }
-                      Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["annual_fluxes_per_pft"]][[fpft_i]][["key2"]] <- key2
+                      Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["annual_fluxes_per_pft"]][[paste0("pft", fpft_i)]][["key2"]] <- key2
+                      key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Fluxes", "annual_fluxes_per_pft", paste0("pft", fpft_i), "key2", fsep = "/")
+                      pos_list[[key]] <- pos
+                      siz_list[[key]] <- 8
+                      
                       for(flux_i in seq_len(key2)){
                         # is this double?
-                        Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["annual_fluxes_per_pft"]][[fpft_i]][[PerPFTFluxType[flux_i]]] <- readBin(zz, "double", 1, 8)
+                        pos <- seek(zz)
+                        Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["annual_fluxes_per_pft"]][[paste0("pft", fpft_i)]][[PerPFTFluxType[flux_i]]] <- readBin(zz, "double", 1, 8)
+                        key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Fluxes", "annual_fluxes_per_pft", paste0("pft", fpft_i), PerPFTFluxType[flux_i], fsep = "/")
+                        pos_list[[key]] <- pos
+                        siz_list[[key]] <- 8
                       }
                     }
                     
@@ -864,16 +910,25 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                     # double monthly_fluxes_patch[12][NPERPATCHFLUXTYPES];
                     # maybe read this as a matrix?
                     n_monthly_fluxes_patch <- 12 * LPJ_GUESS_CONST_INTS$val[LPJ_GUESS_CONST_INTS$var =="PerPatchFluxType"]
+                    pos <- seek(zz)
                     Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["monthly_fluxes_patch"]] <- readBin(zz, "double", n_monthly_fluxes_patch, 8)
+                    key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Fluxes", "monthly_fluxes_patch", fsep = "/")
+                    pos_list[[key]] <- pos
+                    siz_list[[key]] <- 8
                     
                     # monthly_fluxes_pft read as a vector at once
                     # double monthly_fluxes_pft[12][NPERPFTFLUXTYPES];
                     # maybe read this as a matrix?
                     n_monthly_fluxes_pft <- 12 * LPJ_GUESS_CONST_INTS$val[LPJ_GUESS_CONST_INTS$var =="PerPFTFluxType"]
+                    pos <- seek(zz)
                     Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Fluxes"]][["monthly_fluxes_pft"]] <- readBin(zz, "double", n_monthly_fluxes_pft, 8)
+                    key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Fluxes", "monthly_fluxes_pft", fsep = "/")
+                    pos_list[[key]] <- pos
+                    siz_list[[key]] <- 8
                     
                   }else{
-                    # NOT VEGETATION OR FLUX
+                    # NOT VEGETATION OR FLUX. 
+                    # Patchpft or Soil in this case
                     streamed_vars <- find_stream_var(file_in = guesscpp_in, line_nos = beg_end)
                     # NO CROPS, NATURAL VEG
                     if("*cropphen" %in% streamed_vars) streamed_vars <- streamed_vars[!(streamed_vars == "*cropphen")]
@@ -883,6 +938,12 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                       Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][[current_stream_type$name]][[varname]] <- vector("list", num_pft) 
                     }
                     
+                    if (current_stream == "soil"){
+                      past_stream <-  tools::toTitleCase(current_stream)
+                    } else{
+                      past_stream <-  current_stream
+                    }
+                    
                     # maybe try modifying this bit later to make it a function
                     for(pft_i in seq_len(num_pft)){
                       for(sv_i in seq_along(streamed_vars)){ 
@@ -921,10 +982,17 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                               current_stream_specs <- find_stream_size(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS)
                               
                               if(current_stream_specs$single){
+                                pos <- seek(zz)
                                 Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][["Soil"]][["Sompool"]][[current_stream_type$name]][[som_i]] <- readBin(con  = zz,
                                                                                                                                                           what = current_stream_specs$what, 
                                                                                                                                                           n    = current_stream_specs$n, 
                                                                                                                                                           size = current_stream_specs$size)
+                                
+
+                                key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, "Soil", "Sompool", current_stream_type$name, som_i, fsep = "/")
+                                pos_list[[key]] <- pos
+                                siz_list[[key]] <- current_stream_specs$size 
+                                
                               }else{
                                 PEcAn.logger::logger.severe("Historic under sompool.") # Not expecting any
                               }
@@ -935,16 +1003,24 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                           current_stream_specs <- find_stream_size(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS)
                           # and read!
                           if(current_stream_specs$single){ # maybe use current_stream in sublist names to find correct place
-                            Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][[length( Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]])]][[current_stream_type$name]][[pft_i]] <- readBin(con  = zz, 
+                            pos <- seek(zz)
+                            Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][[length(Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]])]][[current_stream_type$name]][[pft_i]] <- readBin(con  = zz, 
                                                                                                                                                                                               what = current_stream_specs$what, 
                                                                                                                                                                                               n    = current_stream_specs$n, 
                                                                                                                                                                                               size = current_stream_specs$size)
+                            key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, past_stream, current_stream_type$name, pft_i, fsep = "/")
+                            pos_list[[key]] <- pos
+                            siz_list[[key]] <- current_stream_specs$size 
                           }else{ # only for historic type?
                             for(css.i in seq_along(current_stream_specs$what)){ # maybe use current_stream in sublist names to find correct place
-                              Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][[length( Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]])]][[current_stream_specs$names[css.i]]]<- readBin(con  = zz, 
+                              pos <- seek(zz)
+                              Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][[length(Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]])]][[current_stream_specs$names[css.i]]]<- readBin(con  = zz, 
                                                                                                                                                                                                what = current_stream_specs$what[css.i], 
                                                                                                                                                                                                n    = current_stream_specs$n[css.i], 
                                                                                                                                                                                                size = current_stream_specs$size[css.i])
+                              key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, past_stream, current_stream_type$names[css.i], pft_i, fsep = "/")
+                              pos_list[[key]] <- pos
+                              siz_list[[key]] <- current_stream_specs$size[css.i]
                             }
                           }
                         }
@@ -958,18 +1034,25 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                   current_stream_specs <- find_stream_size(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS)
                   # and read!
                   if(current_stream_specs$single){
-                    
+                    pos <- seek(zz)
                     Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][[current_stream_type$name]] <- readBin(con  = zz, 
                                                                                                               what = current_stream_specs$what, 
                                                                                                               n    = current_stream_specs$n, 
                                                                                                               size = current_stream_specs$size)
+                    key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, current_stream_type$name, fsep = "/")
+                    pos_list[[key]] <- pos
+                    siz_list[[key]] <- current_stream_specs$size
                   }else{ # probably don't need this but let's keep
                     for(css_i in seq_along(current_stream_specs$what)){
-                      # CHANGE ALL THESE HISTORIC TYPES SO THAT cirrent_index and full goes together with the variable
+                      # CHANGE ALL THESE HISTORIC TYPES SO THAT current_index and full goes together with the variable
+                      pos <- seek(zz)
                       Gridcell[["Stand"]][[stnd_i]][["Patch"]][[ptch_i]][[current_stream_specs$names[css_i]]] <- readBin(con  = zz, 
                                                                                                                          what = current_stream_specs$what[css_i], 
                                                                                                                          n    = current_stream_specs$n[css_i], 
                                                                                                                          size = current_stream_specs$size[css_i])
+                      key <- file.path("Gridcell", "Stand", stnd_i, "Patch", ptch_i, current_stream_specs$names[css_i], fsep = "/")
+                      pos_list[[key]] <- pos
+                      siz_list[[key]] <- current_stream_specs$size[css_i] 
                     }
                   }
                 }# end if-class within Patch
@@ -988,6 +1071,9 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
             Gridcell[["Stand"]][[stnd_i]][[length(Gridcell[["Stand"]][[stnd_i]])+1]] <- list()
             names(Gridcell[["Stand"]][[stnd_i]])[length(Gridcell[["Stand"]][[stnd_i]])] <- current_stream_type$name
             
+            # Save the past stream like Standpft
+            past_stream <- current_stream
+            
             if(current_stream_type$type == "class"){
               
               # CLASS
@@ -1019,16 +1105,24 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
                     current_stream_specs <- find_stream_size(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS)
                     # and read!
                     if(current_stream_specs$single){
+                      pos <- seek(zz)
                       Gridcell[["Stand"]][[stnd_i]][[length(Gridcell[["Stand"]][[stnd_i]])]][[current_stream_type$name]][[pft_i]] <- readBin(con  = zz, 
                                                                                                                                              what = current_stream_specs$what, 
                                                                                                                                              n    = current_stream_specs$n, 
                                                                                                                                              size = current_stream_specs$size)
+                      key <- file.path("Gridcell", "Stand", stnd_i, past_stream, current_stream_type$name, pft_i, fsep = "/")
+                      pos_list[[key]] <- pos
+                      siz_list[[key]] <- current_stream_specs$size
                     }else{
                       for(css.i in seq_along(current_stream_specs$what)){
+                        pos <- seek(zz)
                         Gridcell[[length(Gridcell)]][[current_stream_type$name]][[pft_i]][[current_stream_specs$names[css.i]]]<- readBin(con  = zz, 
                                                                                                                                          what = current_stream_specs$what[css.i], 
                                                                                                                                          n    = current_stream_specs$n[css.i], 
                                                                                                                                          size = current_stream_specs$size[css.i])
+                        key <- file.path("Gridcell", "Stand", stnd_i, past_stream, current_stream_type$name[css.i], fsep = "/")
+                        pos_list[[key]] <- pos
+                        siz_list[[key]] <- current_stream_specs$size[css.i]
                       }
                     }
                   }
@@ -1040,16 +1134,24 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
               current_stream_specs <- find_stream_size(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS)
               # and read!
               if(current_stream_specs$single){
+                pos <- seek(zz)
                 Gridcell[["Stand"]][[stnd_i]][[current_stream_type$name]] <- readBin(con  = zz, 
                                                                                      what = current_stream_specs$what, 
                                                                                      n    = current_stream_specs$n, 
                                                                                      size = current_stream_specs$size)
+                key <- file.path("Gridcell", "Stand", stnd_i, current_stream_type$name, fsep = "/")
+                pos_list[[key]] <- pos
+                siz_list[[key]] <- current_stream_specs$size
               }else{ # probably don't need this but let's keep
                 for(css_i in seq_along(current_stream_specs$what)){
+                  pos <- seek(zz)
                   Gridcell[[length(Gridcell)]][[current_stream_specs$names[css_i]]] <- readBin(con  = zz, 
                                                                                                what = current_stream_specs$what[css_i], 
                                                                                                n    = current_stream_specs$n[css_i], 
                                                                                                size = current_stream_specs$size[css_i])
+                  key <- file.path("Gridcell", "Stand", stnd_i, current_stream_type$name, fsep = "/")
+                  pos_list[[key]] <- pos
+                  siz_list[[key]] <- current_stream_specs$size[css_i] 
                 }
               }
             }# end if-class within Stand
@@ -1062,6 +1164,7 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
     }else{ #not reading in Stand variables
       
       # NOT STAND
+      past_stream <- current_stream
       
       current_stream_type <- find_stream_type(NULL, current_stream, LPJ_GUESS_CLASSES, LPJ_GUESS_TYPES, guessh_in)
       
@@ -1092,27 +1195,40 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
             current_stream_specs <- find_stream_size(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS)
             # and read!
             if(current_stream_specs$single){
+              pos <- seek(zz)
               Gridcell[[length(Gridcell)]][[current_stream_type$name]][[pft_i]] <- readBin(con  = zz, 
                                                                                            what = current_stream_specs$what, 
                                                                                            n    = current_stream_specs$n, 
                                                                                            size = current_stream_specs$size)
-            }else if(current_stream_specs$name %in% c("hmtemp_20", "hmprec_20", "hmeet_20")){
+              key <- file.path("Gridcell", past_stream, current_stream_type$name, pft_i, fsep = "/")
+              pos_list[[key]] <- pos
+              siz_list[[key]] <- current_stream_specs$size 
+              
+            }else if(current_stream_type$name %in% c("hmtemp_20", "hmprec_20", "hmeet_20")){
               # these three are just too different, maybe extract their names in the beginning
               # be careful while writing back to the binary
               # Gridcell[[length(Gridcell)]][[current_stream_type$name]] <- readBin(con = zz, double(), 264, 8)
               Gridcell[[length(Gridcell)]][[current_stream_type$name]] <- vector("list", length(current_stream_specs) - 2)
               for(css.i in seq_len(length(current_stream_specs) - 2)){
+                 pos <- seek(zz)
                  Gridcell[[length(Gridcell)]][[current_stream_type$name]][[css.i]] <- readBin(con  = zz, 
                                                                  what = current_stream_specs[[css.i]]$what, 
                                                                  n    = current_stream_specs[[css.i]]$n, 
                                                                  size = current_stream_specs[[css.i]]$size)
+                 key <- file.path("Gridcell", past_stream, current_stream_type$name, css.i, fsep = "/")
+                 pos_list[[key]] <- pos
+                 siz_list[[key]] <- current_stream_specs[[css.i]]$size 
               }
             }else{
               for(css.i in seq_along(current_stream_specs$what)){
+                pos <- seek(zz)
                 Gridcell[[length(Gridcell)]][[current_stream_type$name]][[pft_i]][[current_stream_specs$names[css.i]]]<- readBin(con  = zz, 
                                                                                                                                  what = current_stream_specs$what[css.i], 
                                                                                                                                  n    = current_stream_specs$n[css.i], 
                                                                                                                                  size = current_stream_specs$size[css.i])
+                key <- file.path("Gridcell", past_stream, current_stream_type$name, pft_i, current_stream_specs$names[css.i], fsep = "/")
+                pos_list[[key]] <- pos
+                siz_list[[key]] <- current_stream_specs$size[css.i] 
               }
             }
             
@@ -1124,16 +1240,24 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
         current_stream_specs <- find_stream_size(current_stream_type, guessh_in, LPJ_GUESS_TYPES, LPJ_GUESS_CONST_INTS)
         # and read!
         if(current_stream_specs$single){
+          pos <- seek(zz)
           Gridcell[[length(Gridcell)]][[current_stream_type$name]] <- readBin(con  = zz, 
                                                                               what = current_stream_specs$what, 
                                                                               n    = current_stream_specs$n, 
                                                                               size = current_stream_specs$size)
-        }else{ # probably don't need this but let's keep
+          key <- file.path("Gridcell", past_stream, current_stream_type$name, fsep = "/")
+          pos_list[[key]] <- pos
+          siz_list[[key]] <- current_stream_specs$size
+        }else{
           for(css_i in seq_along(current_stream_specs$what)){
+            pos <- seek(zz)
             Gridcell[[length(Gridcell)]][[current_stream_specs$names[css_i]]] <- readBin(con  = zz, 
                                                                                          what = current_stream_specs$what[css_i], 
                                                                                          n    = current_stream_specs$n[css_i], 
                                                                                          size = current_stream_specs$size[css_i])
+            key <- file.path("Gridcell", past_stream, current_stream_type$name[css_i], fsep = "/")
+            pos_list[[key]] <- pos
+            siz_list[[key]] <- current_stream_specs$size[css_i] 
           }
         }
       }# end if-class within Gridcell
@@ -1145,7 +1269,12 @@ read_binary_LPJGUESS <- function(outdir, version = "PalEON"){
   
   Gridcell$meta_data <- meta_data
   
-  return(Gridcell)
+  # return(Gridcell)
+  return(list(
+    state = Gridcell,
+    pos_list = pos_list,
+    siz_list = siz_list
+  ))
 } # read_binary_LPJGUESS end
 
 
diff --git a/models/lpjguess/R/split_inputs.LPJGUESS.R b/models/lpjguess/R/split_inputs.LPJGUESS.R
index 9a598e50e8b..8ba5ef99491 100644
--- a/models/lpjguess/R/split_inputs.LPJGUESS.R
+++ b/models/lpjguess/R/split_inputs.LPJGUESS.R
@@ -14,120 +14,125 @@
 ##' @importFrom PEcAn.utils days_in_year
 ##' @export
 split_inputs.LPJGUESS <- function(settings, start.time, stop.time, inputs, overwrite = FALSE, outpath = NULL){
-  
-  #### Lubridate start and end times
-  start.day  <- lubridate::yday(start.time)
-  start.year <- lubridate::year(start.time)
-  end.day    <- lubridate::yday(stop.time)
-  end.year   <- lubridate::year(stop.time)
-  
-  # Whole run period
-  run.start <- lubridate::year(settings$run$start.date)
-  run.end   <- lubridate::year(settings$run$end.date)
-  
-  #### Get met paths
-  met    <- inputs
-  path   <- dirname(met)
-  prefix <- substr(basename(met), 1, nchar(basename(met))-16) #assuming we'll always have "PREFIX.1920.2010.tmp"
-  if(is.null(outpath)){
-    outpath <- path
-  }
-  if(!dir.exists(outpath)) dir.create(outpath)
-  
-  var.names  <- c("tmp", "pre", "cld")
-  long.names <- c("air_temperature",
-                  "precipitation_flux",
-                  "surface_downwelling_shortwave_flux_in_air")
-  
-  # !!! always full years with LPJ-GUESS !!!
-  files.in  <- file.path(outpath, paste0(prefix, run.start,  ".", run.end,  ".", var.names, ".nc"))
-  files.out <- file.path(outpath, paste0(prefix, start.year, ".", end.year, ".", var.names, ".nc"))
-  
-  if(file.exists(files.out[1]) & !overwrite){
-    return(files.out[1])
-  }
-  
-  ## open netcdf files
-  fnc.tmp <- ncdf4::nc_open(files.in[1])
-  fnc.pre <- ncdf4::nc_open(files.in[2])
-  fnc.cld <- ncdf4::nc_open(files.in[3])
-  
-  ## read climate data
-  nc.tmp <- ncdf4::ncvar_get(fnc.tmp, var.names[1])
-  nc.pre <- ncdf4::ncvar_get(fnc.pre, var.names[2])
-  nc.cld <- ncdf4::ncvar_get(fnc.cld, var.names[3])
-
-  # cut where
-  if(start.year == run.start){
-    years <- start.year:end.year
-    inds  <- 1:sum(PEcAn.utils::days_in_year(years))
+  #### If using CRU input, return directly
+  if (grepl("\\.cru(\\.bin)?$", inputs, ignore.case = TRUE)) {
+    PEcAn.logger::logger.info(paste("Input is a CRU file:", inputs, "- returning path directly without splitting."))
+    return(inputs)   # Without cropping, use the original file directly
   }else{
-    ### come back
-  }
-
-  # split
-  nc.tmp <- nc.tmp[1,1,inds]
-  nc.pre <- nc.pre[1,1,inds]
-  nc.cld <- nc.cld[1,1,inds]
-
-  var.list <- list(nc.tmp, nc.pre, nc.cld)
-  
-  # not that these will be different than "K", "kg m-2 s-1", "W m-2"
-  var.units <- c(fnc.tmp$var$tmp$units,
-                 fnc.pre$var$pre$units,
-                 fnc.cld$var$cld$units)
-    
-  # get other stuff to be written to ncdf
-  
-  ## retrieve lat/lon
-  lon <- ncdf4::ncvar_get(fnc.tmp, "lon")
-  lat <- ncdf4::ncvar_get(fnc.tmp, "lat")
-
-  # write back
-  ## write climate data define dimensions
-  
-  latdim  <- ncdf4::ncdim_def(name = "lat", "degrees_north", as.double(lat))
-  londim  <- ncdf4::ncdim_def(name = "lon", "degrees_east", as.double(lon))
-  timedim <- ncdf4::ncdim_def("time", paste0("days since ", start.year - 1, "-12-31", sep = ""), as.double(c(1:length(nc.tmp))))
-  
-  fillvalue <- 9.96920996838687e+36
-  
-  for (n in seq_along(var.names)) {
-    # define variable
-    var.def <- ncdf4::ncvar_def(name = var.names[n],
-                                units = var.units[n],
-                                dim = (list(londim, latdim, timedim)),
-                                fillvalue, long.names[n],
-                                verbose = FALSE,
-                                prec = "float")
-    
-    # create netCD file for LPJ-GUESS
-    ncfile <- ncdf4::nc_create(files.out[[n]], vars = var.def, force_v4 = TRUE)
-    
-
-    # put variable, rep(...,each=4) is a hack to write the same data for all grids (which all are the
-    # same)
-    ncdf4::ncvar_put(ncfile, var.def, rep(var.list[[n]], each = 4))
-
-    
-    # additional attributes for LPJ-GUESS
-    ncdf4::ncatt_put(nc = ncfile, varid = var.names[n], attname = "standard_name", long.names[n])
-    
-    ncdf4::ncatt_put(nc = ncfile, varid = "lon", attname = "axis", "X")
-    ncdf4::ncatt_put(nc = ncfile, varid = "lon", attname = "standard_name", "longitude")
-    
-    ncdf4::ncatt_put(nc = ncfile, varid = "lat", attname = "axis", "Y")
-    ncdf4::ncatt_put(nc = ncfile, varid = "lat", attname = "standard_name", "latitude")
-    
-    ncdf4::ncatt_put(nc = ncfile, varid = "time", attname = "calendar", "gregorian")
-    
-    ncdf4::nc_close(ncfile)
+    #### Lubridate start and end times
+    start.day  <- lubridate::yday(start.time)
+    start.year <- lubridate::year(start.time)
+    end.day    <- lubridate::yday(stop.time)
+    end.year   <- lubridate::year(stop.time)
+    
+    # Whole run period
+    run.start <- lubridate::year(settings$run$start.date)
+    run.end   <- lubridate::year(settings$run$end.date)
+    
+    #### Get met paths
+    met    <- inputs
+    path   <- dirname(met)
+    prefix <- substr(basename(met), 1, nchar(basename(met))-16) #assuming we'll always have "PREFIX.1920.2010.tmp"
+    if(is.null(outpath)){
+      outpath <- path
+    }
+    if(!dir.exists(outpath)) dir.create(outpath)
+    
+    var.names  <- c("tmp", "pre", "cld")
+    long.names <- c("air_temperature",
+                    "precipitation_flux",
+                    "surface_downwelling_shortwave_flux_in_air")
+    
+    # !!! always full years with LPJ-GUESS !!!
+    files.in  <- file.path(outpath, paste0(prefix, run.start,  ".", run.end,  ".", var.names, ".nc"))
+    files.out <- file.path(outpath, paste0(prefix, start.year, ".", end.year, ".", var.names, ".nc"))
+    
+    if(file.exists(files.out[1]) & !overwrite){
+      return(files.out[1])
+    }
+    
+    ## open netcdf files
+    fnc.tmp <- ncdf4::nc_open(files.in[1])
+    fnc.pre <- ncdf4::nc_open(files.in[2])
+    fnc.cld <- ncdf4::nc_open(files.in[3])
+    
+    ## read climate data
+    nc.tmp <- ncdf4::ncvar_get(fnc.tmp, var.names[1])
+    nc.pre <- ncdf4::ncvar_get(fnc.pre, var.names[2])
+    nc.cld <- ncdf4::ncvar_get(fnc.cld, var.names[3])
+    
+    # cut where
+    if(start.year == run.start){
+      years <- start.year:end.year
+      inds  <- 1:sum(PEcAn.utils::days_in_year(years))
+    }else{
+      ### come back
+    }
+    
+    # split
+    nc.tmp <- nc.tmp[1,1,inds]
+    nc.pre <- nc.pre[1,1,inds]
+    nc.cld <- nc.cld[1,1,inds]
+    
+    var.list <- list(nc.tmp, nc.pre, nc.cld)
+    
+    # not that these will be different than "K", "kg m-2 s-1", "W m-2"
+    var.units <- c(fnc.tmp$var$tmp$units,
+                   fnc.pre$var$pre$units,
+                   fnc.cld$var$cld$units)
+    
+    # get other stuff to be written to ncdf
+    
+    ## retrieve lat/lon
+    lon <- ncdf4::ncvar_get(fnc.tmp, "lon")
+    lat <- ncdf4::ncvar_get(fnc.tmp, "lat")
+    
+    # write back
+    ## write climate data define dimensions
+    
+    latdim  <- ncdf4::ncdim_def(name = "lat", "degrees_north", as.double(lat))
+    londim  <- ncdf4::ncdim_def(name = "lon", "degrees_east", as.double(lon))
+    timedim <- ncdf4::ncdim_def("time", paste0("days since ", start.year - 1, "-12-31", sep = ""), as.double(c(1:length(nc.tmp))))
+    
+    fillvalue <- 9.96920996838687e+36
+    
+    for (n in seq_along(var.names)) {
+      # define variable
+      var.def <- ncdf4::ncvar_def(name = var.names[n],
+                                  units = var.units[n],
+                                  dim = (list(londim, latdim, timedim)),
+                                  fillvalue, long.names[n],
+                                  verbose = FALSE,
+                                  prec = "float")
+      
+      # create netCD file for LPJ-GUESS
+      ncfile <- ncdf4::nc_create(files.out[[n]], vars = var.def, force_v4 = TRUE)
+      
+      
+      # put variable, rep(...,each=4) is a hack to write the same data for all grids (which all are the
+      # same)
+      ncdf4::ncvar_put(ncfile, var.def, rep(var.list[[n]], each = 4))
+      
+      
+      # additional attributes for LPJ-GUESS
+      ncdf4::ncatt_put(nc = ncfile, varid = var.names[n], attname = "standard_name", long.names[n])
+      
+      ncdf4::ncatt_put(nc = ncfile, varid = "lon", attname = "axis", "X")
+      ncdf4::ncatt_put(nc = ncfile, varid = "lon", attname = "standard_name", "longitude")
+      
+      ncdf4::ncatt_put(nc = ncfile, varid = "lat", attname = "axis", "Y")
+      ncdf4::ncatt_put(nc = ncfile, varid = "lat", attname = "standard_name", "latitude")
+      
+      ncdf4::ncatt_put(nc = ncfile, varid = "time", attname = "calendar", "gregorian")
+      
+      ncdf4::nc_close(ncfile)
+    }
+    
+    # close nc
+    ncdf4::nc_close(fnc.tmp)
+    ncdf4::nc_close(fnc.pre)
+    ncdf4::nc_close(fnc.cld)
+    
+    return(files.out[1])
   }
-  
-  # close nc
-  ncdf4::nc_close(fnc.tmp)
-  ncdf4::nc_close(fnc.pre)
-  ncdf4::nc_close(fnc.cld)
-  
-  return(files.out[1])
 } # split_inputs.LPJGUESS
\ No newline at end of file
diff --git a/models/lpjguess/R/write.config.LPJGUESS.R b/models/lpjguess/R/write.config.LPJGUESS.R
index 6f2ea4f5d7f..b1ac07b5c45 100644
--- a/models/lpjguess/R/write.config.LPJGUESS.R
+++ b/models/lpjguess/R/write.config.LPJGUESS.R
@@ -14,7 +14,6 @@
 ##' @export
 ##' @author Istem Fer, Tony Gardella
 write.config.LPJGUESS <- function(defaults, trait.values, settings, run.id, restart = NULL) {
-  
   # find out where to write run/ouput
   rundir <- file.path(settings$host$rundir, run.id)
   if (!file.exists(rundir)) {
@@ -194,16 +193,20 @@ write.insfile.LPJGUESS <- function(settings, trait.values, rundir, outdir, run.i
   paramsins <- paramsins[-pftindx] 
   paramsins <- c(paramsins, unlist(write2pftblock))
   
-  
-  # write clim file names
-  
-  tmp.file <- settings$run$inputs$met$path
-  pre.file <- gsub(".tmp.nc", ".pre.nc", tmp.file)
-  cld.file <- gsub(".tmp.nc", ".cld.nc", tmp.file)
-  
-  guessins <- gsub("@TEMP_FILE@", tmp.file, guessins)
-  guessins <- gsub("@PREC_FILE@", pre.file, guessins)
-  guessins <- gsub("@INSOL_FILE@", cld.file, guessins)
+  # # Past version: write clim file names (cf input)
+  # tmp.file <- settings$run$inputs$met$path
+  # pre.file <- gsub(".tmp.nc", ".pre.nc", tmp.file)
+  # cld.file <- gsub(".tmp.nc", ".cld.nc", tmp.file)
+  # 
+  # guessins <- gsub("@TEMP_FILE@", tmp.file, guessins)
+  # guessins <- gsub("@PREC_FILE@", pre.file, guessins)
+  # guessins <- gsub("@INSOL_FILE@", cld.file, guessins)
+  
+  # when using cru input, lpjguess will not use these clim files
+  cru.file <- settings$run$inputs$met$path
+  misc.file <- sub("\\.bin$", "misc.bin", cru.file)
+  guessins <- gsub("@MET_AND_SOIL_FILE@", cru.file, guessins)
+  guessins <- gsub("@MISC_FILE@", misc.file, guessins)
   
   # create and write CO2 file
   start.year <- lubridate::year(settings$run$start.date)
@@ -232,9 +235,12 @@ write.insfile.LPJGUESS <- function(settings, trait.values, rundir, outdir, run.i
   utils::write.table(CO2, file = co2.file, row.names = FALSE, col.names = FALSE, sep = "\t", eol = "\n")
   guessins <- gsub("@CO2_FILE@", co2.file, guessins)
   
-  # write soil file path
-  soil.file <- settings$run$inputs$soil$path
-  guessins <- gsub("@SOIL_FILE@", soil.file, guessins)
+  # # write soil file path
+  # # when using cru input, it's also climate file
+  # soil.file <- settings$run$inputs$soil$path
+  # misc.file <- sub("\\.bin$", "misc.bin", soil.file)
+  # guessins <- gsub("@SOIL_FILE@", soil.file, guessins)
+  # guessins <- gsub("@MISC_FILE@", misc.file, guessins)
   
   settings$model$insfile <- file.path(settings$rundir, run.id, "guess.ins")
   
diff --git a/models/lpjguess/R/write_restart.LPJGUESS.R b/models/lpjguess/R/write_restart.LPJGUESS.R
new file mode 100644
index 00000000000..49271ef0e55
--- /dev/null
+++ b/models/lpjguess/R/write_restart.LPJGUESS.R
@@ -0,0 +1,114 @@
+##' write_restart.LPJGUESS
+##'
+##' Write restart files for LPJGUESS
+##' new.state includes X (AGB.pft) from Analysis
+##' new.params includes LPJGUESS_state
+##' 
+##' @param outdir output directory
+##' @param runid run ID
+##' @param start.time start date and time for each SDA ensemble
+##' @param stop.time stop date and time for each SDA ensemble
+##' @param settings PEcAn settings object
+##' @param new.state analysis state vector
+##' @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 verbose decide if we want to print the runid
+##' 
+##' @return NONE
+##'
+##' @export
+##' @author Yinghao Sun
+write_restart.LPJGUESS <- function(outdir, runid,
+                                   start.time, stop.time, settings,
+                                   new.state, RENAME = TRUE,
+                                   new.params, inputs = NULL, verbose = FALSE){
+  
+  rundir <- settings$host$rundir
+  variables <- colnames(new.state)
+  
+  ## ---- Rename old output, remove old clim ----
+  if (RENAME) {
+    file.rename(file.path(outdir, runid, "lpjguess.out"),
+                file.path(outdir, runid, paste0("lpjguess.", as.Date(start.time), ".out")))
+    system(paste("rm", file.path(rundir, runid, "lpjguess.clim")))
+  } else {
+    PEcAn.logger::logger.severe(paste("rename = FALSE: Restart cannot proceed without output file",
+                                      "lpjguess.out being renamed for", start.time))
+    stop("RENAME flag is FALSE. Must rerun this timestep before continuing.")
+  }
+  
+  settings$run$start.date <- start.time
+  settings$run$end.date <- stop.time
+  
+  ## ---- Pull old state ----
+  if (is.null(new.params$LPJGUESS_state))
+    PEcAn.logger::logger.severe("LPJGUESS_state missing in new.params")
+  # new.params$LPJGUESS_state include state, pos_list, siz_list
+  Gridcell <- new.params$LPJGUESS_state$state
+  pos_list <- new.params$LPJGUESS_state$pos_list
+  siz_list <- new.params$LPJGUESS_state$siz_list
+  
+  ## ---- Build PFT parameter table from new.params ----
+  # TODO: find accurate parameters; read params from settings
+  pft_par_table <- data.frame()
+  # PFTs <- c("Ace_rub","Bet_all","Fag_gra","Que_rub","Tsu_can")
+  PFTs <- names(new.params)
+  for(PFT in PFTs) {
+    this.param.row <- c()
+    this.param.row["sla"] <- new.params[[PFT]]$SLA
+    this.param.row["k_latosa"] <- new.params[[PFT]]$sapwood_ratio
+    this.param.row["wooddens"] <- 200   #kg/m-3
+    # this.param.row["wooddens"] <- 0.2 #g/cm-3
+    this.param.row["lifeform"] <- 1
+    this.param.row["k_rp"] <- 1.6
+    this.param.row["k_allom1"] <- 250
+    this.param.row["k_allom2"] <- 60
+    this.param.row["k_allom3"] <- 0.67
+    this.param.row["crownarea_max"] <- 50
+    # conifer special case
+    if(PFT == "Tsu_can") {
+      this.param.row["k_allom1"] <- 150
+    }
+    pft_par_table <- rbind(pft_par_table , this.param.row)
+  }
+  names(pft_par_table) <- c("sla", "k_latosa", "wooddens", "lifeform", "k_rp", "k_allom1", "k_allom2", "k_allom3", "crownarea_max") 
+  rownames(pft_par_table) <- PFTs
+  
+  ## --- Build initial & target AGB vectors (kg m-2) ---
+  agb.init <- calculateGridcellVariablePerPFT(Gridcell, "AbvGrndWood", min.diam=min.diam, pft.params=pft_par_table)
+  if (any(grepl("^AGB.pft", variables))) {          # column names were set in read.restart
+    agb.targ <- PEcAn.utils::ud_convert(
+      unlist(new.state[, grepl("^AGB.pft", variables), drop=TRUE]),
+      "Mg/ha","kg/m^2")
+  }
+  
+  ### dens will not change because we wont do dens SDA temporarily
+  dens.init <- calculateGridcellVariablePerPFT(Gridcell, "densindiv", min.diam=min.diam, pft.params=pft_par_table)
+  dens.targ <- dens.init 
+  
+  ## --- Update state ---
+  # choose a minimum diameter
+  min.diam = 0.5
+    Gridcell_updated <- update_state_LPJGUESS(Gridcell, pft_par_table, 
+                                            dens.init, dens.targ,
+                                            agb.init, agb.targ,
+                                            AbvGrndWood.epsilon = 0.05,
+                                            trace = FALSE, min.diam)
+  
+  State_updated <- list(state = Gridcell_updated,
+                        pos_list = pos_list,
+                        siz_list = siz_list)
+  
+  write_binary_LPJGUESS(State_updated, file.path(outdir, runid))
+  
+  ## --- Regenerate config for next run ---
+  do.call(write.config.LPJGUESS,
+          list(defaults = NULL,
+               trait.values = new.params,
+               settings = settings,
+               run.id   = runid)
+  )
+  
+  if(verbose) PEcAn.logger::logger.info("restart written for", runid)
+}
diff --git a/models/lpjguess/R/write_state.R b/models/lpjguess/R/write_state.R
new file mode 100644
index 00000000000..a65c27e58bf
--- /dev/null
+++ b/models/lpjguess/R/write_state.R
@@ -0,0 +1,103 @@
+#' Extract nested value from a state list using flat key
+#'
+#' @param state A nested list (usually the model.state$state)
+#' @param key A flat string like "Gridcell/Stand/1/Patch/1/Vegetation/Individuals/3/cmass_leaf"
+#' @return The value stored at that nested position
+#' @keywords internal
+#' @author Yinghao Sun
+extract_from_state_by_key <- function(state, key) {
+  # Optional: remove "Gridcell/" prefix
+  key <- sub("^Gridcell/", "", key)
+  
+  parts <- strsplit(key, "/")[[1]]
+  val <- state
+  
+  for (p in parts) {
+    if (is.null(val)) {
+      warning("NULL reached prematurely at: ", p)
+      return(NULL)
+    }
+    
+    # Case 1: numeric index
+    if (grepl("^[0-9]+$", p)) {
+      idx <- as.integer(p)
+      if (idx > length(val)) {
+        warning("Index out of bounds: ", idx)
+        return(NULL)
+      }
+      val <- val[[idx]]
+      
+      # Case 2: named element (case-insensitive match)
+    } else {
+      val_names <- names(val)
+      match_idx <- which(tolower(val_names) == tolower(p))
+      
+      if (length(match_idx) == 0) {
+        warning("Name not found (case-insensitive): ", p)
+        return(NULL)
+      }
+      
+      val <- val[[match_idx[1]]]  # use first match
+    }
+  }
+  
+  return(val)
+}
+
+
+#' Write updated variables into a copy of the original LPJ-GUESS .state file
+#'
+#' @param State_updated    A list containing updated state variables, position list and size list (get from read_binary)
+#' @param outdir           Path to a directory containing the `0.state` and `meta.bin` files.
+#' 
+#' @return No return value. Writes files to disk as side effect.
+#' @author Yinghao Sun
+#' @export
+write_binary_LPJGUESS <- function(State_updated, outdir) {
+  
+  # Build full paths to source files
+  src_state <- file.path(outdir, "0.state")
+  meta_file <- file.path(outdir, "meta.bin")
+  
+  # back-up
+  bak_state <- file.path(outdir, "bak.state")
+  file.copy(src_state, bak_state, overwrite = TRUE)
+  
+  # a copy to the temporary file
+  new_state <- file.path(outdir, "new.state")
+  file.copy(src_state, new_state, overwrite = TRUE)
+  
+  # # Ensure output directory exists
+  # dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
+  # 
+  # # Copy template files to output directory so we don't overwrite it
+  # file.copy(c(meta_file, original_state), to = output_dir, overwrite = TRUE)
+  # 
+  # # Open copied 0.state file for binary modification
+  # state_path <- file.path(outdir, "0.state")
+  # con <- file(state_path, open = "r+b")
+  
+  # Open temporary new.state file for binary modification
+  con <- file(new_state, open = "r+b")
+  
+  # A named list of byte positions for each variable (generated during reading)
+  pos_list <- State_updated$pos_list
+  # A named list of writeBin sizes for each variable (same keys as pos_list)
+  siz_list <- State_updated$siz_list
+  
+  # Loop over all keys
+  for (key in names(pos_list)) {
+    value <- extract_from_state_by_key(State_updated$state, key)
+    pos <- pos_list[[key]]
+    size <- siz_list[[key]]
+    
+    # Seek and write
+    seek(con, where = pos, origin = "start")
+    writeBin(object = value, con = con, size = size)
+  }
+  
+  close(con)
+  
+  # Atomic substitution
+  file.rename(new_state, src_state)   # After success, bak is still there and can be manually deleted
+}
diff --git a/models/lpjguess/inst/pecan.ins b/models/lpjguess/inst/pecan.ins
index d8d52274477..9d9cdb06ad1 100755
--- a/models/lpjguess/inst/pecan.ins
+++ b/models/lpjguess/inst/pecan.ins
@@ -73,7 +73,7 @@ title 'LPJ-GUESS cohort mode - global pfts'
 
 vegmode "cohort"		! "cohort", "individual" or "population"
 nyear_spinup 500		! number of years to spin up the simulation for
-spinup_lifeform "nolifeform"
+! spinup_lifeform "tree"
 ifcalcsla 0				! whether to calculate SLA from leaf longevity
 ! (PFT-specific value can be specified in this file instead)
 ifcalccton 1			! whether to calculate leaf C:N min from leaf longevity
@@ -84,7 +84,7 @@ patcharea 1000			! patch area (m2)
 estinterval 5			! years between establishment events in cohort mode
 ifdisturb 1				! whether generic patch-destroying disturbances enabled
 distinterval 500		! average return time for generic patch-destroying disturbances
-disturb_year -1
+! disturb_year -1
 ifbgestab 1				! whether background establishment enabled
 ifsme 1					! whether spatial mass effect enabled
 ifstochestab 1			! whether establishment stochastic
diff --git a/models/lpjguess/inst/template.ins b/models/lpjguess/inst/template.ins
index aa5520e66c4..25094f1b2e7 100755
--- a/models/lpjguess/inst/template.ins
+++ b/models/lpjguess/inst/template.ins
@@ -14,10 +14,11 @@ coordinates_precision 2
 ! Forcing Data & gridlists
 !
 
-param "file_gridlist_cf" (str "@GRID_FILE@")
+param "file_gridlist" (str "@GRID_FILE@")
 
 param "file_co2"      (str "@CO2_FILE@")
-param "file_cru"      (str "@SOIL_FILE@")
+param "file_cru"      (str "@MET_AND_SOIL_FILE@")
+param "file_cru_misc"  (str "@MISC_FILE@")
 
 ! N deposition (blank string to use constant pre-industrial level of 2 kgN/ha/year)
 param "file_ndep"     (str "")
diff --git a/models/lpjguess/man/extract_from_state_by_key.Rd b/models/lpjguess/man/extract_from_state_by_key.Rd
new file mode 100644
index 00000000000..4d72aface7a
--- /dev/null
+++ b/models/lpjguess/man/extract_from_state_by_key.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/write_state.R
+\name{extract_from_state_by_key}
+\alias{extract_from_state_by_key}
+\title{Extract nested value from a state list using flat key}
+\usage{
+extract_from_state_by_key(state, key)
+}
+\arguments{
+\item{state}{A nested list (usually the model.state$state)}
+
+\item{key}{A flat string like "Gridcell/Stand/1/Patch/1/Vegetation/Individuals/3/cmass_leaf"}
+}
+\value{
+The value stored at that nested position
+}
+\description{
+Extract nested value from a state list using flat key
+}
+\author{
+Yinghao Sun
+}
+\keyword{internal}
diff --git a/models/lpjguess/man/find_stream_size.Rd b/models/lpjguess/man/find_stream_size.Rd
index bc57bf17949..2aebe539c95 100644
--- a/models/lpjguess/man/find_stream_size.Rd
+++ b/models/lpjguess/man/find_stream_size.Rd
@@ -26,3 +26,4 @@ A numeric value representing the size (number of streamed variables).
 \description{
 Determines the size (number of variables) in a stream based on the file content.
 }
+\keyword{internal}
diff --git a/models/lpjguess/man/find_stream_type.Rd b/models/lpjguess/man/find_stream_type.Rd
index d89e8c40136..df03f662ef2 100644
--- a/models/lpjguess/man/find_stream_type.Rd
+++ b/models/lpjguess/man/find_stream_type.Rd
@@ -29,3 +29,4 @@ A character string indicating the stream type.
 \description{
 Determines the type of a given stream variable in an LPJ-GUESS file.
 }
+\keyword{internal}
diff --git a/models/lpjguess/man/find_stream_var.Rd b/models/lpjguess/man/find_stream_var.Rd
index 2de82ee1895..95c672f2c88 100644
--- a/models/lpjguess/man/find_stream_var.Rd
+++ b/models/lpjguess/man/find_stream_var.Rd
@@ -17,3 +17,4 @@ A character vector of streamed variable names.
 \description{
 A helper function that lists streamed variables. It returns the names of streamed variables.
 }
+\keyword{internal}
diff --git a/models/lpjguess/man/read_binary_LPJGUESS.Rd b/models/lpjguess/man/read_binary_LPJGUESS.Rd
index 5ff044916b9..81a68fa47e8 100644
--- a/models/lpjguess/man/read_binary_LPJGUESS.Rd
+++ b/models/lpjguess/man/read_binary_LPJGUESS.Rd
@@ -7,7 +7,7 @@
 read_binary_LPJGUESS(outdir, version = "PalEON")
 }
 \arguments{
-\item{outdir}{A character string specifying the output directory containing the binary state files.}
+\item{outdir}{The output directory where ".state" and "meta.bin" will be written}
 
 \item{version}{A character string specifying the LPJ-GUESS version (default is "PalEON").}
 }
@@ -17,3 +17,6 @@ A matrix or list containing the extracted data.
 \description{
 Reads a binary file formatted for LPJ-GUESS and extracts relevant data.
 }
+\author{
+Istem Fer, Yinghao Sun
+}
diff --git a/models/lpjguess/man/read_restart.LPJGUESS.Rd b/models/lpjguess/man/read_restart.LPJGUESS.Rd
new file mode 100644
index 00000000000..0bacaa4070b
--- /dev/null
+++ b/models/lpjguess/man/read_restart.LPJGUESS.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/read_restart.LPJGUESS.R
+\name{read_restart.LPJGUESS}
+\alias{read_restart.LPJGUESS}
+\title{Read Restart for LPJGUESS}
+\usage{
+read_restart.LPJGUESS(outdir, runid, stop.time, settings, var.names, params)
+}
+\arguments{
+\item{outdir}{output directory}
+
+\item{runid}{run ID}
+
+\item{stop.time}{year that is being read}
+
+\item{settings}{PEcAn settings object}
+
+\item{var.names}{var.names to be extracted}
+
+\item{params}{passed on to return value}
+}
+\value{
+X_tmp      vector of forecasts
+}
+\description{
+Read Restart for LPJGUESS
+}
+\examples{
+\dontrun{
+  rx <- read_restart.LPJGUESS(
+           outdir   = "/projectnb/…/LPJ_output",
+           runid    = "123456",
+           stop.time = as.POSIXct("2001-12-31 23:59:59", tz = "UTC"),
+           settings = settings,
+           var.names = c("AGB.pft"),
+           params = params)
+}
+}
+\author{
+Istem Fer, Yinghao Sun
+}
diff --git a/models/lpjguess/man/serialize_starts_ends.Rd b/models/lpjguess/man/serialize_starts_ends.Rd
index 5c743458b55..3b3bfa242fd 100644
--- a/models/lpjguess/man/serialize_starts_ends.Rd
+++ b/models/lpjguess/man/serialize_starts_ends.Rd
@@ -17,3 +17,4 @@ A numeric vector of length 2, giving the start and end line numbers.
 \description{
 Finds the start and end lines for serialization.
 }
+\keyword{internal}
diff --git a/models/lpjguess/man/write_binary_LPJGUESS.Rd b/models/lpjguess/man/write_binary_LPJGUESS.Rd
new file mode 100644
index 00000000000..e6718a5155e
--- /dev/null
+++ b/models/lpjguess/man/write_binary_LPJGUESS.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/write_state.R
+\name{write_binary_LPJGUESS}
+\alias{write_binary_LPJGUESS}
+\title{Write updated variables into a copy of the original LPJ-GUESS .state file}
+\usage{
+write_binary_LPJGUESS(State_updated, outdir)
+}
+\arguments{
+\item{State_updated}{A list containing updated state variables, position list and size list (get from read_binary)}
+
+\item{outdir}{Path to a directory containing the `0.state` and `meta.bin` files.}
+}
+\value{
+No return value. Writes files to disk as side effect.
+}
+\description{
+Write updated variables into a copy of the original LPJ-GUESS .state file
+}
+\author{
+Yinghao Sun
+}
diff --git a/models/lpjguess/man/write_restart.LPJGUESS.Rd b/models/lpjguess/man/write_restart.LPJGUESS.Rd
new file mode 100644
index 00000000000..2b06018b468
--- /dev/null
+++ b/models/lpjguess/man/write_restart.LPJGUESS.Rd
@@ -0,0 +1,51 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/write_restart.LPJGUESS.R
+\name{write_restart.LPJGUESS}
+\alias{write_restart.LPJGUESS}
+\title{write_restart.LPJGUESS}
+\usage{
+write_restart.LPJGUESS(
+  outdir,
+  runid,
+  start.time,
+  stop.time,
+  settings,
+  new.state,
+  RENAME = TRUE,
+  new.params,
+  inputs = NULL,
+  verbose = FALSE
+)
+}
+\arguments{
+\item{outdir}{output directory}
+
+\item{runid}{run ID}
+
+\item{start.time}{start date and time for each SDA ensemble}
+
+\item{stop.time}{stop date and time for each SDA ensemble}
+
+\item{settings}{PEcAn settings object}
+
+\item{new.state}{analysis state vector}
+
+\item{RENAME}{flag to either rename output file or not}
+
+\item{new.params}{list of parameters to convert between different states}
+
+\item{inputs}{list of model inputs to use in write.configs.SIPNET}
+
+\item{verbose}{decide if we want to print the runid}
+}
+\value{
+NONE
+}
+\description{
+Write restart files for LPJGUESS
+new.state includes X (AGB.pft) from Analysis
+new.params includes LPJGUESS_state
+}
+\author{
+Yinghao Sun
+}