diff --git a/build/CMakeLists.txt b/build/CMakeLists.txt
index 00d0e9add17356c6dda8c921b5a397031fd814ec..337979ed7b01c81d88349558190a0e869f1e8198 100644
--- a/build/CMakeLists.txt
+++ b/build/CMakeLists.txt
@@ -46,6 +46,7 @@ if (CMAKE_BUILD_TYPE MATCHES Cluster)
         ${PARENT_DIR}/build/includes/summa_actor 
         ${PARENT_DIR}/build/includes/job_actor 
         ${PARENT_DIR}/build/includes/file_access_actor 
+        ${PARENT_DIR}/build/includes/gru_actor
         ${PARENT_DIR}/build/includes/hru_actor)
     set(LIB_ACTORS 
         ${CAF_LIBRARIES} 
@@ -63,6 +64,7 @@ else()
         "${PARENT_DIR}/build/includes/summa_actor"
         "${PARENT_DIR}/build/includes/job_actor"
         "${PARENT_DIR}/build/includes/file_access_actor"
+        "${PARENT_DIR}/build/includes/gru_actor"
         "${PARENT_DIR}/build/includes/hru_actor")
     link_directories("/usr/local/lib")
     set(LIB_ACTORS CAF::core CAF::io -lopenblas -lnetcdff -lnetcdf)
@@ -81,6 +83,7 @@ set(ACTORS_DIR      ${PARENT_DIR}/build/source)
 set(SYS_INIT_DIR    ${ACTORS_DIR}/system_initialization)
 set(FILE_ACCESS_DIR ${ACTORS_DIR}/file_access_actor)
 set(JOB_ACTOR_DIR   ${ACTORS_DIR}/job_actor)
+set(GRU_ACTOR_DIR   ${ACTORS_DIR}/gru_actor)
 set(HRU_ACTOR_DIR   ${ACTORS_DIR}/hru_actor)
 
 # NOAHMP modules
@@ -212,11 +215,14 @@ set(FILE_ACCESS_INTERFACE
     ${FILE_ACCESS_DIR}/fileAccess_writeOutput.f90)
 set(JOB_INTERFACE
     ${JOB_ACTOR_DIR}/gru_struc.f90)
