After chatting with Ahmad, we determined that it is best if the adapt…

…ive 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

@@ -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

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

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

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) {