From 31bee3d3fb2e272e5d8c13712f74d385986af7d4 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:15:24 -0600
Subject: [PATCH 01/19] feat: add listDelete for linked list

---
 include/all.hxx     |  2 +-
 src/linked_list.cxx | 14 ++++++++++++++
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/include/all.hxx b/include/all.hxx
index 6b7526c..8df67f3 100755
--- a/include/all.hxx
+++ b/include/all.hxx
@@ -239,7 +239,7 @@ extern void                     listReverseOrder(struct wetting_front** head_ref
 extern bool                     listFindLayer(struct wetting_front* link, int num_layers, double *cum_layer_thickness_cm,
 					      int *lives_in_layer, bool *extends_to_bottom_flag);
 extern struct wetting_front*    listCopy(struct wetting_front* current, struct wetting_front* state_previous=NULL);
-
+extern void listDelete(struct wetting_front* head);
 
 
 
diff --git a/src/linked_list.cxx b/src/linked_list.cxx
index 926f68c..c63dbdd 100755
--- a/src/linked_list.cxx
+++ b/src/linked_list.cxx
@@ -26,6 +26,20 @@
 //___________________________________________________________
 
 
+/*###########################################################*/
+/* listDelete() - deletes memory allocated to a linked list  */
+/* This function must be called on any list to deallocate    */
+/* the dynamic memory used in creating and manipultating the */
+/* list. (added by NJF)                                                     */
+/*###########################################################*/
+extern void listDelete(struct wetting_front* head)
+{
+  if( head != NULL ){
+    struct wetting_front *next = head->next;
+    free( head );
+    listDelete(next);
+  }
+}
 
 /*#########################################################*/
 /* listPrint() - prints a linked list to screen             */

From e099fa78539c404998cb702092f099346578f9ba Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:22:36 -0600
Subject: [PATCH 02/19] fix: deallocate memory with listDelete before
 overwriting with listCopy

---
 src/bmi_lgar.cxx | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index c11a8aa..bff2bb4 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -143,7 +143,10 @@ Update()
       std::cerr<<"BMI Update |Timesteps = "<< state->lgar_bmi_params.timesteps<<", Time [h] = "<<this->state->lgar_bmi_params.time_s / 3600.<<", Subcycle = "<< cycle <<" of "<<subcycles<<std::endl;
     }
 
-    state->state_previous = NULL;
+    if( state->state_previous != NULL ){
+      listDelete(state->state_previous);
+      state->state_previous = NULL;
+    }
     state->state_previous = listCopy(state->head);
 
     // ensure precip and PET are non-negative
@@ -270,7 +273,10 @@ Update()
         listPrint(state->head);
       }
 
-      state->state_previous = NULL;
+      if(state->state_previous != NULL ){
+        listDelete(state->state_previous);
+        state->state_previous = NULL;
+      }
       state->state_previous = listCopy(state->head);
 
       volin_timestep_cm += volin_subtimestep_cm;

From 093a8a1013789174a1c5cb2c1131df79e8020c74 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:18:09 -0600
Subject: [PATCH 03/19] fix: add destructor to deallocate class held pointers

---
 include/bmi_lgar.hxx | 1 +
 src/bmi_lgar.cxx     | 8 ++++++++
 2 files changed, 9 insertions(+)

diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx
index 606f59d..abcec7c 100644
--- a/include/bmi_lgar.hxx
+++ b/include/bmi_lgar.hxx
@@ -29,6 +29,7 @@ class NotImplemented : public std::logic_error {
 
 class BmiLGAR : public bmi::Bmi {
 public:
+  ~BmiLGAR();
   BmiLGAR() {
     this->input_var_names[0] = "precipitation_rate";
     this->input_var_names[1] = "potential_evapotranspiration_rate";
diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index bff2bb4..805a798 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -17,6 +17,14 @@
 // default verbosity is set to 'none' other option 'high' or 'low' needs to be specified in the config file
 string verbosity="none";
 
+/**
+ * @brief Delete dynamic arrays allocated in Initialize() and held by this object
+ * 
+ */
+BmiLGAR::~BmiLGAR(){
+  if( giuh_ordinates != nullptr ) delete [] giuh_ordinates;
+  if( giuh_runoff_queue != nullptr ) delete [] giuh_runoff_queue;
+}
 
 /* The `head` pointer stores the address in memory of the first member of the linked list containing
    all the wetting fronts. The contents of struct wetting_front are defined in "all.h" */

From 997a2b79a46187d8f9362a8e04ca9e7a52daaf77 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:20:20 -0600
Subject: [PATCH 04/19] fix: ensure cleanup of old memory when allocating soil
 params

---
 include/bmi_lgar.hxx |  1 +
 src/bmi_lgar.cxx     | 17 +++++++++++++++--
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx
index abcec7c..741b419 100644
--- a/include/bmi_lgar.hxx
+++ b/include/bmi_lgar.hxx
@@ -132,6 +132,7 @@ public:
   struct model_state* get_model();
   
 private:
+  void realloc_soil();
   struct model_state* state;
   static const int input_var_name_count  = 3;
   static const int output_var_name_count = 15;
diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index 805a798..6533791 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -56,6 +56,20 @@ Initialize (std::string config_file)
 
 }
 
+/**
+ * @brief Allocate (or reallocate) storage for soil parameters
+ * 
+ */
+void BmiLGAR::realloc_soil(){
+  if(state->lgar_bmi_params.soil_depth_wetting_fronts != nullptr)
+    delete [] state->lgar_bmi_params.soil_depth_wetting_fronts;
+  if(state->lgar_bmi_params.soil_moisture_wetting_fronts != nullptr)
+    delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts;
+
+  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];
+}
+
 /*
   This is the main function calling lgar subroutines for creating, moving, and merging wetting fronts.
   Calls to AET and mass balance module are also happening here
@@ -436,8 +450,7 @@ Update()
   state->lgar_bmi_params.num_wetting_fronts = listLength(state->head);
 
   // allocate new memory based on updated wetting fronts; we could make it conditional i.e. create only if no. of wf are changed
-  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];
+  realloc_soil();
 
   // update thickness/depth and soil moisture of wetting fronts (used for state coupling)
   struct wetting_front *current = state->head;

From f9f302b0eba72e32d965abce0c6e6bc74a34bd14 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:25:09 -0600
Subject: [PATCH 05/19] fix: delete memory when deleting front from list

---
 src/linked_list.cxx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/linked_list.cxx b/src/linked_list.cxx
index c63dbdd..3782d12 100755
--- a/src/linked_list.cxx
+++ b/src/linked_list.cxx
@@ -254,7 +254,7 @@ extern struct wetting_front* listDeleteFront(int front_num, struct wetting_front
     previous = current->next;
 
   }
-
+  if( current != NULL ) free( current );
   current = previous;
 
   while(previous != NULL) { // decrement all front numbers

From 0d1f4938a8a80004bb4eb7d8a34b9a9308a3d29d Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:27:04 -0600
Subject: [PATCH 06/19] fix: delete tempory heap allocated variables in
 bmi_main

---
 src/bmi_main_lgar.cxx | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/bmi_main_lgar.cxx b/src/bmi_main_lgar.cxx
index 1c64c93..37b0ac9 100644
--- a/src/bmi_main_lgar.cxx
+++ b/src/bmi_main_lgar.cxx
@@ -155,6 +155,8 @@ int main(int argc, char *argv[])
       // write layers data to file
       fprintf(outlayer_fptr,"# Timestep = %d, %s \n", i, time[i].c_str());
       write_state(outlayer_fptr, model_state.get_model()->head);
+      delete [] soil_moisture_wetting_front;
+      delete [] soil_thickness_wetting_front;
     }
 
   }

From 335332253f4ee5c66e66f797639f1575901b21fd Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:33:07 -0600
Subject: [PATCH 07/19] fix: delete temporary dynamic allocated arrays

---
 src/lgar.cxx | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/lgar.cxx b/src/lgar.cxx
index fd412a3..55d5330 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -1134,7 +1134,9 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
       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;
-      
+      //done with delta_thetas and delta_thickness, cleanup memory
+      delete delta_thetas;
+      delete delta_thickness;
       current->theta = fmax(theta_r, fmin(theta_new, theta_e));
 
       double Se = calc_Se_from_theta(current->theta,theta_e,theta_r);

From bb5684e08c94a922fe6f2ec4a2ae037962dbb585 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:28:59 -0600
Subject: [PATCH 08/19] fix: clean up list memory prior to listCopy in main.cxx

---
 src/main.cxx | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/main.cxx b/src/main.cxx
index 5c9e586..8f15113 100755
--- a/src/main.cxx
+++ b/src/main.cxx
@@ -435,6 +435,7 @@ for(time_step_num=0;time_step_num<num_time_steps;time_step_num++)
       printf(" ************ time step %d ******** \n",time_step_num);
     }
     
+    if( state_previous != NULL ) listDelete(state_previous);
     state_previous = NULL;
     state_previous = listCopy(head);
     
@@ -542,6 +543,7 @@ for(time_step_num=0;time_step_num<num_time_steps;time_step_num++)
 	lgar_create_surfacial_front(&ponded_depth_cm, &volin_timestep_cm, dry_depth, theta1, soil_type_by_layer, soil_properties, cum_layer_thickness_cm, nint, time_step_h);
 	
 	//listPrint();
+  if( state_previous != NULL ) listDelete(state_previous);
 	state_previous = NULL;
 	state_previous = listCopy(head);
 	

From 7f633857ee3aec8dd23dfc6b8a53c086589951f0 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:29:55 -0600
Subject: [PATCH 09/19] fix: cleanup dynamic memory in bmi state during
 Finalize

---
 src/bmi_lgar.cxx | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index 6533791..317a14e 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -610,6 +610,30 @@ void BmiLGAR::
 Finalize()
 {
   global_mass_balance();
+  if( state->head != NULL ) listDelete(state->head);
+  if( state->state_previous != NULL ) listDelete(state->state_previous);
+
+  if( state->soil_properties != nullptr ) delete [] state->soil_properties;
+
+  if( state->lgar_bmi_params.soil_depth_wetting_fronts != nullptr ) delete [] state->lgar_bmi_params.soil_depth_wetting_fronts;
+  if( state->lgar_bmi_params.soil_moisture_wetting_fronts != nullptr) delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts;
+
+  if( state->lgar_bmi_params.soil_temperature != nullptr ) delete [] state->lgar_bmi_params.soil_temperature;
+  if( state->lgar_bmi_params.soil_temperature_z != nullptr)  delete [] state->lgar_bmi_params.soil_temperature_z;
+  if( state->lgar_bmi_params.layer_soil_type != nullptr ) delete [] state->lgar_bmi_params.layer_soil_type;
+
+  if( state->lgar_calib_params.theta_e != nullptr ) delete [] state->lgar_calib_params.theta_e;
+  if( state->lgar_calib_params.theta_r != nullptr ) delete [] state->lgar_calib_params.theta_r;
+  if( state->lgar_calib_params.vg_n != nullptr ) delete [] state->lgar_calib_params.vg_n;
+  if( state->lgar_calib_params.vg_alpha != nullptr ) delete [] state->lgar_calib_params.vg_alpha;
+  if( state->lgar_calib_params.Ksat != nullptr ) delete [] state->lgar_calib_params.Ksat;
+
+  if( state->lgar_bmi_params.layer_thickness_cm != nullptr ) delete [] state->lgar_bmi_params.layer_thickness_cm;
+  if( state->lgar_bmi_params.cum_layer_thickness_cm != nullptr ) delete [] state->lgar_bmi_params.cum_layer_thickness_cm;
+  if( state->lgar_bmi_params.giuh_ordinates != nullptr ) delete [] state->lgar_bmi_params.giuh_ordinates;
+  if( state->lgar_bmi_params.frozen_factor != nullptr ) delete [] state->lgar_bmi_params.frozen_factor;
+  if( state->lgar_bmi_input_params != nullptr ) delete state->lgar_bmi_input_params;
+  if( state != nullptr ) delete state;
 }
 
 

From 06a3c2e7160f8459f9888c9e540e0dbc61b759ec Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:27:36 -0600
Subject: [PATCH 10/19] fix: call Finalize on bmi model

---
 src/bmi_main_lgar.cxx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/bmi_main_lgar.cxx b/src/bmi_main_lgar.cxx
index 37b0ac9..26c0ed8 100644
--- a/src/bmi_main_lgar.cxx
+++ b/src/bmi_main_lgar.cxx
@@ -163,7 +163,7 @@ int main(int argc, char *argv[])
 
   // do final mass balance
   model_state.global_mass_balance();
-
+  model_state.Finalize();
   if (outdata_fptr) {
     fclose(outdata_fptr);
     fclose(outlayer_fptr);

From 18580d76e749f44fc653830c373526c88dde944a Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 19:32:03 -0600
Subject: [PATCH 11/19] chore: inline comments on things that don't make sense

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

diff --git a/src/lgar.cxx b/src/lgar.cxx
index 55d5330..900a3f8 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -643,6 +643,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
     }
   }
   else {
+    //NJF FIXME these arrays should be allocated based on num_cells_temp...
     state->lgar_bmi_params.soil_temperature   = new double[1]();
     state->lgar_bmi_params.soil_temperature_z = new double[1]();
     state->lgar_bmi_params.num_cells_temp     = 1;
@@ -1117,6 +1118,8 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
 	double theta_below = 0.0;
 
 	new_mass += layer_thickness * (theta - theta_below);
+  //NJF theta_below is always 0, so all delta_thetas are always 0...
+  //does this really need a dynamic array in this case???
 	delta_thetas[k] = theta_below;
 	delta_thickness[k] = layer_thickness;
       }

From 1914cd80f66b747db45e88ecbbbcb68e9f8a62b1 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 12 Jun 2024 20:51:20 -0600
Subject: [PATCH 12/19] fix: free instead of delete for malloc'd mem

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

diff --git a/src/lgar.cxx b/src/lgar.cxx
index 900a3f8..491a8fb 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -1138,8 +1138,8 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
 						 delta_thetas, delta_thickness, soil_type, soil_properties);
       actual_ET_demand = *AET_demand_cm;
       //done with delta_thetas and delta_thickness, cleanup memory
-      delete delta_thetas;
-      delete delta_thickness;
+      free(delta_thetas);
+      free(delta_thickness);
       current->theta = fmax(theta_r, fmin(theta_new, theta_e));
 
       double Se = calc_Se_from_theta(current->theta,theta_e,theta_r);

From e652bb2582e391adee7e4147a2de6c674c664fc3 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Thu, 13 Jun 2024 08:39:34 -0600
Subject: [PATCH 13/19] fix: free another set of tempory allocated pointers

---
 src/lgar.cxx | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/lgar.cxx b/src/lgar.cxx
index 491a8fb..e4a899d 100755
--- a/src/lgar.cxx
+++ b/src/lgar.cxx
@@ -1281,7 +1281,9 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
 	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;
-
+  //done with delta_thetas and delta_thickness, cleanup memory
+  free(delta_thetas);
+  free(delta_thickness);
 	current->theta = fmax(theta_r, fmin(theta_new, theta_e));
 
       }

From 8605e5e6bb10182695c42bf0d01c74c9bed619a6 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Thu, 13 Jun 2024 12:27:33 -0600
Subject: [PATCH 14/19] fix: explicitly intialize class pointers to nullptr

---
 include/bmi_lgar.hxx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx
index 741b419..b29339a 100644
--- a/include/bmi_lgar.hxx
+++ b/include/bmi_lgar.hxx
@@ -30,7 +30,7 @@ class NotImplemented : public std::logic_error {
 class BmiLGAR : public bmi::Bmi {
 public:
   ~BmiLGAR();
-  BmiLGAR() {
+  BmiLGAR():giuh_ordinates(nullptr), giuh_runoff_queue(nullptr) {
     this->input_var_names[0] = "precipitation_rate";
     this->input_var_names[1] = "potential_evapotranspiration_rate";
     this->input_var_names[2] = "soil_temperature_profile";

From a046a340de3468a11b0e8885f8be7a46600a7563 Mon Sep 17 00:00:00 2001
From: Nels <hellkite500@gmail.com>
Date: Wed, 10 Jul 2024 13:16:21 -0700
Subject: [PATCH 15/19] chore: change recursive call to loop in listDelete

Co-authored-by: Phil Miller - NOAA <philip.miller@noaa.gov>
---
 src/linked_list.cxx | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/linked_list.cxx b/src/linked_list.cxx
index 3782d12..0176d41 100755
--- a/src/linked_list.cxx
+++ b/src/linked_list.cxx
@@ -34,10 +34,10 @@
 /*###########################################################*/
 extern void listDelete(struct wetting_front* head)
 {
-  if( head != NULL ){
+  while (head != NULL) {
     struct wetting_front *next = head->next;
     free( head );
-    listDelete(next);
+    head = next;
   }
 }
 

From 025b4c70e8d6f9390ec0b0f91e315b91b9ec2dfe Mon Sep 17 00:00:00 2001
From: Nels <hellkite500@gmail.com>
Date: Wed, 10 Jul 2024 13:20:13 -0700
Subject: [PATCH 16/19] chore: remove redundant NULL check

Co-authored-by: Phil Miller - NOAA <philip.miller@noaa.gov>
---
 src/main.cxx | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/main.cxx b/src/main.cxx
index 8f15113..c3e51dd 100755
--- a/src/main.cxx
+++ b/src/main.cxx
@@ -435,8 +435,7 @@ for(time_step_num=0;time_step_num<num_time_steps;time_step_num++)
       printf(" ************ time step %d ******** \n",time_step_num);
     }
     
-    if( state_previous != NULL ) listDelete(state_previous);
-    state_previous = NULL;
+    listDelete(state_previous);
     state_previous = listCopy(head);
     
     if(debug_flag) {

From fb08aeac22b52c6636e78b19cbfc1412ae0c4fe6 Mon Sep 17 00:00:00 2001
From: PeterLaFollette-NOAA <59893312+peterlafollette@users.noreply.github.com>
Date: Wed, 26 Jun 2024 11:31:29 -0700
Subject: [PATCH 17/19] Ptl add adaptive timestep (#25)

* adding the capability for an adaptive time step. If not provided in the config file, it defaults to false.

* updating readme

* updating readme

* updating giuh fxn with adaptive time step

* 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.

* updating readme in configs

* cleaned up code that was used for debugging

* updating Phillipsburg config file

* cleaned up code that was used for debugging

* 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.

* updating in response to comments, mostly focused on removing hardcoding from min and max time step values

* 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

* removing code that was commented out

---------

Co-authored-by: Peter La Follette <peterlafollette@Peters-MBP-2.lan>
Co-authored-by: Peter La Follette <peterlafollette@Peters-MacBook-Pro-2.local>
---
 configs/README.md                     |  3 +-
 configs/config_lasam_Bushland.txt     |  1 +
 configs/config_lasam_Phillipsburg.txt |  1 +
 configs/config_lasam_sft_ngen.txt     |  1 +
 include/all.hxx                       |  2 ++
 src/bmi_lgar.cxx                      | 49 +++++++++++++++++++--------
 src/lgar.cxx                          | 37 +++++++++++++++-----
 src/soil_funcs.cxx                    |  2 +-
 8 files changed, 72 insertions(+), 24 deletions(-)

diff --git a/configs/README.md b/configs/README.md
index 60ea2c4..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,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, 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/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/include/all.hxx b/include/all.hxx
index 8df67f3..e829422 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,6 +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. 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 317a14e..e3fae68 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -41,18 +41,19 @@ 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*/
-  
+
   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;
+  }
 
 }
 
@@ -98,7 +99,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)
@@ -131,7 +132,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;
@@ -139,6 +139,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!)
@@ -153,6 +154,26 @@ 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) {
+    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 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 
+    }
+    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;
+  }
+
+  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) {
+    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++) {
@@ -386,16 +407,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(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;
@@ -443,6 +457,13 @@ 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 e4a899d..f5a497d 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") || (param_value == "0")) {
+        state->lgar_bmi_params.adaptive_timestep = false;
+      }
+      else if ( (param_value == "true") || (param_value == "1")) {
+        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);
 
@@ -411,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";
@@ -605,29 +622,28 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
     throw runtime_error(errMsg.str());
   }
 