+set(GRU_INTERFACE
+    ${GRU_ACTOR_DIR}/gru_interface.f90)
 set(HRU_INTERFACE
     ${HRU_ACTOR_DIR}/hru_init.f90
     ${HRU_ACTOR_DIR}/hru_read.f90
     ${HRU_ACTOR_DIR}/hru_modelRun.f90
-    ${HRU_ACTOR_DIR}/hru_writeOutput.f90)
+    ${HRU_ACTOR_DIR}/hru_writeOutput.f90
+    ${HRU_ACTOR_DIR}/hru_interface.f90)
 
 # Actors actual actor modules
 set(ACTORS_GLOBAL
@@ -250,6 +256,8 @@ set(JOB_ACTOR
     ${JOB_ACTOR_DIR}/job_utils.cpp
     ${JOB_ACTOR_DIR}/distributed_job_actor.cpp
     ${JOB_ACTOR_DIR}/node_actor.cpp)
+set(GRU_ACTOR
+    ${GRU_ACTOR_DIR}/gru_actor.cpp)
 set(HRU_ACTOR
     ${HRU_ACTOR_DIR}/hru_utils.cpp
     ${HRU_ACTOR_DIR}/hru_actor.cpp
@@ -278,6 +286,7 @@ set(SUMMA_ALL
     ${SYS_INIT_INTERFACE}
     ${FILE_ACCESS_INTERFACE}
     ${JOB_INTERFACE}
+    ${GRU_INTERFACE}
     ${HRU_INTERFACE})
 
 set(MAIN_ACTOR ${ACTORS_DIR}/main.cpp)
@@ -318,6 +327,7 @@ add_executable(${EXEC_NAME}
                ${ACTORS_GLOBAL}
                ${FILE_ACCESS_ACTOR}
                ${JOB_ACTOR}
+               ${GRU_ACTOR}
                ${HRU_ACTOR}
                ${SYS_INIT}
                ${SUMMA_CLIENT}
diff --git a/build/includes/file_access_actor/file_access_actor.hpp b/build/includes/file_access_actor/file_access_actor.hpp
index 95b259728e8486751d9349dcbcaba4b064a484cc..bb13ff968f2f350633a5fdce435bcbc32a1e41e3 100644
--- a/build/includes/file_access_actor/file_access_actor.hpp
+++ b/build/includes/file_access_actor/file_access_actor.hpp
@@ -41,6 +41,7 @@ struct file_access_state {
   caf::actor parent; 
   int start_gru;
   int num_gru;
+  int num_hru;
 
   NumGRUInfo num_gru_info;
 
diff --git a/build/includes/global/auxilary.hpp b/build/includes/global/auxilary.hpp
index 9c29e84c3dbb4c9144e871c45b926af0d0dc6c47..a95069aecc51112269f77de499c40cfdccb03d60 100644
--- a/build/includes/global/auxilary.hpp
+++ b/build/includes/global/auxilary.hpp
@@ -42,7 +42,8 @@ std::vector<long int> get_var_i8_by_indx(void* handle, int struct_indx);
 std::vector<int> get_flagVec_by_indx(void* handle, int struct_indx);
 
 void get_scalar_data(void* handle, double fracJulDay, 
-    double tmZoneOffsetFracDay, int year_length, int computeVegFlux);
+                     double tmZoneOffsetFracDay, int year_length, 
+                     int computeVegFlux, double dt_init, double upArea);
 
 // Setters by index
 void set_var_dlength_by_indx(void* handle, 
@@ -63,4 +64,5 @@ void set_flagVec_by_indx(void* handle, std::vector<int>& summa_struct,
     int struct_indx);
 
 void set_scalar_data(void* handle, double fracJulDay, 
-    double tmZoneOffsetFracDay, int year_length, int computeVegFlux);
\ No newline at end of file
+                     double tmZoneOffsetFracDay, int year_length, 
+                     int computeVegFlux, double dt_init, double upArea);
\ No newline at end of file
diff --git a/build/includes/global/fortran_data_types.hpp b/build/includes/global/fortran_data_types.hpp
index 2eb562e72afb3ef6ba3e3b053d3bc354504dac40..6c0dfd9a9c27e8beaaf4ce1b62a1168445c75784 100644
--- a/build/includes/global/fortran_data_types.hpp
+++ b/build/includes/global/fortran_data_types.hpp
@@ -117,6 +117,9 @@ extern "C" {
     // hru_type
     void* new_handle_hru_type();
     void delete_handle_hru_type(void* handle);
+    // gru_type
+    void* new_handle_gru_type(int& num_hru);
+    void delete_handle_gru_type(void* handle);
 
     // var_dlength_by_indx
     void get_size_var_dlength_by_indx(void* handle, int* indx, int* num_var);
@@ -161,7 +164,11 @@ extern "C" {
 
     // scalar types
     void get_scalar_data_fortran(void* handle, double* fracJulDay, 
-        double* tmZoneOffsetFracDay, int* year_length, int* computeVegFlux);
+                                 double* tmZoneOffsetFracDay, int* year_length, 
+                                 int* computeVegFlux, double* dt_init,
+                                 double* upArea);
     void set_scalar_data_fortran(void* handle, double* fracJulDay,
-        double* tmZoneOffsetFracDay, int* year_length, int* computeVegFlux);
+                                 double* tmZoneOffsetFracDay, int* year_length, 
+                                 int* computeVegFlux, double* dt_init,
+                                 double* upArea);
 }
\ No newline at end of file
diff --git a/build/includes/gru_actor/gru_actor.hpp b/build/includes/gru_actor/gru_actor.hpp
index 8111d5b484112e85c389de39c5157c0cffe82b2f..a4eb0818cb7992c39ab83c54a4274a8fdca2642b 100644
--- a/build/includes/gru_actor/gru_actor.hpp
+++ b/build/includes/gru_actor/gru_actor.hpp
@@ -1,51 +1,58 @@
 #pragma once
-
 #include "caf/all.hpp"
 #include "fortran_data_types.hpp"
-#include "timing_info.hpp"
-#include "global.hpp"
+#include "settings_functions.hpp"
+#include "hru_actor.hpp"
 #include <vector>
-#include <array>
-#include <chrono>
-#include <string>
-#include "var_lookup.hpp"
-
-namespace caf{
-struct gru_state {
-    // GRU Initalization
-    int ref_gru;
-    int indx_gru;
-    std::string config_path;
-    caf::actor file_access_actor;
-    int output_struc_size;
-    caf::actor parent;
-        
-    // HRU information
-    std::vector<caf::actor> hru_list;
-    int num_hrus;
-    int hrus_complete = 0;
-
-    // Global Data
-    int nTimeDelay = 2000; // number of hours in the time delay histogram (default: ~1 season = 24*365/4)
-    struct iLookVarType var_type_lookup;
-    
-    // Data Structure local to the GRU
-    int num_bpar_vars;                               // number of variables in the fortran structure for bpar_struct
-    std::vector<double> bpar_struct;   
-    std::vector<int> bpar_struct_var_type_list;  // The types to the variable at each element in bpar_struct
-
-    int num_bvar_vars;                               // number of variables in the fortran structure for bvar_struct
-    std::vector<std::vector<double>> bvar_struct;
-    std::vector<int> bvar_struct_var_type_list;
-
-    int num_var_types;
-    std::vector<int> i_look_var_type_list; // The types to the variable at each element in bvar_struct
+
+extern "C" {
+  void getNumHRU(int& index_gru, int& num_hru);
+
+  void initGRU_fortran(int& index_gru, void* gru_data, int& err, void* message);
+  void setupGRU_fortran(int& index_gru, void* gru_data, int& err, 
+                        void* message);
+  void readGRURestart_fortran(int& index_gru, void* gru_data, int& err, 
+                              void* message);
+  void setTimeZoneOffsetGRU_fortran(int& iFile, void* gru_data, int& err, 
+                                    void* message);
+  void readGRUForcing_fortran(int& index_gru, int& iStep, int& iRead, 
+                              int& iFile, void* gru_data, int& err, 
+                              void* message);
+  void runGRU_fortran(int& index_gru, int& timestep, void* gru_data, 
+                      int& dt_init_factor, int& err, void* message);
+
+  void writeGRUOutput_fortran(int& index_gru, int& timestep, int& output_step, 
+                              void* gru_data, int& err, void* message);
+
+}
+
+struct gru_actor_state {
+  int netcdf_index;
+  int gru_job_index;
+  HRU_Actor_Settings hru_actor_settings;
+  caf::actor file_access_actor;
+  caf::actor parent;
+
+  int num_hrus;
+  std::vector<void*> hrus;
+  void* bvar_stat = new_handle_var_dlength();
+  void* bvar_struct = new_handle_var_dlength();
+
+  void* gru_data;
+
+  double dt_init = 0.0;
+  int dt_init_factor = 1;
+  int num_steps_until_write;
+  int num_steps = 0;                    // number of time steps
+  int timestep = 1;	                   // Current Timestep of HRU simulation
+  int iFile = 1;
+	int stepsInCurrentFFile;             // number of time steps in current forcing file
+  int forcingStep = 1;                 // index of current time step in current forcing file
+  int output_structure_step_index = 1; // index of current time step in output structure
+
 };
 
-behavior gru_actor(stateful_actor<gru_state>* self, 
-    int ref_gru, 
-    int indx_gru, 
-    std::string config_path, 
-    caf::actor file_access_actor, 
-    caf::actor parent);
-}
\ No newline at end of file
+caf::behavior gru_actor(caf::stateful_actor<gru_actor_state>* self, 
+                        int netcdf_index, int gru_job_index, int num_steps,
+                        HRU_Actor_Settings hru_actor_settings,
+                        caf::actor file_access_actor, caf::actor parent);
\ No newline at end of file
diff --git a/build/includes/hru_actor/hru.hpp b/build/includes/hru_actor/hru.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c121d43f852ed45ad3b5245d78d9bc37853c899e
--- /dev/null
+++ b/build/includes/hru_actor/hru.hpp
@@ -0,0 +1,17 @@
+#pragma once
+#include "fortran_data_types.hpp"
+
+class HRU {
+
+  private:
+    void* hru_data_;
+    int indx_hru_;
+    int indx_gru_;
+    int netcdf_index_;
+    int timestep_;
+    int forcing_step_;
+    int num_steps_; // TODO: May not need this
+    
+};
+
+
diff --git a/build/includes/hru_actor/hru_actor.hpp b/build/includes/hru_actor/hru_actor.hpp
index 408c603194a5114a2c1e612b10e6fb29b5024191..5581ef99f8f303d9b49c20c83c334d4bbcd37033 100644
--- a/build/includes/hru_actor/hru_actor.hpp
+++ b/build/includes/hru_actor/hru_actor.hpp
@@ -23,30 +23,35 @@ struct Date {
  *********************************************/
 extern "C" {
   // Initialize HRU data_structures
-  void initHRU(int* indxGRU, int* num_steps, void* hru_data, int* err);
+  void initHRU_fortran(int& indxGRU, int& indx_hru, int& num_steps, 
+                       void* hru_data, int& err, void* message);
   
-  void setupHRUParam(int* indxGRU, int* indxHRU, void* hru_data, 
-                     double* upArea, int* err);
+  void setupHRU_fortran(int& indxGRU, int& indxHRU, void* hru_data, int& err, 
+                        void* message);
   
   // Setup summa_readRestart File if this option has been chosen 
-  void summa_readRestart(int* indxGRU, int* indxHRU, void* hru_data, 
-                         double* dtInit, int* err);
-  
+  void readHRURestart_fortran(int& indxGRU, int& indxHRU, void* hru_data, 
+                              int& err, void* message);
+                         
+  void readHRUForcing_fortran(int& index_gru, int& indx_hru, int& iStep, 
+                              int& iRead, int& iFile, void* hru_data,  
+                              int& err, void* message);
   // Run the model for one timestep
-  void RunPhysics(int* id, int* stepIndex, void* hru_data, double* dt, 
-                  int* dt_int_factor, double* walltime_timestep, int* err);
+  void runHRU_fortran(int& indx_gru, int& indx_hru, int& timestep, void* hru_data, 
+                      int& dt_int_factor, double& walltime_timestep, 
+                      int& err, void* message);
   
-  void hru_writeOutput(int* index_hru, int* index_gru, int* timestep, 
-                       int* output_step, void* hru_data, int* y, int* m, 
-                       int* d, int* h, int* err);
+  void writeHRUOutput_fortran(int& index_gru, int& index_hru, int& timestep, 
+                              int& output_step, void* hru_data, int& y, int& m, 
+                              int& d, int& h, int& err, void* message);
     
   int  hru_writeRestart(int* index_hru, int* index_gru, int* timestep, 
                         int* output_step, void* hru_data, int* err);
 
-  void setTimeZoneOffset(int* iFile, void* hru_data, int* err);
+  void setTimeZoneOffset_fortran(int& iFile, void* hru_data, int& err, 
+                                 void* message);
+
 
-  void HRU_readForcing(int* index_gru, int* iStep, int* iRead, int* iFile, 
-                       void* hru_data,  int* err);
   
   // hru_writeOutput.f90
   void setFinalizeStatsFalse(int* indx_gru);
@@ -93,9 +98,7 @@ struct hru_state {
   int     iFile = 1;              // index of current forcing file from forcing file list
   int     dt_init_factor = 1; // factor of dt_init (coupled_em)
   int     output_structure_step_index = 1; // index of current time step in output structure
-  double  dt_init;            // used to initialize the length of the sub-step for each HRU
-  double	upArea;             // area upslope of each HRU
-
+  
   // Sundials variables
   double rtol = -9999; // -9999 uses default
   double atol = -9999; // -9999 uses default
diff --git a/build/includes/hru_actor/hru_utils.hpp b/build/includes/hru_actor/hru_utils.hpp
index 0275d60d5e9d405fdbe55eac2325588468896c5f..1ec8926c642b406924a02f2eafc0bcb41f59ba1c 100644
--- a/build/includes/hru_actor/hru_utils.hpp
+++ b/build/includes/hru_actor/hru_utils.hpp
@@ -14,8 +14,6 @@ struct hru {
   int iFile;
   int dt_init_factor;
   int output_structure_step_index;
-  double dt_init;
-  double upArea;
 
   // Sundials variables
   double rtol;
@@ -61,6 +59,8 @@ struct hru {
   double tm_zone_offset_frac_day;
   int year_length;
   int compute_veg_flux;
+  double dt_init;
+  double upArea;
 };
 
 template <class Inspector>
@@ -76,8 +76,6 @@ bool inspect(Inspector& inspector, hru& hru_data) {
          inspector.field("dt_init_factor", hru_data.dt_init_factor),
          inspector.field("output_structure_step_index", 
                          hru_data.output_structure_step_index),
-         inspector.field("dt_init", hru_data.dt_init),
-         inspector.field("upArea", hru_data.upArea),
          inspector.field("rtol", hru_data.rtol),
          inspector.field("atol", hru_data.atol),
          inspector.field("forc_stat", hru_data.forc_stat),
@@ -111,5 +109,7 @@ bool inspect(Inspector& inspector, hru& hru_data) {
          inspector.field("tm_zone_offset_frac_day", 
                          hru_data.tm_zone_offset_frac_day),
          inspector.field("year_length", hru_data.year_length),
-         inspector.field("compute_veg_flux", hru_data.compute_veg_flux));
+         inspector.field("compute_veg_flux", hru_data.compute_veg_flux),
+         inspector.field("dt_init", hru_data.dt_init),
+         inspector.field("upArea", hru_data.upArea));
 }
\ No newline at end of file
diff --git a/build/includes/job_actor/GRU.hpp b/build/includes/job_actor/GRU.hpp
index 69dc87275da50f368443772996922c4d086c8485..6ada68925e4904351c0ec7b1c3f8115fdbe491b8 100644
--- a/build/includes/job_actor/GRU.hpp
+++ b/build/includes/job_actor/GRU.hpp
@@ -11,6 +11,8 @@ class GRU {
     int index_job_;          // The index of the GRU within this job
     caf::actor actor_ref_;   // The actor for the GRU
 
+    int num_hrus_;           // The number of HRUs in the GRU
+
     // Modifyable Parameters
     int dt_init_factor_;     // The initial dt for the GRU
     double rel_tol_;         // The relative tolerance for the GRU
diff --git a/build/includes/job_actor/gru_struc.hpp b/build/includes/job_actor/gru_struc.hpp
index 1d9daf49f14f0eed97aa41deecb27e0c0ad7a7a4..1966ec02b8089453ef1bbc0776b59f29c842f3e2 100644
--- a/build/includes/job_actor/gru_struc.hpp
+++ b/build/includes/job_actor/gru_struc.hpp
@@ -9,6 +9,8 @@ extern "C" {
 
   void read_icond_nlayers_fortran(int& num_gru, int& err, void* message);
 
+  void get_num_hru_per_gru_fortran(int& arr_size, int& num_hru_per_gru_array);
+
   void deallocate_gru_struc_fortran();
 }
 
@@ -22,6 +24,7 @@ class GruStruc {
     inline int getStartGru() const { return start_gru_; }
     inline int getNumGrus() const { return num_gru_; }
     inline int get_file_gru() const { return file_gru_; }
+    inline int getNumHrus() const { return num_hru_; }
     inline int get_gru_info_size() const { return gru_info_.size(); }
     inline int getNumGrusDone() const { return num_gru_done_; }
     inline int getNumGRUFailed() const { return num_gru_failed_; }
@@ -49,6 +52,9 @@ class GruStruc {
       return -1;
     }
 
+    void getNumHrusPerGru();
+    inline int getNumHruPerGru(int index) { return num_hru_per_gru_[index]; }
+
   private:
     // Inital Information about the GRUs
     int start_gru_;
@@ -59,6 +65,7 @@ class GruStruc {
     
     // GRU specific Information
     std::vector<std::unique_ptr<GRU>> gru_info_;
+    std::vector<int> num_hru_per_gru_;
 
     // Runtime status of the GRUs
     int num_gru_done_ = 0;
diff --git a/build/includes/job_actor/job_actor.hpp b/build/includes/job_actor/job_actor.hpp
index f69871c4cbbe4fa1ea7aafa77325cad7112d3608..a6c96c22d78986e0872e199ad45156d88d6af29b 100644
--- a/build/includes/job_actor/job_actor.hpp
+++ b/build/includes/job_actor/job_actor.hpp
@@ -18,6 +18,7 @@
 #include <vector>
 #include <tuple>
 #include "summa_init_struc.hpp"
+#include "gru_actor.hpp"
 
 
 
diff --git a/build/source/file_access_actor/fileAccess_writeOutput.f90 b/build/source/file_access_actor/fileAccess_writeOutput.f90
index bff869754962029a0eaaf7212453c0e42ad2e51b..2054183f9b0b3bb00fcd940010c986964de96a96 100644
--- a/build/source/file_access_actor/fileAccess_writeOutput.f90
+++ b/build/source/file_access_actor/fileAccess_writeOutput.f90
@@ -107,7 +107,7 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, write
   integer(c_int),intent(out)           :: err               ! Error code
   ! local variables
   type(var_i),pointer                  :: ncid
-  integer(i4b)                         :: iGRU              ! loop through GRUs
+  integer(i4b)                         :: iGRU,iHRU         ! loop through GRUs
   integer(i4b)                         :: iStep             ! loop through time steps
   integer(i4b)                         :: iFreq             ! loop through output frequencies
   integer(i4b)                         :: indxHRU=1         ! index of HRU to write
@@ -125,17 +125,19 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, write
   if (write_parm_flag)then
     do iStruct=1,size(structInfo)
       do iGRU=start_gru, max_gru
-        select case(trim(structInfo(iStruct)%structName))
-        case('attr'); call writeParm(ncid,gru_struc(iGRU)%hruInfo(indxHRU)%hru_ix, &
-          summa_struct(1)%attrStruct%gru(iGRU)%hru(indxHRU),attr_meta,err,cmessage)
-        case('type'); call writeParm(ncid,gru_struc(iGRU)%hruInfo(indxHRU)%hru_ix, &
-          summa_struct(1)%typeStruct%gru(iGRU)%hru(indxHRU),type_meta,err,cmessage)
-        case('mpar'); call writeParm(ncid,gru_struc(iGRU)%hruInfo(indxHRU)%hru_ix, &
-          summa_struct(1)%mparStruct%gru(iGRU)%hru(indxHRU),mpar_meta,err,cmessage)
-        end select
-        if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
-        call writeParm(ncid,iGRU,summa_struct(1)%bparStruct%gru(iGRU),bpar_meta,err,cmessage)
-        if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
+        do iHRU=1,size(gru_struc(iGRU)%hruInfo)
+          select case(trim(structInfo(iStruct)%structName))
+          case('attr'); call writeParm(ncid,gru_struc(iGRU)%hruInfo(iHRU)%hru_ix, &
+            summa_struct(1)%attrStruct%gru(iGRU)%hru(iHRU),attr_meta,err,cmessage)
+          case('type'); call writeParm(ncid,gru_struc(iGRU)%hruInfo(iHRU)%hru_ix, &
+            summa_struct(1)%typeStruct%gru(iGRU)%hru(iHRU),type_meta,err,cmessage)
+          case('mpar'); call writeParm(ncid,gru_struc(iGRU)%hruInfo(iHRU)%hru_ix, &
+            summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU),mpar_meta,err,cmessage)
+          end select
+          if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
+          call writeParm(ncid,iGRU,summa_struct(1)%bparStruct%gru(iGRU),bpar_meta,err,cmessage)
+          if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
+        end do ! HRU
       end do ! GRU
     end do ! structInfo
   end if
@@ -160,8 +162,8 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, write
   ! *** write basin data
   ! ****************************************************************************
   call writeBasin(ncid,outputTimeStep(start_gru)%dat(:),outputTimeStepUpdate,num_steps,&
-                  start_gru, max_gru, numGRU, &
-                  bvar_meta,summa_struct(1)%bvarStat,summa_struct(1)%bvarStruct, &
+                  start_gru, max_gru, numGRU, bvar_meta, &
+                  summa_struct(1)%bvarStat,summa_struct(1)%bvarStruct, &
                   bvarChild_map,err,cmessage)
 
   ! ****************************************************************************
@@ -399,6 +401,7 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,nSteps,
   integer(i4b)                     :: iStep
   integer(i4b)                     :: iGRU
   real(rkind)                      :: val
+  integer(i4b)                     :: nHRUrun
   ! initialize error control
   err=0;message="writeData/"
   ! loop through output frequencies
@@ -413,14 +416,10 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,nSteps,
         ! get variable index
         err = nf90_inq_varid(ncid%var(iFreq),trim(meta(iVar)%varName),ncVarID)
         call netcdf_err(err,message); if (err/=0) return
-        
         do iStep = 1, nSteps
-
           if(.not.summa_struct(1)%finalizeStats%gru(minGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
-
           stepCounter = stepCounter+1
           timeVec(stepCounter) = summa_struct(1)%forcStruct%gru(minGRU)%hru(1)%var(iVar)%tim(iStep)
-
         end do ! iStep
         err = nf90_put_var(ncid%var(iFreq),ncVarID,timeVec(1:stepCounter),start=(/outputTimestep(iFreq)/),count=(/stepCounter/))
         call netcdf_err(err,message); if (err/=0)then; print*, "err"; return; endif
@@ -429,6 +428,13 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,nSteps,
         cycle
       end if  ! id time
 
+      ! Calculate the number of HRUs to write
+      nHRUrun = 0
+      do iGRU=minGRU, maxGRU
+        nHRUrun = nHRUrun + size(gru_struc(iGRU)%hruInfo)
+      end do ! iGRU
+
+
       ! define the statistics index
       iStat = meta(iVar)%statIndex(iFreq)
       ! check that the variable is desired
@@ -437,11 +443,11 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,nSteps,
         ! stats output: only scalar variable type
         if(meta(iVar)%varType==iLookVarType%scalarv) then
           call writeScalar(ncid, outputTimeStep, outputTimeStepUpdate, nSteps, &
-                           minGRU, maxGRU, numGRU, iFreq, iVar, meta, stat,    &
+                           minGRU, maxGRU, nHRUrun, iFreq, iVar, meta, stat,   &
                            map, err, message)
         else ! non-scalar variables: regular data structures
           call writeVector(ncid, outputTimeStep, maxLayers, nSteps, minGRU, &
-                           maxGRU, numGRU, iFreq, iVar, meta, dat, indx,    &
+                           maxGRU, nHRUrun, iFreq, iVar, meta, dat, indx,   &
                            err, message)
         end if ! not scalarv
 
@@ -455,7 +461,7 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,nSteps,
 end subroutine writeData
 
 subroutine writeScalar(ncid, outputTimestep, outputTimestepUpdate, nSteps, minGRU, maxGRU, &
-  numGRU, iFreq, iVar, meta, stat, map, err, message)
+  nHRUrun, iFreq, iVar, meta, stat, map, err, message)
   USE data_types,only:var_info                       ! metadata type
   USE, intrinsic :: ieee_arithmetic
   implicit none
@@ -466,7 +472,7 @@ subroutine writeScalar(ncid, outputTimestep, outputTimestepUpdate, nSteps, minGR
   integer(i4b)  ,intent(in)         :: nSteps                  ! number of timeSteps
   integer(i4b)  ,intent(in)         :: minGRU                  ! minGRU index to write
   integer(i4b)  ,intent(in)         :: maxGRU                  ! maxGRU index to write - probably not needed
-  integer(i4b)  ,intent(in)         :: numGRU
+  integer(i4b)  ,intent(in)         :: nHRUrun
   integer(i4b)  ,intent(in)         :: iFreq                   ! output file index (year, month, day, timesteps)
   integer(i4b)  ,intent(in)         :: iVar                    ! netcdf variable we are writing data for
   type(var_info),intent(in)         :: meta(:)                 ! meta data
@@ -477,11 +483,12 @@ subroutine writeScalar(ncid, outputTimestep, outputTimestepUpdate, nSteps, minGR
 
   ! local variables
   integer(i4b)                      :: gruCounter=0             ! counter for the realVecs
+  integer(i4b)                      :: hru_counter=0
   integer(i4b)                      :: iStep=1                  ! counter for looping over nSteps
   integer(i4b)                      :: stepCounter=0            ! counter for the realVec
-  integer(i4b)                      :: iGRU
+  integer(i4b)                      :: iGRU,iHRU
   ! output array
-  real(rkind)                       :: realVec(numGRU, nSteps)! real vector for all HRUs in the run domain
+  real(rkind)                       :: realVec(nHRUrun, nSteps)! real vector for all HRUs in the run domain
   real(rkind)                       :: val
 
   err=0; message="writeOutput.f90-writeScalar/"
@@ -489,20 +496,19 @@ subroutine writeScalar(ncid, outputTimestep, outputTimestepUpdate, nSteps, minGR
   select type(stat)
     class is (gru_hru_time_doubleVec)
       gruCounter=0
+      hru_counter=0
       do iGRU = minGRU, maxGRU
-        stepCounter = 0
-        gruCounter = gruCounter + 1
-        do iStep = 1, nSteps
-          
-          if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
-
-
-          stepCounter = stepCounter + 1
-
-          realVec(gruCounter, stepCounter) = stat%gru(iGRU)%hru(1)%var(map(iVar))%tim(iStep)%dat(iFreq)
-
-          outputTimeStepUpdate(iFreq) = stepCounter
-        end do ! iStep
+          ! gruCounter = gruCounter + 1
+        do iHRU = 1, size(gru_struc(iGRU)%hruInfo)
+          hru_counter = hru_counter + 1
+          stepCounter = 0
+          do iStep = 1, nSteps
+            if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(iHRU)%tim(iStep)%dat(iFreq)) cycle
+            stepCounter = stepCounter + 1
+            realVec(hru_counter, stepCounter) = stat%gru(iGRU)%hru(iHRU)%var(map(iVar))%tim(iStep)%dat(iFreq)
+            outputTimeStepUpdate(iFreq) = stepCounter
+          end do ! iStep
+        end do ! iHRU
       end do ! iGRU 
 
       if (outputTimeStepUpdate(iFreq) /= stepCounter ) then
@@ -515,19 +521,22 @@ subroutine writeScalar(ncid, outputTimestep, outputTimestepUpdate, nSteps, minGR
         print*, "   maxGRU = ", maxGRU
         print*, "   nSteps = ", nSteps
         print*, "   gruCounter = ", gruCounter
+        print*, "   hru_counter = ", hru_counter
         print*, "   iStep = ", iStep
         err = 20
         return
       endif
-
-      err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realVec(1:gruCounter, 1:stepCounter),start=(/minGRU,outputTimestep(iFreq)/),count=(/numGRU,stepCounter/))
+      err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),&
+                         realVec(1:hru_counter, 1:stepCounter),    &
+                         start=(/minGRU,outputTimestep(iFreq)/),   & 
+                         count=(/nHRUrun,stepCounter/))
     class default; err=20; message=trim(message)//'stats must be scalarv and of type gru_hru_doubleVec'; return
   end select  ! stat
 
 end subroutine
 
 subroutine writeVector(ncid, outputTimestep, maxLayers, nSteps, minGRU, maxGRU, &
-  numGRU, iFreq, iVar, meta, dat, indx, err, message)
+  nHRUrun, iFreq, iVar, meta, dat, indx, err, message)
   USE data_types,only:var_info                       ! metadata type
   USE var_lookup,only:iLookIndex                     ! index into index structure
   USE var_lookup,only:iLookVarType                   ! index into type structure
@@ -538,7 +547,7 @@ subroutine writeVector(ncid, outputTimestep, maxLayers, nSteps, minGRU, maxGRU,
   integer(i4b)  ,intent(in)             :: nSteps                  ! number of timeSteps
   integer(i4b)  ,intent(in)             :: minGRU                  ! minGRU index to write
   integer(i4b)  ,intent(in)             :: maxGRU                  ! maxGRU index to write - probably not needed
-  integer(i4b)  ,intent(in)             :: numGRU
+  integer(i4b)  ,intent(in)             :: nHRUrun
   integer(i4b)  ,intent(in)             :: iFreq                   ! output file index (year, month, day, timesteps)
   integer(i4b)  ,intent(in)             :: iVar                    ! netcdf variable we are writing data for
   type(var_info),intent(in)             :: meta(:)                 ! meta data
@@ -549,9 +558,10 @@ subroutine writeVector(ncid, outputTimestep, maxLayers, nSteps, minGRU, maxGRU,
 
   ! local variables
   integer(i4b)                          :: gruCounter             ! counter for the realVecs
+  integer(i4b)                          :: hru_counter
   integer(i4b)                          :: iStep                  ! counter for looping over nSteps
   integer(i4b)                          :: stepCounter            ! counter for the realVec
-  integer(i4b)                          :: iGRU
+  integer(i4b)                          :: iGRU,iHRU
   integer(i4b)                          :: nSoil
   integer(i4b)                          :: nSnow
   integer(i4b)                          :: nLayers
@@ -561,8 +571,8 @@ subroutine writeVector(ncid, outputTimestep, maxLayers, nSteps, minGRU, maxGRU,
   integer(i4b)                          :: dataType          ! type of data
   integer(i4b),parameter                :: ixInteger=1001    ! named variable for integer
   integer(i4b),parameter                :: ixReal=1002       ! named variable for real
-  real(rkind)                           :: realArray(numGRU,maxLayers+1)  ! real array for all HRUs in the run domain
-  integer(i4b)                          :: intArray(numGRU,maxLayers+1)   ! integer array for all HRUs in the run domain
+  real(rkind)                           :: realArray(nHRUrun,maxLayers+1)  ! real array for all HRUs in the run domain
+  integer(i4b)                          :: intArray(nHRUrun,maxLayers+1)   ! integer array for all HRUs in the run domain
   err=0; message="writeOutput.f90-writeVector/"
 
   ! initialize the data vectors
@@ -571,65 +581,66 @@ subroutine writeVector(ncid, outputTimestep, maxLayers, nSteps, minGRU, maxGRU,
     class is (gru_hru_time_intVec);     intArray(:,:) = integerMissing; dataType=ixInteger
     class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return
   end select
-
   ! Loop over GRUs
-  
   stepCounter = outputTimeStep(iFreq)
   do iStep = 1, nSteps
     gruCounter = 1
+    hru_counter = 1
     do iGRU = minGRU, maxGRU
-      ! get the model layers
-      nSoil   = indx%gru(iGRU)%hru(1)%var(iLookIndex%nSoil)%tim(iStep)%dat(1)
-      nSnow   = indx%gru(iGRU)%hru(1)%var(iLookIndex%nSnow)%tim(iStep)%dat(1)
-      nLayers = indx%gru(iGRU)%hru(1)%var(iLookIndex%nLayers)%tim(iStep)%dat(1)
+      do iHRU=1, size(gru_struc(iGRU)%hruInfo)
+        ! get the model layers
+        nSoil   = indx%gru(iGRU)%hru(iHRU)%var(iLookIndex%nSoil)%tim(iStep)%dat(1)
+        nSnow   = indx%gru(iGRU)%hru(iHRU)%var(iLookIndex%nSnow)%tim(iStep)%dat(1)
+        nLayers = indx%gru(iGRU)%hru(iHRU)%var(iLookIndex%nLayers)%tim(iStep)%dat(1)
+
+        ! get the length of each data vector
+        select case (meta(iVar)%varType)
+            case(iLookVarType%wLength); datLength = maxSpectral
+            case(iLookVarType%midToto); datLength = nLayers
+            case(iLookVarType%midSnow); datLength = nSnow
+            case(iLookVarType%midSoil); datLength = nSoil
+            case(iLookVarType%ifcToto); datLength = nLayers+1
+            case(iLookVarType%ifcSnow); datLength = nSnow+1
+            case(iLookVarType%ifcSoil); datLength = nSoil+1
+            case default; cycle
+        end select ! vartype
+
+        ! get the data vectors
+        select type (dat)
+            class is (gru_hru_time_doubleVec)
+                if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(iHRU)%tim(iStep)%dat(iFreq)) cycle
+                realArray(hru_counter,1:datLength) = dat%gru(iGRU)%hru(iHRU)%var(iVar)%tim(iStep)%dat(1:datLength)
+
+            class is (gru_hru_time_intVec)
+                if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(iHRU)%tim(iStep)%dat(iFreq)) cycle
+                intArray(hru_counter,1:datLength) = dat%gru(iGRU)%hru(iHRU)%var(iVar)%tim(iStep)%dat(1:datLength)
+            class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return
+        end select
 
-      ! get the length of each data vector
-      select case (meta(iVar)%varType)
-          case(iLookVarType%wLength); datLength = maxSpectral
-          case(iLookVarType%midToto); datLength = nLayers
-          case(iLookVarType%midSnow); datLength = nSnow
-          case(iLookVarType%midSoil); datLength = nSoil
-          case(iLookVarType%ifcToto); datLength = nLayers+1
-          case(iLookVarType%ifcSnow); datLength = nSnow+1
-          case(iLookVarType%ifcSoil); datLength = nSoil+1
+        ! get the maximum length of each data vector
+        select case (meta(iVar)%varType)
+          case(iLookVarType%wLength); maxLength = maxSpectral
+          case(iLookVarType%midToto); maxLength = maxLayers
+          case(iLookVarType%midSnow); maxLength = maxLayers-nSoil
+          case(iLookVarType%midSoil); maxLength = nSoil
+          case(iLookVarType%ifcToto); maxLength = maxLayers+1
+          case(iLookVarType%ifcSnow); maxLength = (maxLayers-nSoil)+1
+          case(iLookVarType%ifcSoil); maxLength = nSoil+1
           case default; cycle
-      end select ! vartype
-
-      ! get the data vectors
-      select type (dat)
-          class is (gru_hru_time_doubleVec)
-              if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
-              realArray(gruCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(1:datLength)
-
-          class is (gru_hru_time_intVec)
-              if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
-              intArray(gruCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(1:datLength)
-          class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return
-      end select
-
-      ! get the maximum length of each data vector
-      select case (meta(iVar)%varType)
-        case(iLookVarType%wLength); maxLength = maxSpectral
-        case(iLookVarType%midToto); maxLength = maxLayers
-        case(iLookVarType%midSnow); maxLength = maxLayers-nSoil
-        case(iLookVarType%midSoil); maxLength = nSoil
-        case(iLookVarType%ifcToto); maxLength = maxLayers+1
-        case(iLookVarType%ifcSnow); maxLength = (maxLayers-nSoil)+1
-        case(iLookVarType%ifcSoil); maxLength = nSoil+1
-        case default; cycle
-      end select ! vartype
-      gruCounter = gruCounter + 1
+        end select ! vartype
+        hru_counter = hru_counter + 1
+      end do ! iHRU
     end do ! iGRU
 
    ! write the data vectors
     select case(dataType)
 
       case(ixReal)
-        err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realArray(1:numGRU,1:maxLength),start=(/minGRU,1,stepCounter/),count=(/numGRU,maxLength,1/))
+        err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realArray(1:nHRUrun,1:maxLength),start=(/minGRU,1,stepCounter/),count=(/nHRUrun,maxLength,1/))
         if(err/=0)then; print*, "ERROR: with nf90_put_var in data vector (ixReal)"; return; endif
         realArray(:,:) = realMissing ! reset the realArray
       case(ixInteger)
-        err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),intArray(1:numGRU,1:maxLength),start=(/minGRU,1,stepCounter/),count=(/numGRU,maxLength,1/))
+        err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),intArray(1:nHRUrun,1:maxLength),start=(/minGRU,1,stepCounter/),count=(/nHRUrun,maxLength,1/))
         if(err/=0)then; print*, "ERROR: with nf90_put_var in data vector (ixInteger)"; return; endif
         intArray(:,:) = integerMissing ! reset the intArray
       case default; err=20; message=trim(message)//'data must be of type integer or real'; return
@@ -670,37 +681,57 @@ subroutine writeBasin(ncid,outputTimestep,outputTimestepUpdate,nSteps,&
   integer(i4b)                  :: iVar              ! variable index
   integer(i4b)                  :: iStat             ! statistics index
   integer(i4b)                  :: iFreq             ! frequency index
+  integer(i4b)                  :: step_counter
+  integer(i4b)                  :: gru_counter
+  integer(i4b)                  :: iGRU, iStep
+  real(rkind)                   :: realVec(numGRU, nSteps)! real vector for all HRUs in the run domain
+
   ! initialize error control
   err=0;message="f-writeBasin/"
 
-
   ! loop through output frequencies
   do iFreq=1,maxvarFreq
-
     ! skip frequencies that are not needed
     if(.not.outFreq(iFreq)) cycle
-
     ! loop through model variables
     do iVar = 1,size(meta)
-
       ! define the statistics index
       iStat = meta(iVar)%statIndex(iFreq)
-
       ! check that the variable is desired
       if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle
-
-      ! stats/data output - select data type
       select case (meta(iVar)%varType)
-
         case (iLookVarType%scalarv)
-          call writeScalar(ncid, outputTimeStep, outputTimeStepUpdate, nSteps, &
-                           minGRU, maxGRU, numGRU, iFreq, iVar, meta, stat, map, &
-                           err, message)
-
+          select type (stat)
+            class is (gru_hru_time_doubleVec)
+              gru_counter = 0
+              do iGRU = minGRU, maxGRU
+                step_counter = 0
+                gru_counter = gru_counter + 1
+                do iStep = 1, nSteps
+                  step_counter = step_counter + 1
+                  if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+                  realVec(gru_counter, step_counter) = stat%gru(iGRU)%hru(1)%var(map(iVar))%tim(iStep)%dat(iFreq)
+                  outputTimeStepUpdate(iFreq) = step_counter
+                end do ! iStep
+              end do ! iGRU
+              err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq), &
+                                 realVec(1:numGRU,1:step_counter),           &
+                                 start=(/minGRU,outputTimestep(iFreq)/),     & 
+                                 count=(/numGRU,step_counter/))
+            class default; err=20; message=trim(message)//'stats must be scalarv and of type gru_hru_doubleVec'; return
+          end select ! stat
         case (iLookVarType%routing)
-          if (iFreq==1 .and. outputTimestep(iFreq)==1) then
-            ! err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),(/dat(iVar)%tim(iStep)%dat/),start=(/1/),count=(/1000/))
-          end if
+          select type (dat)
+            class is (gru_hru_time_doubleVec)
+              if (iFreq==1 .and. outputTimestep(iFreq)==1) then
+                do iGRU = minGRU, maxGRU
+                  err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq), &
+                                     dat%gru(iGRU)%hru(1)%var(iVar)%tim(1)%dat,        &
+                                     start=(/1/), count=(/1000/))
+                end do
+              end if
+            class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return
+          end select
         case default
           err=40; message=trim(message)//"unknownVariableType[name='"//trim(meta(iVar)%varName)//"';type='"//trim(get_varTypeName(meta(iVar)%varType))//    "']"; return
       end select ! variable type
diff --git a/build/source/file_access_actor/file_access_actor.cpp b/build/source/file_access_actor/file_access_actor.cpp
index eca8711f8860f8f3ee2f5d217111b5346026771a..ca932df713c04448710d5f1fb7d6b2051d358cf3 100644
--- a/build/source/file_access_actor/file_access_actor.cpp
+++ b/build/source/file_access_actor/file_access_actor.cpp
@@ -32,10 +32,10 @@ behavior file_access_actor(stateful_actor<file_access_state>* self,
   self->state.num_output_steps = fa_settings.num_timesteps_in_output_buffer;
   
   return {
-    [=](init_file_access_actor, int file_gru) {
+    [=](init_file_access_actor, int file_gru, int num_hru) {
       aout(self) << "File Access Actor: Intializing\n";
       auto& fa_settings = self->state.file_access_actor_settings;
-      int num_hru = self->state.num_gru;
+      self->state.num_hru = num_hru;
       
       // Get the information about the forcing files
       self->state.forcing_files = std::make_unique<forcingFileContainer>();
@@ -72,7 +72,7 @@ behavior file_access_actor(stateful_actor<file_access_state>* self,
         actor_address = "_" + to_string(self->address());
       }
       defOutputFortran(self->state.handle_ncid, &self->state.start_gru, 
-                       &self->state.num_gru, &num_hru, &file_gru, 
+                       &self->state.num_gru, &self->state.num_hru, &file_gru, 
                        &self->state.num_gru_info.use_global_for_data_structures,
                        actor_address.c_str(), &err);
       if (err != 0) return -1;
diff --git a/build/source/file_access_actor/forcing_file_info.f90 b/build/source/file_access_actor/forcing_file_info.f90
index 9cc76e9721d3c8bfff220e522e6eae7310c0d248..c957d78c18d1c5078874968209c78049421627ec 100644
--- a/build/source/file_access_actor/forcing_file_info.f90
+++ b/build/source/file_access_actor/forcing_file_info.f90
@@ -2,7 +2,6 @@ module forcing_file_info
   USE, intrinsic :: iso_c_binding
   USE nrtype
   USE globalData,only:integerMissing   ! integer missing value
-  USE actor_data_types,only:var_forc   ! global data structure for forcing data
   USE data_types,only:dlength          ! global data structure for forcing data
   implicit none
   public::getNumFrocingFiles_fortran
@@ -12,6 +11,23 @@ module forcing_file_info
   public::openForcingFile
   public::freeForcingFiles_fortran
 
+    
+  type,public :: forcingFileData
+    real(rkind), dimension (:,:), allocatable   :: dataFromFile  ! (hru, time)
+  end type forcingFileData
+
+  type,public :: var_forc
+    type(forcingFileData), allocatable   :: var(:)       ! var(:)%dataFromFile(:,:)
+    character(len=256)                   :: refTimeString
+    real(rkind)                          :: convTime2Days
+    integer(i4b)                         :: nVars
+    integer(i4b),allocatable             :: var_ix(:)
+    real(rkind)                          :: tmZoneOffsetFracDay
+    real(rkind)                          :: refJulDay_data 
+    integer(i4b)                         :: nTimeSteps    ! Number of Timesteps in the file
+  end type var_forc
+
+
   ! Module Data Structures (global to hrus that import this module)
   type(var_forc),allocatable,save,public  :: forcingDataStruct(:) ! forcingDataStruct(:)%var(:)%dataFromFile(:,:)
   type(dlength),allocatable,save,public   :: vecTime(:)
@@ -102,6 +118,7 @@ subroutine read_forcingFile(iFile, startGRU, numGRU, err, message_r) &
   USE globalData,only:forc_meta 
   USE var_lookup,only:iLookTIME,iLookFORCE      ! named variables to define structure elements
   USE summaFileManager,only:FORCING_PATH        ! path of the forcing data file
+  USE globalData,only:ixHRUfile_min,ixHRUfile_max
 
   implicit none
   integer(c_int),intent(in)               :: iFile
@@ -109,7 +126,8 @@ subroutine read_forcingFile(iFile, startGRU, numGRU, err, message_r) &
   integer(c_int),intent(in)               :: numGRU
   integer(c_int),intent(out)              :: err
   type(c_ptr), intent(out)                :: message_r
-  ! local varibles            
+  ! local varibles
+  integer(i4b)                            :: nHRUlocal            
   integer(i4b)                            :: iHRU_Global
   integer(i4b)                            :: varId
   integer(i4b)                            :: ncid
@@ -125,7 +143,8 @@ subroutine read_forcingFile(iFile, startGRU, numGRU, err, message_r) &
   character(len = nf90_max_name)          :: varName          ! dimenison name
   logical(lgt),dimension(size(forc_meta)) :: checkForce       ! flags to check forcing data variables exist
   character(len=256)                      :: message          ! error message 
-  
+
+  nHRUlocal = sum(gru_struc(:)%hruCount)
 
   ! Start Procedure here
   err=0; message="read_forcingFile/"
@@ -185,7 +204,7 @@ subroutine read_forcingFile(iFile, startGRU, numGRU, err, message_r) &
     iVar = forcFileInfo(iFile)%var_ix(iNC)
     checkForce(iVar) = .true.
     if (.not.allocated(forcingDataStruct(iFile)%var(iVar)%dataFromFile))then
-      allocate(forcingDataStruct(iFile)%var(iVar)%dataFromFile(numGRU,nTimeSteps))
+      allocate(forcingDataStruct(iFile)%var(iVar)%dataFromFile(nHRUlocal,nTimeSteps))
     endif
 
     ! Get Forcing Data
@@ -198,7 +217,9 @@ subroutine read_forcingFile(iFile, startGRU, numGRU, err, message_r) &
     numHRU = sum(gru_struc(:)%hruCount)
     
 
-    err=nf90_get_var(ncid,forcFileInfo(iFile)%data_id(ivar),forcingDataStruct(iFile)%var(iVar)%dataFromFile, start=(/startGRU,1/),count=(/numHRU, nTimeSteps/))
+    err=nf90_get_var(ncid,forcFileInfo(iFile)%data_id(ivar), &
+                     forcingDataStruct(iFile)%var(iVar)%dataFromFile, &
+                     start=(/ixHRUfile_min,1/),count=(/nHRUlocal, nTimeSteps/))
     if(err/=nf90_noerr)then; message=trim(message)//'problem reading forcing data: '//trim(varName)//'/'//trim(nf90_strerror(err)); call f_c_string_ptr(trim(message),message_r); return; endif
   end do
 
diff --git a/build/source/file_access_actor/output_structure.f90 b/build/source/file_access_actor/output_structure.f90
index d3321f63dcf2be51a99a6d549b137410d3f7b8bd..8544685c7f7c9368b807ec95e64d3e37c827de0c 100644
--- a/build/source/file_access_actor/output_structure.f90
+++ b/build/source/file_access_actor/output_structure.f90
@@ -296,7 +296,7 @@ subroutine initOutputStructure(maxSteps, num_gru, err)
               call alloc_outputStruc(indx_meta,summa_struct(1)%indxStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,str_name='indx',message=message);    ! model variables
               ! Statistics
-              call alloc_outputStruc(statIndx_meta(:)%var_info,summa_struct(1)%indxStat%gru(iGRU)%hru(1), &
+              call alloc_outputStruc(statIndx_meta(:)%var_info,summa_struct(1)%indxStat%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! index vars
             case('prog')
               ! Structure
@@ -347,9 +347,6 @@ subroutine initOutputStructure(maxSteps, num_gru, err)
       end do ! timeSteps
     end do ! Looping through GRUs
   end do
-
-
-
 end subroutine initOutputStructure
 
 subroutine deallocateOutputStructure(err) bind(C, name="deallocateOutputStructure")
@@ -650,7 +647,10 @@ subroutine alloc_outputStruc(metaStruct,dataStruct,nSteps,nSnow,nSoil,str_name,e
       do iVar=1, nVars
         ! Check if this variable is desired within any timeframe
         if(is_var_desired(metaStruct,iVar) .or. allocAllFlag)then
-          allocate(dataStruct%var(iVar)%tim(nSteps))
+          if (allocated(dataStruct%var(iVar)%tim)) then
+            print*, "Already Allocated"; return;
+          end if
+          allocate(dataStruct%var(iVar)%tim(nSteps), stat=err)
           call allocateDat_rkind_nSteps(metaStruct,dataStruct,nSnow,nSoil,nSteps,iVar,err,cmessage)
         end if
       end do
diff --git a/build/source/global/actor_data_types.f90 b/build/source/global/actor_data_types.f90
index ffe1132a9fac9feb9b5065aa0c542ef5868b3f6b..f73e365a55b49d63977d5de86fefd8482acabcee 100644
--- a/build/source/global/actor_data_types.f90
+++ b/build/source/global/actor_data_types.f90
@@ -5,21 +5,6 @@ module actor_data_types
   implicit none
   private
 
-  
-  type,public :: forcingFileData
-    real(rkind), dimension (:,:), allocatable   :: dataFromFile
-  end type forcingFileData
-
-  type,public :: var_forc
-    type(forcingFileData), allocatable   :: var(:)       ! var(:)%dataFromFile(:,:)
-    character(len=256)                   :: refTimeString
-    real(rkind)                          :: convTime2Days
-    integer(i4b)                         :: nVars
-    integer(i4b),allocatable             :: var_ix(:)
-    real(rkind)                          :: tmZoneOffsetFracDay
-    real(rkind)                          :: refJulDay_data 
-    integer(i4b)                         :: nTimeSteps    ! Number of Timesteps in the file
-  end type var_forc
 
 
   ! ** double precision type of for time series
@@ -172,5 +157,13 @@ module actor_data_types
     integer(c_int)                             :: yearLength                 ! number of days in the current year
     ! Misc Variables
     integer(c_int)                             :: computeVegFlux             ! flag to indicate if we are computing fluxes over vegetation
+    real(c_double)                             :: dt_init
+    real(c_double)                             :: upArea
   end type hru_type
+
+  type, public :: gru_type
+    type(hru_type),allocatable :: hru(:)
+    type(var_dlength),pointer  :: bvarStat
+    type(var_dlength),pointer  :: bvarStruct
+  end type gru_type
 end module
diff --git a/build/source/global/auxiliary.cpp b/build/source/global/auxiliary.cpp
index 2afef1e0e432ef3ae390cd4bbc2a8aa83cd94d75..c12978545a99d4ad4239a48e64e347f4a6fe6716 100644
--- a/build/source/global/auxiliary.cpp
+++ b/build/source/global/auxiliary.cpp
@@ -523,15 +523,17 @@ std::vector<std::vector<std::vector<double>>> get_lookup_struct(void *handle) {
 #endif
 
 void get_scalar_data(void* handle, double fracJulDay, 
-    double tmZoneOffsetFracDay, int year_length, int computeVegFlux) {
+                     double tmZoneOffsetFracDay, int year_length, 
+                     int computeVegFlux, double dt_init, double upArea) {
   get_scalar_data_fortran(handle, &fracJulDay, &tmZoneOffsetFracDay, 
-      &year_length, &computeVegFlux);
+      &year_length, &computeVegFlux, &dt_init, &upArea);
 }
 
 void set_scalar_data(void* handle, double fracJulDay, 
-    double tmZoneOffsetFracDay, int year_length, int computeVegFlux) {
+                     double tmZoneOffsetFracDay, int year_length, 
+                     int computeVegFlux, double dt_init, double upArea) {
   set_scalar_data_fortran(handle, &fracJulDay, &tmZoneOffsetFracDay, 
-      &year_length, &computeVegFlux);
+      &year_length, &computeVegFlux, &dt_init, &upArea);
 }
 
 
diff --git a/build/source/global/cppwrap_datatypes.f90 b/build/source/global/cppwrap_datatypes.f90
index 0cc81f77579cb7e2c1d946951186824fd9475db8..bb65c59ffc0d737af37e2e64ca9a2bf31d0cbb4f 100644
--- a/build/source/global/cppwrap_datatypes.f90
+++ b/build/source/global/cppwrap_datatypes.f90
@@ -2480,12 +2480,14 @@ end subroutine
 ! ****************************** flag_vec ****************************
 
 subroutine get_scalar_data_fortran(handle, fracJulDay, tmZoneOffsetFracDay, &
-    year_length, computeVegFlux) bind(C, name='get_scalar_data_fortran')
+    year_length, computeVegFlux, dt_init, upArea) bind(C, name='get_scalar_data_fortran')
   type(c_ptr), intent(in), value :: handle
   real(c_double), intent(out) :: fracJulDay
   real(c_double), intent(out) :: tmZoneOffsetFracDay
   integer(c_int), intent(out) :: year_length
   integer(c_int), intent(out) :: computeVegFlux
+  real(c_double), intent(out) :: dt_init
+  real(c_double), intent(out) :: upArea
   type(hru_type), pointer :: hru_data
 
   call c_f_pointer(handle, hru_data)
@@ -2494,16 +2496,20 @@ subroutine get_scalar_data_fortran(handle, fracJulDay, tmZoneOffsetFracDay, &
   tmZoneOffsetFracDay = hru_data%tmZoneOffsetFracDay
   year_length = hru_data%yearLength
   computeVegFlux = hru_data%computeVegFlux
+  dt_init = hru_data%dt_init
+  upArea = hru_data%upArea
 
 end subroutine get_scalar_data_fortran
 
 subroutine set_scalar_data_fortran(handle, fracJulDay, tmZoneOffsetFracDay, &
-    year_length, computeVegFlux) bind(C, name='set_scalar_data_fortran')
+    year_length, computeVegFlux, dt_init, upArea) bind(C, name='set_scalar_data_fortran')
   type(c_ptr), intent(in), value :: handle
   real(c_double), intent(in) :: fracJulDay
   real(c_double), intent(in) :: tmZoneOffsetFracDay
   integer(c_int), intent(in) :: year_length
   integer(c_int), intent(in) :: computeVegFlux
+  real(c_double), intent(in) :: dt_init
+  real(c_double), intent(in) :: upArea
   type(hru_type), pointer :: hru_data
 
   call c_f_pointer(handle, hru_data)
@@ -2512,6 +2518,8 @@ subroutine set_scalar_data_fortran(handle, fracJulDay, tmZoneOffsetFracDay, &
   hru_data%tmZoneOffsetFracDay = tmZoneOffsetFracDay
   hru_data%yearLength = year_length
   hru_data%computeVegFlux = computeVegFlux
+  hru_data%dt_init = dt_init
+  hru_data%upArea = upArea
 
 end subroutine set_scalar_data_fortran
 
@@ -2594,6 +2602,65 @@ subroutine delete_handle_hru_type(handle) bind(C, name="delete_handle_hru_type")
 
 end subroutine
 
+function new_handle_gru_type(num_hru) result(handle) bind(C, name="new_handle_gru_type")
+  type(c_ptr)                :: handle
+  integer(c_int), intent(in) :: num_hru
+  type(gru_type), pointer    :: p
+  integer(c_int)             :: i
+
+  allocate(p)
+  allocate(p%hru(num_hru))
+  allocate(p%bvarStat)
+  allocate(p%bvarStruct)
+
+  do i=1,num_hru
+#ifdef V4_ACTIVE
+  allocate(p%hru(i)%lookupStruct)
+#endif
+  allocate(p%hru(i)%forcStat)
+  allocate(p%hru(i)%progStat)
+  allocate(p%hru(i)%diagStat)
+  allocate(p%hru(i)%fluxStat)
+  allocate(p%hru(i)%indxStat)
+  allocate(p%hru(i)%bvarStat)
+  allocate(p%hru(i)%timeStruct)
+  allocate(p%hru(i)%forcStruct)
+  allocate(p%hru(i)%attrStruct)
+  allocate(p%hru(i)%typeStruct)
+  allocate(p%hru(i)%idStruct)
+  allocate(p%hru(i)%indxStruct)
+  allocate(p%hru(i)%mparStruct)
+  allocate(p%hru(i)%progStruct)
+  allocate(p%hru(i)%diagStruct)
+  allocate(p%hru(i)%fluxStruct)
+  allocate(p%hru(i)%bparStruct)
+  allocate(p%hru(i)%bvarStruct)
+  allocate(p%hru(i)%dparStruct)
+  allocate(p%hru(i)%startTime_hru)
+  allocate(p%hru(i)%finishTime_hru)
+  allocate(p%hru(i)%refTime_hru)
+  allocate(p%hru(i)%oldTime_hru)
+  allocate(p%hru(i)%statCounter)
+  allocate(p%hru(i)%outputTimeStep)
+  allocate(p%hru(i)%resetStats)
+  allocate(p%hru(i)%finalizeStats)
+  end do
+
+
+  handle = c_loc(p)
+end function new_handle_gru_type
+
+subroutine delete_handle_gru_type(handle) bind(C, name="delete_handle_gru_type")
+  type(c_ptr), intent(in), value :: handle
+  type(gru_type), pointer :: p
+
+  call c_f_pointer(handle, p)
+  deallocate(p%hru)
+  deallocate(p%bvarStat)
+  deallocate(p%bvarStruct)
+  deallocate(p)
+end subroutine delete_handle_gru_type
+
 end module cppwrap_datatypes
 
 
diff --git a/build/source/gru_actor/gru_actor.cpp b/build/source/gru_actor/gru_actor.cpp
index 81720a66bb273ac52be504d32c4606c5cc1c907e..efd395313fb6741bee6f371fc29efcb085b211ba 100644
--- a/build/source/gru_actor/gru_actor.cpp
+++ b/build/source/gru_actor/gru_actor.cpp
@@ -1,151 +1,115 @@
-#include "caf/all.hpp"
 #include "gru_actor.hpp"
-#include "global.hpp"
-#include "hru_actor.hpp"
-#include "message_atoms.hpp"
-#include <vector>
-#include "gru_actor_subroutine_wrappers.hpp"
 
-namespace caf {
 
-behavior gru_actor(stateful_actor<gru_state>* self, 
-    int ref_gru, int indx_gru, std::string config_path, caf::actor file_access_actor,
-    caf::actor parent) {
-
-    aout(self) << "GRU Actor Has Started\n";
-    self->state.parent = parent;
-    self->state.ref_gru = ref_gru;
-    self->state.indx_gru = indx_gru;
-    self->state.config_path = config_path;
-    self->state.file_access_actor = file_access_actor;
-    self->state.parent = parent;
-
-    self->state.num_hrus = getNumHRU(&self->state.indx_gru);
-
-    self->send(self, init_hru_v);
-
-    return {
-
-        [=](init_hru) {
-            for (int i = 0; i < self->state.num_hrus; i++) {
-                // auto hru = self->spawn(hru_actor,
-                //         self->state.ref_gru, self->state.indx_gru, 
-                //         self->state.config_path, self->state.file_access_actor,
-                //         self->state.output_struc_size,
-                //         self);
-                // self->state.hru_list.push_back(hru);
-            }
-        },
-
-        [=](done_hru, int indx_gru, double total_duration, double init_duration, 
-            double forcing_duration, double run_physics_duration, double write_output_duration) {
-            aout(self) << "GRU Received HRU is Done\n";
-
-            self->state.hrus_complete++;
-            if (self->state.hrus_complete >= self->state.num_hrus) {
-                aout(self) << "All HRUs have finished";
-
-                self->send(self->state.parent,
-                    done_hru_v,
-                    indx_gru, 
-                    total_duration,
-                    init_duration, 
-                    forcing_duration,
-                    run_physics_duration,
-                    write_output_duration);
-            }
-        
-        },
-
-        // What does a GRU need to assemble its data structure?
-        [=](init_gru) {
-            // Get the variable data length, we also need the type information
-            aout(self) << "init GRU \n";
-            int integer_missing_value = -9999;
-
-            // Get the sizes we need for the lists from Fortran
-            getVarSizes(&self->state.num_var_types,
-                        &self->state.num_bpar_vars,
-                        &self->state.num_bvar_vars);
-            
-            aout(self) << "GRU: GOT VAR SIZES\n";
-            aout(self) << "NUM VAR Type = " << self->state.num_var_types << "\n";
-            aout(self) << "NUM BPAR = " << self->state.num_bpar_vars << "\n";
-            aout(self) << "NUM BVAR = " << self->state.num_bvar_vars << "\n";
-
-            for (int i = 0; i < self->state.num_var_types; i++) {
-                self->state.i_look_var_type_list.push_back(integer_missing_value);
-            }
-            for (int i = 0; i < self->state.num_bpar_vars; i++) {
-                self->state.bpar_struct_var_type_list.push_back(integer_missing_value);
-            }
-            for (int i = 0; i < self->state.num_bvar_vars; i++) {
-                self->state.bvar_struct_var_type_list.push_back(integer_missing_value);
-            }
-
-            initVarType(&self->state.var_type_lookup);
-                // aout(self) << "************C++************\n";
-                // aout(self) << self->state.var_type_lookup.scalarv << "\n";
-                // aout(self) << self->state.var_type_lookup.wLength << "\n";
-                // aout(self) << self->state.var_type_lookup.midSnow << "\n";
-                // aout(self) << self->state.var_type_lookup.midSoil << "\n";
-                // aout(self) << self->state.var_type_lookup.midToto << "\n";
-                // aout(self) << self->state.var_type_lookup.ifcSnow << "\n";
-                // aout(self) << self->state.var_type_lookup.ifcSoil << "\n";
-                // aout(self) << self->state.var_type_lookup.ifcToto << "\n";
-                // aout(self) << self->state.var_type_lookup.parSoil << "\n";
-                // aout(self) << self->state.var_type_lookup.routing << "\n";
-                // aout(self) << self->state.var_type_lookup.outstat << "\n";
-                // aout(self) << self->state.var_type_lookup.unknown << "\n";
-                // aout(self) << "************C++************\n";
-
-            // Fill The lists with the values from the fortran lists
-            int err = 0;
-            fillVarTypeLists(&self->state.num_bpar_vars,
-                             &self->state.num_bvar_vars,
-                             &self->state.bpar_struct_var_type_list[0],
-                             &self->state.bvar_struct_var_type_list[0],
-                             &err);
-            // aout(self) << "Printing BPAR\n";
-            // for(int i = 0; i < self->state.num_bpar_vars; i++) {
-            //     aout(self) << i << ": " << self->state.bpar_struct_var_type_list[i] << "\n";
-            // }
-            // aout(self) << "Printing BVAR\n";
-            // for(int i = 0; i < self->state.num_bvar_vars; i++) {
-            //     aout(self) << i << ": " << self->state.bvar_struct_var_type_list[i] << "\n";
-            // }
-
-            // Now we can allocate space for the structures
-            for(int i = 0; i < self->state.num_bpar_vars; i++) {
-                if (self->state.bpar_struct_var_type_list[i] == self->state.var_type_lookup.scalarv) {
-                    self->state.bpar_struct.push_back(0.0);
-                } else {
-                    aout(self) << "ERROR: GRU - bpar_struct contains type that is not a scalar\n";
-                }
-            }
-            for(int i = 0; i < self->state.num_bvar_vars; i++) {
-                if (self->state.bvar_struct_var_type_list[i] == self->state.var_type_lookup.scalarv) {
-                    std::vector<double> temp;
-                    self->state.bvar_struct.push_back(temp);
-                    self->state.bvar_struct[i].push_back(0.0);
-
-                } else if(self->state.bvar_struct_var_type_list[i] == self->state.var_type_lookup.routing) {
-                    std::vector<double> temp;
-                    self->state.bvar_struct.push_back(temp);
-                    for (int x = 0; x < self->state.nTimeDelay; x++) {
-                        self->state.bvar_struct[i].push_back(0.0);
-                    }
-                } else {
-                    aout(self) << "ERROR: GRU - bvar_struct contains type that is not a scalar or routing\n";
-
-                }
-            }
-
-            self->send(self->state.parent, done_init_gru_v);
+using namespace caf;
+
+behavior gru_actor(stateful_actor<gru_actor_state>* self, int netcdf_index, 
+                   int gru_job_index, int num_steps, 
+                   HRU_Actor_Settings hru_actor_settings,
+                   actor file_access_actor, actor parent) {
+  self->state.netcdf_index = netcdf_index;
+  self->state.gru_job_index = gru_job_index;
+  self->state.num_steps = num_steps;
+  self->state.hru_actor_settings = hru_actor_settings;
+  self->state.file_access_actor = file_access_actor;
+  self->state.parent = parent;
+
+  // Check for lateral flows
+  getNumHRU(self->state.gru_job_index, self->state.num_hrus);
+  aout(self) << "NUM HRUS: " << self->state.num_hrus << "\n";
+  self->state.hrus.resize(self->state.num_hrus);
+  self->state.gru_data = new_handle_gru_type(self->state.num_hrus);
+  int err = 0;
+  std::unique_ptr<char[]> message(new char[256]);
+  initGRU_fortran(self->state.gru_job_index, self->state.gru_data, err, 
+                  &message);
+  std::fill(message.get(), message.get() + 256, '\0'); // Clear message
+  setupGRU_fortran(self->state.gru_job_index, self->state.gru_data, err, 
+                   &message);
+  std::fill(message.get(), message.get() + 256, '\0'); // Clear message
+  readGRURestart_fortran(self->state.gru_job_index, self->state.gru_data, err, 
+                         &message);
+
+  aout(self) << "GRU Actor: HRUs Initialized\n";
+  self->send(self, update_hru_async_v);
+
+  return {
+    [=](update_hru_async) {
+      self->request(self->state.file_access_actor, caf::infinite,
+                    get_num_output_steps_v).await([=](int num_steps) {
+        self->state.num_steps_until_write = num_steps;
+        self->send(self->state.file_access_actor, access_forcing_v, 
+                   self->state.iFile, self);
+      });
+    },
+    [=](new_forcing_file, int num_forcing_steps_in_iFile, int iFile) {
+      int err;
+      std::unique_ptr<char[]> message(new char[256]);
+      self->state.iFile = iFile;
+      self->state.stepsInCurrentFFile = num_forcing_steps_in_iFile;
+      setTimeZoneOffsetGRU_fortran(self->state.iFile, self->state.gru_data, 
+                                   err, &message);
+      if (err != 0) {
+        aout(self) << "GRU_Actor: Error setting time zone offset\n";
+        self->quit();
+        return;
+      }
+      self->state.forcingStep = 1;
+      self->send(self, run_hru_v);
+    },
+
+    [=](num_steps_before_write, int num_steps) {
+      self->state.num_steps_until_write = num_steps;
+      self->state.output_structure_step_index = 1;
+    },
+    
+    [=](run_hru) {
+      int err = 0;
+      std::unique_ptr<char[]> message(new char[256]);
+      while(self->state.num_steps_until_write > 0) {
+        if (self->state.forcingStep > self->state.stepsInCurrentFFile) {
+          aout(self) << "GRU Actor: New Forcing File\n";
+          self->send(self->state.file_access_actor, access_forcing_v, 
+                     self->state.iFile+1, self);
+          break;
         }
-
-
-    };
-}
-
+        self->state.num_steps_until_write--;
+        aout(self) << "GRU Actor: timestep=" << self->state.timestep << "\n";
+        readGRUForcing_fortran(self->state.gru_job_index, 
+                               self->state.forcingStep, 
+                               self->state.timestep, self->state.iFile, 
+                               self->state.gru_data, err, &message);
+        std::fill(message.get(), message.get() + 256, '\0'); // Clear message
+        runGRU_fortran(self->state.gru_job_index, self->state.timestep, 
+                       self->state.gru_data, self->state.dt_init_factor,
+                       err, &message);
+        std::fill(message.get(), message.get() + 256, '\0'); // Clear message
+        writeGRUOutput_fortran(self->state.gru_job_index, self->state.timestep, 
+                               self->state.output_structure_step_index, 
+                               self->state.gru_data, err, &message);
+
+        self->state.timestep++;
+        self->state.forcingStep++;
+        self->state.output_structure_step_index++;
+
+        if (self->state.timestep > self->state.num_steps) {
+          aout(self) << "GRU Actor: Done\n";
+          self->send(self, done_hru_v);
+          break;
+        }
+      }
+      // Our output structure is full
+      if (self->state.num_steps_until_write <= 0) {
+        aout(self) << "GRU Actor: Writing Output\n";
+        self->send(self->state.file_access_actor, write_output_v,
+                   self->state.gru_job_index, 1, self);
+      }
+    },
+    
+    [=](done_hru) {
+      self->send(self->state.parent,done_hru_v,self->state.gru_job_index);
+      self->quit();
+      return;
+    },
+  };
 }
\ No newline at end of file
diff --git a/build/source/gru_actor/gru_actor.f90 b/build/source/gru_actor/gru_actor.f90
deleted file mode 100644
index cf52f9b8f59d6516f223c0005b5f70aa6ec1eb4c..0000000000000000000000000000000000000000
--- a/build/source/gru_actor/gru_actor.f90
+++ /dev/null
@@ -1,222 +0,0 @@
-!!! Lets try and build this for only the lateral flows case.
-!!! If lateral flows exits use the code
-module gru_actor
-USE,intrinsic :: iso_c_binding
-USE nrtype
-
-implicit none
-
-! public::run_gru
-public::getVarSizes
-public::fillvarTypeLists
-public::getNumHRU
-
-contains
-
-subroutine getVarSizes(num_var_types, &
-                       num_bpar_vars, &
-                       num_bvar_vars) bind(C,name="getVarSizes")
-    USE var_lookup,only:maxvarBpar, maxvarBvar, maxvarVarType
-    implicit none
-    
-    integer(c_int), intent(out)                   :: num_var_types
-    integer(c_int), intent(out)                   :: num_bpar_vars
-    integer(c_int), intent(out)                   :: num_bvar_vars
-
-    num_var_types = maxvarVarType
-    num_bpar_vars = maxvarBpar
-    num_bvar_vars = maxvarBvar
-
-end subroutine getVarSizes
-
-subroutine initVarType(var_type_lookup) bind(C, name="initVarType")
-    USE var_lookup,only:iLookVarType
-    USE var_lookup,only:iLook_VarType
-
-    implicit none
-    type(iLook_VarType) :: var_type_lookup
-
-    ! Get the indexes to match C++ offset starting at 0
-    var_type_lookup%scalarv = iLookVarType%scalarv - 1
-    var_type_lookup%wLength = iLookVarType%wLength - 1
-    var_type_lookup%midSnow = iLookVarType%midSnow - 1
-    var_type_lookup%midSoil = iLookVarType%midSoil - 1
-    var_type_lookup%midToto = iLookVarType%midToto - 1
-    var_type_lookup%ifcSnow = iLookVarType%ifcSnow - 1
-    var_type_lookup%ifcSoil = iLookVarType%ifcSoil - 1
-    var_type_lookup%ifcToto = iLookVarType%ifcToto - 1
-    var_type_lookup%parSoil = iLookVarType%parSoil - 1
-    var_type_lookup%routing = iLookVarType%routing - 1
-    var_type_lookup%outstat = iLookVarType%outstat - 1
-    var_type_lookup%unknown = iLookVarType%unknown - 1
-
-    ! check the values
-    ! print*, "************FORTRAN************"
-    ! print*, "iLookVarType%scalarv", iLookVarType%scalarv
-    ! print*, "iLookVarType%wLength", iLookVarType%wLength
-    ! print*, "iLookVarType%midSnow", iLookVarType%midSnow
-    ! print*, "iLookVarType%midSoil", iLookVarType%midSoil
-    ! print*, "iLookVarType%midToto", iLookVarType%midToto
-    ! print*, "iLookVarType%ifcSnow", iLookVarType%ifcSnow
-    ! print*, "iLookVarType%ifcSoil", iLookVarType%ifcSoil
-    ! print*, "iLookVarType%ifcToto", iLookVarType%ifcToto
-    ! print*, "iLookVarType%parSoil", iLookVarType%parSoil
-    ! print*, "iLookVarType%routing", iLookVarType%routing
-    ! print*, "iLookVarType%outstat", iLookVarType%outstat
-    ! print*, "iLookVarType%unknown", iLookVarType%unknown
-    ! print*, "************FORTRAN************"
-
-
-end subroutine
-subroutine fillVarTypeLists(num_bpar_vars, &
-                            num_bvar_vars, &
-                            bpar_struct_var_type_list, &
-                            bvar_struct_var_type_list, &
-                            err) bind(C, name="fillVarTypeLists")
-    
-    USE globalData,only:bpar_meta,bvar_meta
-    USE var_lookup,only:iLookBVAR,iLookBPAR,iLookVarType
-    implicit none
-    integer(c_int), intent(in)                              :: num_bpar_vars
-    integer(c_int), intent(in)                              :: num_bvar_vars
-    integer(c_int), intent(out), dimension(num_bpar_vars)   :: bpar_struct_var_type_list
-    integer(c_int), intent(out), dimension(num_bvar_vars)   :: bvar_struct_var_type_list
-    integer(c_int), intent(out)                             :: err  
-    integer(i4b)                                            :: iVar
-
-
-    ! Get the types of the variables for bparStruct
-    ! Index in bpar_struct_var_type_list will match bpar_Struct
-    do iVar=1, num_bpar_vars
-        select case(bpar_meta(iVar)%vartype)
-            case(iLookVarType%scalarv); bpar_struct_var_type_list(iVar) = iLookVarType%scalarv - 1
-            case(iLookVarType%wLength); bpar_struct_var_type_list(iVar) = iLookVarType%wLength - 1
-            case(iLookVarType%midSnow); bpar_struct_var_type_list(iVar) = iLookVarType%midSnow - 1
-            case(iLookVarType%midSoil); bpar_struct_var_type_list(iVar) = iLookVarType%midSoil - 1
-            case(iLookVarType%midToto); bpar_struct_var_type_list(iVar) = iLookVarType%midToto - 1
-            case(iLookVarType%ifcSnow); bpar_struct_var_type_list(iVar) = iLookVarType%ifcSnow - 1
-            case(iLookVarType%ifcSoil); bpar_struct_var_type_list(iVar) = iLookVarType%ifcSoil - 1
-            case(iLookVarType%ifcToto); bpar_struct_var_type_list(iVar) = iLookVarType%ifcToto - 1
-            case(iLookVarType%parSoil); bpar_struct_var_type_list(iVar) = iLookVarType%parSoil - 1
-            case(iLookVarType%routing); bpar_struct_var_type_list(iVar) = iLookVarType%routing - 1
-            case(iLookVarType%outstat); bpar_struct_var_type_list(iVar) = iLookVarType%outstat - 1
-            case(iLookVarType%unknown); bpar_struct_var_type_list(iVar) = iLookVarType%unknown - 1
-            case default
-                err = 40;
-                return
-        end select
-    end do
-
-    do iVar=1, num_bvar_vars
-        select case(bvar_meta(iVar)%vartype)
-            case(iLookVarType%scalarv); bvar_struct_var_type_list(iVar) = iLookVarType%scalarv - 1
-            case(iLookVarType%wLength); bvar_struct_var_type_list(iVar) = iLookVarType%wLength - 1
-            case(iLookVarType%midSnow); bvar_struct_var_type_list(iVar) = iLookVarType%midSnow - 1
-            case(iLookVarType%midSoil); bvar_struct_var_type_list(iVar) = iLookVarType%midSoil - 1
-            case(iLookVarType%midToto); bvar_struct_var_type_list(iVar) = iLookVarType%midToto - 1
-            case(iLookVarType%ifcSnow); bvar_struct_var_type_list(iVar) = iLookVarType%ifcSnow - 1
-            case(iLookVarType%ifcSoil); bvar_struct_var_type_list(iVar) = iLookVarType%ifcSoil - 1
-            case(iLookVarType%ifcToto); bvar_struct_var_type_list(iVar) = iLookVarType%ifcToto - 1
-            case(iLookVarType%parSoil); bvar_struct_var_type_list(iVar) = iLookVarType%parSoil - 1
-            case(iLookVarType%routing); bvar_struct_var_type_list(iVar) = iLookVarType%routing - 1
-            case(iLookVarType%outstat); bvar_struct_var_type_list(iVar) = iLookVarType%outstat - 1
-            case(iLookVarType%unknown); bvar_struct_var_type_list(iVar) = iLookVarType%unknown - 1
-            case default
-                err = 40;
-                return
-        end select
-    end do
-end subroutine fillVarTypeLists
-
-
-
-! subroutine run_gru_for_timestep()
-!     implicit none
-
-!     ! initialize runoff variables
-!     bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1)    = 0._dp  ! surface runoff (m s-1)
-!     bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)     = 0._dp 
-!     bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)    = 0._dp  ! outflow from all "outlet" HRUs (those with no downstream HRU)
-!     bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1)      = 0._dp 
-  
-!     ! initialize baseflow variables
-!     bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)  = 0._dp ! recharge to the aquifer (m s-1)
-!     bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  = 0._dp ! baseflow from the aquifer (m s-1)
-!     bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = 0._dp ! transpiration loss from the aquifer (m s-1)
-
-!     ! ----- calculate weighted basin (GRU) fluxes --------------------------------------------------------------------------------------
- 
-!     ! increment basin surface runoff (m s-1)
-!     bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) = bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + fluxStruct%var(iLookFLUX%scalarSurfaceRunoff)%dat(1) * fracHRU
-    
-!     !increment basin soil drainage (m s-1)
-!     bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)   = bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)  + fluxStruct%var(iLookFLUX%scalarSoilDrainage)%dat(1)  * fracHRU
-    
-!     ! increment aquifer variables -- ONLY if aquifer baseflow is computed individually for each HRU and aquifer is run
-!     ! NOTE: groundwater computed later for singleBasin
-!     if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == localColumn .and. model_decisions(iLookDECISIONS%groundwatr)%iDecision == bigBucket) then
-   
-!       bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)  = bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)   + fluxStruct%var(iLookFLUX%scalarSoilDrainage)%dat(1)     * fracHRU
-!       bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1)  + fluxStruct%var(iLookFLUX%scalarAquiferTranspire)%dat(1) * fracHRU
-!       bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  =  bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  &
-!               +  fluxStruct%var(iLookFLUX%scalarAquiferBaseflow)%dat(1) * fracHRU
-!      end if
-   
-!      ! perform the routing
-!      associate(totalArea => bvarStruct%var(iLookBVAR%basin__totalArea)%dat(1) )
-   
-!     ! compute water balance for the basin aquifer
-!     if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == singleBasin)then
-!      message=trim(message)//'multi_driver/bigBucket groundwater code not transferred from old code base yet'
-!      err=20; return
-!     end if
-   
-!     ! calculate total runoff depending on whether aquifer is connected
-!     if(model_decisions(iLookDECISIONS%groundwatr)%iDecision == bigBucket) then
-!      ! aquifer
-!      bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)
-!     else
-!      ! no aquifer
-!      bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)
-!     endif
-   
-!     call qOverland(&
-!                    ! input
-!                    model_decisions(iLookDECISIONS%subRouting)%iDecision,            &  ! intent(in): index for routing method
-!                    bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1),            &  ! intent(in): total runoff to the channel from all active components (m s-1)
-!                    bvarStruct%var(iLookBVAR%routingFractionFuture)%dat,             &  ! intent(in): fraction of runoff in future time steps (m s-1)
-!                    bvarStruct%var(iLookBVAR%routingRunoffFuture)%dat,               &  ! intent(in): runoff in future time steps (m s-1)
-!                    ! output
-!                    bvarStruct%var(iLookBVAR%averageInstantRunoff)%dat(1),           &  ! intent(out): instantaneous runoff (m s-1)
-!                    bvarStruct%var(iLookBVAR%averageRoutedRunoff)%dat(1),            &  ! intent(out): routed runoff (m s-1)
-!                    err,message)                                                                  ! intent(out): error control
-!     if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
-!     end associate
-    
-
-
-! end subroutine
-
-integer(c_int) function getNumHRU(indx_gru) bind(C, name="getNumHRU")
-    USE globalData,only:gru_struc
-    implicit none
-    integer(c_int), intent(in)        :: indx_gru
-    integer(c_int)                    :: num_hru
-
-    num_hru = gru_struc(indx_gru)%hruCount
-
-    getNumHRU = num_hru
-
-end function
-
-
-end module gru_actor
-
-! the way the lateral flow interface should work seems to be at the level of the hru. 
-! The HRU would have an interface where it has a donwslope compontent. It would then ask
-! for this downslope component. This all does have to be setup before hand so the gru needs to 
-! set these up because the hrus all need to comptue 
-! a grus timestep depends on all the hrus
-
-! So the Gru needs to be controlling the hrus. It needs to know which ones have lateral flows and which ones do not. 
-! This is so the other hrus know they do not need to contact another actor. Otherwise they will have to.
\ No newline at end of file
diff --git a/build/source/gru_actor/gru_interface.f90 b/build/source/gru_actor/gru_interface.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ca0a9cddf11b5c05ad8dd94443c824e671cf5a56
--- /dev/null
+++ b/build/source/gru_actor/gru_interface.f90
@@ -0,0 +1,406 @@
+module gru_interface
+USE,intrinsic :: iso_c_binding
+USE nrtype
+
+
+implicit none
+public :: getNumHRU
+public :: initGRU_fortran
+public :: setupGRU_fortran
+public :: readGRURestart_fortran
+public :: setTimeZoneOffsetGRU_fortran
+public :: readGRUForcing_fortran
+public :: runGRU_fortran
+public :: writeGRUOutput_fortran
+
+
+contains
+
+subroutine getNumHRU(indx_gru, num_hru) bind(C, name="getNumHRU")
+  USE globalData,only:gru_struc
+  implicit none
+  integer(c_int), intent(in)  :: indx_gru
+  integer(c_int), intent(out) :: num_hru
+
+  num_hru = gru_struc(indx_gru)%hruCount
+end subroutine getNumHRU
+
+subroutine initGRU_fortran(indx_gru, handle_gru_data, err, message_r) &
+    bind(C, name="initGRU_fortran")
+  USE actor_data_types,only:gru_type             
+  USE data_types,only:var_dlength
+  USE globalData,only:statBvar_meta                           ! child metadata for stats
+  USE globalData,only:bvar_meta                     ! metadata structures
+  USE allocspace_module,only:allocLocal
+  USE INIT_HRU_ACTOR,only:initHRU
+
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  implicit none
+  ! Dummy variables
+  integer(c_int), intent(in)          :: indx_gru
+  type(c_ptr),    intent(in),value    :: handle_gru_data
+  integer(c_int), intent(out)         :: err
+  type(c_ptr),   intent(out)          :: message_r
+
+  ! local variables
+  type(gru_type),pointer              :: gru_data
+  integer(i4b)                        :: iHRU
+  character(len=256)                  :: message = ""
+  character(len=256)                  :: cmessage
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_gru_data, gru_data)
+
+  call allocLocal(bvar_meta,gru_data%bvarStruct,nSnow=0,nSoil=0,err=err,message=cmessage);
+  if(err /= 0) then; message=trim(message)//cmessage; call f_c_string_ptr(trim(message), message_r);return;end if 
+  call allocLocal(statBvar_meta(:)%var_info,gru_data%bvarStat,nSnow=0,nSoil=0,err=err,message=cmessage);
+  if(err /= 0) then; message=trim(message)//cmessage; call f_c_string_ptr(trim(message), message_r);return;end if 
+
+
+  do iHRU = 1, size(gru_data%hru)
+    call initHRU(indx_gru, iHRU, gru_data%hru(iHRU), err, message)
+    if(err /= 0) then; call f_c_string_ptr(trim(message), message_r);return; end if
+  end do
+end subroutine initGru_fortran
+
+subroutine setupGRU_fortran(indx_gru, handle_gru_data, err, message_r) & 
+    bind(C, name="setupGRU_fortran")
+  USE summa_init_struc,only:init_struc
+  USE actor_data_types,only:gru_type
+  USE INIT_HRU_ACTOR,only:setupHRU
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  type(c_ptr),    intent(in),value :: handle_gru_data
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  integer(i4b)                     :: iHRU
+  integer(i4b)                     :: iVar
+  type(gru_type),pointer           :: gru_data
+  character(len=256)               :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_gru_data, gru_data)
+
+  do iHRU = 1, size(gru_data%hru)
+    call setupHRU(indx_gru, iHRU, gru_data%hru(iHRU), err, message)
+    if(err /= 0) then; call f_c_string_ptr(trim(message), message_r);return; end if
+  end do
+
+  do ivar=1, size(init_struc%bvarStruct%gru(indx_gru)%var(:))
+    gru_data%bvarStruct%var(ivar)%dat(:) = init_struc%bvarStruct%gru(indx_gru)%var(ivar)%dat(:)
+  enddo
+end subroutine setupGRU_fortran
+
+subroutine readGRURestart_fortran(indx_gru, handle_gru_data, err, message_r) &
+    bind(C, name="readGRURestart_fortran")
+  USE actor_data_types,only:gru_type
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  USE INIT_HRU_ACTOR,only:readHRURestart
+
+  USE var_lookup,only:iLookDECISIONS                          ! look-up values for model decisions
+  USE var_lookup,only:iLookBVAR                               ! look-up values for basin-average model variables
+  USE globalData,only:model_decisions                         ! model decision structure
+  USE mDecisions_module,only:localColumn, & ! separate groundwater representation in each local soil column
+                             singleBasin    ! single groundwater store over the entire basin
+  implicit none
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  type(c_ptr),    intent(in),value :: handle_gru_data
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  integer(i4b)                     :: iHRU
+  type(gru_type),pointer           :: gru_data
+  character(len=256)               :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_gru_data, gru_data)
+
+  do iHRU = 1, size(gru_data%hru)
+    call readHRURestart(indx_gru, iHRU, gru_data%hru(iHRU), err, message)
+    if(err /= 0) then; call f_c_string_ptr(trim(message), message_r);return; end if
+  end do
+
+  ! Set the basin variables that pertain to the GRU
+  select case(model_decisions(iLookDECISIONS%spatial_gw)%iDecision)
+    case(localColumn) 
+      gru_data%bvarStruct%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 0._dp
+    case(singleBasin)
+      gru_data%bvarStruct%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 1._dp
+    case default
+      message=trim(message)//'unable to identify decision for regional representation of groundwater'
+      call f_c_string_ptr(trim(message), message_r)
+      err = 1
+      return
+  end select
+
+end subroutine readGRURestart_fortran
+
+subroutine setTimeZoneOffsetGRU_fortran(iFile, handle_gru_data, err, message_r) & 
+    bind(C, name="setTimeZoneOffsetGRU_fortran")
+  USE actor_data_types,only:gru_type
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  USE hru_read,only:setTimeZoneOffset
+  implicit none
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: iFile
+  type(c_ptr),    intent(in),value :: handle_gru_data
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  integer(i4b)                     :: iHRU
+  type(gru_type),pointer           :: gru_data
+  character(len=256)               :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_gru_data, gru_data)
+
+  do iHRU = 1, size(gru_data%hru)
+    call setTimeZoneOffset(iFile, gru_data%hru(iHRU), err, message)
+    if(err /= 0) then; call f_c_string_ptr(trim(message), message_r);return; end if
+  end do
+
+end subroutine setTimeZoneOffsetGRU_fortran
+
+subroutine readGRUForcing_fortran(indx_gru, iStep, iRead, iFile, & 
+    handle_gru_data, err, message_r) bind(C, name="readGRUForcing_fortran")
+  USE actor_data_types,only:gru_type
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  USE hru_read,only:readHRUForcing
+  implicit none
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  integer(c_int), intent(in)       :: iStep
+  integer(c_int), intent(inout)    :: iRead
+  integer(c_int), intent(in)       :: iFile
+  type(c_ptr),    intent(in),value :: handle_gru_data
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  integer(i4b)                     :: iHRU
+  type(gru_type),pointer           :: gru_data
+  character(len=256)               :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_gru_data, gru_data)
+
+  do iHRU = 1, size(gru_data%hru)
+    call readHRUForcing(indx_gru, iHRU, iStep, iRead, iFile, gru_data%hru(iHRU), err, message)
+    if(err /= 0) then; call f_c_string_ptr(trim(message), message_r);return; end if
+  end do
+
+end subroutine readGRUForcing_fortran
+
+subroutine runGRU_fortran(indx_gru, modelTimeStep, handle_gru_data, &
+    dt_init_factor, err, message_r) bind(C, name="runGRU_fortran")
+  USE actor_data_types,only:gru_type
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  USE summa_modelRun,only:runPhysics
+  
+  USE globalData,only:model_decisions          ! model decision structure
+  USE globalData,only:gru_struc
+  USE qTimeDelay_module,only:qOverland         ! module to route water through an "unresolved" river network
+  
+  USE mDecisions_module,only:&               ! look-up values for LAI decisions
+      monthlyTable,& ! LAI/SAI taken directly from a monthly table for different vegetation classes
+      specified,&    ! LAI/SAI computed from green vegetation fraction and winterSAI and summerLAI parameters   
+      localColumn, & ! separate groundwater representation in each local soil column
+      singleBasin, & ! single groundwater store over the entire basin
+      bigBucket
+
+  USE var_lookup,only:iLookBVAR              ! look-up values for basin-average model variables
+  USE var_lookup,only:iLookFLUX              ! look-up values for local column model fluxes
+  USE var_lookup,only:iLookATTR              ! look-up values for local attributes
+  USE var_lookup,only:iLookDECISIONS         ! look-up values for model decisions
+  USE var_lookup,only:iLookTYPE              ! look-up values for HRU types
+  USE var_lookup,only:iLookID                ! look-up values for HRU IDs
+  implicit none
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  integer(c_int), intent(in)       :: modelTimeStep
+  type(c_ptr),    intent(in),value :: handle_gru_data
+  integer(c_int), intent(in)       :: dt_init_factor
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  integer(i4b)                     :: iHRU, kHRU, jHRU
+  integer(i4b)                     :: iVar
+  type(gru_type),pointer           :: gru_data
+  character(len=256)               :: message = ""
+  character(len=256)               :: cmessage
+  real(rkind)                      :: fracHRU                ! fractional area of a given HRU (-)
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_gru_data, gru_data)
+
+  ! ----- basin initialization --------------------------------------------------------------------------------------------
+  ! initialize runoff variables
+  gru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1)    = 0._dp  ! surface runoff (m s-1)
+  gru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)     = 0._dp 
+  gru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)    = 0._dp  ! outflow from all "outlet" HRUs (those with no downstream HRU)
+  gru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1)      = 0._dp 
+
+  ! initialize baseflow variables
+  gru_data%bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)  = 0._dp ! recharge to the aquifer (m s-1)
+  gru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  = 0._dp ! baseflow from the aquifer (m s-1)
+  gru_data%bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = 0._dp ! transpiration loss from the aquifer (m s-1)
+
+  do iHRU = 1, size(gru_data%hru) 
+    gru_data%hru(iHRU)%fluxStruct%var(iLookFLUX%mLayerColumnInflow)%dat(:) = 0._rkind
+  end do
+
+
+  do iHRU = 1, size(gru_data%hru)
+    ! Give the HRU the up to date basin variables
+    do iVar=1, size(gru_data%bvarStruct%var(:))
+      gru_data%hru(iHRU)%bvarStruct%var(iVar)%dat(:) = gru_data%bvarStruct%var(iVar)%dat(:)
+    end do
+    
+    call runPhysics(indx_gru, iHRU, modelTimeStep, gru_data%hru(iHRU), &
+                    dt_init_factor, err, message)
+    if(err /= 0) then; call f_c_string_ptr(trim(message), message_r);return; end if
+
+    fracHRU = gru_data%hru(iHRU)%attrStruct%var(iLookATTR%HRUarea) / &
+              gru_data%hru(iHRU)%bvarStruct%var(iLookBVAR%basin__totalArea)%dat(1)
+
+    ! Compute Fluxes Across HRUs
+    ! identify lateral connectivity
+    ! (Note:  for efficiency, this could this be done as a setup task, not every timestep)
+    kHRU = 0
+    ! identify the downslope HRU
+    dsHRU: do jHRU=1,gru_struc(indx_gru)%hruCount
+      if(gru_data%hru(iHRU)%typeStruct%var(iLookTYPE%downHRUindex) == gru_data%hru(jHRU)%idStruct%var(iLookID%hruId))then
+        if(kHRU==0)then  ! check there is a unique match
+          kHRU=jHRU
+          exit dsHRU
+        end if  ! (check there is a unique match)
+      end if  ! (if identified a downslope HRU)
+    end do dsHRU
+
+    ! if lateral flows are active, add inflow to the downslope HRU
+    if(kHRU > 0)then  ! if there is a downslope HRU
+      gru_data%hru(kHRU)%fluxStruct%var(iLookFLUX%mLayerColumnInflow)%dat(:) = &
+          gru_data%hru(kHRU)%fluxStruct%var(iLookFLUX%mLayerColumnInflow)%dat(:) + &
+          gru_data%hru(iHRU)%fluxStruct%var(iLookFLUX%mLayerColumnOutflow)%dat(:)
+
+    ! otherwise just increment basin (GRU) column outflow (m3 s-1) with the hru fraction
+    else
+      gru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1) = & 
+          gru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1) + &
+          sum( gru_data%hru(iHRU)%fluxStruct%var(iLookFLUX%mLayerColumnOutflow)%dat(:))
+    end if
+
+
+    ! ----- calculate weighted basin (GRU) fluxes --------------------------------------------------------------------------------------
+    
+    ! increment basin surface runoff (m s-1)
+    gru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) = &
+        gru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + &
+        gru_data%hru(iHRU)%fluxStruct%var(iLookFLUX%scalarSurfaceRunoff)%dat(1) * &
+        fracHRU
+    
+    !increment basin soil drainage (m s-1)
+    gru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1) = &
+        gru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1) + & 
+        gru_data%hru(iHRU)%fluxStruct%var(iLookFLUX%scalarSoilDrainage)%dat(1) * &
+        fracHRU
+    
+    ! increment aquifer variables -- ONLY if aquifer baseflow is computed individually for each HRU and aquifer is run
+    ! NOTE: groundwater computed later for singleBasin
+    if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == localColumn .and. &
+       model_decisions(iLookDECISIONS%groundwatr)%iDecision == bigBucket) then
+
+      gru_data%bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)  = &
+          gru_data%bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1) + &
+          gru_data%hru(iHRU)%fluxStruct%var(iLookFLUX%scalarSoilDrainage)%dat(1) * &
+          fracHRU
+      gru_data%bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = &
+          gru_data%bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1) +& 
+          gru_data%hru(iHRU)%fluxStruct%var(iLookFLUX%scalarAquiferTranspire)%dat(1) * &
+          fracHRU
+      gru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1) = & 
+          gru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1) + &
+          gru_data%hru(iHRU)%fluxStruct%var(iLookFLUX%scalarAquiferBaseflow)%dat(1) &
+          * fracHRU
+    end if
+  end do
+  ! ***********************************************************************************************************************
+  ! ********** END LOOP THROUGH HRUS **************************************************************************************
+  ! ***********************************************************************************************************************
+  ! perform the routing
+  associate(totalArea => gru_data%bvarStruct%var(iLookBVAR%basin__totalArea)%dat(1) )
+
+  ! compute water balance for the basin aquifer
+  if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == singleBasin)then
+    message=trim(message)//'multi_driver/bigBucket groundwater code not transferred from old code base yet'
+    err=20; call f_c_string_ptr(trim(message), message_r); return
+  end if
+
+  ! calculate total runoff depending on whether aquifer is connected
+  if(model_decisions(iLookDECISIONS%groundwatr)%iDecision == bigBucket) then
+    ! aquifer
+    gru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = &
+        gru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + &
+        gru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + &
+        gru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)
+  else
+    ! no aquifer
+    gru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = &
+        gru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + &
+        gru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + &
+        gru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)
+  endif
+
+  call qOverland(&! input
+                  model_decisions(iLookDECISIONS%subRouting)%iDecision,            &  ! intent(in): index for routing method
+                  gru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1),             &  ! intent(in): total runoff to the channel from all active components (m s-1)
+                  gru_data%bvarStruct%var(iLookBVAR%routingFractionFuture)%dat,             &  ! intent(in): fraction of runoff in future time steps (m s-1)
+                  gru_data%bvarStruct%var(iLookBVAR%routingRunoffFuture)%dat,               &  ! intent(in): runoff in future time steps (m s-1)
+                  ! output
+                  gru_data%bvarStruct%var(iLookBVAR%averageInstantRunoff)%dat(1),           &  ! intent(out): instantaneous runoff (m s-1)
+                  gru_data%bvarStruct%var(iLookBVAR%averageRoutedRunoff)%dat(1),            &  ! intent(out): routed runoff (m s-1)
+                  err,message)                                                                  ! intent(out): error control
+  if(err/=0)then; err=20; message=trim(message)//trim(cmessage); print*, message; return; endif;
+  end associate
+
+  ! update hru's bvarStruct with the basin's bvarStruct
+  do iHRU = 1, size(gru_data%hru)
+    do iVar=1, size(gru_data%bvarStruct%var(:))
+      gru_data%hru(iHRU)%bvarStruct%var(iVar)%dat(:) = gru_data%bvarStruct%var(iVar)%dat(:)
+    end do
+  end do
+
+end subroutine runGRU_fortran
+
+subroutine writeGRUOutput_fortran(indx_gru, timestep, outputstep, &
+    handle_gru_data, err, message_r) bind(C, name="writeGRUOutput_fortran")
+  USE actor_data_types,only:gru_type
+  USE HRUwriteoOutput_module,only:writeHRUOutput
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  implicit none
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  integer(c_int), intent(in)       :: timestep
+  integer(c_int), intent(in)       :: outputstep
+  type(c_ptr),    intent(in),value :: handle_gru_data
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  integer(i4b)                     :: iHRU
+  type(gru_type),pointer           :: gru_data
+  character(len=256)               :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_gru_data, gru_data)
+
+  do iHRU = 1, size(gru_data%hru)
+    call writeHRUOutput(indx_gru, iHRU, timestep, outputstep, gru_data%hru(iHRU), & 
+                        err, message)
+    if(err /= 0) then; call f_c_string_ptr(trim(message), message_r);return; end if
+  end do
+
+end subroutine writeGRUOutput_fortran
+
+end module gru_interface
diff --git a/build/source/hru_actor/hru.cpp b/build/source/hru_actor/hru.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/build/source/hru_actor/hru_actor.cpp b/build/source/hru_actor/hru_actor.cpp
index d99b698e8eac3d91578de7a80c804c3adfa101f1..25d101ff13a384a74664ae0303483c86d047031e 100644
--- a/build/source/hru_actor/hru_actor.cpp
+++ b/build/source/hru_actor/hru_actor.cpp
@@ -82,7 +82,7 @@ behavior hru_actor(stateful_actor<hru_state>* self, int ref_gru, int indx_gru,
       // Our output structure is full
       if (self->state.num_steps_until_write <= 0) {
           self->send(self->state.file_access_actor, write_output_v, 
-          self->state.indxGRU, self->state.indxHRU, self);
+                     self->state.indxGRU, self->state.indxHRU, self);
       }
     },
 
