From 2a281c6de81517dbb2bdfff62adaffadc076ab31 Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@MacBook-Air.local>
Date: Wed, 13 May 2026 19:58:16 -0700
Subject: [PATCH 1/2] Add standalone restart support

Write standalone restart files for wetting fronts, non-vadose state, and GIUH queue. Load restart state from the last data row of configured restart files. Save wetting front state in model-native units with round-trip double precision, preserving depth, theta, layer/front numbers, to_bottom, psi, and dzdt for each wetting front. Restore conceptual reservoir storage, ponded/surface memory, flux cache state, and GIUH queue. Ignore restart files for invalid soil types and ignore saved cache fields when flux caching is disabled. Document restart config keys, saved variables, standalone-only writing, and restart behavior.
---
 CMakeLists.txt                                |   6 +-
 configs/README.md                             |  47 +-
 configs/config_lasam_Bushland.txt             |  11 +-
 ...nfig_lasam_Phillipsburg_with_reservoir.txt |   3 +
 include/all.hxx                               |  23 +
 include/bmi_lgar.hxx                          |   2 +
 src/bmi_lgar.cxx                              |  46 +-
 src/bmi_main_lgar.cxx                         |  75 +-
 src/init_from_state.cxx                       | 672 ++++++++++++++++++
 src/lgar.cxx                                  | 272 +++++--
 10 files changed, 1104 insertions(+), 53 deletions(-)
 create mode 100644 src/init_from_state.cxx

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 226881f..739d7c2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -20,13 +20,13 @@ message(CMAKE_C_COMPILER " ${CMAKE_C_COMPILER}")
 message("CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}")
 
 # standalone
