From 8082e41bc3d2a37b1634b337871269959c7dba82 Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MBP-2.lan>
Date: Mon, 6 May 2024 15:10:27 -0700
Subject: [PATCH 01/13] adding the capability for an adaptive time step. If not
 provided in the config file, it defaults to false.

---
 configs/config_lasam_Bushland.txt     |  1 +
 configs/config_lasam_Phillipsburg.txt |  1 +
 configs/config_lasam_sft_ngen.txt     |  1 +
 giuh/giuh.c                           | 50 ++++++++++++++++++++-------
 giuh/giuh.h                           |  2 +-
 include/all.hxx                       |  1 +
 src/bmi_lgar.cxx                      | 41 ++++++++++++++++++++--
 src/lgar.cxx                          | 15 ++++++++
 8 files changed, 95 insertions(+), 17 deletions(-)

diff --git a/configs/config_lasam_Bushland.txt b/configs/config_lasam_Bushland.txt
index 06b6e89..cd141ef 100644
--- a/configs/config_lasam_Bushland.txt
+++ b/configs/config_lasam_Bushland.txt
@@ -13,3 +13,4 @@ max_soil_types=25
 wilting_point_psi=15495.0[cm]
 field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03
+adaptive_timestep=true
\ No newline at end of file
diff --git a/configs/config_lasam_Phillipsburg.txt b/configs/config_lasam_Phillipsburg.txt
index 1d230cf..d9f3041 100644
--- a/configs/config_lasam_Phillipsburg.txt
+++ b/configs/config_lasam_Phillipsburg.txt
@@ -13,3 +13,4 @@ max_soil_types=25
 wilting_point_psi=15495.0[cm]
 field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03
+adaptive_timestep=true
diff --git a/configs/config_lasam_sft_ngen.txt b/configs/config_lasam_sft_ngen.txt
index e9ff131..6c22a0f 100644
--- a/configs/config_lasam_sft_ngen.txt
+++ b/configs/config_lasam_sft_ngen.txt
@@ -15,3 +15,4 @@ field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03
 sft_coupled=true
 soil_z=10,20,30,40,50,60,70,80,90,100.0,110.,120,130.,140.,150.,160.,170.,180.,190.,200.0[cm]
+adaptive_timestep=true
diff --git a/giuh/giuh.c b/giuh/giuh.c
index 6a0c3ce..efe303c 100644
--- a/giuh/giuh.c
+++ b/giuh/giuh.c
@@ -4,36 +4,60 @@
 
 #include "giuh.h"
 #include <stdio.h>
+#include <math.h>
 
 
 //##############################################################
 //############### GIUH CONVOLUTION INTEGRAL   ##################
 //##############################################################