@@ -96,7 +96,8 @@ behavior hru_actor(stateful_actor<hru_state>* self, int ref_gru, int indx_gru,
       int err;
       self->state.iFile = iFile;
       self->state.stepsInCurrentFFile = num_forcing_steps_in_iFile;
-      setTimeZoneOffset(&self->state.iFile, self->state.hru_data, &err);
+      std::unique_ptr<char[]> message(new char[256]);
+      setTimeZoneOffset_fortran(iFile, self->state.hru_data, err, &message);
       if (err != 0) {
         aout(self) << "Error: HRU_Actor - setTimeZoneOffset - HRU = " 
                    << self->state.indxHRU << " - indxGRU = " 
@@ -126,7 +127,8 @@ behavior hru_actor(stateful_actor<hru_state>* self, int ref_gru, int indx_gru,
                      << " Updating timeZoneOffset\n";
       int err;
       self->state.iFile = iFile;
-      setTimeZoneOffset(&iFile, self->state.hru_data, &err);
+      std::unique_ptr<char[]> message(new char[256]);
+      setTimeZoneOffset_fortran(iFile, self->state.hru_data, err, &message);
     },
 
     // BMI - Functions
@@ -169,8 +171,9 @@ behavior hru_actor(stateful_actor<hru_state>* self, int ref_gru, int indx_gru,
 
 void Initialize_HRU(stateful_actor<hru_state>* self) {
   int err = 0;
-  initHRU(&self->state.indxGRU, &self->state.num_steps, self->state.hru_data, 
-          &err);
+  std::unique_ptr<char[]> message(new char[256]);
+  initHRU_fortran(self->state.indxGRU, self->state.indxHRU, 
+                  self->state.num_steps, self->state.hru_data, err, &message);
   if (err != 0) {
     aout(self) << "Error: HRU_Actor - Initialize - HRU = " 
                << self->state.indxHRU  
@@ -179,9 +182,10 @@ void Initialize_HRU(stateful_actor<hru_state>* self) {
                << "\nError Code = " << err << "\n";
     self->quit();
   }
-
-  setupHRUParam(&self->state.indxGRU, &self->state.indxHRU, 
-                self->state.hru_data, &self->state.upArea, &err);
+  
+  std::fill(message.get(), message.get() + 256, '\0'); // Clear message
+  setupHRU_fortran(self->state.indxGRU, self->state.indxHRU, 
+                   self->state.hru_data, err, &message);
   if (err != 0) {
     aout(self) << "Error: HRU_Actor - SetupHRUParam - HRU = " 
                 << self->state.indxHRU
@@ -190,9 +194,10 @@ void Initialize_HRU(stateful_actor<hru_state>* self) {
     self->quit();
     return;
   }
-          
-  summa_readRestart(&self->state.indxGRU, &self->state.indxHRU,
-                    self->state.hru_data, &self->state.dt_init, &err);
+
+  std::fill(message.get(), message.get() + 256, '\0'); // Clear message
+  readHRURestart_fortran(self->state.indxGRU, self->state.indxHRU,
+                        self->state.hru_data, err, &message);
   if (err != 0) {
     aout(self) << "Error: HRU_Actor - summa_readRestart - HRU = " 
                << self->state.indxHRU
@@ -213,9 +218,11 @@ void Initialize_HRU(stateful_actor<hru_state>* self) {
 
 int Run_HRU(stateful_actor<hru_state>* self) {
   int err = 0;
-  HRU_readForcing(&self->state.indxGRU, &self->state.timestep, 
-                  &self->state.forcingStep, &self->state.iFile, 
-                  self->state.hru_data, &err);
+  std::unique_ptr<char[]> message(new char[256]);
+  readHRUForcing_fortran(self->state.indxGRU, self->state.indxHRU, 
+                         self->state.timestep, self->state.forcingStep, 
+                         self->state.iFile, self->state.hru_data, err,
+                         &message);
   if (err != 0) {
     aout(self) << "Error---HRU_Actor: ReadForcingHRU\n" 
                << "\tIndxGRU = " << self->state.indxGRU << "\n"
@@ -235,10 +242,11 @@ int Run_HRU(stateful_actor<hru_state>* self) {
     aout(self) << self->state.ref_gru << " - Timestep = " 
                << self->state.timestep << "\n";
   }
-    
-  RunPhysics(&self->state.indxHRU, &self->state.timestep, self->state.hru_data, 
-             &self->state.dt_init,  &self->state.dt_init_factor, 
-             &self->state.walltime_timestep, &err);
+
+  std::fill(message.get(), message.get() + 256, '\0'); // Clear message
+  runHRU_fortran(self->state.indxGRU, self->state.indxHRU, self->state.timestep, 
+                 self->state.hru_data, self->state.dt_init_factor, 
+                 self->state.walltime_timestep, err, &message);
   if (err != 0) {
     aout(self) << "Error---RunPhysics:\n"
                << "\tIndxGRU = "  << self->state.indxGRU 
@@ -252,11 +260,11 @@ int Run_HRU(stateful_actor<hru_state>* self) {
     // write variables for output to summa_struct, extract y,m,d,h from 
     // fortran side ald save it to the hru's state
     int y,m,d,h;
-
-    hru_writeOutput(&self->state.indxHRU, &self->state.indxGRU,
-                    &self->state.timestep, 
-                    &self->state.output_structure_step_index,
-                    self->state.hru_data, &y, &m, &d, &h, &err);
+    std::fill(message.get(), message.get() + 256, '\0'); // Clear message
+    writeHRUOutput_fortran(self->state.indxGRU, self->state.indxHRU,
+                           self->state.timestep, 
+                           self->state.output_structure_step_index,
+                           self->state.hru_data, y, m, d, h, err, &message);
     if (err != 0) {
       aout(self) << "Error: HRU_Actor - writeHRUToOutputStructure - HRU = " 
                  << self->state.indxHRU << " - indxGRU = " 
@@ -301,7 +309,8 @@ int Run_HRU(stateful_actor<hru_state>* self) {
   return 0;      
 }
 
-// given a hru_actor with a state, compared current date with starting date to deterimine if hru is on a checkpoint
+// given a hru_actor with a state, compared current date with starting date
+// to deterimine if hru is on a checkpoint
 bool isCheckpoint(stateful_actor<hru_state>* self){
   switch(self->state.restartFrequency){
     case 0: // restart not enabled
diff --git a/build/source/hru_actor/hru_init.f90 b/build/source/hru_actor/hru_init.f90
index d0be85341077a36e8f54e750212cb7551cf0079a..2a7c7a73b0c64d1052027d7fa36276f925aba158 100755
--- a/build/source/hru_actor/hru_init.f90
+++ b/build/source/hru_actor/hru_init.f90
@@ -41,13 +41,11 @@ USE var_lookup,only:maxVarFreq                               ! # of available ou
 USE var_lookup,only:iLookATTR                               ! look-up values for local attributes
 USE var_lookup,only:iLookTYPE                               ! look-up values for classification of veg, soils etc.
 USE var_lookup,only:iLookPARAM                              ! look-up values for local column model parameters
-USE var_lookup,only:iLookID                                 ! look-up values for local column model parameters
+USE var_lookup,only:iLookID                                   ! look-up values for local column model parameters
 
 USE var_lookup,only:iLookPROG                               ! look-up values for local column model prognostic (state) variables
 USE var_lookup,only:iLookDIAG                               ! look-up values for local column model diagnostic variables
 USE var_lookup,only:iLookFLUX                               ! look-up values for local column model fluxes
-USE var_lookup,only:iLookBVAR                               ! look-up values for basin-average model variables
-USE var_lookup,only:iLookDECISIONS                          ! look-up values for model decisions
 USE globalData,only:urbanVegCategory                        ! vegetation category for urban areas
 
 ! named variables to define LAI decisions
@@ -59,16 +57,13 @@ USE mDecisions_module,only:&
 implicit none
 private
 public::initHRU
-public::setupHRUParam
-public::summa_readRestart
+public::setupHRU
+public::readHRURestart
 contains
 ! **************************************************************************************************
 ! public subroutine initHRU: ! used to declare and allocate summa data structures and initialize model state to known values
 ! **************************************************************************************************
-subroutine initHRU(indxGRU,            & !  Index of HRU's GRU parent
-                   num_steps,          &
-                   handle_hru_data,    &
-                   err) bind(C,name='initHRU')
+subroutine initHRU(indx_gru, indx_hru, hru_data, err, message)
   ! ---------------------------------------------------------------------------------------
   ! * desired modules
   ! ---------------------------------------------------------------------------------------
@@ -84,35 +79,20 @@ subroutine initHRU(indxGRU,            & !  Index of HRU's GRU parent
   ! miscellaneous global data
   USE globalData,only:gru_struc                               ! gru-hru mapping structures
   USE globalData,only:structInfo                              ! information on the data structures
-  USE globalData,only:numtim
   USE globalData,only:startTime,finshTime,refTime,oldTime
 
   USE var_lookup,only:maxvarFreq                              ! maximum number of output files
   USE var_lookup,only:iLookFreq                               ! output frequency lookup table
   implicit none
-  
-  ! ---------------------------------------------------------------------------------------
-  ! * variables from C++
-  ! ---------------------------------------------------------------------------------------
-  integer(c_int),intent(in)                  :: indxGRU                    ! indx of the parent GRU
-  integer(c_int),intent(out)                 :: num_steps                  ! number of steps in model, local to the HRU                 
-  type(c_ptr), intent(in), value             :: handle_hru_data            ! hru data structure (hru_type
-  ! ancillary data structures
-  integer(c_int),intent(inout)               :: err  
-  ! ---------------------------------------------------------------------------------------
-  ! * Fortran Variables For Conversion
-  ! ---------------------------------------------------------------------------------------
-  type(hru_type),pointer                     :: hru_data                   ! hru data structure (hru_type
-  ! ---------------------------------------------------------------------------------------
-  ! * Local Subroutine Variables
-  ! ---------------------------------------------------------------------------------------
-  character(LEN=256)                         :: message                    ! error message
-  character(LEN=256)                         :: cmessage                   ! error message of downwind routine
-  integer(i4b)                               :: iStruct                    ! looping variables
-  ! ---------------------------------------------------------------------------------------
-  ! * Convert From C++ to Fortran
-  ! ---------------------------------------------------------------------------------------
-  call c_f_pointer(handle_hru_data,   hru_data)
+  ! Dummy Variables
+  integer(c_int),intent(in)                  :: indx_gru      ! indx of the parent GRU
+  integer(c_int),intent(in)                  :: indx_hru      ! indx of the HRU
+  type(hru_type),intent(out)                 :: hru_data      ! hru data structure (hru_type
+  integer(c_int),intent(out)                 :: err  
+  character(len=256),intent(out)             :: message       ! error message
+  ! Local Variables
+  character(LEN=256)                         :: cmessage      ! error message of downwind routine
+  integer(i4b)                               :: iStruct       ! looping variables
   ! ---------------------------------------------------------------------------------------
   ! initialize error control
   err=0; message='hru_init/'
@@ -125,13 +105,9 @@ subroutine initHRU(indxGRU,            & !  Index of HRU's GRU parent
   elapsedWrite=0._dp
   elapsedPhysics=0._dp
 
-  ! copy the number of the steps for the hru
-  num_steps = numtim
-
   ! *****************************************************************************
   ! *** allocate space for data structures
-  ! *****************************************************************************
-
+  ! ****************************************************************************
   ! allocate time structures
   do iStruct=1,4
   select case(iStruct)
@@ -152,8 +128,8 @@ subroutine initHRU(indxGRU,            & !  Index of HRU's GRU parent
 
   ! get the number of snow and soil layers
   associate(&
-  nSnow => gru_struc(indxGRU)%hruInfo(1)%nSnow, & ! number of snow layers for each HRU
-  nSoil => gru_struc(indxGRU)%hruInfo(1)%nSoil  ) ! number of soil layers for each HRU
+  nSnow => gru_struc(indx_gru)%hruInfo(indx_hru)%nSnow, & ! number of snow layers for each HRU
+  nSoil => gru_struc(indx_gru)%hruInfo(indx_hru)%nSoil  ) ! number of soil layers for each HRU
 
   ! allocate other data structures
   do iStruct=1,size(structInfo)
@@ -236,20 +212,14 @@ subroutine initHRU(indxGRU,            & !  Index of HRU's GRU parent
 end subroutine initHRU
 
 
-
 ! **************************************************************************************************
 ! public subroutine setupHRUParam: initializes parameter data structures (e.g. vegetation and soil parameters).
 ! **************************************************************************************************
-subroutine setupHRUParam(indxGRU,                 & ! ID of hru
-                         indxHRU,                 & ! Index of the parent GRU of the HRU 
-                         handle_hru_data,         &
-                         upArea,                  & ! area upslope of each HRU,
-                         err) bind(C, name='setupHRUParam')
+subroutine setupHRU(indxGRU, indxHRU, hru_data, err, message)
   ! ---------------------------------------------------------------------------------------
   ! * desired modules
   ! ---------------------------------------------------------------------------------------
   USE nrtype                                                  ! variable types, etc.
-  USE output_structure_module,only:summa_struct
   USE summa_init_struc,only:init_struc
   ! subroutines and functions
   use time_utils_module,only:elapsedSec                       ! calculate the elapsed time
@@ -281,26 +251,19 @@ subroutine setupHRUParam(indxGRU,                 & ! ID of hru
   ! calling variables
   integer(c_int),intent(in)                :: indxGRU              ! Index of the parent GRU of the HRU
   integer(c_int),intent(in)                :: indxHRU              ! ID to locate correct HRU from netcdf file 
-  type(c_ptr), intent(in), value           :: handle_hru_data      ! pointer to the hru data structure (for error messages
-  real(c_double),intent(inout)             :: upArea
+  type(hru_type),intent(out)               :: hru_data             ! local hru data structure
   integer(c_int),intent(inout)             :: err
+  character(len=256),intent(out)           :: message
 
   ! local variables
-  type(hru_type),pointer                   :: hru_data             ! local hru data structure
 
   integer(i4b)                             :: ivar                 ! loop counter
   integer(i4b)                             :: i_z                  ! loop counter
-  character(len=256)                       :: message              ! error message
   character(len=256)                       :: cmessage             ! error message of downwind routine
 
   ! ---------------------------------------------------------------------------------------
   ! initialize error control
-  err=0; message='setupHRUParam/'
-
-  ! convert to fortran pointer from C++ pointer
-  call c_f_pointer(handle_hru_data, hru_data)
-
-  ! ffile_info and mDecisions moved to their own seperate subroutine call
+  err=0; message='setupHRU'
 
   hru_data%oldTime_hru%var(:) = hru_data%startTime_hru%var(:)
   hru_data%attrStruct%var(:) = init_struc%attrStruct%gru(indxGRU)%hru(indxHRU)%var(:)
@@ -327,26 +290,16 @@ subroutine setupHRUParam(indxGRU,                 & ! ID of hru
   do ivar=1, size(init_struc%indxStruct%gru(indxGRU)%hru(indxHRU)%var(:))
     hru_data%indxStruct%var(ivar)%dat(:) = init_struc%indxStruct%gru(indxGRU)%hru(indxHRU)%var(ivar)%dat(:)
   enddo
-end subroutine setupHRUParam
+end subroutine setupHRU
 
 
 ! **************************************************************************************************
 ! public subroutine summa_readRestart: read restart data and reset the model state
 ! **************************************************************************************************
-subroutine summa_readRestart(indxGRU,         & ! index of GRU in gru_struc
-                            indxHRU,          & ! index of HRU in gru_struc
-                            handle_hru_data,  & ! data structure for the HRU
-                            dt_init,          & ! used to initialize the length of the sub-step for each HRU
-                            err) bind(C,name='summa_readRestart')
-  ! ---------------------------------------------------------------------------------------
-  ! * desired modules
-  ! ---------------------------------------------------------------------------------------
-  ! data types
+subroutine readHRURestart(indxGRU, indxHRU, hru_data, err, message)
   USE nrtype                                                  ! variable types, etc.
   ! functions and subroutines
   USE time_utils_module,only:elapsedSec                       ! calculate the elapsed time
-  ! USE read_icond_module,only:read_icond               ! module to read initial conditions
-  ! USE check_icond4chm_module,only:check_icond4chm             ! module to check initial conditions
   USE var_derive_module,only:calcHeight                       ! module to calculate height at layer interfaces and layer mid-point
   USE var_derive_module,only:v_shortcut                       ! module to calculate "short-cut" variables
   USE var_derive_module,only:rootDensty                       ! module to calculate the vertical distribution of roots
@@ -356,40 +309,29 @@ subroutine summa_readRestart(indxGRU,         & ! index of GRU in gru_struc
   ! timing variables
   USE globalData,only:startRestart,endRestart                 ! date/time for the start and end of reading model restart files
   USE globalData,only:elapsedRestart                          ! elapsed time to read model restart files
+  ! Lookup values
+  USE var_lookup,only:iLookDECISIONS                          ! look-up values for model decisions
+  USE var_lookup,only:iLookBVAR                               ! look-up values for basin-average model variables
   ! model decisions
   USE mDecisions_module,only:&                                ! look-up values for the choice of method for the spatial representation of groundwater
   localColumn, & ! separate groundwater representation in each local soil column
   singleBasin    ! single groundwater store over the entire basin
   implicit none
-  ! ---------------------------------------------------------------------------------------
   ! Dummy variables
-  ! ---------------------------------------------------------------------------------------
   integer(c_int),intent(in)               :: indxGRU            !  index of GRU in gru_struc
   integer(c_int),intent(in)               :: indxHRU            !  index of HRU in gru_struc
-  type(c_ptr), intent(in), value          :: handle_hru_data   !  data structure for the HRU
-
-  real(c_double), intent(inout)           :: dt_init
-  integer(c_int), intent(inout)           :: err
-  ! ---------------------------------------------------------------------------------------
-  ! Fortran Pointers
-  ! ---------------------------------------------------------------------------------------
-  type(hru_type),pointer                  :: hru_data
-  ! ---------------------------------------------------------------------------------------
+  type(hru_type),intent(out)              :: hru_data
+  integer(c_int), intent(out)             :: err
+  character(len=256),intent(out)          :: message
   ! local variables
-  ! ---------------------------------------------------------------------------------------
   integer(i4b)                            :: ivar               ! index of variable
-  character(len=256)                      :: message            ! error message
   character(LEN=256)                      :: cmessage           ! error message of downwind routine
   character(LEN=256)                      :: restartFile        ! restart file name
   integer(i4b)                            :: nGRU
   ! ---------------------------------------------------------------------------------------
-
-  call c_f_pointer(handle_hru_data, hru_data)
-
   ! initialize error control
   err=0; message='hru_actor_readRestart/'
 
-
   ! *****************************************************************************
   ! *** compute ancillary variables
   ! *****************************************************************************
@@ -465,9 +407,9 @@ subroutine summa_readRestart(indxGRU,         & ! index of GRU in gru_struc
   ! *****************************************************************************
 
   ! initialize time step length
-  dt_init = hru_data%progStruct%var(iLookPROG%dt_init)%dat(1) ! seconds
+  hru_data%dt_init = hru_data%progStruct%var(iLookPROG%dt_init)%dat(1) ! seconds
 
-end subroutine summa_readRestart
+end subroutine readHRURestart
 
 ! Set the HRU's relative and absolute tolerances
 subroutine setIDATolerances(handle_hru_data,    &
diff --git a/build/source/hru_actor/hru_interface.f90 b/build/source/hru_actor/hru_interface.f90
new file mode 100644
index 0000000000000000000000000000000000000000..1dc3585d485dab7596fd9a82ed3b8426bfe80c55
--- /dev/null
+++ b/build/source/hru_actor/hru_interface.f90
@@ -0,0 +1,290 @@
+module hru_interface
+  USE,intrinsic :: iso_c_binding
+  USE nrtype
+
+  implicit none
+  private
+  public::initHRU_fortran
+  public::setupHRU_fortran
+  public::readHRURestart_fortran
+  public::setTimeZoneOffset_fortran
+  public::readHRUForcing_fortran
+  public::runHRU_fortran
+  public::writeHRUOutput_fortran
+
+  contains
+
+subroutine initHRU_fortran(indx_gru, indx_hru, num_steps, handle_hru_data, &
+    err, message_r) bind(C, name="initHRU_fortran")
+  USE globalData,only:numtim
+  USE actor_data_types,only:hru_type             
+  USE INIT_HRU_ACTOR,only:initHRU
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  implicit none
+  ! Dummy Variables
+  integer(c_int),intent(in)       :: indx_gru
+  integer(c_int),intent(in)       :: indx_hru
+  integer(c_int),intent(out)      :: num_steps
+  type(c_ptr),   intent(in),value :: handle_hru_data
+  integer(c_int),intent(out)      :: err
+  type(c_ptr),   intent(out)      :: message_r
+  ! Local Variables
+  type(hru_type),pointer         :: hru_data
+  character(len=256)             :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_hru_data, hru_data)
+
+  num_steps = numtim
+  call initHRU(indx_gru, indx_hru, hru_data, err, message)
+  if(err /= 0) then; call f_c_string_ptr(trim(message), message_r); end if
+end subroutine initHRU_fortran
+
+subroutine setupHRU_fortran(indx_gru, indx_hru, handle_hru_data, err, message_r) & 
+    bind(C, name="setupHRU_fortran")
+  USE actor_data_types,only:hru_type
+  USE INIT_HRU_ACTOR,only:setupHRU
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  integer(c_int), intent(in)       :: indx_hru
+  type(c_ptr),    intent(in),value :: handle_hru_data
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  type(hru_type),pointer         :: hru_data
+  character(len=256)             :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_hru_data, hru_data)
+
+  call setupHRU(indx_gru, indx_hru, hru_data, err, message)
+  if(err /= 0) then; call f_c_string_ptr(trim(message), message_r); end if
+end subroutine setupHRU_fortran
+
+subroutine readHRURestart_fortran(indx_gru, indx_hru, handle_hru_data, &
+    err, message_r) bind(C, name="readHRURestart_fortran")
+  USE actor_data_types,only:hru_type
+  USE INIT_HRU_ACTOR,only:readHRURestart
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  integer(c_int), intent(in)       :: indx_hru
+  type(c_ptr),    intent(in),value :: handle_hru_data
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  type(hru_type),pointer           :: hru_data
+  character(len=256)               :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_hru_data, hru_data)
+
+  call readHRURestart(indx_gru, indx_hru, hru_data, err, message)
+  if(err /= 0) then; call f_c_string_ptr(trim(message), message_r); end if
+end subroutine readHRURestart_fortran
+
+subroutine setTimeZoneOffset_fortran(iFile, handle_hru_data, err, message_r) &
+    bind(C, name="setTimeZoneOffset_fortran")
+  USE actor_data_types,only:hru_type
+  USE hru_read,only:setTimeZoneOffset
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  ! Dummy Variables
+  integer(c_int), intent(in)        :: iFile
+  type(c_ptr),    intent(in),value  :: handle_hru_data
+  integer(c_int), intent(out)       :: err
+  type(c_ptr),    intent(out)       :: message_r
+  ! Local Variables
+  type(hru_type),pointer            :: hru_data
+  character(len=256)                :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_hru_data, hru_data)
+
+  call setTimeZoneOffset(iFile, hru_data, err, message)
+  if(err /= 0) then; call f_c_string_ptr(trim(message), message_r); end if
+end subroutine setTimeZoneOffset_fortran
+
+subroutine readHRUForcing_fortran(indx_gru, indx_hru, iStep, iRead, iFile, &
+    handle_hru_data, err, message_r) bind(C, name="readHRUForcing_fortran")
+  USE actor_data_types,only:hru_type
+  USE hru_read,only:readHRUForcing
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  integer(c_int), intent(in)       :: indx_hru
+  integer(c_int), intent(in)       :: iStep
+  integer(c_int), intent(inout)    :: iRead
+  integer(c_int), intent(in)       :: iFile
+  type(c_ptr),    intent(in),value :: handle_hru_data
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  type(hru_type),pointer           :: hru_data
+  character(len=256)               :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_hru_data, hru_data)
+
+  call readHRUForcing(indx_gru, indx_hru, iStep, iRead, iFile, hru_data, &
+                      err, message)
+  if(err /= 0) then; call f_c_string_ptr(trim(message), message_r); end if
+end subroutine readHRUForcing_fortran
+
+subroutine runHRU_fortran(indx_gru, indx_hru, modelTimeStep, handle_hru_data, &
+    dt_init_factor, wallTimeTimeStep, err, message_r) bind(C, name="runHRU_fortran")
+  USE actor_data_types,only:hru_type
+  USE summa_modelRun,only:runPhysics
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  USE globalData,only:model_decisions          ! model decision structure
+  USE qTimeDelay_module,only:qOverland         ! module to route water through an "unresolved" river network
+  
+  USE mDecisions_module,only:&               ! look-up values for LAI decisions
+      monthlyTable,& ! LAI/SAI taken directly from a monthly table for different vegetation classes
+      specified,&    ! LAI/SAI computed from green vegetation fraction and winterSAI and summerLAI parameters   
+      localColumn, & ! separate groundwater representation in each local soil column
+      singleBasin, & ! single groundwater store over the entire basin
+      bigBucket
+
+  USE var_lookup,only:iLookFLUX              ! look-up values for local column model fluxes
+  USE var_lookup,only:iLookBVAR              ! look-up values for basin-average model variables
+  USE var_lookup,only:iLookDIAG              ! look-up values for local column model diagnostic variables
+  USE var_lookup,only:iLookDECISIONS         ! look-up values for model decisions
+  USE var_lookup,only:iLookATTR              ! look-up values for local attributes
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  integer(c_int), intent(in)       :: indx_hru
+  integer(c_int), intent(in)       :: modelTimeStep
+  type(c_ptr),    intent(in),value :: handle_hru_data
+  integer(c_int), intent(in)       :: dt_init_factor
+  real(c_double), intent(out)      :: wallTimeTimeStep
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  type(hru_type),pointer           :: hru_data
+  character(len=256)               :: message = ""
+  character(len=256)               :: cmessage
+  real(rkind)                      :: fracHRU                ! fractional area of a given HRU (-)
+
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_hru_data, hru_data)
+
+  !****************************************************************************** 
+  !****************************** From run_oneGRU *******************************
+  !******************************************************************************
+  ! ----- basin initialization --------------------------------------------------------------------------------------------
+  ! initialize runoff variables
+  hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1)    = 0._dp  ! surface runoff (m s-1)
+  hru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)     = 0._dp 
+  hru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)    = 0._dp  ! outflow from all "outlet" HRUs (those with no downstream HRU)
+  hru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1)      = 0._dp 
+
+  ! initialize baseflow variables
+  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)  = 0._dp ! recharge to the aquifer (m s-1)
+  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  = 0._dp ! baseflow from the aquifer (m s-1)
+  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = 0._dp ! transpiration loss from the aquifer (m s-1)
+
+  ! initialize total inflow for each layer in a soil column
+  ! if (modelTimeStep == 0 .and. indxHRU == 1)then
+  hru_data%fluxStruct%var(iLookFLUX%mLayerColumnInflow)%dat(:) = 0._dp
+  ! end if
+
+
+  call runPhysics(indx_gru, indx_hru, modelTimeStep, hru_data, dt_init_factor, err, message)
+  if(err /= 0) then; call f_c_string_ptr(trim(message), message_r); return; end if
+
+  fracHRU = hru_data%attrStruct%var(iLookATTR%HRUarea) / hru_data%bvarStruct%var(iLookBVAR%basin__totalArea)%dat(1)
+
+
+
+  ! ----- calculate weighted basin (GRU) fluxes --------------------------------------------------------------------------------------
+  
+  ! increment basin surface runoff (m s-1)
+  hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) = hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + hru_data%fluxStruct%var(iLookFLUX%scalarSurfaceRunoff)%dat(1) * fracHRU
+  
+  !increment basin soil drainage (m s-1)
+  hru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)   = hru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)  + hru_data%fluxStruct%var(iLookFLUX%scalarSoilDrainage)%dat(1)  * fracHRU
+  
+  ! increment aquifer variables -- ONLY if aquifer baseflow is computed individually for each HRU and aquifer is run
+  ! NOTE: groundwater computed later for singleBasin
+  if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == localColumn .and. model_decisions(iLookDECISIONS%groundwatr)%iDecision == bigBucket) then
+
+    hru_data%bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)  = hru_data%bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)   + hru_data%fluxStruct%var(iLookFLUX%scalarSoilDrainage)%dat(1)     * fracHRU
+    hru_data%bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = hru_data%bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1)  + hru_data%fluxStruct%var(iLookFLUX%scalarAquiferTranspire)%dat(1) * fracHRU
+    hru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  =  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  &
+            +  hru_data%fluxStruct%var(iLookFLUX%scalarAquiferBaseflow)%dat(1) * fracHRU
+    end if
+
+  ! perform the routing
+  associate(totalArea => hru_data%bvarStruct%var(iLookBVAR%basin__totalArea)%dat(1) )
+
+  ! compute water balance for the basin aquifer
+  if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == singleBasin)then
+    message=trim(message)//'multi_driver/bigBucket groundwater code not transferred from old code base yet'
+    err=20; call f_c_string_ptr(trim(message), message_r); return
+  end if
+
+  ! calculate total runoff depending on whether aquifer is connected
+  if(model_decisions(iLookDECISIONS%groundwatr)%iDecision == bigBucket) then
+    ! aquifer
+    hru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + hru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + hru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)
+  else
+    ! no aquifer
+    hru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + hru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + hru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)
+  endif
+
+  call qOverland(&! input
+                  model_decisions(iLookDECISIONS%subRouting)%iDecision,            &  ! intent(in): index for routing method
+                  hru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1),             &  ! intent(in): total runoff to the channel from all active components (m s-1)
+                  hru_data%bvarStruct%var(iLookBVAR%routingFractionFuture)%dat,             &  ! intent(in): fraction of runoff in future time steps (m s-1)
+                  hru_data%bvarStruct%var(iLookBVAR%routingRunoffFuture)%dat,               &  ! intent(in): runoff in future time steps (m s-1)
+                  ! output
+                  hru_data%bvarStruct%var(iLookBVAR%averageInstantRunoff)%dat(1),           &  ! intent(out): instantaneous runoff (m s-1)
+                  hru_data%bvarStruct%var(iLookBVAR%averageRoutedRunoff)%dat(1),            &  ! intent(out): routed runoff (m s-1)
+                  err,message)                                                                  ! intent(out): error control
+  if(err/=0)then; err=20; message=trim(message)//trim(cmessage); call f_c_string_ptr(trim(message), message_r); return; endif;
+  end associate
+  wallTimeTimeStep = hru_data%diagStruct%var(iLookDIAG%wallClockTime)%dat(1)
+
+
+end subroutine runHRU_fortran
+
+subroutine writeHRUOutput_fortran(indx_gru, indx_hru, timestep, outputstep, &
+    handle_hru_data, y, m, d, h, err, message_r) bind(C, name="writeHRUOutput_fortran")
+  USE actor_data_types,only:hru_type
+  USE HRUwriteoOutput_module,only:writeHRUOutput
+  USE var_lookup,only:iLookTIME                 ! named variables for time data structure
+  USE C_interface_module,only:f_c_string_ptr  ! convert fortran string to c string
+  ! Dummy Variables
+  integer(c_int), intent(in)       :: indx_gru
+  integer(c_int), intent(in)       :: indx_hru
+  integer(c_int), intent(in)       :: timestep
+  integer(c_int), intent(in)       :: outputstep
+  type(c_ptr),    intent(in),value :: handle_hru_data
+  integer(c_int), intent(out)      :: y
+  integer(c_int), intent(out)      :: m
+  integer(c_int), intent(out)      :: d
+  integer(c_int), intent(out)      :: h
+  integer(c_int), intent(out)      :: err
+  type(c_ptr),    intent(out)      :: message_r
+  ! Local Variables
+  type(hru_type),pointer           :: hru_data
+  character(len=256)               :: message = ""
+
+  call f_c_string_ptr(trim(message), message_r)
+  call c_f_pointer(handle_hru_data, hru_data)
+
+  ! updating date variables to be passed back to the actors
+  y = hru_data%timeStruct%var(iLookTIME%iyyy)
+  m = hru_data%timeStruct%var(iLookTIME%im)
+  d = hru_data%timeStruct%var(iLookTIME%id)
+  h = hru_data%timeStruct%var(iLookTIME%ih)
+
+  call writeHRUOutput(indx_gru, indx_hru, timestep, outputstep, hru_data, &
+                      err, message)
+  if(err /= 0) then; call f_c_string_ptr(trim(message), message_r); end if
+
+
+end subroutine writeHRUOutput_fortran
+end module hru_interface
\ No newline at end of file
diff --git a/build/source/hru_actor/hru_modelRun.f90 b/build/source/hru_actor/hru_modelRun.f90
index b72f6b445d16e4dd25f5a4f1223eac9bd5a38fb9..97d962c38148b0571ba5b43ec422b64a3c2517e0 100644
--- a/build/source/hru_actor/hru_modelRun.f90
+++ b/build/source/hru_actor/hru_modelRun.f90
@@ -76,14 +76,8 @@ public::set_sundials_tolerances
 contains
 
 ! Runs the model physics for an HRU
-subroutine runPhysics(&
-              indxHRU,             &
-              modelTimeStep,       &
-              handle_hru_data,     &
-              dt_init,             & ! used to initialize the length of the sub-step for each HRU
-              dt_init_factor,      & ! used to adjust the length of the timestep in the event of a failure
-              wallTimeTimeStep,    &
-              err) bind(C, name='RunPhysics')
+subroutine runPhysics(indxGRU, indxHRU, modelTimeStep, hru_data, &
+    dt_init_factor, err, message)
   ! ---------------------------------------------------------------------------------------
   ! * desired modules
   ! ---------------------------------------------------------------------------------------
@@ -98,30 +92,21 @@ subroutine runPhysics(&
   USE coupled_em_module,only:coupled_em        ! module to run the coupled energy and mass model
   USE qTimeDelay_module,only:qOverland         ! module to route water through an "unresolved" river network
   ! global data
+  USE globalData,only:gru_struc
   USE globalData,only:model_decisions          ! model decision structure
   USE globalData,only:startPhysics,endPhysics  ! date/time for the start and end of the initialization
   USE globalData,only:elapsedPhysics           ! elapsed time for the initialization
   implicit none
- 
-  ! ---------------------------------------------------------------------------------------
   ! Dummy Variables
-  ! ---------------------------------------------------------------------------------------
-  integer(c_long),intent(in)                :: indxHRU                ! id of HRU                   
+  integer(c_int),intent(in)                 :: indxGRU                ! id of GRU
+  integer(c_int),intent(in)                 :: indxHRU                ! id of HRU                   
   integer(c_int), intent(in)                :: modelTimeStep          ! time step index
-  type(c_ptr),    intent(in), value         :: handle_hru_data        ! c_ptr to -- hru data
-  real(c_double), intent(inout)             :: dt_init                ! used to initialize the length of the sub-step for each HRU
+  type(hru_type), intent(inout)             :: hru_data               ! c_ptr to -- hru data
   integer(c_int), intent(in)                :: dt_init_factor         ! used to adjust the length of the timestep in the event of a failure
-  real(c_double), intent(out)               :: wallTimeTimeStep       ! wall time for the time step
   integer(c_int), intent(inout)             :: err                    ! error code
-  ! ---------------------------------------------------------------------------------------
-  ! FORTRAN POINTERS
-  ! ---------------------------------------------------------------------------------------
-  type(hru_type),pointer                    :: hru_data               ! hru data
-
-  ! ---------------------------------------------------------------------------------------
+  character(len=256), intent(out)           :: message                ! error message
   ! local variables: general
-  ! ---------------------------------------------------------------------------------------
-  real(dp)                                  :: fracHRU                ! fractional area of a given HRU (-)
+  integer(8)                                :: hruId                  ! hruId
   character(LEN=256)                        :: cmessage               ! error message of downwind routine
   ! local variables: veg phenology
   logical(lgt)                              :: computeVegFluxFlag     ! flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
@@ -131,10 +116,8 @@ subroutine runPhysics(&
   integer(i4b)                              :: nSoil                  ! number of soil layers
   integer(i4b)                              :: nLayers                ! total number of layers
   real(dp), allocatable                     :: zSoilReverseSign(:)    ! height at bottom of each soil layer, negative downwards (m)
-  character(len=256)                        :: message                ! error message
   ! ---------------------------------------------------------------------------------------
-
-  call c_f_pointer(handle_hru_data, hru_data)
+  hruId = gru_struc(indxGRU)%hruInfo(indxHRU)%hru_id
 
   ! ---------------------------------------------------------------------------------------
   ! initialize error control
@@ -178,28 +161,6 @@ subroutine runPhysics(&
   ! ****************************************************************************
   ! *** model simulation
   ! ****************************************************************************
-
-  
-  !****************************************************************************** 
-  !****************************** From run_oneGRU *******************************
-  !******************************************************************************
-  ! ----- basin initialization --------------------------------------------------------------------------------------------
-  ! initialize runoff variables
-  hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1)    = 0._dp  ! surface runoff (m s-1)
-  hru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)     = 0._dp 
-  hru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)    = 0._dp  ! outflow from all "outlet" HRUs (those with no downstream HRU)
-  hru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1)      = 0._dp 
-
-  ! initialize baseflow variables
-  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)  = 0._dp ! recharge to the aquifer (m s-1)
-  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  = 0._dp ! baseflow from the aquifer (m s-1)
-  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = 0._dp ! transpiration loss from the aquifer (m s-1)
-
-  ! initialize total inflow for each layer in a soil column
-  if (modelTimeStep == 0 .and. indxHRU == 1)then
-    hru_data%fluxStruct%var(iLookFLUX%mLayerColumnInflow)%dat(:) = 0._dp
-  end if
- 
   ! update the number of layers
   nSnow   = hru_data%indxStruct%var(iLookINDEX%nSnow)%dat(1)    ! number of snow layers
   nSoil   = hru_data%indxStruct%var(iLookINDEX%nSoil)%dat(1)    ! number of soil layers
@@ -213,7 +174,6 @@ subroutine runPhysics(&
   ! water pixel: do nothing
   if (hru_data%typeStruct%var(iLookTYPE%vegTypeIndex) == isWater) then
       ! Set wall_clock time to zero so it does not get a random value
-    wallTimeTimeStep = 0._dp
     hru_data%diagStruct%var(iLookDIAG%wallClockTime)%dat(1) = 0._dp 
     return
   endif
@@ -271,12 +231,12 @@ subroutine runPhysics(&
   ! run the model for a single HRU
   call coupled_em(&
                   ! model control
-                  indxHRU,                     & ! intent(in):    hruID
-                  dt_init,                     & ! intent(inout): initial time step
+                  hruId,                       & ! intent(in):    hruID
+                  hru_data%dt_init,            & ! intent(inout): initial time step
                   dt_init_factor,              & ! Used to adjust the length of the timestep in the event of a failure
                   computeVegFluxFlag,          & ! intent(inout): flag to indicate if we are computing fluxes over vegetation
-                  hru_data%fracJulDay,                  & ! intent(in):    fractional julian days since the start of year
-                  hru_data%yearLength,                  & ! intent(in):    number of days in the current year
+                  hru_data%fracJulDay,         & ! intent(in):    fractional julian days since the start of year
+                  hru_data%yearLength,         & ! intent(in):    number of days in the current year
                   ! data structures (input)
                   hru_data%typeStruct,         & ! intent(in):    local classification of soil veg etc. for each HRU
                   hru_data%attrStruct,         & ! intent(in):    local attributes for each HRU
@@ -301,59 +261,33 @@ subroutine runPhysics(&
   if(computeVegFluxFlag)      hru_data%ComputeVegFlux = yes
   if(.not.computeVegFluxFlag) hru_data%ComputeVegFlux = no
 
-  fracHRU = hru_data%attrStruct%var(iLookATTR%HRUarea) / hru_data%bvarStruct%var(iLookBVAR%basin__totalArea)%dat(1)
 
-  ! ----- calculate weighted basin (GRU) fluxes --------------------------------------------------------------------------------------
-  
-  ! increment basin surface runoff (m s-1)
-  hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) = hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + hru_data%fluxStruct%var(iLookFLUX%scalarSurfaceRunoff)%dat(1) * fracHRU
-  
-  !increment basin soil drainage (m s-1)
-  hru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)   = hru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)  + hru_data%fluxStruct%var(iLookFLUX%scalarSoilDrainage)%dat(1)  * fracHRU
-  
-  ! increment aquifer variables -- ONLY if aquifer baseflow is computed individually for each HRU and aquifer is run
-  ! NOTE: groundwater computed later for singleBasin
-  if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == localColumn .and. model_decisions(iLookDECISIONS%groundwatr)%iDecision == bigBucket) then
-
-    hru_data%bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)  = hru_data%bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1)   + hru_data%fluxStruct%var(iLookFLUX%scalarSoilDrainage)%dat(1)     * fracHRU
-    hru_data%bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = hru_data%bvarStruct%var(iLookBVAR%basin__AquiferTranspire)%dat(1)  + hru_data%fluxStruct%var(iLookFLUX%scalarAquiferTranspire)%dat(1) * fracHRU
-    hru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  =  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)  &
-            +  hru_data%fluxStruct%var(iLookFLUX%scalarAquiferBaseflow)%dat(1) * fracHRU
-    end if
-
-  ! perform the routing
-  associate(totalArea => hru_data%bvarStruct%var(iLookBVAR%basin__totalArea)%dat(1) )
-
-  ! compute water balance for the basin aquifer
-  if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == singleBasin)then
-    message=trim(message)//'multi_driver/bigBucket groundwater code not transferred from old code base yet'
-    err=20; print*, message; return
-  end if
 