-add_executable(lasam_standalone ./src/bmi_main_lgar.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx ./src/conceptual_reservoir.cxx
+add_executable(lasam_standalone ./src/bmi_main_lgar.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx ./src/conceptual_reservoir.cxx ./src/init_from_state.cxx
              ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/util_funcs.cxx ./src/aet.cxx
              ./giuh/giuh.h ./giuh/giuh.c)
 target_link_libraries(lasam_standalone PRIVATE m)
 
 # unittest
-add_executable(lasam_unitest ./tests/main_unit_test_bmi.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx ./src/conceptual_reservoir.cxx
+add_executable(lasam_unitest ./tests/main_unit_test_bmi.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx ./src/conceptual_reservoir.cxx ./src/init_from_state.cxx
              ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/util_funcs.cxx ./src/aet.cxx ./giuh/giuh.h
              ./giuh/giuh.c)
 target_link_libraries(lasam_unitest PRIVATE m)
@@ -34,7 +34,7 @@ target_link_libraries(lasam_unitest PRIVATE m)
 # Make sure these are compiled with this directive
 add_compile_definitions(BMI_ACTIVE)
 
-add_library(lasambmi SHARED src/bmi_lgar.cxx src/lgar.cxx ./src/soil_funcs.cxx ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/conceptual_reservoir.cxx
+add_library(lasambmi SHARED src/bmi_lgar.cxx src/lgar.cxx ./src/soil_funcs.cxx ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/conceptual_reservoir.cxx ./src/init_from_state.cxx
         ./src/util_funcs.cxx ./src/aet.cxx ./giuh/giuh.c include/all.hxx ./giuh/giuh.h)
 
 target_compile_definitions(lasambmi PRIVATE NGEN)
diff --git a/configs/README.md b/configs/README.md
index 5bdae7c..69934af 100644
--- a/configs/README.md
+++ b/configs/README.md
@@ -35,6 +35,51 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | spf_factor | double (scalar) | 0.1 <= spf_factor <= 1 | - | parameter for fluxes to nonlinear reservoir | storage that contributes directly to streamflow | Simple bypass of surface water to the nonlinear reservoir will occur when the most superficial wetting front achieves the theta_e value of its layer times spf_factor. When this occurs, the amount of water sent to the nonlinear reservoir is equal to the precipitation plus any ponded water times frac_to_CR. This is a rather simple representation of preferential flow that intends to simulate the episodic nature of streamflow events in arid or semi arid environments. Defaults to 0.98. |
 | allow_flux_caching | Boolean | true, false | - | trades a small amount of accuracy for a lot of speed | flux caching | During dry periods, it is often the case that wetting fronts will move very slowly and AET will be significantly less than PET. In these cases, in the context of streamflow simulation, it is not efficient to recompute fluxes and soil moisture dynamics for each time step. If this is set to true, then fluxes and wetting front movement will only be recomputed once every 24 hours, or when the conditions resulting in dry and slow wetting fronts and low AET cease. During the times for which fluxes are not recomputed, instead they are stored in a cache and fluxes for subsequent time steps are set using this cache. Sligtly different strategies are used for fluxes through the lower boundary and AET. Also note that because NextGen models should ideally provide output for each hour, simply setting an adaptive time step to be larger than one hour is not a preferred runtime reduction method here. Note that this can cause small mass balance errors when the lower boundary condition is set to free drainage. Defaults to false. |
 | log_mode | Boolean | true, false | - | helps calibration search space exploration | log transform of parameters | When this is set to true, then all inputs for the van Genuchten parameter alpha, saturated hydraulic conductivity, and the nonlinear reservoir parameter a must be input as their log values rather than the normal values. For example, if an saturated hydraulic conductivity of 0.1 cm/h is desired, then the input value must be -1 because 10^-1 = 0.1. The reasoning for this is that these parameters are not distributed normally in nature but rather are distributed log normally, such that simply sampling the parameter space normally during calibration will vastly undersample a big region of the parameter space in which we expect useful parameter sets to be. Defaults to false. |
+| init_state_path | string | - | - | filename | optional initialization | If set, CASAM will initialize by loading the wetting fronts from the last data row in this file. The expected format is the standalone-mode `data_layers.csv` restart file. |
+| init_non_vadose_state_path | string | - | - | filename | optional initialization | If set, CASAM will initialize non-vadose restart variables from the last data row in this file, including conceptual reservoir storage, `volon_timestep_cm`, `runoff_in_prev_step`, `precip_previous_timestep_cm`, and flux-caching state. The expected format is the standalone-mode `data_non_vadose_state.csv` restart file. If `allow_flux_caching=false`, saved flux-cache fields are ignored. |
+| init_giuh_state_path | string | - | - | filename | optional initialization | If set, CASAM will initialize the saved GIUH runoff queue from the last data row in this file. The expected format is the standalone-mode `data_giuh_state.csv` restart file. |
 | a_slow | double (scalar) | 1E-8 < a_slow < 1E-1 | cm^(1-b) h^-1 | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This is exactly like the parameter a, except it corresponds to a second nonlinear reservoir, which was added to simulate cases where receding limbs have behaviors that can not easily be captured by one reservoir. Defaults to 0.|
 | b_slow | double (scalar) | 0.01 < b_slow < 5 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This is exactly like the parameter b, except it corresponds to a second nonlinear reservoir, which was added to simulate cases where receding limbs have behaviors that can not easily be captured by one reservoir. Defaults to 0.|
-| frac_slow | double (scalar) | 0.0 < frac_slow <= 1 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This describes the partitioning of water to the two reservoris, where the the input to the slow reservoir is equal to the total input for the nonlinear reservoirs times frac_slow. Note that either all or none of a_slow, b_slow, and frac_slow must be specified. If none are specified then the model will not simulate a second nonlinear reservoir. Defaults to 0.|
\ No newline at end of file
+| frac_slow | double (scalar) | 0.0 < frac_slow <= 1 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This describes the partitioning of water to the two reservoris, where the the input to the slow reservoir is equal to the total input for the nonlinear reservoirs times frac_slow. Note that either all or none of a_slow, b_slow, and frac_slow must be specified. If none are specified then the model will not simulate a second nonlinear reservoir. Defaults to 0.|
+
+### Restart Files
+
+Restart files are currently written only by standalone mode. Restart initialization reads the last data row from each requested restart file. If an invalid soil type is specified, restart files are ignored because CASAM returns input precipitation as output discharge for invalid soil types.
+
+The standalone `data_layers.csv` restart file saves and loads the following wetting front variables:
+
+| Variable | Units | Description |
+| -------- | ----- | ----------- |
+| depth_cm | cm | wetting front depth |
+| theta | - | volumetric water content |
+| layer_num | - | 1-based soil layer number |
+| front_num | - | 1-based wetting front number |
+| to_bottom | - | whether the wetting front is in contact with a layer bottom |
+| psi_cm | cm | capillary pressure head |
+| dzdt_cm_per_h | cm/h | wetting front velocity |
+
+The standalone `data_non_vadose_state.csv` restart file saves and loads the following non-vadose and conceptual reservoir variables:
+
+| Variable | Units | Description |
+| -------- | ----- | ----------- |
+| CR_fast_storage_cm | cm | fast conceptual reservoir storage |
+| CR_slow_storage_cm | cm | slow conceptual reservoir storage |
+| volon_timestep_cm | cm | surface water available at the previous timestep |
+| runoff_in_prev_step | - | whether runoff occurred in the previous timestep |
+| precip_previous_timestep_cm | cm | precipitation from the previous timestep |
+| cache_fluxes | - | whether cached fluxes were active |
+| cache_count | - | number of timesteps using cached fluxes |
+| previous_AET | cm | cached actual evapotranspiration |
+| previous_PET | cm | cached potential evapotranspiration |
+| previous_recharge | cm | cached recharge |
+| accumulated_PET | cm | accumulated PET while using cached fluxes |
+| accumulated_free_drainage | cm | accumulated free drainage while using cached fluxes |
+
+If `allow_flux_caching=false`, saved values for `cache_fluxes`, `cache_count`, `previous_AET`, `previous_PET`, `previous_recharge`, `accumulated_PET`, and `accumulated_free_drainage` are ignored and initialized as inactive cache state.
+
+The standalone `data_giuh_state.csv` restart file saves and loads the following GIUH variables:
+
+| Variable | Units | Description |
+| -------- | ----- | ----------- |
+| num_giuh_ordinates | - | number of GIUH ordinates expected by the saved queue |
+| queue | cm | GIUH runoff queue values from index 0 through `num_giuh_ordinates` |
diff --git a/configs/config_lasam_Bushland.txt b/configs/config_lasam_Bushland.txt
index 5d2e89d..d2f3844 100644
--- a/configs/config_lasam_Bushland.txt
+++ b/configs/config_lasam_Bushland.txt
@@ -1,13 +1,13 @@
-verbosity=none
+verbosity=high
 forcing_file=./forcing/forcing_data_resampled_uniform_Bushland.csv
 soil_params_file=./data/vG_default_params.dat
 layer_thickness=18.0,76.0,135.0[cm]
 initial_psi=2000.0[cm]
-timestep=3600[sec]
+timestep=300[sec]
 endtime=8333.0[hr]
 forcing_resolution=3600[sec]
 ponded_depth_max=0[cm]
-use_closed_form_G=true
+use_closed_form_G=false
 layer_soil_type=16,17,18
 max_valid_soil_types=25
 wilting_point_psi=15495.0[cm]
@@ -18,4 +18,7 @@ a=0.001
 b=3.0
 frac_to_CR=0.1
 spf_factor=0.9
-allow_flux_caching=true
\ No newline at end of file
+allow_flux_caching=true
+init_state_path=./data_layers.csv
+init_non_vadose_state_path=./data_non_vadose_state.csv
+init_giuh_state_path=./data_giuh_state.csv
diff --git a/configs/config_lasam_Phillipsburg_with_reservoir.txt b/configs/config_lasam_Phillipsburg_with_reservoir.txt
index 83251a0..a91ed3e 100644
--- a/configs/config_lasam_Phillipsburg_with_reservoir.txt
+++ b/configs/config_lasam_Phillipsburg_with_reservoir.txt
@@ -19,3 +19,6 @@ a=0.001
 b=2.5
 spf_factor=0.9
 frac_to_CR=0.2
+#init_state_path=./data_layers.csv
+#init_non_vadose_state_path=./data_non_vadose_state.csv
+#init_giuh_state_path=./data_giuh_state.csv
diff --git a/include/all.hxx b/include/all.hxx
index aacbf6e..a81ab56 100755
--- a/include/all.hxx
+++ b/include/all.hxx
@@ -132,6 +132,9 @@ struct lgar_bmi_parameters
   double  field_capacity_psi_cm;          // field capacity represented as a capillary head. Note that both wilting point and field capacity are specified for the whole model domain with single values
   bool   use_closed_form_G = false;      /* true if closed form of capillary drive calculation is desired, false if numeric integral
 					    for capillary drive calculation is desired */
+  string init_state_path;                // optional path for saved wetting front state
+  string init_non_vadose_state_path;      // optional path for conceptual reservoir / restart state
+  string init_giuh_state_path;            // optional path for saved GIUH queue state
   bool   PET_affects_precip = false;     // set to true in config file if you want PET to be taken from precip 
   bool   adaptive_timestep = false;      // if set to true, model uses adaptive timestep. In this case, the minimum timestep is the timestep specified in the config file. The maximum time step will be equal to the forcing resolution.
   bool   free_drainage_enabled = false;  // free_drainage_enabled will specify whether the lower boundary condition is no flow (false), or free drainage (true). Defaults to false.
@@ -403,6 +406,21 @@ extern void InitFromConfigFile(string config_file, struct model_state *state);
 extern vector<double> ReadVectorData(string key);
 extern void InitializeWettingFronts(bool is_invalid_soil_type, int num_layers, double initial_psi_cm, int *layer_soil_type, double *cum_layer_thickness_cm,
 				    double *frozen_factor, struct wetting_front** head, struct soil_properties_ *soil_properties);
+extern void InitializeWettingFrontsFromCSV(
+    int num_layers,
+    const char *data_layers_csv_path,
+    int *layer_soil_type,
+    double *cum_layer_thickness_cm,
+    double *frozen_factor,
+    struct wetting_front **head,
+    struct soil_properties_ *soil_properties);
+extern void InitializenonvadoseStateFromCSV(
+    const char *non_vadose_state_csv_path,
+    struct model_state *state);
+extern void InitializeGIUHRunoffQueueFromCSV(
+    const char *giuh_state_csv_path,
+    double *giuh_runoff_queue,
+    int num_giuh_ordinates);
 
 /********************************************************************/
 /*Other function prototypes for doing hydrology calculations, etc.  */
@@ -423,6 +441,11 @@ extern void lgar_global_mass_balance(struct model_state *state, double *giuh_run
 
 // writes full state of wetting fronts (depth, theta, no. of wetting front, no. of layer, dz/dt, psi) to a file at each time step
 extern void write_state(FILE *out, struct wetting_front* head);
+extern void write_non_vadose_state(FILE *out, struct model_state *state);
+extern void write_giuh_runoff_queue_state(
+    FILE *out,
+    const double *giuh_runoff_queue,
+    int num_giuh_ordinates);
 
 
 /********************************************************************/
diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx
index 4136791..950fb1c 100644
--- a/include/bmi_lgar.hxx
+++ b/include/bmi_lgar.hxx
@@ -147,6 +147,8 @@ public:
   void global_mass_balance();
   double update_calibratable_parameters();
   struct model_state* get_model();
+  double* get_giuh_runoff_queue();
+  int get_num_giuh_ordinates();
   
 private:
   void realloc_soil();
diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index e47930b..d73530f 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -9,6 +9,7 @@
 #include <vector>
 #include <algorithm>
 #include <iostream>
+#include <iomanip>
 #include "../bmi/bmi.hxx"
 #include "../include/bmi_lgar.hxx"
 #include "../include/all.hxx"
@@ -106,8 +107,35 @@ Initialize (std::string config_file)
     giuh_ordinates[i] = state->lgar_bmi_params.giuh_ordinates[i+1]; // note lgar uses 1-indexing
   }
 
-  for (int i=0; i<=num_giuh_ordinates;i++){
-    giuh_runoff_queue[i] = 0.0;
+  if (!state->lgar_bmi_params.is_invalid_soil_type &&
+      !state->lgar_bmi_params.init_giuh_state_path.empty()) {
+    InitializeGIUHRunoffQueueFromCSV(
+        state->lgar_bmi_params.init_giuh_state_path.c_str(),
+        giuh_runoff_queue,
+        num_giuh_ordinates);
+
+    if (verbosity.compare("high") == 0) {
+      streamsize old_precision = std::cerr.precision();
+      ios::fmtflags old_flags = std::cerr.flags();
+
+      std::cerr << "GIUH runoff queue initialized from saved state file: "
+                << state->lgar_bmi_params.init_giuh_state_path << "\n";
+      std::cerr << "GIUH runoff queue values: [";
+      std::cerr << std::fixed << std::setprecision(12);
+      for (int i = 0; i <= num_giuh_ordinates; i++) {
+        std::cerr << giuh_runoff_queue[i];
+        if (i < num_giuh_ordinates) std::cerr << ",";
+      }
+      std::cerr << "]\n";
+      std::cerr.flags(old_flags);
+      std::cerr.precision(old_precision);
+      std::cerr << "          *****         \n";
+    }
+  }
+  else {
+    for (int i=0; i<=num_giuh_ordinates;i++){
+      giuh_runoff_queue[i] = 0.0;
+    }
   }
 
 }
@@ -269,6 +297,10 @@ Update()
     state->lgar_bmi_params.timestep_h = subtimestep_h;
   }
 
+  if (!state->lgar_bmi_params.allow_flux_caching) {
+    state->lgar_mass_balance.cache_fluxes = FALSE;
+  }
+
   bool caching_at_start = state->lgar_mass_balance.cache_fluxes;
   bool switch_caching = FALSE;
 
@@ -1663,4 +1695,14 @@ GetGridNodesPerFace(const int grid, int *nodes_per_face)
   throw bmi_lgar::NotImplemented();
 }
 
+double* BmiLGAR::get_giuh_runoff_queue()
+{
+  return giuh_runoff_queue;
+}
+
+int BmiLGAR::get_num_giuh_ordinates()
+{
+  return num_giuh_ordinates;
+}
+
 #endif
diff --git a/src/bmi_main_lgar.cxx b/src/bmi_main_lgar.cxx
index 01c28b6..75bb63c 100644
--- a/src/bmi_main_lgar.cxx
+++ b/src/bmi_main_lgar.cxx
@@ -94,13 +94,22 @@ int main(int argc, char *argv[])
 
   FILE *outdata_fptr = NULL;
   FILE *outlayer_fptr = NULL;
+  FILE *outnonvadose_fptr = NULL;
+  FILE *outgiuh_fptr = NULL;
 
   if (!is_IO_supress) {
     outdata_fptr = fopen("data_variables.csv", "w");  // write output variables (e.g. infiltration, storage etc.) to this file pointer
     outlayer_fptr = fopen("data_layers.csv", "w");    // write output layers to this file pointer
+    outnonvadose_fptr = fopen("data_non_vadose_state.csv", "w");
+    outgiuh_fptr = fopen("data_giuh_state.csv", "w");
 
     // write heading (variable names)
     fprintf(outdata_fptr,"Time,");
+    fprintf(outnonvadose_fptr,
+        "Time,CR_fast_storage_cm,CR_slow_storage_cm,volon_timestep_cm,"
+        "runoff_in_prev_step,precip_previous_timestep_cm,cache_fluxes,"
+        "cache_count,previous_AET,previous_PET,previous_recharge,"
+        "accumulated_PET,accumulated_free_drainage\n");
     for (int j = 0; j < num_output_var; j++) {
       fprintf(outdata_fptr,"%s",output_var_names[j].c_str());
       if (j == num_output_var-1)
@@ -144,13 +153,21 @@ int main(int argc, char *argv[])
 	std::string name = output_var_names[j];
 	double value = 0.0;
 	model_state.GetValue(name,&value);
-	fprintf(outdata_fptr,"%6.15f",value);
+	fprintf(outdata_fptr,"%.12f",value);
 	if (j == num_output_var-1)
 	  fprintf(outdata_fptr,"\n");
 	else
 	  fprintf(outdata_fptr,",");
       }
 
+      fprintf(outnonvadose_fptr, "%s,", time[i].c_str());
+      write_non_vadose_state(outnonvadose_fptr, model_state.get_model());
+
+      fprintf(outgiuh_fptr, "%s,", time[i].c_str());
+      write_giuh_runoff_queue_state(
+          outgiuh_fptr,
+          model_state.get_giuh_runoff_queue(),
+          model_state.get_num_giuh_ordinates());
 
       // write layers data to file
       fprintf(outlayer_fptr,"# Timestep = %d, %s \n", i, time[i].c_str());
@@ -166,6 +183,8 @@ int main(int argc, char *argv[])
   if (outdata_fptr) {
     fclose(outdata_fptr);
     fclose(outlayer_fptr);
+    fclose(outnonvadose_fptr);
+    fclose(outgiuh_fptr);
   }
 
   end_time = clock();
@@ -270,11 +289,61 @@ extern void write_state(FILE *out, struct wetting_front* head){
   while(current != NULL)
   {
     if (current == head)
-      fprintf(out,"(%lf,%lf,%d,%d,%lf)",current->depth_cm*10., current->theta, current->layer_num,current->front_num, current->psi_cm*10.);
+      fprintf(out,"(%.17g,%.17g,%d,%d,%d,%.17g,%.17g)",
+              current->depth_cm, current->theta, current->layer_num,
+              current->front_num, current->to_bottom ? 1 : 0,
+              current->psi_cm, current->dzdt_cm_per_h);
     else
-      fprintf(out,"|(%lf,%lf,%d,%d,%lf)",current->depth_cm*10., current->theta, current->layer_num,current->front_num, current->psi_cm*10.);
+      fprintf(out,"|(%.17g,%.17g,%d,%d,%d,%.17g,%.17g)",
+              current->depth_cm, current->theta, current->layer_num,
+              current->front_num, current->to_bottom ? 1 : 0,
+              current->psi_cm, current->dzdt_cm_per_h);
   current = current->next;
   }
   fprintf(out, "]\n");
 
 }
+
+extern void write_non_vadose_state(FILE *out, struct model_state* state)
+{
+  fprintf(out,
+          "CR_fast_storage_cm=%.17g,"
+          "CR_slow_storage_cm=%.17g,"
+          "volon_timestep_cm=%.17g,"
+          "runoff_in_prev_step=%d,"
+          "precip_previous_timestep_cm=%.17g,"
+          "cache_fluxes=%d,"
+          "cache_count=%d,"
+          "previous_AET=%.17g,"
+          "previous_PET=%.17g,"
+          "previous_recharge=%.17g,"
+          "accumulated_PET=%.17g,"
+          "accumulated_free_drainage=%.17g\n",
+          state->lgar_mass_balance.CR_fast_storage_cm,
+          state->lgar_mass_balance.CR_slow_storage_cm,
+          state->lgar_mass_balance.volon_timestep_cm,
+          state->lgar_bmi_params.runoff_in_prev_step ? 1 : 0,
+          state->lgar_bmi_params.precip_previous_timestep_cm,
+          state->lgar_mass_balance.cache_fluxes ? 1 : 0,
+          state->lgar_bmi_params.cache_count,
+          state->lgar_mass_balance.previous_AET,
+          state->lgar_mass_balance.previous_PET,
+          state->lgar_mass_balance.previous_recharge,
+          state->lgar_mass_balance.accumulated_PET,
+          state->lgar_mass_balance.accumulated_free_drainage);
+}
+
+extern void write_giuh_runoff_queue_state(
+    FILE *out,
+    const double *giuh_runoff_queue,
+    int num_giuh_ordinates)
+{
+  fprintf(out, "num_giuh_ordinates=%d,queue=[", num_giuh_ordinates);
+
+  for (int i = 0; i <= num_giuh_ordinates; i++) {
+    if (i > 0) fprintf(out, ",");
+    fprintf(out, "%.17g", giuh_runoff_queue[i]);
+  }
+
+  fprintf(out, "]\n");
+}
diff --git a/src/init_from_state.cxx b/src/init_from_state.cxx
new file mode 100644
index 0000000..0aba673
--- /dev/null
+++ b/src/init_from_state.cxx
@@ -0,0 +1,672 @@
+#include "../include/all.hxx"
+
+#include <cstddef>
+#include <iostream>
+#include <sstream>
+#include <string.h>
+
+using namespace std;
+
+/*
+  Restart-state readers for standalone LASAM output.
+
+  The standalone driver writes one restart row per model output time.  These
+  readers intentionally use the last non-header data row from each file, which
+  lets a later standalone run restart from the final state saved by an earlier
+  run.  Values are read in the model's natural units, which are cm for depths,
+  pressure heads, storage, and flux-like accumulated depths.
+
+  The public functions in this file exit with an error message if a requested
+  restart file is missing or malformed.  That is deliberate: a partial or
+  silently misread restart state would be much harder to diagnose than a hard
+  initialization failure.
+*/
+
+/* One parsed wetting front tuple from data_layers.csv. */
+typedef struct {
+  double depth_cm;        // Wetting front depth in model-native cm.
+  double theta;           // Volumetric water content for the front.
+  int layer_num;          // 1-based soil layer containing this front.
+  int front_num;          // 1-based wetting front ordering index.
+  bool to_bottom;         // Whether this front is in contact with a layer bottom.
+  double psi_cm;          // Capillary pressure head in model-native cm.
+  double dzdt_cm_per_h;   // Saved front velocity in cm/h, restored for continuity.
+} WFRecord;
+
+/* qsort callback used to rebuild the linked list in front_num order. */
+static int cmp_record_by_front_num(const void *a, const void *b)
+{
+  const WFRecord *ra = (const WFRecord*)a;
+  const WFRecord *rb = (const WFRecord*)b;
+  return (ra->front_num - rb->front_num);
+}
+
+/*
+  Return true for rows that can contain restart data.
+
+  The standalone files include human-readable headers and, for wetting fronts,
+  timestep comment lines.  The restart reader skips those rows and keeps only
+  actual saved-state records.
+*/
+static bool is_restart_data_line(const char *buf)
+{
+  const char *p = buf;
+  while (*p == ' ' || *p == '\t' || *p == '\r' || *p == '\n') p++;
+
+  if (*p == '\0' || *p == '#') return false;
+  if (strncmp(p, "Time,", 5) == 0) return false;
+
+  return true;
+}
+
+/*
+  Read through the whole file and keep the final restart data row.
+
+  This is intentionally not "first data row" behavior: data_layers.csv,
+  data_non_vadose_state.csv, and data_giuh_state.csv are time series of saved
+  states, so the final valid row is the state at the end of the previous run.
+*/
+static bool read_last_data_line(FILE *fp, char *buf, size_t buflen)
+{
+  char line[65536];
+  bool found = false;
+
+  while (fgets(line, sizeof(line), fp)) {
+    if (!is_restart_data_line(line)) continue;
+
+    snprintf(buf, buflen, "%s", line);
+    found = true;
+  }
+
+  return found;
+}
+
+/*
+  Parses state lines like:
+    [(18,0.197519,1,1,1,2000,1e-05)|(...)]
+
+  Each tuple is:
+    (depth_cm,theta,layer_num,front_num,to_bottom,psi_cm,dzdt_cm_per_h)
+
+  The parser only translates text into WFRecord values.  Physical consistency
+  checks, sorting, linked-list insertion, and hydraulic-property reconstruction
+  happen later in InitializeWettingFrontsFromCSV().
+*/
+static int parse_state_line_to_records(const char *line_in, WFRecord **recs_out, int *n_out)
+{
+  /* Give callers predictable outputs even if parsing fails early. */
+  *recs_out = NULL;
+  *n_out = 0;
+
+  /*
+    strdup lets strtok_r and bracket replacement modify a private copy of the
+    row.  The input line may point into a caller-owned buffer.
+  */
+  char *line = strdup(line_in);
+  if (!line) return 1;
+
+  /* Wetting front records are enclosed in square brackets. */
+  char *lbr = strchr(line, '[');
+  char *rbr = strrchr(line, ']');
+  if (!lbr || !rbr || rbr <= lbr) {
+    free(line);
+    return 2;
+  }
+  *rbr = '\0';
+  char *p = lbr + 1;
+
+  /*
+    The number of wetting fronts changes during a simulation, so grow the
+    record array dynamically rather than assuming one front per soil layer.
+  */
+  int cap = 16;
+  int n = 0;
+  WFRecord *recs = (WFRecord*)malloc(sizeof(WFRecord) * cap);
+  if (!recs) {
+    free(line);
+    return 3;
+  }
+
+  char *saveptr = NULL;
+  for (char *tok = strtok_r(p, "|", &saveptr); tok != NULL; tok = strtok_r(NULL, "|", &saveptr)) {
+    /* Be tolerant of spaces and a leading comma before each tuple. */
+    while (*tok == ' ' || *tok == '\t' || *tok == ',') tok++;
+
+    double depth_cm, theta, psi_cm, dzdt_cm_per_h;
+    int layer_num, front_num, to_bottom_int;
+
+    int matched = sscanf(tok, " ( %lf , %lf , %d , %d , %d , %lf , %lf ) ",
+                         &depth_cm, &theta, &layer_num, &front_num,
+                         &to_bottom_int, &psi_cm, &dzdt_cm_per_h);
+
+    /* A restart written by the current writer must contain all seven fields. */
+    if (matched != 7) {
+      free(recs);
+      free(line);
+      return 4;
+    }
+
+    if (to_bottom_int != 0 && to_bottom_int != 1) {
+      free(recs);
+      free(line);
+      return 6;
+    }
+
+    if (n == cap) {
+      cap *= 2;
+      WFRecord *tmp = (WFRecord*)realloc(recs, sizeof(WFRecord) * cap);
+      if (!tmp) {
+        free(recs);
+        free(line);
+        return 7;
+      }
+      recs = tmp;
+    }
+
+    recs[n].depth_cm = depth_cm;
+    recs[n].theta = theta;
+    recs[n].layer_num = layer_num;
+    recs[n].front_num = front_num;
+    recs[n].to_bottom = (to_bottom_int == 1);
+    recs[n].psi_cm = psi_cm;
+    recs[n].dzdt_cm_per_h = dzdt_cm_per_h;
+    n++;
+  }
+
+  free(line);
+  *recs_out = recs;
+  *n_out = n;
+  return 0;
+}
+
+extern void InitializeWettingFrontsFromCSV(
+    int num_layers,
+    const char *data_layers_csv_path,
+    int *layer_soil_type,
+    double *cum_layer_thickness_cm,
+    double *frozen_factor,
+    struct wetting_front **head,
+    struct soil_properties_ *soil_properties)
+{
+  /*
+    Restart initialization owns the wetting front linked list from scratch.
+    If anything fails after this point, the function exits rather than leaving
+    a partly initialized model.
+  */
+  *head = NULL;
+
+  FILE *fp = fopen(data_layers_csv_path, "r");
+  if (!fp) {
+    fprintf(stderr, "ERROR: could not open data_layers file: %s\n", data_layers_csv_path);
+    exit(1);
+  }
+
+  /* Use the final saved state, not the first state in the file. */
+  char state_line[65536];
+  if (!read_last_data_line(fp, state_line, sizeof(state_line))) {
+    fclose(fp);
+    fprintf(stderr, "ERROR: no wetting front restart state lines found in %s\n",
+            data_layers_csv_path);
+    exit(1);
+  }
+  fclose(fp);
+
+  WFRecord *recs = NULL;
+  int nrecs = 0;
+  int perr = parse_state_line_to_records(state_line, &recs, &nrecs);
+  if (perr != 0 || !recs || nrecs <= 0) {
+    fprintf(stderr,
+            "ERROR: failed to parse state line in %s (err=%d). Expected "
+            "(depth_cm,theta,layer_num,front_num,to_bottom,psi_cm,dzdt_cm_per_h) tuples.\n",
+            data_layers_csv_path, perr);
+    exit(1);
+  }
+
+  /*
+    The linked list order is hydrologically meaningful.  Sort by the saved
+    front_num before rebuilding it, then check that no front numbers are
+    missing or duplicated.
+  */
+  qsort(recs, nrecs, sizeof(WFRecord), cmp_record_by_front_num);
+
+  if (recs[0].front_num != 1) {
+    fprintf(stderr, "ERROR: first front_num in file is %d, but restart requires numbering from 1\n",
+            recs[0].front_num);
+    free(recs);
+    exit(1);
+  }
+
+  for (int i = 1; i < nrecs; i++) {
+    if (recs[i].front_num != recs[i - 1].front_num + 1) {
+      fprintf(stderr, "ERROR: front_num values are not contiguous at %d -> %d\n",
+              recs[i - 1].front_num, recs[i].front_num);
+      free(recs);
+      exit(1);
+    }
+  }
+
+  for (int i = 0; i < nrecs; i++) {
+    WFRecord r = recs[i];
+
+    /* layer_num is 1-based everywhere in this model. */
+    if (r.layer_num < 1 || r.layer_num > num_layers) {
+      fprintf(stderr, "ERROR: record layer_num=%d out of range 1..%d\n", r.layer_num, num_layers);
+      free(recs);
+      exit(1);
+    }
+
+    /*
+      A front must fall inside the layer it says it belongs to.  This also
+      catches old restart files written in mm/x10 units before this reader and
+      writer were changed to use cm directly.
+    */
+    double layer_top_cm = cum_layer_thickness_cm[r.layer_num - 1];
+    double layer_bottom_cm = cum_layer_thickness_cm[r.layer_num];
+    if (r.depth_cm <= 0.0 ||
+        r.depth_cm < layer_top_cm - 1.0e-8 ||
+        r.depth_cm > layer_bottom_cm + 1.0e-8) {
+      fprintf(stderr,
+              "ERROR: restart front_num=%d has depth_cm=%.17g outside layer %d "
+              "bounds [%.17g, %.17g] cm. Regenerate old x10/mm restart files "
+              "with the current cm-unit writer before loading.\n",
+              r.front_num, r.depth_cm, r.layer_num,
+              layer_top_cm, layer_bottom_cm);
+      free(recs);
+      exit(1);
+    }
+
+    /*
+      listInsertFront restores depth, theta, front number, layer number, and
+      to_bottom.  to_bottom is restart state and is intentionally loaded from
+      the file rather than inferred from layer ordering.
+    */
+    struct wetting_front *current =
+        listInsertFront(r.depth_cm, r.theta, r.front_num, r.layer_num, r.to_bottom, head);
+
+    if (current == NULL) {
+      fprintf(stderr, "ERROR: listInsertFront returned NULL inserting front_num=%d\n", r.front_num);
+      free(recs);
+      exit(1);
+    }
+
+    /* These state variables are not set by listInsertFront, so restore them. */
+    current->psi_cm = r.psi_cm;
+    current->dzdt_cm_per_h = r.dzdt_cm_per_h;
+
+    /*
+      Hydraulic conductivity is derived state.  Recompute it from the restored
+      theta and the soil properties for the front's current layer instead of
+      trusting an additional saved value.
+    */
+    int soil = layer_soil_type[r.layer_num];
+    double Se = calc_Se_from_theta(current->theta,
+                                   soil_properties[soil].theta_e,
+                                   soil_properties[soil].theta_r);
+    double Ksat_cm_per_h = frozen_factor[r.layer_num] * soil_properties[soil].Ksat_cm_per_h;
+    current->K_cm_per_h = calc_K_from_Se(Se, Ksat_cm_per_h, soil_properties[soil].vg_m);
+  }
+
+  free(recs);
+}
+
+/*
+  Restart variables that are not part of the wetting front linked list.
+
+  These values carry conceptual-reservoir storage, the surface-runoff memory
+  needed for wetting front creation, and flux-caching state used during dry
+  timesteps.  The writer stores them as key=value pairs so this parser can be
+  insensitive to column position after the leading time column.
+*/
+typedef struct {
+  double CR_fast_storage_cm;
+  double CR_slow_storage_cm;
+  double volon_timestep_cm;
+  bool runoff_in_prev_step;
+  double precip_previous_timestep_cm;
+  bool cache_fluxes;
+  int cache_count;
+  double previous_AET;
+  double previous_PET;
+  double previous_recharge;
+  double accumulated_PET;
+  double accumulated_free_drainage;
+} nonvadoseRestartState;
+
+/* Booleans are written as 0/1 in the restart CSV files. */
+static bool parse_bool_01(const char *s, bool *out)
+{
+  if (strcmp(s, "1") == 0) {
+    *out = true;
+    return true;
+  }
+  if (strcmp(s, "0") == 0) {
+    *out = false;
+    return true;
+  }
+  return false;
+}
+
+static int parse_non_vadose_state_kv_line(const char *line_in, nonvadoseRestartState *rst)
+{
+  /*
+    Older restart rows may not have all cache-history fields.  Defaults here
+    keep those optional cache terms cold-started unless explicit values are
+    found below.
+  */
+  rst->cache_fluxes = false;
+  rst->cache_count = 1;
+  rst->previous_AET = 0.0;
+  rst->previous_PET = 0.0;
+  rst->previous_recharge = 0.0;
+  rst->accumulated_PET = 0.0;
+  rst->accumulated_free_drainage = 0.0;
+
+  char *line = strdup(line_in);
+  if (!line) return 1;
+
+  /*
+    These fields are required because they affect immediate restart behavior:
+    reservoir storage and the previous-step surface-water state.
+  */
+  bool got_fast = false;
+  bool got_slow = false;
+  bool got_volon = false;
+  bool got_runoff = false;
+  bool got_precip_prev = false;
+
+  char *saveptr = NULL;
+  for (char *tok = strtok_r(line, ",\r\n", &saveptr);
+       tok != NULL;
+       tok = strtok_r(NULL, ",\r\n", &saveptr)) {
+    while (*tok == ' ' || *tok == '\t') tok++;
+
+    /* Each non-vadose state token is key=value. */
+    char *eq = strchr(tok, '=');
+    if (!eq) {
+      free(line);
+      return 2;
+    }
+
+    *eq = '\0';
+    const char *key = tok;
+    const char *val = eq + 1;
+
+    /*
+      Unknown keys are ignored.  That keeps the reader forward-tolerant if a
+      future writer adds more restart metadata that older code can safely skip.
+    */
+    if (strcmp(key, "CR_fast_storage_cm") == 0) {
+      rst->CR_fast_storage_cm = strtod(val, NULL);
+      got_fast = true;
+    }
+    else if (strcmp(key, "CR_slow_storage_cm") == 0) {
+      rst->CR_slow_storage_cm = strtod(val, NULL);
+      got_slow = true;
+    }
+    else if (strcmp(key, "volon_timestep_cm") == 0) {
+      rst->volon_timestep_cm = strtod(val, NULL);
+      got_volon = true;
+    }
+    else if (strcmp(key, "runoff_in_prev_step") == 0) {
+      if (!parse_bool_01(val, &rst->runoff_in_prev_step)) {
+        free(line);
+        return 3;
+      }
+      got_runoff = true;
+    }
+    else if (strcmp(key, "precip_previous_timestep_cm") == 0) {
+      rst->precip_previous_timestep_cm = strtod(val, NULL);
+      got_precip_prev = true;
+    }
+    else if (strcmp(key, "cache_fluxes") == 0) {
+      if (!parse_bool_01(val, &rst->cache_fluxes)) {
+        free(line);
+        return 5;
+      }
+    }
+    else if (strcmp(key, "cache_count") == 0) {
+      rst->cache_count = (int)strtol(val, NULL, 10);
+      if (rst->cache_count < 1) rst->cache_count = 1;
+    }
+    else if (strcmp(key, "previous_AET") == 0) {
+      rst->previous_AET = strtod(val, NULL);
+    }
+    else if (strcmp(key, "previous_PET") == 0) {
+      rst->previous_PET = strtod(val, NULL);
+    }
+    else if (strcmp(key, "previous_recharge") == 0) {
+      rst->previous_recharge = strtod(val, NULL);
+    }
+    else if (strcmp(key, "accumulated_PET") == 0) {
+      rst->accumulated_PET = strtod(val, NULL);
+    }
+    else if (strcmp(key, "accumulated_free_drainage") == 0) {
+      rst->accumulated_free_drainage = strtod(val, NULL);
+    }
+  }
+
+  free(line);
+
+  /* Refuse restart if any required non-vadose state is absent. */
+  if (!(got_fast && got_slow && got_volon && got_runoff && got_precip_prev)) {
+    return 4;
+  }
+
+  return 0;
+}
+
+extern void InitializenonvadoseStateFromCSV(
+    const char *non_vadose_state_csv_path,
+    struct model_state *state)
+{
+  FILE *fp = fopen(non_vadose_state_csv_path, "r");
+  if (!fp) {
+    fprintf(stderr, "ERROR: could not open non-vadose state file: %s\n",
+            non_vadose_state_csv_path);
+    exit(1);
+  }
+
+  /* Use the final saved non-vadose state from the previous run. */
+  char line[65536];
+  if (!read_last_data_line(fp, line, sizeof(line))) {
+    fclose(fp);
+    fprintf(stderr, "ERROR: no data lines found in non-vadose state file: %s\n",
+            non_vadose_state_csv_path);
+    exit(1);
+  }
+
+  nonvadoseRestartState rst;
+  char *parse_start = line;
+
+  /*
+    data_non_vadose_state.csv starts with a Time column.  If a comma appears
+    before the first key=value pair, skip that leading time field and parse only
+    the restart-state tokens.
+  */
+  char *comma = strchr(line, ',');
+  char *first_equal = strchr(line, '=');
+  if (comma && first_equal && comma < first_equal) parse_start = comma + 1;
+
+  int perr = parse_non_vadose_state_kv_line(parse_start, &rst);
+
+  fclose(fp);
+
+  if (perr != 0) {
+    fprintf(stderr, "ERROR: failed to parse non-vadose restart state in %s (err=%d)\n",
+            non_vadose_state_csv_path, perr);
+    exit(1);
+  }
+
+  /*
+    Restore reservoir storage and the surface-water memory used on the next
+    timestep.  lgar_initialize() will use these reservoir values to initialize
+    volCRstart/volCRend consistently for the restart mass balance.
+  */
+  state->lgar_mass_balance.CR_fast_storage_cm = rst.CR_fast_storage_cm;
+  state->lgar_mass_balance.CR_slow_storage_cm = rst.CR_slow_storage_cm;
+  state->lgar_mass_balance.volon_timestep_cm = rst.volon_timestep_cm;
+  state->lgar_bmi_params.runoff_in_prev_step = rst.runoff_in_prev_step;
+  state->lgar_bmi_params.precip_previous_timestep_cm = rst.precip_previous_timestep_cm;
+
+  /*
+    Flux-cache values are only meaningful when the current config allows flux
+    caching.  If the restart file contains cached fluxes but the current config
+    disables caching, intentionally reset those fields.
+  */
+  if (state->lgar_bmi_params.allow_flux_caching) {
+    state->lgar_mass_balance.cache_fluxes = rst.cache_fluxes;
+    state->lgar_bmi_params.cache_count = rst.cache_count;
+    state->lgar_mass_balance.previous_AET = rst.previous_AET;
+    state->lgar_mass_balance.previous_PET = rst.previous_PET;
+    state->lgar_mass_balance.previous_recharge = rst.previous_recharge;
+    state->lgar_mass_balance.accumulated_PET = rst.accumulated_PET;
+    state->lgar_mass_balance.accumulated_free_drainage = rst.accumulated_free_drainage;
+  }
+  else {
+    state->lgar_mass_balance.cache_fluxes = false;
+    state->lgar_bmi_params.cache_count = 1;
+    state->lgar_mass_balance.previous_AET = 0.0;
+    state->lgar_mass_balance.previous_PET = 0.0;
+    state->lgar_mass_balance.previous_recharge = 0.0;
+    state->lgar_mass_balance.accumulated_PET = 0.0;
+    state->lgar_mass_balance.accumulated_free_drainage = 0.0;
+  }
+
+  /*
+    These are set here for completeness; lgar_initialize() also recomputes the
+    reservoir start/end bookkeeping immediately after config initialization.
+  */
+  state->lgar_mass_balance.volCRend_cm = rst.CR_fast_storage_cm + rst.CR_slow_storage_cm;
+  state->lgar_mass_balance.volCRend_timestep_cm = rst.CR_fast_storage_cm + rst.CR_slow_storage_cm;
+}
+
+/* Parsed representation of one data_giuh_state.csv row. */
+typedef struct {
+  int num_giuh_ordinates;
+  double *queue_vals;
+} GIUHRestartState;
+
+/*
+  Parse one GIUH restart row, for example:
+    num_giuh_ordinates=5,queue=[0,0.1,0.2,0,0,0]
+
+  The queue has num_giuh_ordinates + 1 entries because the existing GIUH code
+  uses indices 0..num_giuh_ordinates.
+*/
+static int parse_giuh_state_line(const char *line_in, GIUHRestartState *rst)
+{
+  /* Initialize outputs so the caller can safely inspect/free after failure. */
+  rst->num_giuh_ordinates = -1;
+  rst->queue_vals = NULL;
+
+  char *line = strdup(line_in);
+  if (!line) return 1;
+
+  char *num_ptr = strstr(line, "num_giuh_ordinates=");
+  char *queue_ptr = strstr(line, "queue=[");
+  if (!num_ptr || !queue_ptr) {
+    free(line);
+    return 2;
+  }
+
+  if (sscanf(num_ptr, "num_giuh_ordinates=%d", &rst->num_giuh_ordinates) != 1) {
+    free(line);
+    return 3;
+  }
+
+  /* Queue values are enclosed in square brackets and separated by commas. */
+  char *lbr = strchr(queue_ptr, '[');
+  char *rbr = strrchr(queue_ptr, ']');
+  if (!lbr || !rbr || rbr <= lbr) {
+    free(line);
+    return 4;
+  }
+
+  *rbr = '\0';
+  char *vals = lbr + 1;
+
+  int expected_n = rst->num_giuh_ordinates + 1;
+  rst->queue_vals = (double*)malloc(sizeof(double) * expected_n);
+  if (!rst->queue_vals) {
+    free(line);
+    return 5;
+  }
+
+  int nread = 0;
+  char *saveptr = NULL;
+  for (char *tok = strtok_r(vals, ",", &saveptr);
+       tok != NULL;
+       tok = strtok_r(NULL, ",", &saveptr)) {
+    while (*tok == ' ' || *tok == '\t') tok++;
+
+    /* More values than expected means this file does not match the config. */
+    if (nread >= expected_n) {
+      free(rst->queue_vals);
+      free(line);
+      return 6;
+    }
+
+    rst->queue_vals[nread++] = strtod(tok, NULL);
+  }
+
+  free(line);
+
+  /* Fewer values than expected is also a malformed restart row. */
+  if (nread != expected_n) {
+    free(rst->queue_vals);
+    rst->queue_vals = NULL;
+    return 7;
+  }
+
+  return 0;
+}
+
+extern void InitializeGIUHRunoffQueueFromCSV(
+    const char *giuh_state_csv_path,
+    double *giuh_runoff_queue,
+    int num_giuh_ordinates)
+{
+  FILE *fp = fopen(giuh_state_csv_path, "r");
+  if (!fp) {
+    fprintf(stderr, "ERROR: could not open GIUH state file: %s\n", giuh_state_csv_path);
+    exit(1);
+  }
+
+  /* Use the final saved GIUH queue from the previous run. */
+  char line[65536];
+  if (!read_last_data_line(fp, line, sizeof(line))) {
+    fclose(fp);
+    fprintf(stderr, "ERROR: no data lines found in GIUH state file: %s\n", giuh_state_csv_path);
+    exit(1);
+  }
+
+  GIUHRestartState rst;
+  int perr = parse_giuh_state_line(line, &rst);
+
+  fclose(fp);
+
+  if (perr != 0) {
+    fprintf(stderr, "ERROR: failed to parse GIUH restart state in %s (err=%d)\n",
+            giuh_state_csv_path, perr);
+    exit(1);
+  }
+
+  /*
+    The saved queue length must match the current GIUH configuration.  Loading
+    a queue from a different GIUH shape would shift runoff timing.
+  */
+  if (rst.num_giuh_ordinates != num_giuh_ordinates) {
+    fprintf(stderr,
+            "ERROR: GIUH restart file has num_giuh_ordinates=%d but model expects %d\n",
+            rst.num_giuh_ordinates, num_giuh_ordinates);
+    free(rst.queue_vals);
+    exit(1);
+  }
+
+  /* Copy the parsed queue into the BMI-owned queue array. */
+  for (int i = 0; i <= num_giuh_ordinates; i++) {
+    giuh_runoff_queue[i] = rst.queue_vals[i];
+  }
+
+  free(rst.queue_vals);
+}
diff --git a/src/lgar.cxx b/src/lgar.cxx
index 0e8a966..fe910fa 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -6,6 +6,7 @@
 #include <fstream>
 #include <string.h>
 #include <sstream>
+#include <iomanip>
 
 using namespace std;
 
@@ -69,6 +70,87 @@ using namespace std;
 #define DEPTH_AVOIDS_SAME_WF_DEPTH 1.E-6       // in the event that multiple WFs all would cross a layer boundary and would each have their depth in the new layer limited by FACTOR_LIMITS_LAYER_CROSSING_SPEED, this just prevents these WFs from being exactly at the same depth.
 #define PSI_UPPER_LIM 1.E7                     // in loops that close the mass balance by iterating theta and psi, we impose an upper limit on capillary head because some values are just not physically realistic
 
+static void PrintSavedStateInitializationDetails(struct model_state *state,
+    bool is_state_path_set,
+    bool is_non_vadose_state_path_set,
+    bool is_giuh_state_path_set)
+{
+  if (verbosity.compare("high") != 0) return;
+  if (!(is_state_path_set || is_non_vadose_state_path_set || is_giuh_state_path_set)) return;
+
+  streamsize old_precision = std::cerr.precision();
+  ios::fmtflags old_flags = std::cerr.flags();
+
+  std::cerr << "--- Saved-state initialization details --- \n";
+
+  if (state->lgar_bmi_params.is_invalid_soil_type) {
+    std::cerr << "Restart files were provided but are ignored because an "
+              << "invalid soil type makes the model return input precipitation "
+              << "as output Qout. \n";
+    std::cerr << "No wetting fronts, non-vadose state, or GIUH queue were "
+              << "initialized from saved state. \n";
+    std::cerr.flags(old_flags);
+    std::cerr.precision(old_precision);
+    std::cerr << "          *****         \n";
+    return;
+  }
+
+  std::cerr << "Restart files are read from the last data row. \n";
+
+  if (is_state_path_set) {
+    std::cerr << "Wetting fronts initialized from: "
+              << state->lgar_bmi_params.init_state_path << "\n";
+    std::cerr << "No. of wetting fronts initialized: "
+              << state->lgar_bmi_params.num_wetting_fronts << "\n";
+    std::cerr << std::fixed << std::setprecision(12)
+              << "Initial soil water from saved wetting fronts: "
+              << state->lgar_mass_balance.volstart_cm << " cm\n";
+  }
+
+  if (is_non_vadose_state_path_set) {
+    std::cerr << "Non-vadose restart variables initialized from: "
+              << state->lgar_bmi_params.init_non_vadose_state_path << "\n";
+    std::cerr << std::fixed << std::setprecision(12)
+              << "Initial fast reservoir storage: "
+              << state->lgar_mass_balance.CR_fast_storage_cm << " cm\n";
+    std::cerr << "Initial slow reservoir storage: "
+              << state->lgar_mass_balance.CR_slow_storage_cm << " cm\n";
+    std::cerr << "Saved volon_timestep_cm: "
+              << state->lgar_mass_balance.volon_timestep_cm << " cm\n";
+    std::cerr << "Saved precip_previous_timestep_cm: "
+              << state->lgar_bmi_params.precip_previous_timestep_cm << " cm\n";
+    std::cerr << "Saved runoff_in_prev_step: "
+              << (state->lgar_bmi_params.runoff_in_prev_step ? "true" : "false")
+              << "\n";
+    std::cerr << "Saved cache_fluxes: "
+              << (state->lgar_mass_balance.cache_fluxes ? "true" : "false")
+              << "\n";
+    std::cerr << "Saved cache_count: "
+              << state->lgar_bmi_params.cache_count << "\n";
+    std::cerr << "Saved previous_AET: "
+              << state->lgar_mass_balance.previous_AET << " cm\n";
+    std::cerr << "Saved previous_PET: "
+              << state->lgar_mass_balance.previous_PET << " cm\n";
+    std::cerr << "Saved previous_recharge: "
+              << state->lgar_mass_balance.previous_recharge << " cm\n";
+    std::cerr << "Saved accumulated_PET: "
+              << state->lgar_mass_balance.accumulated_PET << " cm\n";
+    std::cerr << "Saved accumulated_free_drainage: "
+              << state->lgar_mass_balance.accumulated_free_drainage << " cm\n";
+  }
+
+  if (is_giuh_state_path_set) {
+    std::cerr << "GIUH runoff queue requested from: "
+              << state->lgar_bmi_params.init_giuh_state_path << "\n";
+    std::cerr << "Expected no. of GIUH ordinates: "
+              << state->lgar_bmi_params.num_giuh_ordinates << "\n";
+  }
+
+  std::cerr.flags(old_flags);
+  std::cerr.precision(old_precision);
+  std::cerr << "          *****         \n";
+}
+
 
 // ############################################################################################
 /*
@@ -84,11 +166,10 @@ extern void lgar_initialize(string config_file, struct model_state *state)
   if (!state->lgar_bmi_params.is_invalid_soil_type){
     int soil;
     
+    state->lgar_bmi_params.num_wetting_fronts = listLength(state->head);
     state->lgar_bmi_params.shape[0] = state->lgar_bmi_params.num_layers;
     state->lgar_bmi_params.shape[1] = state->lgar_bmi_params.num_wetting_fronts;
 
-    // initial number of wetting fronts are same are number of layers
-    state->lgar_bmi_params.num_wetting_fronts           = state->lgar_bmi_params.num_layers;
     state->lgar_bmi_params.soil_depth_wetting_fronts    = new double[state->lgar_bmi_params.num_wetting_fronts];
     state->lgar_bmi_params.soil_moisture_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts];
 
@@ -113,10 +194,10 @@ extern void lgar_initialize(string config_file, struct model_state *state)
     state->lgar_calib_params.spf_factor = state->lgar_bmi_params.spf_factor;
 
     struct wetting_front *current = state->head;
-    for (int i=0; i<state->lgar_bmi_params.num_wetting_fronts; i++) { // note that this only works because at init there is 1 WF per layer, otherwise get soil type from current
+    for (int i=0; i<state->lgar_bmi_params.num_wetting_fronts; i++) {
       assert (current != NULL);
       
-      soil = state->lgar_bmi_params.layer_soil_type[i+1];
+      soil = state->lgar_bmi_params.layer_soil_type[current->layer_num];
 
       state->lgar_bmi_params.soil_moisture_wetting_fronts[i] = current->theta;
       state->lgar_bmi_params.soil_depth_wetting_fronts[i]    = current->depth_cm * state->units.cm_to_m;
@@ -128,7 +209,7 @@ extern void lgar_initialize(string config_file, struct model_state *state)
       // state->lgar_calib_params.vg_alpha[i] = state->soil_properties[soil].vg_alpha_per_cm;
       // state->lgar_calib_params.Ksat[i]     = state->soil_properties[soil].Ksat_cm_per_h;
 
-      if (i==0){
+      if (current->layer_num==1){
         state->lgar_calib_params.theta_e_1  = state->soil_properties[soil].theta_e;
         state->lgar_calib_params.theta_r_1  = state->soil_properties[soil].theta_r;
         state->lgar_calib_params.vg_n_1     = state->soil_properties[soil].vg_n;
@@ -136,7 +217,7 @@ extern void lgar_initialize(string config_file, struct model_state *state)
         state->lgar_calib_params.Ksat_1     = state->soil_properties[soil].Ksat_cm_per_h;
       }
 
-      if (i==1){
+      if (current->layer_num==2){
         state->lgar_calib_params.theta_e_2  = state->soil_properties[soil].theta_e;
         state->lgar_calib_params.theta_r_2  = state->soil_properties[soil].theta_r;
         state->lgar_calib_params.vg_n_2     = state->soil_properties[soil].vg_n;
@@ -144,7 +225,7 @@ extern void lgar_initialize(string config_file, struct model_state *state)
         state->lgar_calib_params.Ksat_2     = state->soil_properties[soil].Ksat_cm_per_h;
       }
 
-      if (i==2){
+      if (current->layer_num==3){
         state->lgar_calib_params.theta_e_3  = state->soil_properties[soil].theta_e;
         state->lgar_calib_params.theta_r_3  = state->soil_properties[soil].theta_r;
         state->lgar_calib_params.vg_n_3     = state->soil_properties[soil].vg_n;
@@ -174,23 +255,46 @@ extern void lgar_initialize(string config_file, struct model_state *state)
   state->lgar_bmi_input_params->precipitation_mm_per_h = -1.0;
   state->lgar_bmi_input_params->PET_mm_per_h           = -1.0;
 
-  // initialize all global mass balance variables to zero
-  state->lgar_mass_balance.volprecip_cm              = 0.0;
-  state->lgar_mass_balance.volin_cm                  = 0.0;
-  state->lgar_mass_balance.volend_cm                 = 0.0;
-  state->lgar_mass_balance.volCRend_cm               = 0.0;
-  state->lgar_mass_balance.volAET_cm                 = 0.0;
-  state->lgar_mass_balance.volrech_cm                = 0.0;
-  state->lgar_mass_balance.volrunoff_cm              = 0.0;
-  state->lgar_mass_balance.volrunoff_giuh_cm         = 0.0;
-  state->lgar_mass_balance.volQ_cm                   = 0.0;
-  state->lgar_mass_balance.volQ_CR_cm                = 0.0;
-  state->lgar_mass_balance.volPET_cm                 = 0.0;
-  state->lgar_mass_balance.volon_cm                  = 0.0;
-  state->lgar_mass_balance.volon_timestep_cm         = 0.0; /* setting volon and precip at the initial time to 0.0
-							       as they determine the creation of surficail wetting front */
-  state->lgar_bmi_params.precip_previous_timestep_cm = 0.0;
-  state->lgar_mass_balance.volchange_calib_cm        = 0.0;
+  bool non_vadose_restart_loaded =
+    !state->lgar_bmi_params.is_invalid_soil_type &&
+    !state->lgar_bmi_params.init_state_path.empty() &&
+    !state->lgar_bmi_params.init_non_vadose_state_path.empty();
+
+  state->lgar_mass_balance.volprecip_cm       = 0.0;
+  state->lgar_mass_balance.volin_cm           = 0.0;
+  state->lgar_mass_balance.volend_cm          = 0.0;
+  state->lgar_mass_balance.volAET_cm          = 0.0;
+  state->lgar_mass_balance.volrech_cm         = 0.0;
+  state->lgar_mass_balance.volrunoff_cm       = 0.0;
+  state->lgar_mass_balance.volrunoff_giuh_cm  = 0.0;
+  state->lgar_mass_balance.volQ_cm            = 0.0;
+  state->lgar_mass_balance.volQ_CR_cm         = 0.0;
+  state->lgar_mass_balance.volPET_cm          = 0.0;
+  state->lgar_mass_balance.volon_cm           = 0.0;
+  state->lgar_mass_balance.volchange_calib_cm = 0.0;
+
+  if (!non_vadose_restart_loaded) {
+    state->lgar_mass_balance.CR_fast_storage_cm = 0.0;
+    state->lgar_mass_balance.CR_slow_storage_cm = 0.0;
+    state->lgar_mass_balance.volCRend_cm = 0.0;
+    state->lgar_mass_balance.volCRend_timestep_cm = 0.0;
+    state->lgar_mass_balance.volCRstart_cm = 0.0;
+    state->lgar_mass_balance.volon_timestep_cm = 0.0; /* setting volon and precip at the initial time to 0.0
+							       as they determine the creation of surficial wetting front */
+    state->lgar_bmi_params.precip_previous_timestep_cm = 0.0;
+    state->lgar_bmi_params.runoff_in_prev_step = false;
+  }
+  else {
+    state->lgar_mass_balance.volCRend_cm =
+        state->lgar_mass_balance.CR_fast_storage_cm +
+        state->lgar_mass_balance.CR_slow_storage_cm;
+    state->lgar_mass_balance.volCRend_timestep_cm =
+        state->lgar_mass_balance.CR_fast_storage_cm +
+        state->lgar_mass_balance.CR_slow_storage_cm;
+    state->lgar_mass_balance.volCRstart_cm =
+        state->lgar_mass_balance.CR_fast_storage_cm +
+        state->lgar_mass_balance.CR_slow_storage_cm;
+  }
 }
 
 
@@ -300,6 +404,9 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
   bool is_giuh_ordinates_set        = false;
   bool is_soil_z_set                = false;
   bool is_ponded_depth_max_cm_set   = false;
+  bool is_state_path_set            = false;
+  bool is_non_vadose_state_path_set = false;
+  bool is_giuh_state_path_set       = false;
 
   string soil_params_file;
 
@@ -412,11 +519,6 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
       state->lgar_bmi_params.initial_psi_cm = stod(param_value);
       is_initial_psi_set = true;
 
-      if (verbosity.compare("high") == 0) {
-	std::cerr<<"Initial Psi : "<<state->lgar_bmi_params.initial_psi_cm<<"\n";
-	std::cerr<<"          *****         \n";
-      }
-
       continue;
     }
     else if (param_key == "max_valid_soil_types") {
@@ -550,6 +652,51 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
 
       continue;
     }
