diff --git a/build/CMakeLists.txt b/build/CMakeLists.txt
index 5777623f321f6a3dc83afc8e3e7b373b7407b94b..337979ed7b01c81d88349558190a0e869f1e8198 100644
--- a/build/CMakeLists.txt
+++ b/build/CMakeLists.txt
@@ -216,12 +216,13 @@ set(FILE_ACCESS_INTERFACE
 set(JOB_INTERFACE
     ${JOB_ACTOR_DIR}/gru_struc.f90)
 set(GRU_INTERFACE
-    ${GRU_ACTOR_DIR}/gru_actor.f90)
+    ${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
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 401476d0fec81675819b1614c02ef528a2c0893e..a4eb0818cb7992c39ab83c54a4274a8fdca2642b 100644
--- a/build/includes/gru_actor/gru_actor.hpp
+++ b/build/includes/gru_actor/gru_actor.hpp
@@ -1,21 +1,58 @@
 #pragma once
 #include "caf/all.hpp"
 #include "fortran_data_types.hpp"
+#include "settings_functions.hpp"
+#include "hru_actor.hpp"
 #include <vector>
 
 extern "C" {
-  void getHruInfo_fortran(int& index_gru, int& num_hru);
+  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 {
-  caf::actor parent;
   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
+
 };
 
 caf::behavior gru_actor(caf::stateful_actor<gru_actor_state>* self, 
-                        int netcdf_index, int gru_job_index, caf::actor parent);
\ No newline at end of file
+                        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_actor.hpp b/build/includes/hru_actor/hru_actor.hpp
index 0fc6322d3a9d340a3a2960f0f0ab4bbe23403916..5581ef99f8f303d9b49c20c83c334d4bbcd37033 100644
--- a/build/includes/hru_actor/hru_actor.hpp
+++ b/build/includes/hru_actor/hru_actor.hpp
@@ -23,30 +23,33 @@ 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 HRU_readForcing(int& index_gru, int& indx_hru, int& iStep, int& iRead, 
-                       int& iFile, void* hru_data,  int& err);
+  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& indx_gru, int& indx_hru, int& timestep, 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);
 
 
   
@@ -95,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/source/file_access_actor/fileAccess_writeOutput.f90 b/build/source/file_access_actor/fileAccess_writeOutput.f90
index bff869754962029a0eaaf7212453c0e42ad2e51b..72929c4b80570803f34843749280e00bc1ab2d12 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,23 @@ 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 +548,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 +559,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 +572,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 +582,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 +682,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/global/actor_data_types.f90 b/build/source/global/actor_data_types.f90
index 0a13a0871f7d9edafe5e6cfaf128eda846033c6e..f73e365a55b49d63977d5de86fefd8482acabcee 100644
--- a/build/source/global/actor_data_types.f90
+++ b/build/source/global/actor_data_types.f90
@@ -157,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 a63d1d21b79ad44dd9dad52a35149689324facfb..efd395313fb6741bee6f371fc29efcb085b211ba 100644
--- a/build/source/gru_actor/gru_actor.cpp
+++ b/build/source/gru_actor/gru_actor.cpp
@@ -4,20 +4,112 @@
 using namespace caf;
 
 behavior gru_actor(stateful_actor<gru_actor_state>* self, int netcdf_index, 
-                   int gru_job_index, caf::actor parent) {
-  self->state.parent = parent;
+                   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
-  getHruInfo_fortran(self->state.gru_job_index, self->state.num_hrus);
-
+  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);
-  for (int i = 0; i < self->state.num_hrus; i++) {
-    self->state.hrus[i] = new_handle_hru_type();
-  }
+  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++;
 
-  return {};
+        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 69cace452a43e3376eff8667e130bab811b4e6d1..0000000000000000000000000000000000000000
--- a/build/source/gru_actor/gru_actor.f90
+++ /dev/null
@@ -1,31 +0,0 @@
-module gru_actor
-USE,intrinsic :: iso_c_binding
-USE nrtype
-
-
-implicit none
-
-public :: getHruInfo_fortran
-
-
-contains
-
-subroutine getHruInfo_fortran(indx_gru, num_hru) bind(C, name="getHruInfo_fortran")
-  USE globalData,only:gru_struc
-  USE summa_init_struc,only:init_struc
-  USE var_lookup,only:iLookTYPE          ! look-up values for classification of veg, soils etc.
-  USE var_lookup,only:iLookID            ! look-up values for hru and gru IDs
-  implicit none
-  ! Dummy variables
-  integer(c_int), intent(in)  :: indx_gru
-  integer(c_int), intent(out) :: num_hru
-  ! local variables
-  integer(i4b)                :: kHRU
-  integer(i4b)                :: jHRU
-  integer(i4b)                :: iHRU
-
-  num_hru = gru_struc(indx_gru)%hruCount
-end subroutine getHruInfo_fortran
-
-
-end module gru_actor
diff --git a/build/source/gru_actor/gru_interface.f90 b/build/source/gru_actor/gru_interface.f90
new file mode 100644
index 0000000000000000000000000000000000000000..357b4120ecd2c34237585dc0e9dda06478d5fdb6
--- /dev/null
+++ b/build/source/gru_actor/gru_interface.f90
@@ -0,0 +1,374 @@
+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 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
+  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
+  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)
+
+    ! ----- 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_actor.cpp b/build/source/hru_actor/hru_actor.cpp
index 2025cb064a438c8363952b5ce77526ecb4137a23..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.indxHRU, 
-                  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.indxGRU, 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 = " 
diff --git a/build/source/hru_actor/hru_init.f90 b/build/source/hru_actor/hru_init.f90
index cdfb8e70599cfc71eaab819dec75c02d7e1d4c3c..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,21 +290,13 @@ 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
@@ -354,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
   ! *****************************************************************************
@@ -463,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..3d080579166cff7da54621939e8c2f190af56d01
--- /dev/null
+++ b/build/source/hru_actor/hru_interface.f90
@@ -0,0 +1,291 @@
+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
+  ! public::initGRU_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 86b98f4f207bcce465e8418eea79ecffa07f0da2..97d962c38148b0571ba5b43ec422b64a3c2517e0 100644
--- a/build/source/hru_actor/hru_modelRun.f90
+++ b/build/source/hru_actor/hru_modelRun.f90
@@ -76,15 +76,8 @@ public::set_sundials_tolerances
 contains
 
 ! Runs the model physics for an HRU
-subroutine runPhysics(&
-              indxGRU,             &
-              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
   ! ---------------------------------------------------------------------------------------
@@ -104,27 +97,16 @@ subroutine runPhysics(&
   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_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
-  ! ---------------------------------------------------------------------------------------
   integer(8)                                :: hruId                  ! hruId
-  real(dp)                                  :: fracHRU                ! fractional area of a given HRU (-)
   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)
@@ -134,9 +116,7 @@ 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
 
   ! ---------------------------------------------------------------------------------------
@@ -181,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
@@ -216,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
@@ -274,8 +231,8 @@ subroutine runPhysics(&
   ! run the model for a single HRU
   call coupled_em(&
                   ! model control
-                  hruId,                      & ! 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
@@ -304,7 +261,6 @@ 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)
 
 
 
@@ -331,65 +287,8 @@ subroutine runPhysics(&
   !   bvarData%var(iLookBVAR%basin__ColumnOutflow)%dat(1) = bvarData%var(iLookBVAR%basin__ColumnOutflow)%dat(1) + sum(fluxHRU%hru(iHRU)%var(iLookFLUX%mLayerColumnOutflow)%dat(:))
   ! end if
 
-
-
-
-
-
-
-
-  ! ----- 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 *****************************************
 
-  ! Get the elapsed time for the GRU
-  wallTimeTimeStep = hru_data%diagStruct%var(iLookDIAG%wallClockTime)%dat(1)
-
 end subroutine runPhysics
 
 ! *******************************************************************************************
diff --git a/build/source/hru_actor/hru_read.f90 b/build/source/hru_actor/hru_read.f90
index e3452240da18abf5245b8a68d35d08e0d706e8b0..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(indx_gru, indx_hru, 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
@@ -91,10 +90,10 @@ subroutine HRU_readForcing(indx_gru, indx_hru, iStep, iRead, iFile, handle_hru_d
   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
@@ -110,17 +109,14 @@ subroutine HRU_readForcing(indx_gru, indx_hru, iStep, iRead, iFile, handle_hru_d
   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
@@ -234,11 +230,7 @@ subroutine HRU_readForcing(indx_gru, indx_hru, iStep, iRead, iFile, handle_hru_d
                                             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 d66a465cfb0afa9ea848b3c647d3996b0f90b808..44bff51209313d337b59d554f55adcdc459b3474 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,18 +108,9 @@ 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, newOutputFile,        &
@@ -285,7 +259,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/job_utils.cpp b/build/source/job_actor/job_utils.cpp
index 16a291aaae0838d2fcc83f9b92b0d556d6c03c11..042c76ed9843b98040322da3df943d8376c7660f 100644
--- a/build/source/job_actor/job_utils.cpp
+++ b/build/source/job_actor/job_utils.cpp
@@ -8,14 +8,17 @@ void spawnHRUActors(stateful_actor<job_state>* self) {
     auto netcdf_index = gru_struc->getStartGru() + i;
     auto job_index = i + 1;
     
-    // auto gru = self->spawn(gru_actor, netcdf_index, job_index, self);
-    
-    // Start GRU
-    auto gru = self->spawn(hru_actor, netcdf_index, job_index, 
+    auto gru = self->spawn(gru_actor, netcdf_index, job_index, 
+                           self->state.num_steps, 
                            self->state.hru_actor_settings, 
                            self->state.file_access_actor, self);
-    self->send(gru, init_hru_v);
-    self->send(gru, update_hru_async_v);
+    
+    // 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);
 
 
     // Save information about the GRU