-extern double giuh_convolution_integral(double runoff_m,int num_giuh_ordinates, 
+extern double giuh_convolution_integral(int adaptive_timestep, double subtimestep_h, double runoff_m,int num_giuh_ordinates, 
 					double *giuh_ordinates, double *runoff_queue_m_per_timestep)
 {
   //##############################################################
   // This function solves the convolution integral involving N
   //  GIUH ordinates.
   //##############################################################
+
   double runoff_m_now;
-  int N,i;
+  int N,i,j;
   
   N=num_giuh_ordinates;
   runoff_queue_m_per_timestep[N]=0.0;
   
-  for(i=0;i<N;i++)
-    {
-      runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
+  if (subtimestep_h==1.0/12.0){
+    for(i=0;i<N;i++)
+      {
+        runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
+      }
+    
+    runoff_m_now=runoff_queue_m_per_timestep[0];
+    
+    for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
+      {
+        runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
+      }
+    runoff_queue_m_per_timestep[N-1]=0.0;
+  }
+  else {//The time step is not equal to the smallest internal time step. In this case, the giuh code must be repeated j times because the time between each entry in runoff_queue_m_per_timestep is fixed.
+    double factor = (1.0/12.0)/subtimestep_h;
+    runoff_m = runoff_m*factor;
+    int resample_limit = fmax(1,(int) 1.0/factor);
+    for (j=0;j<resample_limit;j++){
+      for(i=0;i<N;i++)
+        {
+          runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
+        }
+      
+      runoff_m_now+=runoff_queue_m_per_timestep[0];
+      
+      for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
+        {
+          runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
+        }
+      runoff_queue_m_per_timestep[N-1]=0.0;
     }
-  
-  runoff_m_now=runoff_queue_m_per_timestep[0];
-  
-  for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
-    {
-      runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
-    }
-  runoff_queue_m_per_timestep[N-1]=0.0;
+
+  }
   
   return runoff_m_now;
 }
diff --git a/giuh/giuh.h b/giuh/giuh.h
index de05e17..88a8997 100644
--- a/giuh/giuh.h
+++ b/giuh/giuh.h
@@ -2,7 +2,7 @@
 
 #define _GIUH_H
 
-extern double giuh_convolution_integral(double runoff_m, int num_giuh_ordinates, 
+extern double giuh_convolution_integral(int adaptive_timestep, double subtimestep_h, double runoff_m, int num_giuh_ordinates, 
                                    double *giuh_ordinates, double *runoff_queue_m_per_timestep);
 
 #endif
diff --git a/include/all.hxx b/include/all.hxx
index 6b7526c..c30a96d 100755
--- a/include/all.hxx
+++ b/include/all.hxx
@@ -131,6 +131,7 @@ 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 */
+  bool   adaptive_timestep = false;      // if set to true, model uses adaptive timestep 
   double ponded_depth_cm;                // amount of water on the surface unavailable for surface runoff
   double ponded_depth_max_cm;            // maximum amount of water on the surface unavailable for surface runoff
   double precip_previous_timestep_cm;    // amount of rainfall (previous time step)
diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index c11a8aa..380b397 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -33,7 +33,6 @@ Initialize (std::string config_file)
 
   num_giuh_ordinates = state->lgar_bmi_params.num_giuh_ordinates;
 
-
   /* giuh ordinates are static and read in the lgar.cxx, and we need to have a copy of it to pass to
      giuh.cxx, so allocating/copying here*/
   
@@ -76,7 +75,7 @@ Update()
   double mm_to_cm = 0.1; // unit conversion
 
   // local variables for readibility
-  int subcycles = state->lgar_bmi_params.forcing_interval;
+  int subcycles;
   int num_layers = state->lgar_bmi_params.num_layers;
 
   // local variables for a full timestep (i.e., timestep of the forcing data)
@@ -117,6 +116,7 @@ Update()
   double wilting_point_psi_cm = state->lgar_bmi_params.wilting_point_psi_cm;
   double field_capacity_psi_cm = state->lgar_bmi_params.field_capacity_psi_cm;
   bool use_closed_form_G = state->lgar_bmi_params.use_closed_form_G; 
+  bool adaptive_timestep = state->lgar_bmi_params.adaptive_timestep;
 
   // constant value used in the AET function
   double AET_thresh_Theta = 0.85;    // scaled soil moisture (0-1) above which AET=PET (fix later!)
@@ -131,6 +131,32 @@ Update()
 
   assert (state->lgar_bmi_input_params->precipitation_mm_per_h >= 0.0);
   assert(state->lgar_bmi_input_params->PET_mm_per_h >=0.0);
+
+  // adaptive time step is set 
+  if (adaptive_timestep){
+    if ( (state->lgar_bmi_input_params->precipitation_mm_per_h > 0.0) || (state->lgar_mass_balance.volon_timestep_cm > 0.0) ){
+      if (state->lgar_bmi_input_params->precipitation_mm_per_h > 10.0 || (state->lgar_mass_balance.volon_timestep_cm > 0.0) ){
+        subtimestep_h = 1.0/12.0;//case where precip > 1 cm/h, or there is any ponded head from the last time step
+        state->lgar_bmi_params.timestep_h = 1.0/12.0;
+      }
+      else {
+        subtimestep_h = 1.0/6.0;//case where there is precip > 0, but it is less than 1 cm/h, and there is no ponded head from the last time step
+        state->lgar_bmi_params.timestep_h = 1.0/6.0;
+      }
+    }
+    else {//case where the precip = 0 and there is no ponded head from the last time step 
+      subtimestep_h = 1.0;
+      state->lgar_bmi_params.timestep_h = 1.0;
+    }
+  }
+
+
+  state->lgar_bmi_params.forcing_interval = int(state->lgar_bmi_params.forcing_resolution_h/state->lgar_bmi_params.timestep_h+1.0e-08); // add 1.0e-08 to prevent truncation error
+  subcycles = state->lgar_bmi_params.forcing_interval;
+
+  if (verbosity.compare("high") == 0) {//and adaptive time step is engaged? 
+    printf("time step size in hours: %lf \n", state->lgar_bmi_params.timestep_h);
+  }
   
   // subcycling loop (loop over model's timestep)
   for (int cycle=1; cycle <= subcycles; cycle++) {
@@ -360,7 +386,7 @@ Update()
     /*----------------------------------------------------------------------*/
     // compute giuh runoff for the subtimestep
     surface_runoff_subtimestep_cm = volrunoff_subtimestep_cm;
-    volrunoff_giuh_subtimestep_cm = giuh_convolution_integral(volrunoff_subtimestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue);
+    volrunoff_giuh_subtimestep_cm = giuh_convolution_integral(adaptive_timestep, subtimestep_h, volrunoff_subtimestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue);
 
     surface_runoff_timestep_cm += surface_runoff_subtimestep_cm ;
     volrunoff_giuh_timestep_cm += volrunoff_giuh_subtimestep_cm;
@@ -415,6 +441,15 @@ Update()
 
   } // end of subcycling
 
+  /*----------------------------------------------------------------------*/
+  // compute giuh runoff for the timestep
+  // volrunoff_giuh_timestep_cm = giuh_convolution_integral(volrunoff_timestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue);
+
+  // // total mass of water leaving the system, at this time it is the giuh-only, but later will add groundwater component as well.
+
+  // volQ_timestep_cm += volrunoff_giuh_timestep_cm;
+
+
   /*----------------------------------------------------------------------*/
   // Everything related to lgar state is done at this point, now time to update some dynamic variables
 
diff --git a/src/lgar.cxx b/src/lgar.cxx
index fd412a3..6b74e84 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -213,6 +213,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
   // setting these options to false (defualt) 
   state->lgar_bmi_params.sft_coupled       = false;
   state->lgar_bmi_params.use_closed_form_G = false;
+  state->lgar_bmi_params.adaptive_timestep = false;
   
   bool is_layer_thickness_set       = false;
   bool is_initial_psi_set           = false;
@@ -398,6 +399,20 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
 
       continue;
     }
+    else if (param_key == "adaptive_timestep") { 
+      if (param_value == "false") {
+        state->lgar_bmi_params.adaptive_timestep = false;
+      }
+      else if (param_value == "true") {
+        state->lgar_bmi_params.adaptive_timestep = true;
+      }
+      else {
+	std::cerr<<"Invalid option: adaptive_timestep must be true or false. \n";
+        abort();
+      }
+
+      continue;
+    }
     else if (param_key == "timestep") {
       state->lgar_bmi_params.timestep_h = stod(param_value);
 

From 4f1af7924b5bb6aee437a40fa3221eac28641087 Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MBP-2.lan>
Date: Mon, 6 May 2024 15:15:14 -0700
Subject: [PATCH 02/13] updating readme

---
 configs/README.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/configs/README.md b/configs/README.md
index 60ea2c4..dc79c3f 100644
--- a/configs/README.md
+++ b/configs/README.md
@@ -24,3 +24,4 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | sft_coupled | Boolean | true, false | - | model coupling | impacts hydraulic conductivity | couples LASAM to SFT. Coupling to SFT reduces hydraulic conducitivity, and hence infiltration, when soil is frozen|
 | soil_z | double (1D array) | - | cm | spatial resolution | - | vertical resolution of the soil column (computational domain of the SFT model) |
 | calib_params | Boolean | true, false | - | calibratable params flag | impacts soil properties | If set to true, soil `smcmax`, `smcmin`, `vg_n`, `vg_alpha`, `hydraulic_conductivity`, `field_capacity_psi`, and `ponded_depth_max` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content |
+| adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep. If set to false, LGAR will use the above specified timestep as a fixed timestep. If set to true or false, the internal giuh timestep will be equal to the above specified timestep. |

From 3912c14df35b8399703e0bdfa87155f0b11b2f6b Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MBP-2.lan>
Date: Mon, 6 May 2024 15:17:52 -0700
Subject: [PATCH 03/13] updating readme

---
 configs/README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configs/README.md b/configs/README.md
index dc79c3f..0c687fa 100644
--- a/configs/README.md
+++ b/configs/README.md
@@ -24,4 +24,4 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | sft_coupled | Boolean | true, false | - | model coupling | impacts hydraulic conductivity | couples LASAM to SFT. Coupling to SFT reduces hydraulic conducitivity, and hence infiltration, when soil is frozen|
 | soil_z | double (1D array) | - | cm | spatial resolution | - | vertical resolution of the soil column (computational domain of the SFT model) |
 | calib_params | Boolean | true, false | - | calibratable params flag | impacts soil properties | If set to true, soil `smcmax`, `smcmin`, `vg_n`, `vg_alpha`, `hydraulic_conductivity`, `field_capacity_psi`, and `ponded_depth_max` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content |
-| adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep. If set to false, LGAR will use the above specified timestep as a fixed timestep. If set to true or false, the internal giuh timestep will be equal to the above specified timestep. |
+| adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep. If set to false, LGAR will use the above specified timestep as a fixed timestep. If set to true or false, the internal giuh timestep will be equal to the above specified timestep. Testing indicates that setting this value to true substantially decreases runtime while negligibly changing the simulation. We recommend this to be set to true. |

From f091570a181a56af6b05638b99bee5ec2d8e0fbb Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MBP-2.lan>
Date: Mon, 6 May 2024 16:29:32 -0700
Subject: [PATCH 04/13] updating giuh fxn with adaptive time step

---
 configs/config_lasam_Phillipsburg.txt |  2 +-
 giuh/giuh.c                           | 78 ++++++++++++++++++---------
 src/lgar.cxx                          |  3 +-
 3 files changed, 55 insertions(+), 28 deletions(-)

diff --git a/configs/config_lasam_Phillipsburg.txt b/configs/config_lasam_Phillipsburg.txt
index d9f3041..6a69349 100644
--- a/configs/config_lasam_Phillipsburg.txt
+++ b/configs/config_lasam_Phillipsburg.txt
@@ -3,7 +3,7 @@ forcing_file=./forcing/forcing_data_resampled_uniform_Phillipsburg.csv
 soil_params_file=./data/vG_default_params.dat
 layer_thickness=44.0,131.0,25.0[cm]
 initial_psi=2000.0[cm]
-timestep=300[sec]
+timestep=3600[sec]
 endtime=7500.0[hr]
 forcing_resolution=3600[sec]
 ponded_depth_max=2[cm]
diff --git a/giuh/giuh.c b/giuh/giuh.c
index efe303c..d56323e 100644
--- a/giuh/giuh.c
+++ b/giuh/giuh.c
@@ -18,37 +18,21 @@ extern double giuh_convolution_integral(int adaptive_timestep, double subtimeste
   //  GIUH ordinates.
   //##############################################################
 
-  double runoff_m_now;
-  int N,i,j;
-  
-  N=num_giuh_ordinates;
-  runoff_queue_m_per_timestep[N]=0.0;
-  
-  if (subtimestep_h==1.0/12.0){
-    for(i=0;i<N;i++)
-      {
-        runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
-      }
+  if (adaptive_timestep){//adaptive time step case
+
+    double runoff_m_now;
+    int N,i,j;
     
-    runoff_m_now=runoff_queue_m_per_timestep[0];
+    N=num_giuh_ordinates;
+    runoff_queue_m_per_timestep[N]=0.0;
     
-    for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
-      {
-        runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
-      }
-    runoff_queue_m_per_timestep[N-1]=0.0;
-  }
-  else {//The time step is not equal to the smallest internal time step. In this case, the giuh code must be repeated j times because the time between each entry in runoff_queue_m_per_timestep is fixed.
-    double factor = (1.0/12.0)/subtimestep_h;
-    runoff_m = runoff_m*factor;
-    int resample_limit = fmax(1,(int) 1.0/factor);
-    for (j=0;j<resample_limit;j++){
+    if (subtimestep_h==1.0/12.0){
       for(i=0;i<N;i++)
         {
           runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
         }
       
-      runoff_m_now+=runoff_queue_m_per_timestep[0];
+      runoff_m_now=runoff_queue_m_per_timestep[0];
       
       for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
         {
@@ -56,10 +40,52 @@ extern double giuh_convolution_integral(int adaptive_timestep, double subtimeste
         }
       runoff_queue_m_per_timestep[N-1]=0.0;
     }
+    else {//The time step is not equal to the smallest internal time step. In this case, the giuh code must be repeated j times because the time between each entry in runoff_queue_m_per_timestep is fixed.
+      double factor = (1.0/12.0)/subtimestep_h;
+      runoff_m = runoff_m*factor;
+      int resample_limit = fmax(1,(int) 1.0/factor);
+      for (j=0;j<resample_limit;j++){
+        for(i=0;i<N;i++)
+          {
+            runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
+          }
+        
+        runoff_m_now+=runoff_queue_m_per_timestep[0];
+        
+        for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
+          {
+            runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
+          }
+        runoff_queue_m_per_timestep[N-1]=0.0;
+      }
+
+    }
+    
+    return runoff_m_now;
+  }
 
+  else {//fixed time step case
+    double runoff_m_now;
+    int N,i;
+    
+    N=num_giuh_ordinates;
+    runoff_queue_m_per_timestep[N]=0.0;
+    
+    for(i=0;i<N;i++)
+      {
+        runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
+      }
+    
+    runoff_m_now=runoff_queue_m_per_timestep[0];
+    
+    for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
+      {
+        runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
+      }
+    runoff_queue_m_per_timestep[N-1]=0.0;
+    
+    return runoff_m_now;
   }
-  
-  return runoff_m_now;
 }
 
 
diff --git a/src/lgar.cxx b/src/lgar.cxx
index 6b74e84..ee957b6 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -622,7 +622,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
 
 
   if (is_giuh_ordinates_set) {
-    int factor = int(1.0/state->lgar_bmi_params.timestep_h);
+    // int factor = int(1.0/state->lgar_bmi_params.timestep_h);
+    int factor = int(12.0);
 
     state->lgar_bmi_params.num_giuh_ordinates = factor * (giuh_ordinates_temp.size() - 1);
     state->lgar_bmi_params.giuh_ordinates = new double[state->lgar_bmi_params.num_giuh_ordinates+1];

From 904a6638d705efb0607df00b0cc1cc6fad777a5f Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MBP-2.lan>
Date: Tue, 7 May 2024 10:57:59 -0700
Subject: [PATCH 05/13] giuh now uses a fixed time step that is hardcoded to be
 the the same as the smallest adaptive time step for the rest of the model. If
 the adaptive time step is off, the giuh time step will be the same as the
 timestep input in the config file.

---
 src/lgar.cxx | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/lgar.cxx b/src/lgar.cxx
index ee957b6..8c0e700 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -622,8 +622,13 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
 
 
   if (is_giuh_ordinates_set) {
-    // int factor = int(1.0/state->lgar_bmi_params.timestep_h);
-    int factor = int(12.0);
+    int factor;
+    if (!state->lgar_bmi_params.adaptive_timestep){
+      factor = int(1.0/state->lgar_bmi_params.timestep_h);
+    }
+    else {
+      factor = int(12.0);
+    }
 
     state->lgar_bmi_params.num_giuh_ordinates = factor * (giuh_ordinates_temp.size() - 1);
     state->lgar_bmi_params.giuh_ordinates = new double[state->lgar_bmi_params.num_giuh_ordinates+1];

From bf3293708a402a93ca84454203b714c407e5be20 Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MBP-2.lan>
Date: Tue, 7 May 2024 11:01:25 -0700
Subject: [PATCH 06/13] updating readme in configs

---
 configs/README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configs/README.md b/configs/README.md
index 0c687fa..30286a1 100644
--- a/configs/README.md
+++ b/configs/README.md
@@ -24,4 +24,4 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | sft_coupled | Boolean | true, false | - | model coupling | impacts hydraulic conductivity | couples LASAM to SFT. Coupling to SFT reduces hydraulic conducitivity, and hence infiltration, when soil is frozen|
 | soil_z | double (1D array) | - | cm | spatial resolution | - | vertical resolution of the soil column (computational domain of the SFT model) |
 | calib_params | Boolean | true, false | - | calibratable params flag | impacts soil properties | If set to true, soil `smcmax`, `smcmin`, `vg_n`, `vg_alpha`, `hydraulic_conductivity`, `field_capacity_psi`, and `ponded_depth_max` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content |
-| adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep. If set to false, LGAR will use the above specified timestep as a fixed timestep. If set to true or false, the internal giuh timestep will be equal to the above specified timestep. Testing indicates that setting this value to true substantially decreases runtime while negligibly changing the simulation. We recommend this to be set to true. |
+| adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep, and the above timestep is ignored. If set to false, LGAR will use the above specified timestep as a fixed timestep. Testing indicates that setting this value to true substantially decreases runtime while negligibly changing the simulation. We recommend this to be set to true. |

From 313524f5c52540d3750448a230f1fbf39805cb5c Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MacBook-Pro-2.local>
Date: Thu, 9 May 2024 10:07:32 -0700
Subject: [PATCH 07/13] cleaned up code that was used for debugging

---
 src/bmi_lgar.cxx   | 12 ++----------
 src/lgar.cxx       | 16 ++++++++++------
 src/soil_funcs.cxx |  2 +-
 3 files changed, 13 insertions(+), 17 deletions(-)

diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index 380b397..96d70af 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -136,11 +136,11 @@ Update()
   if (adaptive_timestep){
     if ( (state->lgar_bmi_input_params->precipitation_mm_per_h > 0.0) || (state->lgar_mass_balance.volon_timestep_cm > 0.0) ){
       if (state->lgar_bmi_input_params->precipitation_mm_per_h > 10.0 || (state->lgar_mass_balance.volon_timestep_cm > 0.0) ){
-        subtimestep_h = 1.0/12.0;//case where precip > 1 cm/h, or there is any ponded head from the last time step
+        subtimestep_h = 1.0/12.0;//case where precip > 1 cm/h, or there is ponded head from the last time step
         state->lgar_bmi_params.timestep_h = 1.0/12.0;
       }
       else {
-        subtimestep_h = 1.0/6.0;//case where there is precip > 0, but it is less than 1 cm/h, and there is no ponded head from the last time step
+        subtimestep_h = 1.0/6.0;//case where there is either nonzero precip or ponded head, but both are less than the threshold that requires the smallest internal time step 
         state->lgar_bmi_params.timestep_h = 1.0/6.0;
       }
     }
@@ -441,14 +441,6 @@ Update()
 
   } // end of subcycling
 
-  /*----------------------------------------------------------------------*/
-  // compute giuh runoff for the timestep
-  // volrunoff_giuh_timestep_cm = giuh_convolution_integral(volrunoff_timestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue);
-
-  // // total mass of water leaving the system, at this time it is the giuh-only, but later will add groundwater component as well.
-
-  // volQ_timestep_cm += volrunoff_giuh_timestep_cm;
-
 
   /*----------------------------------------------------------------------*/
   // Everything related to lgar state is done at this point, now time to update some dynamic variables
diff --git a/src/lgar.cxx b/src/lgar.cxx
index 8c0e700..2e69176 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -627,7 +627,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
       factor = int(1.0/state->lgar_bmi_params.timestep_h);
     }
     else {
-      factor = int(12.0);
+      factor = 12;
     }
 
     state->lgar_bmi_params.num_giuh_ordinates = factor * (giuh_ordinates_temp.size() - 1);
@@ -1293,7 +1293,6 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
 	if (wf_free_drainage_demand == wf)
 	  prior_mass += precip_mass_to_add - (free_drainage_demand + actual_ET_demand);
   // theta mass balance computes new theta that conserves the mass; new theta is assigned to the current wetting front
-
 	double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, AET_demand_cm,
 						   delta_thetas, delta_thickness, soil_type, soil_properties);
   actual_ET_demand = *AET_demand_cm;
@@ -2504,7 +2503,7 @@ extern void lgar_dzdt_calc(bool use_closed_form_G, int nint, double h_p, int *so
 // ############################################################################################
 extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm, double new_mass,
 				      double prior_mass, double *AET_demand_cm, double *delta_theta, double *delta_thickness,
-				      int *soil_type, struct soil_properties_ *soil_properties)
+				      int *soil_type, struct soil_properties_ *soil_properties)//TODO: remove head as an argument, was just used for debugging purposes
 {
 
   double psi_cm_loc = psi_cm; // location psi
@@ -2594,13 +2593,13 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm
 
     // stop the loop if the error between the current and previous psi is less than 10^-15
     // 1. enough accuracy, 2. the algorithm can't improve the error further,
-    // 3. avoid infinite loop, 4. handles a corner case when prior mass is tiny (e.g., <1.E-5)
+    // 3. avoid infinite loop, 4. handles case where theta is very close to theta_r and convergence might be possible but would be extremely slow
+    // 5. handles a corner case when prior mass is tiny (e.g., <1.E-5)
     // printf("A1 = %.20f, %.18f %.18f %.18f %.18f \n ",fabs(psi_cm_loc - psi_cm_loc_prev) , psi_cm_loc, psi_cm_loc_prev, factor, delta_mass);
     
     if (fabs(psi_cm_loc - psi_cm_loc_prev) < 1E-15 && factor < 1E-13) break;
 
-    // another condition to avoid infinite loop when the error does not improve
-    if (fabs(delta_mass - delta_mass_prev) < 1E-15)
+    if (fabs(delta_mass - delta_mass_prev) < 1e-15)
       count_no_mass_change++;
     else
       count_no_mass_change = 0;
@@ -2609,6 +2608,11 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm
     if (count_no_mass_change == break_no_mass_change)
       break;
     
+    if (psi_cm_loc > 1e7){//there are rare cases where theta is very close to theta_r, and delta_mass - delta_mass_prev will change extremely slowly. Convergence might be possible but the model will take hours to converge rather than seconds. 
+    //an alternative solution was to change the threshold in if (fabs(delta_mass - delta_mass_prev) < 1e-15) to 1e-11, but that solution is somewhat slow. 
+      break;
+    }
+
     // -ve pressure will return NAN, so terminate the loop if previous psi is way small and current psi is zero
     // the wetting front is almost saturated
     if (psi_cm_loc <= 0 && psi_cm_loc_prev < 0) break;
diff --git a/src/soil_funcs.cxx b/src/soil_funcs.cxx
index bfedb11..7da1adc 100755
--- a/src/soil_funcs.cxx
+++ b/src/soil_funcs.cxx
@@ -154,7 +154,7 @@ double calc_theta_from_h(double h,double alpha, double m, double n, double theta
 /***********************************/
 double calc_Se_from_h(double h,double alpha, double m, double n)
 {
-  if(is_epsilon_less_than(h,1.0e-01)) return 1.0;  // this function doesn't work well ffor tiny h
+  if(is_epsilon_less_than(h,1.0e-10)) return 1.0;  // this function doesn't work well ffor tiny h
   else return(1.0/(pow(1.0+pow(alpha*h,n),m)));
 }
 

From 38e720215515c1e6b58f63fa8148827400f0e598 Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MacBook-Pro-2.local>
Date: Thu, 9 May 2024 10:10:15 -0700
Subject: [PATCH 08/13] updating Phillipsburg config file

---
 configs/config_lasam_Phillipsburg.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configs/config_lasam_Phillipsburg.txt b/configs/config_lasam_Phillipsburg.txt
index 6a69349..d9f3041 100644
--- a/configs/config_lasam_Phillipsburg.txt
+++ b/configs/config_lasam_Phillipsburg.txt
@@ -3,7 +3,7 @@ forcing_file=./forcing/forcing_data_resampled_uniform_Phillipsburg.csv
 soil_params_file=./data/vG_default_params.dat
 layer_thickness=44.0,131.0,25.0[cm]
 initial_psi=2000.0[cm]
-timestep=3600[sec]
+timestep=300[sec]
 endtime=7500.0[hr]
 forcing_resolution=3600[sec]
 ponded_depth_max=2[cm]

From 3a531c8ef69d8f71fdcbd4290bea49db0cb9acd9 Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MacBook-Pro-2.local>
Date: Thu, 9 May 2024 10:17:07 -0700
Subject: [PATCH 09/13] cleaned up code that was used for debugging

---
 src/lgar.cxx | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/lgar.cxx b/src/lgar.cxx
index 2e69176..621a6db 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -2503,7 +2503,7 @@ extern void lgar_dzdt_calc(bool use_closed_form_G, int nint, double h_p, int *so
 // ############################################################################################
 extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm, double new_mass,
 				      double prior_mass, double *AET_demand_cm, double *delta_theta, double *delta_thickness,
-				      int *soil_type, struct soil_properties_ *soil_properties)//TODO: remove head as an argument, was just used for debugging purposes
+				      int *soil_type, struct soil_properties_ *soil_properties)
 {
 
   double psi_cm_loc = psi_cm; // location psi
@@ -2599,7 +2599,8 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm
     
     if (fabs(psi_cm_loc - psi_cm_loc_prev) < 1E-15 && factor < 1E-13) break;
 
-    if (fabs(delta_mass - delta_mass_prev) < 1e-15)
+    // another condition to avoid infinite loop when the error does not improve
+    if (fabs(delta_mass - delta_mass_prev) < 1E-15)
       count_no_mass_change++;
     else
       count_no_mass_change = 0;

From f8569786846569170998a4b51e054aad50e8a469 Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MacBook-Pro-2.local>
Date: Mon, 10 Jun 2024 13:24:45 -0700
Subject: [PATCH 10/13] After chatting with Ahmad, we determined that it is
 best if the adaptive time step strategy does not alter the function in giuh.c
 because this function is used in multiple models. So, the adaptive time step
 method has been altered such that the giuh is computed at the hourly time
 step rather than at the sub timestep. This allows for the giuh fxn to be
 unaltered, and has the added bonus of not reqiring resampling for the giuh
 based on time step.

---
 giuh/giuh.c      | 90 +++++++++++-------------------------------------
 giuh/giuh.h      |  2 +-
 src/bmi_lgar.cxx | 24 ++++++-------
 src/lgar.cxx     | 15 ++------
 4 files changed, 36 insertions(+), 95 deletions(-)

diff --git a/giuh/giuh.c b/giuh/giuh.c
index d56323e..2605cba 100644
--- a/giuh/giuh.c
+++ b/giuh/giuh.c
@@ -4,89 +4,39 @@
 
 #include "giuh.h"
 #include <stdio.h>
-#include <math.h>
 
 
 //##############################################################
 //############### GIUH CONVOLUTION INTEGRAL   ##################
 //##############################################################
-extern double giuh_convolution_integral(int adaptive_timestep, double subtimestep_h, double runoff_m,int num_giuh_ordinates, 
+extern double giuh_convolution_integral(double runoff_m,int num_giuh_ordinates, 
 					double *giuh_ordinates, double *runoff_queue_m_per_timestep)
 {
   //##############################################################
   // This function solves the convolution integral involving N
   //  GIUH ordinates.
   //##############################################################
-
-  if (adaptive_timestep){//adaptive time step case
-
-    double runoff_m_now;
-    int N,i,j;
-    
-    N=num_giuh_ordinates;
-    runoff_queue_m_per_timestep[N]=0.0;
-    
-    if (subtimestep_h==1.0/12.0){
-      for(i=0;i<N;i++)
-        {
-          runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
-        }
-      
-      runoff_m_now=runoff_queue_m_per_timestep[0];
-      
-      for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
-        {
-          runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
-        }
-      runoff_queue_m_per_timestep[N-1]=0.0;
+  double runoff_m_now;
+  int N,i;
+  
+  N=num_giuh_ordinates;
+  runoff_queue_m_per_timestep[N]=0.0;
+  
+  for(i=0;i<N;i++)
+    {
+      runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
     }
-    else {//The time step is not equal to the smallest internal time step. In this case, the giuh code must be repeated j times because the time between each entry in runoff_queue_m_per_timestep is fixed.
-      double factor = (1.0/12.0)/subtimestep_h;
-      runoff_m = runoff_m*factor;
-      int resample_limit = fmax(1,(int) 1.0/factor);
-      for (j=0;j<resample_limit;j++){
-        for(i=0;i<N;i++)
-          {
-            runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
-          }
-        
-        runoff_m_now+=runoff_queue_m_per_timestep[0];
-        
-        for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
-          {
-            runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
-          }
-        runoff_queue_m_per_timestep[N-1]=0.0;
-      }
-
+  
+  runoff_m_now=runoff_queue_m_per_timestep[0];
+  
+  for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
+    {
+      runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
     }
-    
-    return runoff_m_now;
-  }
-
-  else {//fixed time step case
-    double runoff_m_now;
-    int N,i;
-    
-    N=num_giuh_ordinates;
-    runoff_queue_m_per_timestep[N]=0.0;
-    
-    for(i=0;i<N;i++)
-      {
-        runoff_queue_m_per_timestep[i]+=giuh_ordinates[i]*runoff_m;
-      }
-    
-    runoff_m_now=runoff_queue_m_per_timestep[0];
-    
-    for(i=1;i<N;i++)  // shift all the entries in preperation ffor the next timestep
-      {
-        runoff_queue_m_per_timestep[i-1]=runoff_queue_m_per_timestep[i];
-      }
-    runoff_queue_m_per_timestep[N-1]=0.0;
-    
-    return runoff_m_now;
-  }
+  runoff_queue_m_per_timestep[N-1]=0.0;
+  
+  return runoff_m_now;
 }
 
 
-#endif
+#endif
\ No newline at end of file
diff --git a/giuh/giuh.h b/giuh/giuh.h
index 88a8997..de05e17 100644
--- a/giuh/giuh.h
+++ b/giuh/giuh.h
@@ -2,7 +2,7 @@
 
 #define _GIUH_H
 
-extern double giuh_convolution_integral(int adaptive_timestep, double subtimestep_h, double runoff_m, int num_giuh_ordinates, 
+extern double giuh_convolution_integral(double runoff_m, int num_giuh_ordinates, 
                                    double *giuh_ordinates, double *runoff_queue_m_per_timestep);
 
 #endif
diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index 96d70af..d66f838 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -35,15 +35,17 @@ Initialize (std::string config_file)
 
   /* giuh ordinates are static and read in the lgar.cxx, and we need to have a copy of it to pass to
      giuh.cxx, so allocating/copying here*/
-  
+
   giuh_ordinates = new double[num_giuh_ordinates];
   giuh_runoff_queue = new double[num_giuh_ordinates+1];
 
-  for (int i=0; i<num_giuh_ordinates;i++)
+  for (int i=0; i<num_giuh_ordinates;i++){
     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++)
+  for (int i=0; i<=num_giuh_ordinates;i++){
     giuh_runoff_queue[i] = 0.0;
+  }
 
 }
 
@@ -108,7 +110,6 @@ Update()
   double volrech_subtimestep_cm;
   double surface_runoff_subtimestep_cm; // direct surface runoff
   double precip_previous_subtimestep_cm;
-  double volrunoff_giuh_subtimestep_cm;
   double volQ_gw_subtimestep_cm = 0.0; // fix it for non-zero values after adding groundwater reservoir
   
   double subtimestep_h = state->lgar_bmi_params.timestep_h;
@@ -384,16 +385,9 @@ Update()
 
 
     /*----------------------------------------------------------------------*/
-    // compute giuh runoff for the subtimestep
+    // increment runoff for the subtimestep
     surface_runoff_subtimestep_cm = volrunoff_subtimestep_cm;
-    volrunoff_giuh_subtimestep_cm = giuh_convolution_integral(adaptive_timestep, subtimestep_h, volrunoff_subtimestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue);
-
     surface_runoff_timestep_cm += surface_runoff_subtimestep_cm ;
-    volrunoff_giuh_timestep_cm += volrunoff_giuh_subtimestep_cm;
-
-    // total mass of water leaving the system, at this time it is the giuh-only, but later will add groundwater component as well.
-
-    volQ_timestep_cm += volrunoff_giuh_subtimestep_cm;
 
     // adding groundwater flux to stream channel (note: this will be updated/corrected after adding the groundwater reservoir)
     volQ_gw_timestep_cm += volQ_gw_subtimestep_cm;
@@ -441,6 +435,12 @@ Update()
 
   } // end of subcycling
 
+  //update giuh at the time step level (was previously updated at the sub time step level)
+  volrunoff_giuh_timestep_cm = giuh_convolution_integral(volrunoff_timestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue);
+
+  // total mass of water leaving the system, at this time it is the giuh-only, but later will add groundwater component as well.
+  // when groundwater component is added, it should probably happen inside of the subcycling loop.
+  volQ_timestep_cm = volrunoff_giuh_timestep_cm;
 
   /*----------------------------------------------------------------------*/
   // Everything related to lgar state is done at this point, now time to update some dynamic variables
diff --git a/src/lgar.cxx b/src/lgar.cxx
index 621a6db..28d9099 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -622,22 +622,13 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
 
 
   if (is_giuh_ordinates_set) {
-    int factor;
-    if (!state->lgar_bmi_params.adaptive_timestep){
-      factor = int(1.0/state->lgar_bmi_params.timestep_h);
-    }
-    else {
-      factor = 12;
-    }
 
-    state->lgar_bmi_params.num_giuh_ordinates = factor * (giuh_ordinates_temp.size() - 1);
+    state->lgar_bmi_params.num_giuh_ordinates = (giuh_ordinates_temp.size() - 1);
     state->lgar_bmi_params.giuh_ordinates = new double[state->lgar_bmi_params.num_giuh_ordinates+1];
     
     for (int i=0; i<giuh_ordinates_temp.size()-1; i++) {
-      for (int j=0; j<factor; j++) {
-	int index = j + i * factor + 1;
-	state->lgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1]/double(factor);
-      }
+      int index = i + 1;
+      state->lgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1];
     }
     
     if (verbosity.compare("high") == 0) {

From 8d2c5a3651fd5fe7610ced21b6787603163c9a8a Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MacBook-Pro-2.local>
Date: Tue, 11 Jun 2024 11:00:43 -0700
Subject: [PATCH 11/13] updating in response to comments, mostly focused on
 removing hardcoding from min and max time step values

---
 configs/README.md |  4 ++--
 include/all.hxx   |  3 ++-
 src/bmi_lgar.cxx  | 24 ++++++++++--------------
 src/lgar.cxx      |  6 ++++--
 4 files changed, 18 insertions(+), 19 deletions(-)

diff --git a/configs/README.md b/configs/README.md
index 30286a1..8980fd6 100644
--- a/configs/README.md
+++ b/configs/README.md
@@ -12,7 +12,7 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | initial_psi | double (scalar)| >=0 | cm | capillary head | - | used to initialize layers with a constant head |
 | ponded_depth_max | double (scalar)| >=0 | cm | maximum surface ponding | - | the maximum amount of water unavailable for surface drainage, default is set to zero |
 | timestep | double (scalar)| >0 | sec/min/hr | temporal resolution | - | timestep of the model |
-| forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data |
+| forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data. Recommended value of 3600 seconds. |
 | endtime | double (scalar)| >0 | sec, min, hr, d | simulation duration | - | time at which model simulation ends |
 | layer_soil_type | int (1D array) | - | - | state variable | - | layer soil type (read from the database file soil_params_file) |
 | max_soil_types | int | >1 | - | - | - | maximum number of soil types read from the file soil_params_file (default is set to 15) |
@@ -24,4 +24,4 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | sft_coupled | Boolean | true, false | - | model coupling | impacts hydraulic conductivity | couples LASAM to SFT. Coupling to SFT reduces hydraulic conducitivity, and hence infiltration, when soil is frozen|
 | soil_z | double (1D array) | - | cm | spatial resolution | - | vertical resolution of the soil column (computational domain of the SFT model) |
 | calib_params | Boolean | true, false | - | calibratable params flag | impacts soil properties | If set to true, soil `smcmax`, `smcmin`, `vg_n`, `vg_alpha`, `hydraulic_conductivity`, `field_capacity_psi`, and `ponded_depth_max` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content |
-| adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep, and the above timestep is ignored. If set to false, LGAR will use the above specified timestep as a fixed timestep. Testing indicates that setting this value to true substantially decreases runtime while negligibly changing the simulation. We recommend this to be set to true. |
+| adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep, and the above timestep is used as a minimum timestep (recommended value of 300 seconds). The adaptive timestep will never be larger than the forcing resolution. If set to false, LGAR will use the above specified timestep as a fixed timestep. Testing indicates that setting this value to true substantially decreases runtime while negligibly changing the simulation. We recommend this to be set to true. |
diff --git a/include/all.hxx b/include/all.hxx
index c30a96d..4708964 100755
--- a/include/all.hxx
+++ b/include/all.hxx
@@ -114,6 +114,7 @@ struct lgar_bmi_parameters
   double initial_psi_cm;           // model initial (psi) condition
   double timestep_h;               // model timestep in hours
   double forcing_resolution_h;     // forcing resolution in hours
+  double minimum_timestep_h;       // minimum time step in hours, only used if adaptive_timestep is true
   int    forcing_interval;         // = forcing_resolution_h/timestep_h
   int    num_soil_types;           // number of soil types; must be less than or equal to MAX_NUM_SOIL_TYPES
   double AET_cm;                   // actual evapotranspiration in cm
@@ -131,7 +132,7 @@ 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 */
-  bool   adaptive_timestep = false;      // if set to true, model uses adaptive timestep 
+  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.
   double ponded_depth_cm;                // amount of water on the surface unavailable for surface runoff
   double ponded_depth_max_cm;            // maximum amount of water on the surface unavailable for surface runoff
   double precip_previous_timestep_cm;    // amount of rainfall (previous time step)
diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index d66f838..d110995 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -134,24 +134,20 @@ Update()
   assert(state->lgar_bmi_input_params->PET_mm_per_h >=0.0);
 
   // adaptive time step is set 
-  if (adaptive_timestep){
-    if ( (state->lgar_bmi_input_params->precipitation_mm_per_h > 0.0) || (state->lgar_mass_balance.volon_timestep_cm > 0.0) ){
-      if (state->lgar_bmi_input_params->precipitation_mm_per_h > 10.0 || (state->lgar_mass_balance.volon_timestep_cm > 0.0) ){
-        subtimestep_h = 1.0/12.0;//case where precip > 1 cm/h, or there is ponded head from the last time step
-        state->lgar_bmi_params.timestep_h = 1.0/12.0;
-      }
-      else {
-        subtimestep_h = 1.0/6.0;//case where there is either nonzero precip or ponded head, but both are less than the threshold that requires the smallest internal time step 
-        state->lgar_bmi_params.timestep_h = 1.0/6.0;
-      }
+  if (adaptive_timestep) {
+    subtimestep_h = state->lgar_bmi_params.forcing_resolution_h;
+    if (state->lgar_bmi_input_params->precipitation_mm_per_h > 10.0 || volon_timestep_cm > 0.0 ) {
+      subtimestep_h = state->lgar_bmi_params.minimum_timestep_h;  //case where precip > 1 cm/h, or there is ponded head from the last time step
     }
-    else {//case where the precip = 0 and there is no ponded head from the last time step 
-      subtimestep_h = 1.0;
-      state->lgar_bmi_params.timestep_h = 1.0;
+    else if (state->lgar_bmi_input_params->precipitation_mm_per_h > 0.0) {
+      subtimestep_h = state->lgar_bmi_params.minimum_timestep_h * 2.0;  //case where precip is less than 1 cm/h but greater than 0, and there is no ponded head 
     }
+    if (subtimestep_h>state->lgar_bmi_params.forcing_resolution_h){//just in case the user has specified a minimum time step that would make the subtimestep_h greater than the forcing resolution 
+      subtimestep_h = state->lgar_bmi_params.forcing_resolution_h;
+    }
+    state->lgar_bmi_params.timestep_h = subtimestep_h;
   }
 
-
   state->lgar_bmi_params.forcing_interval = int(state->lgar_bmi_params.forcing_resolution_h/state->lgar_bmi_params.timestep_h+1.0e-08); // add 1.0e-08 to prevent truncation error
   subcycles = state->lgar_bmi_params.forcing_interval;
 
diff --git a/src/lgar.cxx b/src/lgar.cxx
index 28d9099..67b9397 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -400,10 +400,10 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
       continue;
     }
     else if (param_key == "adaptive_timestep") { 
-      if (param_value == "false") {
+      if ((param_value == "false") || (param_value == "0")) {
         state->lgar_bmi_params.adaptive_timestep = false;
       }
-      else if (param_value == "true") {
+      else if ( (param_value == "true") || (param_value == "1")) {
         state->lgar_bmi_params.adaptive_timestep = true;
       }
       else {
@@ -426,6 +426,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
       assert (state->lgar_bmi_params.timestep_h > 0);
       is_timestep_set = true;
 
+      state->lgar_bmi_params.minimum_timestep_h = state->lgar_bmi_params.timestep_h;
+
       if (verbosity.compare("high") == 0) {
 	std::cerr<<"Model timestep [hours,seconds]: "<<state->lgar_bmi_params.timestep_h<<" , "
 		 <<state->lgar_bmi_params.timestep_h*3600<<"\n";

From 38e8037555e7a1065946cf41ddcaac57dbe6208b Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MacBook-Pro-2.local>
Date: Tue, 11 Jun 2024 12:38:29 -0700
Subject: [PATCH 12/13] more changes in response to comments, including
 formatting of giuh.c, removing hardcoding for internal time steps and
 replacing with user specified minimum time step, and establishing a maximum
 time step equal to the forcing resolution

---
 giuh/giuh.c      |  2 +-
 src/bmi_lgar.cxx |  4 +---
 src/lgar.cxx     | 32 +++++++++++++++++++++++++++-----
 3 files changed, 29 insertions(+), 9 deletions(-)

diff --git a/giuh/giuh.c b/giuh/giuh.c
index 2605cba..6a0c3ce 100644
--- a/giuh/giuh.c
+++ b/giuh/giuh.c
@@ -39,4 +39,4 @@ extern double giuh_convolution_integral(double runoff_m,int num_giuh_ordinates,
 }
 
 
-#endif
\ No newline at end of file
+#endif
diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index d110995..0389ee4 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -142,9 +142,7 @@ Update()
     else if (state->lgar_bmi_input_params->precipitation_mm_per_h > 0.0) {
       subtimestep_h = state->lgar_bmi_params.minimum_timestep_h * 2.0;  //case where precip is less than 1 cm/h but greater than 0, and there is no ponded head 
     }
-    if (subtimestep_h>state->lgar_bmi_params.forcing_resolution_h){//just in case the user has specified a minimum time step that would make the subtimestep_h greater than the forcing resolution 
-      subtimestep_h = state->lgar_bmi_params.forcing_resolution_h;
-    }
+    subtimestep_h = fmin(subtimestep_h, state->lgar_bmi_params.forcing_resolution_h);  //just in case the user has specified a minimum time step that would make the subtimestep_h greater than the forcing resolution 
     state->lgar_bmi_params.timestep_h = subtimestep_h;
   }
 
diff --git a/src/lgar.cxx b/src/lgar.cxx
index 67b9397..5c0ca8c 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -623,25 +623,47 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
   }
 
 
+  // if (is_giuh_ordinates_set) {
+
+  //   state->lgar_bmi_params.num_giuh_ordinates = (giuh_ordinates_temp.size() - 1);
+  //   state->lgar_bmi_params.giuh_ordinates = new double[state->lgar_bmi_params.num_giuh_ordinates+1];
+    
+  //   for (int i=0; i<giuh_ordinates_temp.size()-1; i++) {
+  //     int index = i + 1;
+  //     state->lgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1];
+  //   }
+    
+  //   if (verbosity.compare("high") == 0) {
+  //     for (int i=1; i<=state->lgar_bmi_params.num_giuh_ordinates; i++)
+	// std::cerr<<"GIUH ordinates (scaled) : "<<state->lgar_bmi_params.giuh_ordinates[i]<<"\n";
+      
+  //     std::cerr<<"          *****         \n";
+  //   }
+    
+  //   giuh_ordinates_temp.clear();
+  // }
   if (is_giuh_ordinates_set) {
+    int factor = int(1.0/state->lgar_bmi_params.forcing_resolution_h);
 
-    state->lgar_bmi_params.num_giuh_ordinates = (giuh_ordinates_temp.size() - 1);
+    state->lgar_bmi_params.num_giuh_ordinates = factor * (giuh_ordinates_temp.size() - 1);
     state->lgar_bmi_params.giuh_ordinates = new double[state->lgar_bmi_params.num_giuh_ordinates+1];
     
     for (int i=0; i<giuh_ordinates_temp.size()-1; i++) {
-      int index = i + 1;
-      state->lgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1];
+      for (int j=0; j<factor; j++) {
+        int index = j + i * factor + 1;
+        state->lgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1]/double(factor);
+      }
     }
     
     if (verbosity.compare("high") == 0) {
       for (int i=1; i<=state->lgar_bmi_params.num_giuh_ordinates; i++)
-	std::cerr<<"GIUH ordinates (scaled) : "<<state->lgar_bmi_params.giuh_ordinates[i]<<"\n";
+	      std::cerr<<"GIUH ordinates (scaled) : "<<state->lgar_bmi_params.giuh_ordinates[i]<<"\n";
       
       std::cerr<<"          *****         \n";
     }
-    
     giuh_ordinates_temp.clear();
   }
+
   else if (!is_giuh_ordinates_set) {
     stringstream errMsg;
     errMsg << "The configuration file \'" << config_file <<"\' does not set giuh_ordinates. \n";

From b01969c0a51cc06a2fd6d9b52dec629e7ba98f2a Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MacBook-Pro-2.local>
Date: Tue, 11 Jun 2024 13:53:58 -0700
Subject: [PATCH 13/13] removing code that was commented out

---
 src/bmi_lgar.cxx |  2 +-
 src/lgar.cxx     | 20 --------------------
 2 files changed, 1 insertion(+), 21 deletions(-)

diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index 0389ee4..b9f92ba 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -149,7 +149,7 @@ Update()
   state->lgar_bmi_params.forcing_interval = int(state->lgar_bmi_params.forcing_resolution_h/state->lgar_bmi_params.timestep_h+1.0e-08); // add 1.0e-08 to prevent truncation error
   subcycles = state->lgar_bmi_params.forcing_interval;
 
-  if (verbosity.compare("high") == 0) {//and adaptive time step is engaged? 
+  if (verbosity.compare("high") == 0) {
     printf("time step size in hours: %lf \n", state->lgar_bmi_params.timestep_h);
   }
   
diff --git a/src/lgar.cxx b/src/lgar.cxx
index 5c0ca8c..815f37f 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -622,26 +622,6 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
     throw runtime_error(errMsg.str());
   }
 
-
-  // if (is_giuh_ordinates_set) {
-
-  //   state->lgar_bmi_params.num_giuh_ordinates = (giuh_ordinates_temp.size() - 1);
-  //   state->lgar_bmi_params.giuh_ordinates = new double[state->lgar_bmi_params.num_giuh_ordinates+1];
-    
-  //   for (int i=0; i<giuh_ordinates_temp.size()-1; i++) {
-  //     int index = i + 1;
-  //     state->lgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1];
-  //   }
-    
-  //   if (verbosity.compare("high") == 0) {
-  //     for (int i=1; i<=state->lgar_bmi_params.num_giuh_ordinates; i++)
-	// std::cerr<<"GIUH ordinates (scaled) : "<<state->lgar_bmi_params.giuh_ordinates[i]<<"\n";
-      
-  //     std::cerr<<"          *****         \n";
-  //   }
-    
-  //   giuh_ordinates_temp.clear();
-  // }
   if (is_giuh_ordinates_set) {
     int factor = int(1.0/state->lgar_bmi_params.forcing_resolution_h);