-  ! calculate total runoff depending on whether aquifer is connected
-  if(model_decisions(iLookDECISIONS%groundwatr)%iDecision == bigBucket) then
-    ! aquifer
-    hru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + hru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + hru_data%bvarStruct%var(iLookBVAR%basin__AquiferBaseflow)%dat(1)
-  else
-    ! no aquifer
-    hru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = hru_data%bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + hru_data%bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + hru_data%bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1)
-  endif
 
-  call qOverland(&! input
-                  model_decisions(iLookDECISIONS%subRouting)%iDecision,            &  ! intent(in): index for routing method
-                  hru_data%bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1),             &  ! intent(in): total runoff to the channel from all active components (m s-1)
-                  hru_data%bvarStruct%var(iLookBVAR%routingFractionFuture)%dat,             &  ! intent(in): fraction of runoff in future time steps (m s-1)
-                  hru_data%bvarStruct%var(iLookBVAR%routingRunoffFuture)%dat,               &  ! intent(in): runoff in future time steps (m s-1)
-                  ! output
-                  hru_data%bvarStruct%var(iLookBVAR%averageInstantRunoff)%dat(1),           &  ! intent(out): instantaneous runoff (m s-1)
-                  hru_data%bvarStruct%var(iLookBVAR%averageRoutedRunoff)%dat(1),            &  ! intent(out): routed runoff (m s-1)
-                  err,message)                                                                  ! intent(out): error control
-  if(err/=0)then; err=20; message=trim(message)//trim(cmessage); print*, message; return; endif;
-  end associate
-  !************************************* End of run_oneGRU *****************************************
+  ! ----- compute fluxes across HRUs --------------------------------------------------------------------------------------------------
 