+    else if (param_key == "init_state_path") {
+      size_t first = param_value.find_first_not_of(" \t");
+      size_t last = param_value.find_last_not_of(" \t");
+      param_value = (first == string::npos) ? "" : param_value.substr(first, last - first + 1);
+
+      state->lgar_bmi_params.init_state_path = param_value;
+      is_state_path_set = true;
+
+      if (verbosity.compare("high") == 0) {
+        std::cerr << "init_state_path set to: "
+                  << state->lgar_bmi_params.init_state_path << "\n";
+        std::cerr << "          *****         \n";
+      }
+      continue;
+    }
+    else if (param_key == "init_non_vadose_state_path") {
+      size_t first = param_value.find_first_not_of(" \t");
+      size_t last = param_value.find_last_not_of(" \t");
+      param_value = (first == string::npos) ? "" : param_value.substr(first, last - first + 1);
+
+      state->lgar_bmi_params.init_non_vadose_state_path = param_value;
+      is_non_vadose_state_path_set = true;
+
+      if (verbosity.compare("high") == 0) {
+        std::cerr << "init_non_vadose_state_path set to: "
+                  << state->lgar_bmi_params.init_non_vadose_state_path << "\n";
+        std::cerr << "          *****         \n";
+      }
+      continue;
+    }
+    else if (param_key == "init_giuh_state_path") {
+      size_t first = param_value.find_first_not_of(" \t");
+      size_t last = param_value.find_last_not_of(" \t");
+      param_value = (first == string::npos) ? "" : param_value.substr(first, last - first + 1);
+
+      state->lgar_bmi_params.init_giuh_state_path = param_value;
+      is_giuh_state_path_set = true;
+
+      if (verbosity.compare("high") == 0) {
+        std::cerr << "init_giuh_state_path set to: "
+                  << state->lgar_bmi_params.init_giuh_state_path << "\n";
+        std::cerr << "          *****         \n";
+      }
+      continue;
+    }
     else if (param_key == "free_drainage_enabled") { 
       if (param_value == "false") {
         state->lgar_bmi_params.free_drainage_enabled = false;
@@ -862,6 +1009,11 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
     throw runtime_error(errMsg.str());
   }
 
+  if (verbosity.compare("high") == 0 && !is_state_path_set) {
+    std::cerr<<"Initial Psi : "<<state->lgar_bmi_params.initial_psi_cm<<"\n";
+    std::cerr<<"          *****         \n";
+  }
+
   if (!is_timestep_set) {
     stringstream errMsg;
     errMsg << "The configuration file \'" << config_file <<"\' does not set timestep. \n";
@@ -903,12 +1055,29 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
     throw runtime_error(errMsg.str());
   }
 
+  if (!state->lgar_bmi_params.is_invalid_soil_type &&
+      ((is_non_vadose_state_path_set && !is_state_path_set) || (!is_non_vadose_state_path_set && is_state_path_set)) &&
+      is_a_set) {
+    stringstream errMsg;
+    errMsg << "The configuration file \'" << config_file <<"\' sets one of init_non_vadose_state_path or init_state_path but not both, while a nonlinear reservoir is desired. Either both or neither must be set in this case. \n";
+    throw runtime_error(errMsg.str());
+  }
+
   if (state->lgar_bmi_params.free_drainage_to_CR && !state->lgar_bmi_params.free_drainage_enabled){
     stringstream errMsg;
     errMsg << "The configuration file \'" << config_file <<"\' sets free_drainage_to_CR as true but sets free_drainage_enabled as false. free_drainage_to_CR being true requires free drainage being enabled \n";
     throw runtime_error(errMsg.str());
   }
 
+  if (!state->lgar_bmi_params.is_invalid_soil_type &&
+      is_giuh_ordinates_set &&
+      !is_giuh_state_path_set &&
+      is_state_path_set) {
+    stringstream errMsg;
+    errMsg << "The configuration file \'" << config_file <<"\' sets GIUH ordinates and is loading wetting fronts but is not loading a GIUH queue. Either provide init_giuh_state_path, or do not load a wetting front state. \n";
+    throw runtime_error(errMsg.str());
+  }
+
   if (!is_forcing_resolution_set) {
     stringstream errMsg;
     errMsg << "The configuration file \'" << config_file <<"\' does not set forcing_resolution. \n";
@@ -970,9 +1139,29 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
     state->lgar_bmi_params.frozen_factor[i] = 1.0;
 
   state->head = NULL; //this will be updated if there are only valid soil types, but if there are any invalid soil types, it will remain null
-  InitializeWettingFronts(state->lgar_bmi_params.is_invalid_soil_type, state->lgar_bmi_params.num_layers, state->lgar_bmi_params.initial_psi_cm,
-        state->lgar_bmi_params.layer_soil_type, state->lgar_bmi_params.cum_layer_thickness_cm,
-        state->lgar_bmi_params.frozen_factor, &state->head, state->soil_properties);
+  if (!state->lgar_bmi_params.is_invalid_soil_type) {
+    if (!is_state_path_set) {
+      InitializeWettingFronts(state->lgar_bmi_params.is_invalid_soil_type, state->lgar_bmi_params.num_layers, state->lgar_bmi_params.initial_psi_cm,
+            state->lgar_bmi_params.layer_soil_type, state->lgar_bmi_params.cum_layer_thickness_cm,
+            state->lgar_bmi_params.frozen_factor, &state->head, state->soil_properties);
+    }
+    else {
+      InitializeWettingFrontsFromCSV(
+          state->lgar_bmi_params.num_layers,
+          state->lgar_bmi_params.init_state_path.c_str(),
+          state->lgar_bmi_params.layer_soil_type,
+          state->lgar_bmi_params.cum_layer_thickness_cm,
+          state->lgar_bmi_params.frozen_factor,
+          &state->head,
+          state->soil_properties);
+
+      if (!state->lgar_bmi_params.init_non_vadose_state_path.empty()) {
+        InitializenonvadoseStateFromCSV(
+            state->lgar_bmi_params.init_non_vadose_state_path.c_str(),
+            state);
+      }
+    }
+  }
   
   if (verbosity.compare("none") != 0) {
     std::cerr<<"--- Initial state/conditions --- \n";
@@ -990,19 +1179,21 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
 
   state->lgar_bmi_params.ponded_depth_cm    = 0.0; // initially we start with a dry surface (no surface ponding)
   state->lgar_bmi_params.nint               = 120; // hacked, not needed to be an input option
-  state->lgar_bmi_params.num_wetting_fronts = state->lgar_bmi_params.num_layers;
   if (state->lgar_bmi_params.is_invalid_soil_type){
     state->lgar_bmi_params.num_wetting_fronts = 0;
   }
   else {
-    assert (state->lgar_bmi_params.num_layers == listLength(state->head));
+    state->lgar_bmi_params.num_wetting_fronts = listLength(state->head);
+    if (!is_state_path_set) {
+      assert (state->lgar_bmi_params.num_layers == state->lgar_bmi_params.num_wetting_fronts);
+    }
   }
   
-
-  if (verbosity.compare("high") == 0) {
-    std::cerr<<"Initial ponded depth is set to zero. \n";
-    std::cerr<<"No. of spatial intervals used in trapezoidal integration to compute G : "<<state->lgar_bmi_params.nint<<"\n";
-  }
+  PrintSavedStateInitializationDetails(
+      state,
+      is_state_path_set,
+      is_non_vadose_state_path_set,
+      is_giuh_state_path_set);
 
   state->lgar_bmi_input_params     = new lgar_bmi_input_parameters;
   state->lgar_bmi_params.time_s    = 0.0;
@@ -1172,6 +1363,7 @@ extern void lgar_global_mass_balance(struct model_state *state, double *giuh_run
   double volrech            = state->lgar_mass_balance.volrech_cm;
   double volend             = state->lgar_mass_balance.volend_cm;
   double volCRend           = state->lgar_mass_balance.volCRend_cm;
+  double volCRstart         = state->lgar_mass_balance.volCRstart_cm;
   double volrunoff_giuh     = state->lgar_mass_balance.volrunoff_giuh_cm;
   double volend_giuh_cm     = 0.0;
   double total_Q_cm         = state->lgar_mass_balance.volQ_cm;
@@ -1182,7 +1374,7 @@ extern void lgar_global_mass_balance(struct model_state *state, double *giuh_run
   for(int i=0; i <= state->lgar_bmi_params.num_giuh_ordinates; i++)
     volend_giuh_cm += giuh_runoff_queue_cm[i];
 
-  double global_error_cm = volstart + volprecip - volrunoff - volAET - volon - volrech - volend + volchange_calib_cm - volrunoff_CR - volCRend;
+  double global_error_cm = volstart + volprecip - volrunoff - volAET - volon - volrech - volend + volchange_calib_cm - volrunoff_CR - volCRend + volCRstart;
   
   printf("\n********************************************************* \n");
   printf("-------------------- Simulation Summary ----------------- \n");

From 03eff5d78691240c1b740c8579ee3351f5708914 Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@MacBook-Air.local>
Date: Wed, 13 May 2026 20:25:58 -0700
Subject: [PATCH 2/2] disable restart paths in Bushland config

---
 configs/config_lasam_Bushland.txt | 6 +++---
 src/init_from_state.cxx           | 3 +--
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/configs/config_lasam_Bushland.txt b/configs/config_lasam_Bushland.txt
index d2f3844..0d99ee0 100644
--- a/configs/config_lasam_Bushland.txt
+++ b/configs/config_lasam_Bushland.txt
@@ -19,6 +19,6 @@ b=3.0
 frac_to_CR=0.1
 spf_factor=0.9
 allow_flux_caching=true
-init_state_path=./data_layers.csv
-init_non_vadose_state_path=./data_non_vadose_state.csv
-init_giuh_state_path=./data_giuh_state.csv
+#init_state_path=./data_layers.csv
+#init_non_vadose_state_path=./data_non_vadose_state.csv
+#init_giuh_state_path=./data_giuh_state.csv
diff --git a/src/init_from_state.cxx b/src/init_from_state.cxx
index 0aba673..4cecdf5 100644
--- a/src/init_from_state.cxx
+++ b/src/init_from_state.cxx
@@ -267,8 +267,7 @@ extern void InitializeWettingFrontsFromCSV(
         r.depth_cm > layer_bottom_cm + 1.0e-8) {
       fprintf(stderr,
               "ERROR: restart front_num=%d has depth_cm=%.17g outside layer %d "
-              "bounds [%.17g, %.17g] cm. Regenerate old x10/mm restart files "
-              "with the current cm-unit writer before loading.\n",
+              "bounds [%.17g, %.17g] cm.\n",
               r.front_num, r.depth_cm, r.layer_num,
               layer_top_cm, layer_bottom_cm);
       free(recs);