-
   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.forcing_resolution_h);
 
     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++) {
       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 = 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";
@@ -1277,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;
@@ -2580,7 +2595,8 @@ 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;
@@ -2595,6 +2611,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 3fab3e8c58cfd0c329fbaf1ce504c5a3c21d9d17 Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 10 Jul 2024 14:22:46 -0600
Subject: [PATCH 18/19] fix: remove extra mass balance check at end of bmi_main
 test

---
 src/bmi_main_lgar.cxx | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/bmi_main_lgar.cxx b/src/bmi_main_lgar.cxx
index 26c0ed8..a106c8b 100644
--- a/src/bmi_main_lgar.cxx
+++ b/src/bmi_main_lgar.cxx
@@ -161,8 +161,7 @@ int main(int argc, char *argv[])
 
   }
 
-  // do final mass balance
-  model_state.global_mass_balance();
+  // do final mass balance ( inside Finalize() ) and finish the simulation
   model_state.Finalize();
   if (outdata_fptr) {
     fclose(outdata_fptr);

From d52d6f0dd041b4c5d90004a3ff8285cfc28afd6c Mon Sep 17 00:00:00 2001
From: hellkite500 <nfrazier@lynker.com>
Date: Wed, 10 Jul 2024 14:31:35 -0600
Subject: [PATCH 19/19] chore: remove defensive checks for nullptr

---
 src/bmi_lgar.cxx | 59 ++++++++++++++++++++++++------------------------
 1 file changed, 29 insertions(+), 30 deletions(-)

diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx
index e3fae68..f1b744c 100644
--- a/src/bmi_lgar.cxx
+++ b/src/bmi_lgar.cxx
@@ -22,8 +22,8 @@ string verbosity="none";
  * 
  */
 BmiLGAR::~BmiLGAR(){
-  if( giuh_ordinates != nullptr ) delete [] giuh_ordinates;
-  if( giuh_runoff_queue != nullptr ) delete [] giuh_runoff_queue;
+  delete [] giuh_ordinates;
+  delete [] giuh_runoff_queue;
 }
 
 /* The `head` pointer stores the address in memory of the first member of the linked list containing
@@ -62,10 +62,9 @@ Initialize (std::string config_file)
  * 
  */
 void BmiLGAR::realloc_soil(){
-  if(state->lgar_bmi_params.soil_depth_wetting_fronts != nullptr)
-    delete [] state->lgar_bmi_params.soil_depth_wetting_fronts;
-  if(state->lgar_bmi_params.soil_moisture_wetting_fronts != nullptr)
-    delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts;
+  
+  delete [] state->lgar_bmi_params.soil_depth_wetting_fronts;
+  delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts;
 
   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];
@@ -631,30 +630,30 @@ void BmiLGAR::
 Finalize()
 {
   global_mass_balance();
-  if( state->head != NULL ) listDelete(state->head);
-  if( state->state_previous != NULL ) listDelete(state->state_previous);
-
-  if( state->soil_properties != nullptr ) delete [] state->soil_properties;
-
-  if( state->lgar_bmi_params.soil_depth_wetting_fronts != nullptr ) delete [] state->lgar_bmi_params.soil_depth_wetting_fronts;
-  if( state->lgar_bmi_params.soil_moisture_wetting_fronts != nullptr) delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts;
-
-  if( state->lgar_bmi_params.soil_temperature != nullptr ) delete [] state->lgar_bmi_params.soil_temperature;
-  if( state->lgar_bmi_params.soil_temperature_z != nullptr)  delete [] state->lgar_bmi_params.soil_temperature_z;
-  if( state->lgar_bmi_params.layer_soil_type != nullptr ) delete [] state->lgar_bmi_params.layer_soil_type;
-
-  if( state->lgar_calib_params.theta_e != nullptr ) delete [] state->lgar_calib_params.theta_e;
-  if( state->lgar_calib_params.theta_r != nullptr ) delete [] state->lgar_calib_params.theta_r;
-  if( state->lgar_calib_params.vg_n != nullptr ) delete [] state->lgar_calib_params.vg_n;
-  if( state->lgar_calib_params.vg_alpha != nullptr ) delete [] state->lgar_calib_params.vg_alpha;
-  if( state->lgar_calib_params.Ksat != nullptr ) delete [] state->lgar_calib_params.Ksat;
-
-  if( state->lgar_bmi_params.layer_thickness_cm != nullptr ) delete [] state->lgar_bmi_params.layer_thickness_cm;
-  if( state->lgar_bmi_params.cum_layer_thickness_cm != nullptr ) delete [] state->lgar_bmi_params.cum_layer_thickness_cm;
-  if( state->lgar_bmi_params.giuh_ordinates != nullptr ) delete [] state->lgar_bmi_params.giuh_ordinates;
-  if( state->lgar_bmi_params.frozen_factor != nullptr ) delete [] state->lgar_bmi_params.frozen_factor;
-  if( state->lgar_bmi_input_params != nullptr ) delete state->lgar_bmi_input_params;
-  if( state != nullptr ) delete state;
+  listDelete(state->head);
+  listDelete(state->state_previous);
+
+  delete [] state->soil_properties;
+
+  delete [] state->lgar_bmi_params.soil_depth_wetting_fronts;
+  delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts;
+
+  delete [] state->lgar_bmi_params.soil_temperature;
+  delete [] state->lgar_bmi_params.soil_temperature_z;
+  delete [] state->lgar_bmi_params.layer_soil_type;
+
+  delete [] state->lgar_calib_params.theta_e;
+  delete [] state->lgar_calib_params.theta_r;
+  delete [] state->lgar_calib_params.vg_n;
+  delete [] state->lgar_calib_params.vg_alpha;
+  delete [] state->lgar_calib_params.Ksat;
+
+  delete [] state->lgar_bmi_params.layer_thickness_cm;
+  delete [] state->lgar_bmi_params.cum_layer_thickness_cm;
+  delete [] state->lgar_bmi_params.giuh_ordinates;
+  delete [] state->lgar_bmi_params.frozen_factor;
+  delete state->lgar_bmi_input_params;
+  delete state;
 }