-  ! Get the elapsed time for the GRU
-  wallTimeTimeStep = hru_data%diagStruct%var(iLookDIAG%wallClockTime)%dat(1)
+  ! ! identify lateral connectivity
+  ! ! (Note:  for efficiency, this could this be done as a setup task, not every timestep)
+  ! kHRU = 0
+  ! ! identify the downslope HRU
+  ! dsHRU: do jHRU=1,gruInfo%hruCount
+  !   if(typeHRU%hru(iHRU)%var(iLookTYPE%downHRUindex) == idHRU%hru(jHRU)%var(iLookID%hruId))then
+  !     if(kHRU==0)then  ! check there is a unique match
+  !       kHRU=jHRU
+  !       exit dsHRU
+  !     end if  ! (check there is a unique match)
+  !   end if  ! (if identified a downslope HRU)
+  ! end do dsHRU
+
+  ! ! if lateral flows are active, add inflow to the downslope HRU
+  ! if(kHRU > 0)then  ! if there is a downslope HRU
+  !   fluxHRU%hru(kHRU)%var(iLookFLUX%mLayerColumnInflow)%dat(:) = fluxHRU%hru(kHRU)%var(iLookFLUX%mLayerColumnInflow)%dat(:)  + fluxHRU%hru(iHRU)%var(iLookFLUX%mLayerColumnOutflow)%dat(:)
+  ! ! otherwise just increment basin (GRU) column outflow (m3 s-1) with the hru fraction
+  ! else
+  !   bvarData%var(iLookBVAR%basin__ColumnOutflow)%dat(1) = bvarData%var(iLookBVAR%basin__ColumnOutflow)%dat(1) + sum(fluxHRU%hru(iHRU)%var(iLookFLUX%mLayerColumnOutflow)%dat(:))
+  ! end if
+
+  !************************************* End of run_oneGRU *****************************************
 
 end subroutine runPhysics
 
diff --git a/build/source/hru_actor/hru_read.f90 b/build/source/hru_actor/hru_read.f90
index 98b03a2b5f9167d00fadb5a704ab4c82e8ef3228..97165f132bbe9e08013d56df74baad37f76e7f35 100644
--- a/build/source/hru_actor/hru_read.f90
+++ b/build/source/hru_actor/hru_read.f90
@@ -14,7 +14,7 @@ USE data_types,only:&
 USE actor_data_types,only:hru_type
 implicit none
 public::setTimeZoneOffset
-public::HRU_readForcing
+public::readHRUForcing
 private::getFirstTimeStep
 
 
@@ -24,7 +24,7 @@ real(dp),parameter  :: smallOffset=1.e-8_rkind   ! small offset (units=days) to
 contains
 
 ! set the refTimeString and extract the time to set the tmZonOffsetFracDay
-subroutine setTimeZoneOffset(iFile, handle_hru_data, err) bind(C, name="setTimeZoneOffset")
+subroutine setTimeZoneOffset(iFile, hru_data, err, message)
   USE forcing_file_info,only:forcingDataStruct         ! forcing structure
   USE time_utils_module,only:extractTime        ! extract time info from units string
   USE time_utils_module,only:fracDay            ! compute fractional day
@@ -32,18 +32,16 @@ subroutine setTimeZoneOffset(iFile, handle_hru_data, err) bind(C, name="setTimeZ
   implicit none
 
   integer(c_int),intent(in)             :: iFile
-  type(c_ptr),intent(in),value          :: handle_hru_data  ! vector of time data for a given time step
+  type(hru_type)                        :: hru_data         !  model time data
   integer(c_int),intent(out)            :: err
+  character(len=256),intent(out)        :: message
 
   ! local variables
-  type(hru_type),pointer                :: hru_data         !  model time data
-  character(len=256)                    :: message
   character(len=256)                    :: cmessage
   integer(i4b)                          :: iyyy,im,id,ih,imin ! date
   integer(i4b)                          :: ih_tz,imin_tz      ! time zone information
   real(dp)                              :: dsec,dsec_tz       ! seconds
 
-  call c_f_pointer(handle_hru_data, hru_data)
   err=0; message="hru_actor.f90 - setForcingTimeInfo";
 
   ! define the reference time for the model simulation
@@ -51,7 +49,7 @@ subroutine setTimeZoneOffset(iFile, handle_hru_data, err) bind(C, name="setTimeZ
                    iyyy,im,id,ih,imin,dsec,                & ! output = year, month, day, hour, minute, second
                    ih_tz, imin_tz, dsec_tz,                & ! output = time zone information (hour, minute, second)
                    err,cmessage)                             ! output = error code and error message
-  if(err/=0)then; message=trim(message)//trim(cmessage); print*, "message"; return; end if
+  if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
   
   ! set the timezone offset
   select case(trim(NC_TIME_ZONE))
@@ -65,7 +63,8 @@ subroutine setTimeZoneOffset(iFile, handle_hru_data, err) bind(C, name="setTimeZ
 
 end subroutine setTimeZoneOffset
 
-subroutine HRU_readForcing(indxGRU, iStep, iRead, iFile, handle_hru_data, err) bind(C, name="HRU_readForcing")
+subroutine readHRUForcing(indx_gru, indx_hru, iStep, iRead, iFile, &
+    hru_data, err, message)
   USE multiconst,only:secprday                  ! number of seconds in a day
   USE time_utils_module,only:compJulDay         ! convert calendar date to julian day
   ! global Data
@@ -81,21 +80,25 @@ subroutine HRU_readForcing(indxGRU, iStep, iRead, iFile, handle_hru_data, err) b
   USE netcdf,only:nf90_max_name                                   ! used for nf90_max_name
   USE time_utils_module,only:compcalday                 ! convert julian day to calendar date
   USE globalData,only:refJulDay                 ! reference time (fractional julian days)
-
+  USE globalData,only:ixHRUfile_min             ! minimum index of HRU in the forcing file
+  USE globalData,only:ixHRUfile_max             ! maximum index of HRU in the forcing file
+  USE globalData,only:gru_struc                 ! GRU structure
   implicit none
 
-  integer(c_int),intent(in)               :: indxGRU          ! Index of the GRU in gru_struc
+  integer(c_int),intent(in)               :: indx_gru         ! Index of the GRU in gru_struc
+  integer(c_int),intent(in)               :: indx_hru         ! Index of the HRU in hru_struc
   integer(c_int),intent(in)               :: istep            ! Model Timestep
   integer(c_int),intent(inout)            :: iRead            ! Model Timestep 
   integer(c_int),intent(in)               :: iFile            ! index of current forcing file from forcing file list 
-  type(c_ptr),intent(in),value            :: handle_hru_data  ! vector of time data for a given time step
+  type(hru_type)                          :: hru_data         !  model time data
   integer(c_int),intent(out)              :: err              ! Model Timestep
+  character(len=256),intent(out)          :: message          ! error message
   ! local variables
-  type(hru_type),pointer                  :: hru_data         !  model time data
   real(dp)                                :: currentJulDay    ! Julian day of current time step
   real(dp)                                :: dataJulDay       ! julian day of current forcing data step being read
   real(dp)                                :: startJulDay      ! julian day at the start of the year
-  
+  integer(i4b)                            :: iHRU_global      ! index of HRU in the forcing file
+  integer(i4b)                            :: iHRU_local       ! index of HRU in the forcing file
   ! Counters
   integer(i4b)                            :: iline            ! loop through lines in the file
   integer(i4b)                            :: iVar
@@ -103,17 +106,17 @@ subroutine HRU_readForcing(indxGRU, iStep, iRead, iFile, handle_hru_data, err) b
   ! other
   logical(lgt),dimension(size(forc_meta)) :: checkForce       ! flags to check forcing data variables exist
   logical(lgt),parameter                  :: checkTime=.false.! flag to check the time
-  
   real(dp)                                :: dsec             ! double precision seconds (not used)
   real(dp),parameter                      :: dataMin=-1._dp   ! minimum allowable data value (all forcing variables should be positive)
   character(len = nf90_max_name)          :: varName          ! dimenison name
-
-  character(len=256)                      :: message          ! error message
   character(len=256)                      :: cmessage         ! error message
 
-  call c_f_pointer(handle_hru_data, hru_data)
   err=0;message="hru_actor.f90 - readForcingHRU";
 
+  ! Get index into the forcing structure
+  iHRU_global = gru_struc(indx_gru)%hruInfo(indx_hru)%hru_nc
+  iHRU_local  = (iHRU_global - ixHRUfile_min)+1
+
   if(istep == 1) then
     call getFirstTimestep(iFile, iRead, err)
     if(err/=0)then; message=trim(message)//"getFirstTimestep"; print*,message;return; end if
@@ -162,12 +165,12 @@ subroutine HRU_readForcing(indxGRU, iStep, iRead, iFile, handle_hru_data, err) b
     checkForce(iVar) = .true.
 
     ! check individual data value
-    if(forcingDataStruct(iFile)%var(ivar)%dataFromFile(indxGRU,iRead)<dataMin)then
+    if(forcingDataStruct(iFile)%var(ivar)%dataFromFile(indx_gru,iRead)<dataMin)then
       write(message,'(a,f13.5)') trim(message)//'forcing data for variable '//trim(varname)//' is less than minimum allowable value ', dataMin
       err=20; return
     endif
     ! put the data into structures
-    hru_data%forcStruct%var(ivar) = forcingDataStruct(iFile)%var(ivar)%dataFromFile(indxGRU,iRead)
+    hru_data%forcStruct%var(ivar) = forcingDataStruct(iFile)%var(ivar)%dataFromFile(indx_gru,iRead)
   end do  ! loop through forcing variables
     
   ! check if any forcing data is missing
@@ -227,11 +230,7 @@ subroutine HRU_readForcing(indxGRU, iStep, iRead, iFile, handle_hru_data, err) b
                                             hru_data%yearLength                             ! number of days in the current year
     !pause ' checking time'
   end if
-
-
-
-
-end subroutine HRU_readForcing 
+end subroutine readHRUForcing 
 
  ! Find the first timestep within the forcing file
 subroutine getFirstTimestep(iFile, iRead, err)
diff --git a/build/source/hru_actor/hru_utils.cpp b/build/source/hru_actor/hru_utils.cpp
index 87535b73fe70e60e1674ee748fac1c2a3eafea91..b197bcdddec457ab2d9d7ef794f6c9e0f1936748 100644
--- a/build/source/hru_actor/hru_utils.cpp
+++ b/build/source/hru_actor/hru_utils.cpp
@@ -19,8 +19,6 @@ void serializeHru(stateful_actor<hru_state>* self, hru& serialized_state) {
   serialized_state.dt_init_factor = self->state.dt_init_factor;
   serialized_state.output_structure_step_index = 
       self->state.output_structure_step_index; 
-  serialized_state.dt_init = self->state.dt_init;
-  serialized_state.upArea = self->state.upArea;
   serialized_state.rtol = self->state.rtol;
   serialized_state.atol = self->state.atol;
 
@@ -69,7 +67,9 @@ void serializeHru(stateful_actor<hru_state>* self, hru& serialized_state) {
   get_scalar_data(self->state.hru_data, serialized_state.frac_jul_day,
                   serialized_state.tm_zone_offset_frac_day,
                   serialized_state.year_length,
-                  serialized_state.compute_veg_flux);
+                  serialized_state.compute_veg_flux,
+                  serialized_state.dt_init,
+                  serialized_state.upArea);
 }
 
 void deserializeHru(stateful_actor<hru_state>* self, hru& new_state) {
@@ -89,8 +89,6 @@ void deserializeHru(stateful_actor<hru_state>* self, hru& new_state) {
   self->state.dt_init_factor = new_state.dt_init_factor;
   self->state.output_structure_step_index = 
       new_state.output_structure_step_index;
-  self->state.dt_init = new_state.dt_init;
-  self->state.upArea = new_state.upArea;
   self->state.rtol = new_state.rtol;
   self->state.atol = new_state.atol;
   // Statistic Structures
@@ -129,5 +127,6 @@ void deserializeHru(stateful_actor<hru_state>* self, hru& new_state) {
   // scalar data
   set_scalar_data(self->state.hru_data, new_state.frac_jul_day,
                   new_state.tm_zone_offset_frac_day, new_state.year_length,
-                  new_state.compute_veg_flux);
+                  new_state.compute_veg_flux, new_state.dt_init, 
+                  new_state.upArea);
 }
\ No newline at end of file
diff --git a/build/source/hru_actor/hru_writeOutput.f90 b/build/source/hru_actor/hru_writeOutput.f90
index 4e479d8127a4c3f4efaac10ac583e1d2cef190fd..7742200193b93f896dac1a378b23ac82e7401335 100644
--- a/build/source/hru_actor/hru_writeOutput.f90
+++ b/build/source/hru_actor/hru_writeOutput.f90
@@ -47,9 +47,6 @@ USE globalData,only:indxChild_map             ! index of the child data structur
 USE globalData,only:bvarChild_map             ! index of the child data structure: stats bvar
 USE globalData,only:outFreq                   ! output frequencies
 ! named variables
-USE var_lookup,only:iLookTIME                 ! named variables for time data structure
-USE var_lookup,only:iLookDIAG                 ! named variables for local column model diagnostic variables
-USE var_lookup,only:iLookPROG                 ! named variables for local column model prognostic variables
 USE var_lookup,only:iLookINDEX                ! named variables for local column index variables
 USE var_lookup,only:iLookFreq                 ! named variables for the frequency structure
 USE var_lookup,only:iLookBVAR                 ! named variables for basin parameters
@@ -62,7 +59,7 @@ USE output_structure_module,only:summa_struct
 
 implicit none
 private
-public::hru_writeOutput
+public::writeHRUOutput
 public::writeParm
 public::writeData
 public::writeBasin
@@ -71,13 +68,7 @@ public::writeRestart
 public::setFinalizeStatsFalse
 integer(i4b),parameter      :: maxSpectral=2              ! maximum number of spectral bands
 contains
-subroutine hru_writeOutput(&
-                            indxHRU,                   &
-                            indxGRU,                   &
-                            timestep,                  & ! model timestep
-                            outputStep,                & ! index into the output Struc
-                            handle_hru_data, y,m,d,h,  & ! local HRU data  
-                            err) bind(C, name="hru_writeOutput") 
+subroutine writeHRUOutput(indxGRU, indxHRU, timestep, outputStep, hru_data, err, message)
   USE nrtype
   USE globalData,only:structInfo
   USE globalData,only:startWrite,endWrite
@@ -99,25 +90,17 @@ subroutine hru_writeOutput(&
   USE netcdf_util_module,only:nc_file_close                   ! close netcdf file
   USE netcdf_util_module,only:nc_file_open                    ! open netcdf file
   USE var_lookup,only:maxvarFreq                              ! maximum number of output files
-
   implicit none
-  integer(c_int),intent(in)             :: indxHRU               ! index of hru in GRU
-  integer(c_int),intent(in)             :: indxGRU               ! index of the GRU
-  integer(c_int),intent(in)             :: timestep              ! model timestep
-  integer(c_int),intent(in)             :: outputStep            ! index into the output Struc
-  type(c_ptr),intent(in),value          :: handle_hru_data       ! local HRU data
+  ! Dummy variables
+  integer(c_int),intent(in)             :: indxGRU            ! index of the GRU
+  integer(c_int),intent(in)             :: indxHRU            ! index of hru in GRU
+  integer(c_int),intent(in)             :: timestep           ! model timestep
+  integer(c_int),intent(in)             :: outputStep         ! index into the output Struc
+  type(hru_type),intent(in),value       :: hru_data           ! local HRU data
   integer(c_int),intent(out)            :: err
-  integer(c_int),intent(out)            :: y
-  integer(c_int),intent(out)            :: m
-  integer(c_int),intent(out)            :: d
-  integer(c_int),intent(out)            :: h
-
-
-  ! local pointers
-  type(hru_type), pointer               :: hru_data              ! local HRU data
+  character(len=256),intent(out)        :: message
   ! local variables
   character(len=256)                    :: cmessage
-  character(len=256)                    :: message 
   logical(lgt)                          :: defNewOutputFile=.false.
   logical(lgt)                          :: printRestart=.false.
   logical(lgt)                          :: printProgress=.false.
@@ -125,35 +108,36 @@ subroutine hru_writeOutput(&
   character(len=256)                    :: timeString        ! portion of restart file name that contains the write-out time
   integer(i4b)                          :: iStruct           ! index of model structure
   integer(i4b)                          :: iFreq             ! index of the output frequency
-  ! convert the C pointers to Fortran pointers
-  call c_f_pointer(handle_hru_data, hru_data)
   err=0; message='summa_manageOutputFiles/'
   ! identify the start of the writing
-
-  ! updating date variables to be passed back to the actors
-  y = hru_data%timeStruct%var(iLookTIME%iyyy)
-  m = hru_data%timeStruct%var(iLookTIME%im)
-  d = hru_data%timeStruct%var(iLookTIME%id)
-  h = hru_data%timeStruct%var(iLookTIME%ih)
-
-
   ! Many variables get there values from summa4chm_util.f90:getCommandArguments()
-  call summa_setWriteAlarms(hru_data%oldTime_hru%var, hru_data%timeStruct%var, hru_data%finishTime_hru%var,  &   ! time vectors
-                            newOutputFile,  defNewOutputFile,            &
-                            ixRestart,      printRestart,                &   ! flag to print the restart file
-                            ixProgress,     printProgress,               &   ! flag to print simulation progress
-                            hru_data%resetStats%dat, hru_data%finalizeStats%dat,           &   ! flags to reset and finalize stats
-                            hru_data%statCounter%var,                             &   ! statistics counter
-                            err, cmessage)                                  ! error control
+  call summa_setWriteAlarms(hru_data%oldTime_hru%var, hru_data%timeStruct%var, & 
+                            hru_data%finishTime_hru%var, newOutputFile,        &
+                            defNewOutputFile, ixRestart, printRestart,         &
+                            ixProgress, printProgress, hru_data%resetStats%dat,& 
+                            hru_data%finalizeStats%dat,                        &
+                            hru_data%statCounter%var, err, cmessage)
   if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
 
   ! Write Parameters to output structure if this is the first
   if (timestep == 1)then
     do iStruct=1,size(structInfo)
       select case(trim(structInfo(iStruct)%structName))
-        case('attr'); call writeParm(indxGRU,indxHRU,gru_struc(indxGRU)%hruInfo(indxHRU)%hru_ix,hru_data%attrStruct,attr_meta,structInfo(iStruct)%structName,err,cmessage)
-        case('type'); call writeParm(indxGRU,indxHRU,gru_struc(indxGRU)%hruInfo(indxHRU)%hru_ix,hru_data%typeStruct,type_meta,structInfo(iStruct)%structName,err,cmessage)
-        case('mpar'); call writeParm(indxGRU,indxHRU,gru_struc(indxGRU)%hruInfo(indxHRU)%hru_ix,hru_data%mparStruct,mpar_meta,structInfo(iStruct)%structName,err,cmessage)
+        case('attr'); call writeParm(indxGRU,indxHRU,                           &
+                                     gru_struc(indxGRU)%hruInfo(indxHRU)%hru_ix,&
+                                     hru_data%attrStruct,attr_meta,             &
+                                     structInfo(iStruct)%structName,            &
+                                     err,cmessage)
+        case('type'); call writeParm(indxGRU,indxHRU,                           &
+                                     gru_struc(indxGRU)%hruInfo(indxHRU)%hru_ix,&
+                                     hru_data%typeStruct,type_meta,             &
+                                     structInfo(iStruct)%structName,            &
+                                     err,cmessage)
+        case('mpar'); call writeParm(indxGRU,indxHRU,                           &
+                                     gru_struc(indxGRU)%hruInfo(indxHRU)%hru_ix,&
+                                     hru_data%mparStruct,mpar_meta,             &
+                                     structInfo(iStruct)%structName,            &
+                                     err,cmessage)
       end select
       if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
     end do  ! (looping through structures)
@@ -171,22 +155,45 @@ subroutine hru_writeOutput(&
  ! ****************************************************************************
   do iStruct=1,size(structInfo)
     select case(trim(structInfo(iStruct)%structName))
-      case('forc'); call calcStats(hru_data%forcStat%var, hru_data%forcStruct%var, statForc_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
-      case('prog'); call calcStats(hru_data%progStat%var, hru_data%progStruct%var, statProg_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
-      case('diag'); call calcStats(hru_data%diagStat%var, hru_data%diagStruct%var, statDiag_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
-      case('flux'); call calcStats(hru_data%fluxStat%var, hru_data%fluxStruct%var, statFlux_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
-      case('indx'); call calcStats(hru_data%indxStat%var, hru_data%indxStruct%var, statIndx_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)     
+      case('forc'); call calcStats(hru_data%forcStat%var,                     &
+                                   hru_data%forcStruct%var, statForc_meta,    & 
+                                   hru_data%resetStats%dat,                   &     
+                                   hru_data%finalizeStats%dat,                &
+                                   hru_data%statCounter%var, err, cmessage)
+      case('prog'); call calcStats(hru_data%progStat%var,                     &
+                                   hru_data%progStruct%var, statProg_meta,    &
+                                   hru_data%resetStats%dat,                   &
+                                   hru_data%finalizeStats%dat,                &
+                                   hru_data%statCounter%var, err, cmessage)
+      case('diag'); call calcStats(hru_data%diagStat%var,                     &
+                                   hru_data%diagStruct%var, statDiag_meta,    &
+                                   hru_data%resetStats%dat,                   &
+                                   hru_data%finalizeStats%dat,                &
+                                   hru_data%statCounter%var, err, cmessage)
+      case('flux'); call calcStats(hru_data%fluxStat%var,                     &
+                                   hru_data%fluxStruct%var, statFlux_meta,    &
+                                   hru_data%resetStats%dat,                   &
+                                   hru_data%finalizeStats%dat,                &
+                                   hru_data%statCounter%var, err, cmessage)
+      case('indx'); call calcStats(hru_data%indxStat%var,                     &
+                                   hru_data%indxStruct%var, statIndx_meta,    &
+                                   hru_data%resetStats%dat,                   &
+                                   hru_data%finalizeStats%dat,                &
+                                   hru_data%statCounter%var, err, cmessage)     
     end select
     if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
   end do  ! (looping through structures)
     
   ! calc basin stats
-  call calcStats(hru_data%bvarStat%var(:), hru_data%bvarStruct%var(:), statBvar_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
+  call calcStats(hru_data%bvarStat%var(:), hru_data%bvarStruct%var(:),        &
+                 statBvar_meta, hru_data%resetStats%dat,                      &
+                 hru_data%finalizeStats%dat, hru_data%statCounter%var,        &
+                 err, cmessage)                   
   if(err/=0)then; message=trim(message)//trim(cmessage)//'[bvar stats]'; return; endif
-  
   ! write basin-average variables
-  call writeBasin(indxGRU,indxHRU,outputStep,hru_data%finalizeStats%dat, &
-                  hru_data%outputTimeStep%var,bvar_meta,hru_data%bvarStat%var,hru_data%bvarStruct%var,bvarChild_map,err,cmessage)
+  call writeBasin(indxGRU,indxHRU,outputStep,hru_data%finalizeStats%dat,      &
+                  hru_data%outputTimeStep%var,bvar_meta,hru_data%bvarStat%var,&
+                  hru_data%bvarStruct%var,bvarChild_map,err,cmessage)
   if(err/=0)then; message=trim(message)//trim(cmessage)//'[bvar]'; return; endif
 
   ! ****************************************************************************
@@ -201,16 +208,31 @@ subroutine hru_writeOutput(&
   ! Thus, we must also pass the stats parent->child maps from childStruct.
   do iStruct=1,size(structInfo)
     select case(trim(structInfo(iStruct)%structName))
-      case('forc'); call writeData(indxGRU,indxHRU,outputStep,"forc",hru_data%finalizeStats%dat,&
-                    maxLayers,forc_meta,hru_data%forcStat,hru_data%forcStruct,forcChild_map,hru_data%indxStruct,err,cmessage)
-      case('prog'); call writeData(indxGRU,indxHRU,outputStep,"prog",hru_data%finalizeStats%dat,&
-                    maxLayers,prog_meta,hru_data%progStat,hru_data%progStruct,progChild_map,hru_data%indxStruct,err,cmessage)
-      case('diag'); call writeData(indxGRU,indxHRU,outputStep,"diag",hru_data%finalizeStats%dat,&
-                    maxLayers,diag_meta,hru_data%diagStat,hru_data%diagStruct,diagChild_map,hru_data%indxStruct,err,cmessage)
-      case('flux'); call writeData(indxGRU,indxHRU,outputStep,"flux",hru_data%finalizeStats%dat,&
-                    maxLayers,flux_meta,hru_data%fluxStat,hru_data%fluxStruct,fluxChild_map,hru_data%indxStruct,err,cmessage)
-      case('indx'); call writeData(indxGRU,indxHRU,outputStep,"indx",hru_data%finalizeStats%dat,&
-                    maxLayers,indx_meta,hru_data%indxStat,hru_data%indxStruct,indxChild_map,hru_data%indxStruct,err,cmessage)
+      case('forc'); call writeData(indxGRU,indxHRU,outputStep,"forc",         &
+                                   hru_data%finalizeStats%dat, maxLayers,     &
+                                   forc_meta,hru_data%forcStat,               &
+                                   hru_data%forcStruct,forcChild_map,         &
+                                   hru_data%indxStruct,err,cmessage)
+      case('prog'); call writeData(indxGRU,indxHRU,outputStep,"prog",         &
+                                   hru_data%finalizeStats%dat, maxLayers,     &
+                                   prog_meta,hru_data%progStat,               &
+                                   hru_data%progStruct,progChild_map,         &
+                                   hru_data%indxStruct,err,cmessage)
+      case('diag'); call writeData(indxGRU,indxHRU,outputStep,"diag",         &
+                                   hru_data%finalizeStats%dat, maxLayers,     &
+                                   diag_meta,hru_data%diagStat,               &
+                                   hru_data%diagStruct,diagChild_map,         &
+                                   hru_data%indxStruct,err,cmessage)
+      case('flux'); call writeData(indxGRU,indxHRU,outputStep,"flux",         &
+                                   hru_data%finalizeStats%dat, maxLayers,     &
+                                   flux_meta,hru_data%fluxStat,               &
+                                   hru_data%fluxStruct,fluxChild_map,         &
+                                   hru_data%indxStruct,err,cmessage)
+      case('indx'); call writeData(indxGRU,indxHRU,outputStep,"indx",         &
+                                   hru_data%finalizeStats%dat, maxLayers,     &
+                                   indx_meta,hru_data%indxStat,               &
+                                   hru_data%indxStruct,indxChild_map,         &
+                                   hru_data%indxStruct,err,cmessage)
     end select
     if(err/=0)then 
       message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'
@@ -225,7 +247,9 @@ subroutine hru_writeOutput(&
   ! increment output file timestep
   do iFreq = 1,maxvarFreq
     hru_data%statCounter%var(iFreq) = hru_data%statCounter%var(iFreq)+1
-    if(hru_data%finalizeStats%dat(iFreq)) hru_data%outputTimeStep%var(iFreq) = hru_data%outputTimeStep%var(iFreq) + 1
+    if(hru_data%finalizeStats%dat(iFreq)) then
+      hru_data%outputTimeStep%var(iFreq) = hru_data%outputTimeStep%var(iFreq) + 1
+    end if
   end do
 
   ! if finalized stats, then reset stats on the next time step
@@ -234,7 +258,7 @@ subroutine hru_writeOutput(&
   ! save time vector
   hru_data%oldTime_hru%var(:) = hru_data%timeStruct%var(:)
 
-end subroutine hru_writeOutput
+end subroutine writeHRUOutput
 
 
 subroutine hru_writeRestart(&
diff --git a/build/source/job_actor/async_mode.cpp b/build/source/job_actor/async_mode.cpp
index 20eacb1e93accdafdb7992880af8e1da2a5ecdc8..e8edf150e2c38d7cfb474a3da186aa069487d557 100644
--- a/build/source/job_actor/async_mode.cpp
+++ b/build/source/job_actor/async_mode.cpp
@@ -10,7 +10,7 @@ behavior async_mode(stateful_actor<job_state>* self) {
     [=](file_access_actor_ready, int num_timesteps) {
       aout(self) << "Async Mode: init_file_access_actor\n";
       self->state.num_steps = num_timesteps;
-      spawnHRUActors(self);
+      spawnHRUActors(self);    
     },
 
     [=](done_hru, int gru_job_index) {
diff --git a/build/source/job_actor/gru_struc.cpp b/build/source/job_actor/gru_struc.cpp
index 78ef42eddf4e2a1ccf52627db6b7f55ffc15e335..ec5a7a76bd1c5854ecbf5f12fa2ee12445cb004c 100644
--- a/build/source/job_actor/gru_struc.cpp
+++ b/build/source/job_actor/gru_struc.cpp
@@ -10,6 +10,7 @@ GruStruc::GruStruc(int start_gru, int num_gru, int num_retry_attempts) {
 }
 
 int GruStruc::ReadDimension() {
+  // gru_struc is set up in fortran here
   int err = 0; int num_hru, file_gru, file_hru;
   std::unique_ptr<char[]> err_msg(new char[256]);
   read_dimension_fortran(start_gru_, num_gru_, num_hru, file_gru, file_hru,
@@ -33,5 +34,9 @@ int GruStruc::ReadIcondNlayers() {
   return 0;
 }
 
+void GruStruc::getNumHrusPerGru() {
+  num_hru_per_gru_.resize(num_gru_, 0);
+  get_num_hru_per_gru_fortran(num_gru_, num_hru_per_gru_[0]);
+}
 
 
diff --git a/build/source/job_actor/gru_struc.f90 b/build/source/job_actor/gru_struc.f90
index 884373935390b2af3c0106f0b20e06c32779c139..c7eb55452d9a0915bc31fc9f89a10e92466490da 100644
--- a/build/source/job_actor/gru_struc.f90
+++ b/build/source/job_actor/gru_struc.f90
@@ -5,6 +5,7 @@ module gru_struc_module
 
   public::read_dimension_fortran
   public::read_icond_nlayers_fortran
+  public::get_num_hru_per_gru_fortran
   public::deallocate_gru_struc_fortran
   contains
 
@@ -87,6 +88,22 @@ subroutine read_icond_nlayers_fortran(num_gru, err, message_r)&
 
 end subroutine read_icond_nlayers_fortran
 
+subroutine get_num_hru_per_gru_fortran(num_gru, num_hru_per_gru_array) &
+    bind(C, name="get_num_hru_per_gru_fortran")
+  USE globalData,only:gru_struc           ! gru->hru mapping structure
+  implicit none
+  ! Dummy Variables
+  integer(c_int), intent(in)      :: num_gru
+  integer(c_int), intent(out)     :: num_hru_per_gru_array(num_gru)
+  ! Local Variables
+  integer                         :: iGRU
+
+  do iGRU = 1, num_gru
+    num_hru_per_gru_array(iGRU) = gru_struc(iGRU)%hruCount
+  end do
+  
+end subroutine 
+
 subroutine deallocate_gru_struc_fortran() bind(C, name="deallocate_gru_struc_fortran")
     USE globalData,only:gru_struc           ! gru->hru mapping structure
     USE globalData,only:index_map
diff --git a/build/source/job_actor/job_actor.cpp b/build/source/job_actor/job_actor.cpp
index 14cd9843969425a4235724edfc03e82cdcb8bc2d..cf9ab3d25191f727f06a850d7c089ad555c28140 100644
--- a/build/source/job_actor/job_actor.cpp
+++ b/build/source/job_actor/job_actor.cpp
@@ -55,7 +55,7 @@ behavior job_actor(stateful_actor<job_state>* self, int start_gru, int num_gru,
     aout(self) << "ERROR: Job_Actor - ReadIcondNlayers\n";
     return {};
   }
-
+  gru_struc->getNumHrusPerGru();
 
   self->state.summa_init_struc = std::make_unique<SummaInitStruc>();
   if (self->state.summa_init_struc->allocate(self->state.num_gru) != 0) {
@@ -79,11 +79,13 @@ behavior job_actor(stateful_actor<job_state>* self, int start_gru, int num_gru,
                                         gru_struc->get_file_gru(), 
                                         false);
 
-  self->state.file_access_actor = self->spawn(
-      file_access_actor, self->state.num_gru_info, 
-      self->state.file_access_actor_settings, self);
+  self->state.file_access_actor = self->spawn(file_access_actor, 
+                                              self->state.num_gru_info, 
+                                              self->state.file_access_actor_settings, 
+                                              self);
   self->request(self->state.file_access_actor, caf::infinite, 
-                init_file_access_actor_v, gru_struc->get_file_gru())
+                init_file_access_actor_v, gru_struc->get_file_gru(),
+                gru_struc->getNumHrus())
       .await([=](int num_timesteps){
     
     if (num_timesteps < 0) {
@@ -93,7 +95,6 @@ behavior job_actor(stateful_actor<job_state>* self, int start_gru, int num_gru,
       return;
     }
 
-
     aout(self) << "Job_Actor: File Access Actor Ready\n";  
     self->state.job_timing.updateEndPoint("init_duration");
     aout(self) << "Job Actor Initialized \n";
diff --git a/build/source/job_actor/job_utils.cpp b/build/source/job_actor/job_utils.cpp
index 1be82baf9fe952b4fe9091d9eb03654bcb049b07..24fe4af6565c345a31482360808124197b8c03f3 100644
--- a/build/source/job_actor/job_utils.cpp
+++ b/build/source/job_actor/job_utils.cpp
@@ -7,12 +7,20 @@ void spawnHRUActors(stateful_actor<job_state>* self) {
   for (int i = 0; i < gru_struc->getNumGrus(); i++) {
     auto netcdf_index = gru_struc->getStartGru() + i;
     auto job_index = i + 1;
-    // Start GRU
-    auto gru = self->spawn(hru_actor, netcdf_index, job_index, 
-                           self->state.hru_actor_settings, 
-                           self->state.file_access_actor, self);
-    self->send(gru, init_hru_v);
-    self->send(gru, update_hru_async_v);
+    caf::actor gru;
+    if (gru_struc->getNumHruPerGru(i) > 1) {
+      gru = self->spawn(gru_actor, netcdf_index, job_index, 
+                        self->state.num_steps, 
+                        self->state.hru_actor_settings, 
+                        self->state.file_access_actor, self);
+    } else {
+      gru = self->spawn(hru_actor, netcdf_index, job_index, 
+                        self->state.hru_actor_settings, 
+                        self->state.file_access_actor, self);
+      self->send(gru, init_hru_v);
+      self->send(gru, update_hru_async_v);
+    }
+    
     // Save information about the GRU
     std::unique_ptr<GRU> gru_obj = std::make_unique<GRU>(
         netcdf_index, job_index, gru, self->state.dt_init_start_factor, 
diff --git a/build/source/main.cpp b/build/source/main.cpp
index d3a52e163b8007f5355646a43fe5db79db753a09..16a2092eaecbc7d707878a7bc15ac2cede6496fb 100644
--- a/build/source/main.cpp
+++ b/build/source/main.cpp
@@ -114,6 +114,8 @@ void caf_main(actor_system& sys, const config& cfg) {
     if (!std::filesystem::exists((std::filesystem::path) cfg.config_file)) {
       aout(self) << "\n\n**** Config (-c) or Master File (-m) "
                  << "Does Not Exist or Not Specified!! ****\n\n" 
+                 << "Config File: " << cfg.config_file << "\n"
+                 << "Master File: " << cfg.master_file << "\n\n"
                  << command_line_help << std::endl;
       return;
     }