From 4a3bd4f66848bc9cf302db926aea1e73a49e1e42 Mon Sep 17 00:00:00 2001
From: KyleKlenk <kyle.c.klenk@gmail.com>
Date: Mon, 23 Oct 2023 16:37:13 -0600
Subject: [PATCH] Cleaned up HRU_Actor to encapsulate more functions in fortran
 and hide data in hru_type

---
 build/CMakeLists.txt                          |   1 +
 .../file_access_actor/partition_actor.hpp     |  35 ---
 build/includes/hru_actor/hru_actor.hpp        |  26 +-
 .../hru_actor_subroutine_wrappers.hpp         |  18 +-
 .../cpp_code/partition_actor.cpp              |  18 --
 .../fortran_code/cppwrap_fileAccess.f90       |  14 +-
 .../fortran_code/write_to_netcdf.f90          | 115 --------
 .../actors/global/cppwrap_datatypes.f90       |   8 +
 .../actors/hru_actor/cpp_code/hru_actor.cpp   |  73 ++---
 .../hru_actor/fortran_code/hru_actor.f90      | 269 +----------------
 .../hru_actor/fortran_code/hru_init.f90       |  16 +-
 .../hru_actor/fortran_code/hru_modelRun.f90   |  43 ++-
 .../hru_actor/fortran_code/hru_read.f90       | 277 ++++++++++++++++++
 .../hru_actor/fortran_code/hru_setup.f90      |   4 +-
 .../fortran_code/hru_writeOutput.f90          | 101 ++-----
 build/source/dshare/data_types.f90            |  15 +
 16 files changed, 389 insertions(+), 644 deletions(-)
 delete mode 100644 build/includes/file_access_actor/partition_actor.hpp
 delete mode 100644 build/source/actors/file_access_actor/cpp_code/partition_actor.cpp
 create mode 100644 build/source/actors/hru_actor/fortran_code/hru_read.f90

diff --git a/build/CMakeLists.txt b/build/CMakeLists.txt
index 1959119..1ca4bff 100644
--- a/build/CMakeLists.txt
+++ b/build/CMakeLists.txt
@@ -226,6 +226,7 @@ set(JOB_INTERFACE
 set(HRU_INTERFACE
     ${HRU_ACTOR_DIR}/fortran_code/hru_actor.f90
     ${HRU_ACTOR_DIR}/fortran_code/hru_init.f90
+    ${HRU_ACTOR_DIR}/fortran_code/hru_read.f90
     ${HRU_ACTOR_DIR}/fortran_code/hru_modelRun.f90
     ${HRU_ACTOR_DIR}/fortran_code/hru_modelwrite.f90
     ${HRU_ACTOR_DIR}/fortran_code/hru_restart.f90
diff --git a/build/includes/file_access_actor/partition_actor.hpp b/build/includes/file_access_actor/partition_actor.hpp
deleted file mode 100644
index e4385d8..0000000
--- a/build/includes/file_access_actor/partition_actor.hpp
+++ /dev/null
@@ -1,35 +0,0 @@
-#pragma once
-
-#include "caf/all.hpp"
-namespace caf {
-
-/*
-* This struct stores the information about a specific HRU
-*/
-struct hru_output_info {
-    caf::actor hru_actor;
-    int local_hru_index;
-    int local_gru_index;
-    bool ready_to_write;
-};
-
-
-struct partition_actor_state {
-    int start_local_gru_index; // The index of the first GRU in the partition
-    int num_local_grus; // The number of GRUs in the partition
-    int num_timesteps_simulation; // The number of timesteps in the simulation
-    int num_timesteps_buffer; // The number of timesteps this actor can hold before needing to write
-    int num_non_failed_grus = 0; // The number of GRUs that have not failed
-    bool grus_ready_to_write = false; // Whether or not all GRUs in the partition are ready to write
-
-    std::vector<std::unique_ptr<hru_output_info>> gru_list; // The GRU actors and their data
-
-
-};
-
-
-behavior partition_actor(stateful_actor<partition_actor_state>* self, int start_local_gru_index, 
-    int num_local_grus, int num_timesteps_simulation, int num_timesteps_buffer);
-
-
-}
\ 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 285a79f..e05d2f2 100644
--- a/build/includes/hru_actor/hru_actor.hpp
+++ b/build/includes/hru_actor/hru_actor.hpp
@@ -24,35 +24,24 @@ struct hru_state {
     int refGRU;			// The actual ID of the GRU we are
 
     // Variables for forcing structures
-	int stepsInCurrentFFile;
-    int num_steps_until_write;
+	int stepsInCurrentFFile;        // number of time steps in current forcing file
+    int num_steps_until_write;      // number of time steps until we pause for FA_Actor to write
 
     // HRU data structures (formerly summa_type)
     void *hru_data = new_handle_hru_type();
 
-    // Counter variables
-    void *handle_statCounter = new_handle_var_i();
-    void *handle_outputTimeStep = new_handle_var_i();
-    void *handle_resetStats = new_handle_flagVec();
-    void *handle_finalizeStats = new_handle_flagVec();
-
     // Misc Variables
 	int 		timestep = 1;	    // Current Timestep of HRU simulation
-    int         computeVegFlux;     // flag to indicate if we are computing fluxes over vegetation
     double      dt_init;            // used to initialize the length of the sub-step for each HRU
     double		upArea;             // area upslope of each HRU
     int         num_steps = 0;      // number of time steps
-    int         forcingStep;    // index of current time step in current forcing file
-    int         iFile;          // index of current forcing file from forcing file list
+    int         forcingStep;        // index of current time step in current forcing file
+    int         iFile;              // index of current forcing file from forcing file list
     int         dt_init_factor = 1; // factor of dt_init (coupled_em)
     bool        printOutput;
     int         outputFrequency;
     int         output_structure_step_index; // index of current time step in output structure
 
-    // Julian Day variables
-    double      fracJulDay;
-    double      tmZoneOffsetFracDay;
-    int         yearLength;
 
     // Settings
     HRU_Actor_Settings hru_actor_settings;
@@ -61,13 +50,6 @@ struct hru_state {
     
     ~hru_state() {
         delete_handle_hru_type(hru_data);
-
-        // counter variables
-        delete_handle_var_i(handle_statCounter);
-        delete_handle_var_i(handle_outputTimeStep);
-        delete_handle_flagVec(handle_resetStats);
-        delete_handle_flagVec(handle_finalizeStats);
-
     }
 };
 
diff --git a/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp b/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp
index 7cbe601..396fcbf 100644
--- a/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp
+++ b/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp
@@ -3,10 +3,6 @@
 extern "C" {
   // Initialize HRU data_structures
 	void initHRU(int* indxGRU, int* num_steps, void* hru_data, int* err);
-
-      // Initalize statistics flags
-      void initStatisticsFlags(void* handle_statCounter, void* handle_outputTimeStep, void* handle_resetStats,
-            void* handle_finalizeStats, int* err);
     
       void setupHRUParam( int* indxGRU, int* indxHRU, void* hru_data, double* upArea, int* err);
 
@@ -16,22 +12,16 @@ extern "C" {
 
 
       // Run the model for one timestep
-      void  RunPhysics(int* id, int* stepIndex, void* hru_data, double* fracJulDay, double* tmZoneOffsetFracDay, 
-                       int* yearLength, int* flag, double* dt, int* dt_int_factor, int* err);
+      void RunPhysics(int* id, int* stepIndex, void* hru_data, double* dt, int* dt_int_factor, int* err);
 
-      void writeHRUToOutputStructure(int* index_hru, int* index_gru, int* output_step, void* hru_data, 
-      void* handle_statCounter, void* handle_outputTimeStep, void* handle_resetStats, void* handle_finalizeStats,
-      int* err);
+      void writeHRUToOutputStructure(int* index_hru, int* index_gru, int* output_step, void* hru_data, int* err);
       
       
-      void setTimeZoneOffset(int* iFile, double* tmZoneOffsetFracDay, int* err);
+      void setTimeZoneOffset(int* iFile, void* hru_data, int* err);
 
-      void getFirstTimestep(int* iFile, int* forcing_step, int* err);
 
-      void readForcingHRU(int* index_gru, int* iStep, int* iRead, 
-            void* hru_data, int* iFile, int* err);
+      void HRU_readForcing(int* index_gru, int* iStep, int* iRead, int* iFile, void* hru_data,  int* err);
 
-      void computeTimeForcingHRU(void* hru_data, double* fracJulDay, int* yearLength, int* err);
 
       void setIDATolerances(void* hru_data,
                             double* relTolTempCas,
diff --git a/build/source/actors/file_access_actor/cpp_code/partition_actor.cpp b/build/source/actors/file_access_actor/cpp_code/partition_actor.cpp
deleted file mode 100644
index 4a226c7..0000000
--- a/build/source/actors/file_access_actor/cpp_code/partition_actor.cpp
+++ /dev/null
@@ -1,18 +0,0 @@
-#include "partition_actor.hpp"
-
-
-namespace caf {
-
-behavior partition_actor(stateful_actor<partition_actor_state>* self, int start_local_gru_index, 
-    int num_local_grus, int num_timesteps_simulation, int num_timesteps_buffer) {
-    
-    self->state.start_local_gru_index = start_local_gru_index;
-    self->state.num_local_grus = num_local_grus;
-    self->state.num_timesteps_simulation = num_timesteps_simulation;
-    self->state.num_timesteps_buffer = num_timesteps_buffer;
-    
-    
-    return {
-    };
-}
-}
\ No newline at end of file
diff --git a/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90 b/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90
index 16351f8..5357222 100644
--- a/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90
+++ b/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90
@@ -74,7 +74,7 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   USE globalData,only:checkHRU                                ! index of the HRU for a single HRU run
   
   ! look-up values for the choice of heat capacity computation
-#ifdef SUNDIALS_ACTIVE
+#ifdef V4_ACTIVE
   USE mDecisions_module,only:enthalpyFD                       ! heat capacity using enthalpy
   USE t2enthalpy_module,only:T2E_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
 #endif
@@ -296,7 +296,7 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
       if(err/=0)then; print*,message; return; endif
 
       ! calculate a lookup table to compute enthalpy from temperature, only for enthalpyFD
-#ifdef SUNDIALS_ACTIVE      
+#ifdef V4_ACTIVE      
       if(model_decisions(iLookDECISIONS%howHeatCap)%iDecision == enthalpyFD)then
         call T2E_lookup(gru_struc(iGRU)%hruInfo(iHRU)%nSoil,   &   ! intent(in):    number of soil layers
                         outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU),        &   ! intent(in):    parameter data structure
@@ -367,15 +367,7 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
                    outputStructure(1)%mparStruct,      & ! intent(in):    model parameters
                    outputStructure(1)%indxStruct_init, & ! intent(inout): model indices
                    err,message)                          ! intent(out):   error control
-  if(err/=0)then; print*, message; return; endif
-  
-
-
-
-
-
-
-  
+  if(err/=0)then; print*, message; return; endif  
 end subroutine fileAccessActor_init_fortran
 
 
diff --git a/build/source/actors/file_access_actor/fortran_code/write_to_netcdf.f90 b/build/source/actors/file_access_actor/fortran_code/write_to_netcdf.f90
index 4dc4668..18c8909 100644
--- a/build/source/actors/file_access_actor/fortran_code/write_to_netcdf.f90
+++ b/build/source/actors/file_access_actor/fortran_code/write_to_netcdf.f90
@@ -78,123 +78,8 @@ subroutine writeParamToNetCDF(handle_ncid,    &
     print*, message
     return
   endif
-
 end subroutine writeParamToNetCDF
 
-subroutine writeDataToNetCDF(handle_ncid,           &
-                             index_gru,             &
-                             index_hru,             &
-                             handle_finalize_stats, & 
-                             handle_forc_stat,      &
-                             handle_forc_struct,    &
-                             handle_prog_stat,      &
-                             handle_prog_struct,    &
-                             handle_diag_stat,      &
-                             handle_diag_struct,    &
-                             handle_flux_stat,      &
-                             handle_flux_struct,    &
-                             handle_indx_stat,      &
-                             handle_indx_struct,    &
-                             handle_output_timestep,&
-                             err) bind(C, name="writeDataToNetCDF")
-  USE globalData, only:forc_meta, prog_meta, diag_meta, flux_meta, indx_meta
-  USE globalData, only:forcChild_map, progChild_map, diagChild_map, fluxChild_map, indxChild_map
-  USE globalData,only:maxLayers                               ! maximum number of layers
-  USE globalData,only:structInfo
-  USE modelwrite_module,only:writeParm 
-  USE modelwrite_module,only:writeData
-
-  implicit none
-  ! dummy variables
-  type(c_ptr),    intent(in), value  :: handle_ncid
-  integer(c_int), intent(in)         :: index_gru
-  integer(c_int), intent(in)         :: index_hru
-  type(c_ptr),    intent(in), value  :: handle_finalize_stats
-  type(c_ptr),    intent(in), value  :: handle_forc_stat
-  type(c_ptr),    intent(in), value  :: handle_forc_struct
-  type(c_ptr),    intent(in), value  :: handle_prog_stat
-  type(c_ptr),    intent(in), value  :: handle_prog_struct
-  type(c_ptr),    intent(in), value  :: handle_diag_stat
-  type(c_ptr),    intent(in), value  :: handle_diag_struct
-  type(c_ptr),    intent(in), value  :: handle_flux_stat
-  type(c_ptr),    intent(in), value  :: handle_flux_struct
-  type(c_ptr),    intent(in), value  :: handle_indx_stat
-  type(c_ptr),    intent(in), value  :: handle_indx_struct
-  type(c_ptr),    intent(in), value  :: handle_output_timestep
-  integer(c_int), intent(out)        :: err
-  ! local pointers
-  type(var_i), pointer               :: ncid
-  type(flagVec), pointer             :: finalize_stats
-  type(var_dlength),pointer          :: forc_stat 
-  type(var_d),pointer                :: forc_struct 
-  type(var_dlength),pointer          :: prog_stat
-  type(var_dlength),pointer          :: prog_struct 
-  type(var_dlength),pointer          :: diag_stat
-  type(var_dlength),pointer          :: diag_struct 
-  type(var_dlength),pointer          :: flux_stat 
-  type(var_dlength),pointer          :: flux_struct 
-  type(var_dlength),pointer          :: indx_stat 
-  type(var_ilength),pointer          :: indx_struct 
-  type(var_i),pointer                :: output_timestep
-  ! local_variables
-  integer(i4b)                       :: numGRU=1
-  integer(i4b)                       :: nSteps=1
-  integer(i4b)                       :: iStruct
-  character(LEN=256)                 :: cmessage
-  character(LEN=256)                 :: message
-  ! ---------------------------------------------------------------------------------------
-  ! * Convert From C++ to Fortran
-  ! ---------------------------------------------------------------------------------------
-  call c_f_pointer(handle_ncid, ncid)
-  call c_f_pointer(handle_finalize_stats, finalize_stats)
-  call c_f_pointer(handle_forc_stat, forc_stat)
-  call c_f_pointer(handle_forc_struct, forc_struct)
-  call c_f_pointer(handle_prog_stat, prog_stat)
-  call c_f_pointer(handle_prog_struct, prog_struct)
-  call c_f_pointer(handle_diag_stat, diag_stat)
-  call c_f_pointer(handle_diag_struct, diag_struct)
-  call c_f_pointer(handle_flux_stat, flux_stat)
-  call c_f_pointer(handle_flux_struct, flux_struct)
-  call c_f_pointer(handle_indx_stat, indx_stat)
-  call c_f_pointer(handle_indx_struct, indx_struct)
-  call c_f_pointer(handle_output_timestep, output_timestep)
-
-  message="file_access_actor.f90 - writeDataToNetCDF"
-  do iStruct=1,size(structInfo)
-    select case(trim(structInfo(iStruct)%structName))
-      case('forc')
-        call writeData(ncid, finalize_stats%dat(:), output_timestep%var(:),maxLayers,&
-                      index_gru, numGRU, & 
-                      forc_meta,forc_stat,forc_struct,'forc', &
-                      forcChild_map,indx_struct,err,cmessage)
-      case('prog')
-        call writeData(ncid, finalize_stats%dat(:), output_timestep%var(:),maxLayers,&
-                      index_gru, numGRU, &
-                      prog_meta,prog_stat,prog_struct,'prog', &
-                      progChild_map,indx_struct,err,cmessage)
-      case('diag')
-        call writeData(ncid, finalize_stats%dat(:), output_timestep%var(:),maxLayers,&
-                      index_gru, numGRU, &
-                      diag_meta, diag_stat, diag_struct,'diag', &
-                      diagChild_map,indx_struct,err,cmessage)
-      case('flux')
-        call writeData(ncid, finalize_stats%dat(:), output_timestep%var(:),maxLayers,&
-                      index_gru, numGRU, &
-                      flux_meta,flux_stat,flux_struct,'flux', &
-                      fluxChild_map,indx_struct,err,cmessage)
-      case('indx')
-        call writeData(ncid, finalize_stats%dat(:), output_timestep%var(:),maxLayers,&
-                      index_gru, numGRU, &
-                      indx_meta,indx_stat, indx_struct, 'indx', &
-                      indxChild_map,indx_struct,err,cmessage)
-    end select
-    if(err/=0)then
-      message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'
-      print*, message
-      return
-    endif
-  end do  ! (looping through structures)
-end subroutine writeDataToNetCDF
 
 
 subroutine writeGRUStatistics(handle_ncid,      &
diff --git a/build/source/actors/global/cppwrap_datatypes.f90 b/build/source/actors/global/cppwrap_datatypes.f90
index 44770ad..38d207a 100644
--- a/build/source/actors/global/cppwrap_datatypes.f90
+++ b/build/source/actors/global/cppwrap_datatypes.f90
@@ -1296,6 +1296,10 @@ function new_handle_hru_type() result(handle) bind(C, name="new_handle_hru_type"
   allocate(p%finishTime_hru)
   allocate(p%refTime_hru)
   allocate(p%oldTime_hru)
+  allocate(p%statCounter)
+  allocate(p%outputTimeStep)
+  allocate(p%resetStats)
+  allocate(p%finalizeStats)
   handle = c_loc(p)
 end function
 
@@ -1328,6 +1332,10 @@ subroutine delete_handle_hru_type(handle) bind(C, name="delete_handle_hru_type")
   deallocate(p%finishTime_hru)
   deallocate(p%refTime_hru)
   deallocate(p%oldTime_hru)
+  deallocate(p%statCounter)
+  deallocate(p%outputTimeStep)
+  deallocate(p%resetStats)
+  deallocate(p%finalizeStats)
   deallocate(p)
 
 end subroutine
diff --git a/build/source/actors/hru_actor/cpp_code/hru_actor.cpp b/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
index e317bd3..8da70f1 100644
--- a/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
+++ b/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
@@ -30,7 +30,7 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
     self->state.hru_actor_settings = hru_actor_settings;
     self->state.dt_init_factor = hru_actor_settings.dt_init_factor;
 
-
+    // Initialize HRU data and statistics structures
     initHRU(&self->state.indxGRU, &self->state.num_steps, self->state.hru_data, &self->state.err);
     if (self->state.err != 0) {
         aout(self) << "Error: HRU_Actor - Initialize - HRU = " << self->state.indxHRU  
@@ -40,20 +40,12 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
         self->quit();
     }
 
-    // Initialize flags that are used for the output 
-    initStatisticsFlags(self->state.handle_statCounter, 
-                        self->state.handle_outputTimeStep, 
-                        self->state.handle_resetStats, 
-                        self->state.handle_finalizeStats, 
-                        &self->state.err);
-
     // Get the number of timesteps required until needing to write
     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->state.output_structure_step_index = 1;
                     Initialize_HRU(self);
                     self->send(self, start_hru_v);
                   });
@@ -94,10 +86,10 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
         [=](run_hru) {
             int err = 0;
 
-            if (self->state.timestep == 1) {
-                getFirstTimestep(&self->state.iFile, &self->state.forcingStep, &err);
-                if (self->state.forcingStep == -1) { aout(self) << "HRU - Wrong starting forcing file\n";} 
-            }
+            // if (self->state.timestep == 1) {
+            //     getFirstTimestep(&self->state.iFile, &self->state.forcingStep, &err);
+            //     if (self->state.forcingStep == -1) { aout(self) << "HRU - Wrong starting forcing file\n";} 
+            // }
 
             while(self->state.num_steps_until_write > 0) {
                 if (self->state.forcingStep > self->state.stepsInCurrentFFile) {
@@ -121,10 +113,6 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
                                           &self->state.indxGRU, 
                                           &self->state.output_structure_step_index,
                                           self->state.hru_data,
-                                          self->state.handle_statCounter,
-                                          self->state.handle_outputTimeStep,
-                                          self->state.handle_resetStats,
-                                          self->state.handle_finalizeStats,
                                           &err);
                 if (err != 0) {
                     aout(self) << "Error: HRU_Actor - writeHRUToOutputStructure - HRU = " << self->state.indxHRU
@@ -156,7 +144,7 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
             int err;
             self->state.iFile = iFile;
             self->state.stepsInCurrentFFile = num_forcing_steps_in_iFile;
-            setTimeZoneOffset(&self->state.iFile, &self->state.tmZoneOffsetFracDay, &err);
+            setTimeZoneOffset(&self->state.iFile, self->state.hru_data, &err);
             self->state.forcingStep = 1;
             self->send(self, run_hru_v);
         },
@@ -228,34 +216,20 @@ int Run_HRU(stateful_actor<hru_state>* self) {
     ** READ FORCING
     **********************************************************************/    
 
-    readForcingHRU(&self->state.indxGRU,
-                   &self->state.timestep,
-                   &self->state.forcingStep,
-                   self->state.hru_data,
-                   &self->state.iFile,
-                   &self->state.err);
+    HRU_readForcing(&self->state.indxGRU,
+                    &self->state.timestep,
+                    &self->state.forcingStep,
+                    &self->state.iFile,
+                    self->state.hru_data,
+                    &self->state.err);
     if (self->state.err != 0) {
         aout(self) << "Error---HRU_Actor: ReadForcingHRU\n" 
-                   << "     IndxGRU = " << self->state.indxGRU << "\n"
-                   << "     RefGRU = " << self->state.refGRU << "\n"
-                   << "     Forcing Step = " << self->state.forcingStep << "\n"
-                   << "     Timestep = " << self->state.timestep << "\n"
-                   << "     iFile = " << self->state.iFile << "\n"
-                   << "     Steps in Forcing File = " << self->state.stepsInCurrentFFile << "\n";
-        self->quit();
-        return -1;
-    }
-
-    computeTimeForcingHRU(self->state.hru_data,&self->state.fracJulDay,
-                          &self->state.yearLength, &self->state.err);
-    if (self->state.err != 0) {
-        aout(self) << "Error---HRU_Actor - ComputeTimeForcingHRU\n"
-                   << "     IndxGRU = " << self->state.indxGRU << "\n"
-                   << "     RefGRU = " << self->state.refGRU << "\n"
-                   << "     Forcing Step = " << self->state.forcingStep << "\n"
-                   << "     Timestep = " << self->state.timestep << "\n"
-                   << "     iFile = " << self->state.iFile << "\n"
-                   << "     Steps in Forcing File = " << self->state.stepsInCurrentFFile << "\n";
+                   << "\tIndxGRU = "               << self->state.indxGRU << "\n"
+                   << "\tRefGRU = "                << self->state.refGRU << "\n"
+                   << "\tForcing Step = "          << self->state.forcingStep << "\n"
+                   << "\tTimestep = "              << self->state.timestep << "\n"
+                   << "\tiFile = "                 << self->state.iFile << "\n"
+                   << "\tSteps in Forcing File = " << self->state.stepsInCurrentFFile << "\n";
         self->quit();
         return -1;
     }
@@ -274,19 +248,14 @@ int Run_HRU(stateful_actor<hru_state>* self) {
     RunPhysics(&self->state.indxHRU,
                &self->state.timestep,
                self->state.hru_data,
-               &self->state.fracJulDay,
-               &self->state.tmZoneOffsetFracDay,
-               &self->state.yearLength,
-               &self->state.computeVegFlux,
                &self->state.dt_init, 
                &self->state.dt_init_factor,
                &self->state.err);
-
     if (self->state.err != 0) {
         aout(self) << "Error---RunPhysics:\n"
-                   << "     IndxGRU = " << self->state.indxGRU 
-                   << "     RefGRU = " << self->state.refGRU 
-                   << "     Timestep = " << self->state.timestep <<  "\n";
+                   << "\tIndxGRU = "  << self->state.indxGRU 
+                   << "\tRefGRU = "   << self->state.refGRU 
+                   << "\tTimestep = " << self->state.timestep <<  "\n";
         self->quit();
         return 20;
     }
diff --git a/build/source/actors/hru_actor/fortran_code/hru_actor.f90 b/build/source/actors/hru_actor/fortran_code/hru_actor.f90
index 4e08e85..a022349 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_actor.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_actor.f90
@@ -11,8 +11,7 @@ USE data_types,only:&
                     hru_type
 implicit none
 
-public::getFirstTimestep
-public::setTimeZoneOffset
+
 public::setIDATolerances
 
 real(dp),parameter  :: verySmall=1e-3_rkind      ! tiny number
@@ -21,273 +20,7 @@ real(dp),parameter  :: smallOffset=1.e-8_rkind   ! small offset (units=days) to
 contains
 
 
-! Find the first timestep within the forcing file
-subroutine getFirstTimestep(iFile, iRead, err) bind(C, name="getFirstTimestep")
-  USE globalData,only:forcingDataStruct         ! forcing structure
-  USE globalData,only:vecTime                   ! time structure for forcing 
-  USE globalData,only:dJulianStart              ! julian day of start time of simulation
-  USE globalData,only:data_step                 ! length of the data step (s)
-  USE globalData,only:refJulDay_data            ! reference time for data files (fractional julian days)
-   
-  USE multiconst,only:secprday                  ! number of seconds in a day
-  
-  USE nr_utility_module,only:arth               ! get a sequence of numbers
-
-  implicit none
-
-  integer(c_int),intent(in)                               :: iFile
-  integer(c_int),intent(out)                              :: iRead
-  integer(c_int),intent(out)                              :: err
-  ! local variables
-  character(len=256)                                      :: message
-  real(dp)                                                :: timeVal(1)    ! single time value (restrict time read)
-  real(dp),dimension(forcingDataStruct(iFile)%nTimeSteps) :: fileTime      ! array of time from netcdf file
-  real(dp),dimension(forcingDataStruct(iFile)%nTimeSteps) :: diffTime      ! array of time differences
-
-  err=0; message="hru_actor.f90 - getFirstTimeStep"
-
-  ! get time vector & convert units based on offset and data step
-  timeVal(1) = vecTime(iFile)%dat(1)
-  fileTime = arth(0,1,forcingDataStruct(iFile)%nTimeSteps) * data_step/secprday + refJulDay_data &
-              + timeVal(1)/forcingDataStruct(iFile)%convTime2Days
-
-  ! find difference of fileTime from currentJulDay
-  diffTime=abs(fileTime-dJulianStart)
-  
-  if(any(diffTime < verySmall))then
-    iRead=minloc(diffTime,1)
-  else
-    iRead=-1 ! set to -1 to designinate this forcing file is not the start
-  endif
-  
-end subroutine getFirstTimestep
-
-! set the refTimeString and extract the time to set the tmZonOffsetFracDay
-subroutine setTimeZoneOffset(iFile, tmZoneOffsetFracDay, err) bind(C, name="setTimeZoneOffset")
-  USE globalData,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
-
-  USE summafilemanager,only:NC_TIME_ZONE
-  implicit none
-
-  integer(c_int),intent(in)             :: iFile
-  real(c_double),intent(out)            :: tmZoneOffsetFracDay
-  integer(c_int),intent(out)            :: err
-
-  ! local variables
-  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
-
-  err=0; message="hru_actor.f90 - setForcingTimeInfo";
-
- ! define the reference time for the model simulation
- call extractTime(forcingDataStruct(iFile)%refTimeString,& ! input  = units string for time data
-                 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
-  
- ! set the timezone offset
- select case(trim(NC_TIME_ZONE))
-    case('ncTime'); tmZoneOffsetFracDay = sign(1, ih_tz) * fracDay(ih_tz,   & ! time zone hour
-                                                                  imin_tz, & ! time zone minute
-                                                                  dsec_tz)   ! time zone second
-    case('utcTime');   tmZoneOffsetFracDay = 0._dp
-    case('localTime'); tmZoneOffsetFracDay = 0._dp
-    case default; err=20; message=trim(message)//'unable to identify time zone info option'; return
-  end select ! (option time zone option)
-
-end subroutine setTimeZoneOffset
-
-
-subroutine readForcingHRU(indxGRU, iStep, iRead, handle_hru_data, iFile, err) bind(C, name="readForcingHRU")
-  USE multiconst,only:secprday                  ! number of seconds in a day
-  ! global Data
-  USE globalData,only:data_step                 ! length of the data step (s)
-  USE globalData,only:dJulianStart              ! julian day of start time of simulation
-  USE globalData,only:refJulDay_data            ! reference time for data files (fractional julian days)
-  USE globalData,only:integerMissing            ! integer missing value
-  USE globalData,only:vecTime
-  USE globalData,only:forcingDataStruct
-  USE globalData,only:time_meta,forc_meta
-  USE var_lookup,only:iLookTIME,iLookFORCE
-  USE data_types,only:var_i,var_d
-  USE netcdf                                   ! used for nf90_max_name
-  USE time_utils_module,only:compcalday                 ! convert julian day to calendar date
-
-  implicit none
-
-  integer(c_int),intent(in)               :: indxGRU          ! Index of the GRU in gru_struc
-  integer(c_int),intent(in)               :: istep            ! Model Timestep
-  integer(c_int),intent(in)               :: iRead            ! Model Timestep 
-  type(c_ptr),intent(in),value            :: handle_hru_data  ! vector of time data for a given time step
-  integer(c_int),intent(in)               :: iFile            ! index of current forcing file from forcing file list 
-  integer(c_int),intent(out)              :: err              ! Model Timestep
-  ! 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
-  ! Counters
-  integer(i4b)                            :: iline            ! loop through lines in the file
-  integer(i4b)                            :: iVar
-  integer(i4b)                            :: iNC
-  ! other
-  logical(lgt),dimension(size(forc_meta)) :: checkForce       ! flags to check forcing data variables exist
-  
-  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
-
-  ! ---------------------------------------------------------------------------------------
-  ! * Convert From C++ to Fortran
-  ! ---------------------------------------------------------------------------------------
-  call c_f_pointer(handle_hru_data, hru_data)
-
-  err=0;message="hru_actor.f90 - readForcingHRU";
-
-  ! determine the julDay of current model step (istep) we need to read
-  currentJulDay = dJulianStart + (data_step*real(iStep-1,dp))/secprday
-
-  hru_data%timeStruct%var(:) = integerMissing
-  dataJulDay = vecTime(iFile)%dat(iRead)/forcingDataStruct(iFile)%convTime2Days + refJulDay_data
- if(abs(currentJulDay - dataJulDay) > verySmall)then
-    write(message,'(a,f18.8,a,f18.8)') trim(message)//'date for time step: ',dataJulDay,' differs from the expected date: ',currentJulDay
-    print*, message
-    err=40
-    return
-  end if
-
-  ! convert julian day to time vector
-  ! NOTE: use small offset to force ih=0 at the start of the day
-  call compcalday(dataJulDay+smallOffset,         & ! input  = julian day
-                  hru_data%timeStruct%var(iLookTIME%iyyy),      & ! output = year
-                  hru_data%timeStruct%var(iLookTIME%im),        & ! output = month
-                  hru_data%timeStruct%var(iLookTIME%id),        & ! output = day
-                  hru_data%timeStruct%var(iLookTIME%ih),        & ! output = hour
-                  hru_data%timeStruct%var(iLookTIME%imin),dsec, & ! output = minute/second
-                  err,cmessage)                     ! output = error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
-  
-  ! check to see if any of the time data is missing -- note that it is OK if ih_tz or imin_tz are missing
-  if((hru_data%timeStruct%var(iLookTIME%iyyy)==integerMissing) .or. (hru_data%timeStruct%var(iLookTIME%im)==integerMissing) .or. (hru_data%timeStruct%var(iLookTIME%id)==integerMissing) .or. (hru_data%timeStruct%var(iLookTIME%ih)==integerMissing) .or. (hru_data%timeStruct%var(iLookTIME%imin)==integerMissing))then
-      do iline=1,size(hru_data%timeStruct%var)
-          if(hru_data%timeStruct%var(iline)==integerMissing)then; err=40; message=trim(message)//"variableMissing[var='"//trim(time_meta(iline)%varname)//"']"; return; end if
-      end do
-  end if
-
-  ! initialize flags for forcing data
-  checkForce(:) = .false.
-  checkForce(iLookFORCE%time) = .true.  ! time is handled separately
-
-  do iNC=1,forcingDataStruct(iFile)%nVars
-    ! check variable is desired
-    if(forcingDataStruct(iFile)%var_ix(iNC)==integerMissing) cycle
-
-    ! get index in forcing structure
-    iVar = forcingDataStruct(iFile)%var_ix(iNC)
-    checkForce(iVar) = .true.
-
-    ! check individual data value
-    if(forcingDataStruct(iFile)%var(ivar)%dataFromFile(indxGRU,iRead)<dataMin)then
-      write(message,'(a,f13.5)') trim(message)//'forcing data for variable '//trim(varname)//' is less than minimum allowable value ', dataMin
-      err=20; return
-    endif
-    ! put the data into structures
-    hru_data%forcStruct%var(ivar) = forcingDataStruct(iFile)%var(ivar)%dataFromFile(indxGRU,iRead)
-  end do  ! loop through forcing variables
-  
-  ! check if any forcing data is missing
-  if(count(checkForce)<size(forc_meta))then
-    do iline=1,size(forc_meta)
-    if(.not.checkForce(iline))then
-      message=trim(message)//"checkForce_variableMissing[var='"//trim(forc_meta(iline)%varname)//"']"
-      err=20; return
-    endif    ! if variable is missing
-    end do   ! looping through variables
-  end if   ! if any variables are missing
-end subroutine readForcingHRU 
-
-! This is part 2: from the reading of forcing - separated it so we could call from C++
-subroutine computeTimeForcingHRU(handle_hru_data, fracJulDay, yearLength, err) bind(C, name="computeTimeForcingHRU")
-  USE var_lookup,only:iLookTIME,iLookFORCE
-  USE time_utils_module,only:compJulDay         ! convert calendar date to julian day
-  USE data_types,only:var_i,var_d
-  USE multiconst,only:secprday                  ! number of seconds in a day
-  USE globalData,only:refJulDay                 ! reference time (fractional julian days)
-
-  implicit none
-  type(c_ptr),intent(in),value        :: handle_hru_data      ! vector of time data for a given time step
-  real(c_double),intent(out)          :: fracJulDay
-  integer(c_int),intent(out)          :: yearLength    
-  integer(c_int),intent(out)          :: err   
-  ! local variables
-  type(hru_type),pointer              :: hru_data         !  model time data
-  real(dp)                            :: startJulDay      ! julian day at the start of the year
-  real(dp)                            :: currentJulDay    ! Julian day of current time step
-  character(len=256)                  :: message          ! error message for downwind routine
-  character(len=256)                  :: cmessage         ! error message for downwind routine
-  logical(lgt),parameter              :: checkTime=.false.  ! flag to check the time
-
-  ! ---------------------------------------------------------------------------------------
-  ! * Convert From C++ to Fortran
-  ! ---------------------------------------------------------------------------------------
-  call c_f_pointer(handle_hru_data, hru_data)
-
-  err=0; message="hru_actor.f90 - computTimeForcingHRU/"
-
-  ! compute the julian day at the start of the year
-  call compjulday(hru_data%timeStruct%var(iLookTIME%iyyy),          & ! input  = year
-                  1, 1, 1, 1, 0._dp,                  & ! input  = month, day, hour, minute, second
-                  startJulDay,err,cmessage)             ! output = julian day (fraction of day) + error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
-
-  ! compute the fractional julian day for the current time step
-  call compjulday(hru_data%timeStruct%var(iLookTIME%iyyy),           & ! input  = year
-                  hru_data%timeStruct%var(iLookTIME%im),             & ! input  = month
-                  hru_data%timeStruct%var(iLookTIME%id),             & ! input  = day
-                  hru_data%timeStruct%var(iLookTIME%ih),             & ! input  = hour
-                  hru_data%timeStruct%var(iLookTIME%imin),0._dp,     & ! input  = minute/second
-                  currentJulDay,err,cmessage)            ! output = julian day (fraction of day) + error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
-  ! compute the time since the start of the year (in fractional days)
-  fracJulDay = currentJulDay - startJulDay
-  ! set timing of current forcing vector (in seconds since reference day)
-  ! NOTE: It is a bit silly to have time information for each HRU and GRU
-  hru_data%forcStruct%var(iLookFORCE%time) = (currentJulDay-refJulDay)*secprday
-
-  ! compute the number of days in the current year
-  yearLength = 365
-  if(mod(hru_data%timeStruct%var(iLookTIME%iyyy),4) == 0)then
-    yearLength = 366
-    if(mod(hru_data%timeStruct%var(iLookTIME%iyyy),100) == 0)then
-    yearLength = 365
-    if(mod(hru_data%timeStruct%var(iLookTIME%iyyy),400) == 0)then
-      yearLength = 366
-    end if
-    end if
-  end if
 
-  ! test
-  if(checkTime)then
-    write(*,'(i4,1x,4(i2,1x),f9.3,1x,i4)')  hru_data%timeStruct%var(iLookTIME%iyyy),           & ! year
-                                            hru_data%timeStruct%var(iLookTIME%im),             & ! month
-                                            hru_data%timeStruct%var(iLookTIME%id),             & ! day
-                                            hru_data%timeStruct%var(iLookTIME%ih),             & ! hour
-                                            hru_data%timeStruct%var(iLookTIME%imin),           & ! minute
-                                            fracJulDay,                          & ! fractional julian day for the current time step
-                                            yearLength                             ! number of days in the current year
-    !pause ' checking time'
-  end if
-end subroutine computeTimeForcingHRU
 
 
 ! Set the HRU's relative and absolute tolerances
diff --git a/build/source/actors/hru_actor/fortran_code/hru_init.f90 b/build/source/actors/hru_actor/fortran_code/hru_init.f90
index 7b7ed50..6718e02 100755
--- a/build/source/actors/hru_actor/fortran_code/hru_init.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_init.f90
@@ -65,7 +65,7 @@ contains
   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
   
   ! ---------------------------------------------------------------------------------------
@@ -190,6 +190,20 @@ contains
     endif
   end do ! iStruct
 
+
+  ! Intilaize the statistics data structures
+  allocate(hru_data%statCounter%var(maxVarFreq), stat=err)
+  allocate(hru_data%outputTimeStep%var(maxVarFreq), stat=err)
+  allocate(hru_data%resetStats%dat(maxVarFreq), stat=err)
+  allocate(hru_data%finalizeStats%dat(maxVarFreq), stat=err)
+  hru_data%statCounter%var(1:maxVarFreq) = 1
+  hru_data%outputTimeStep%var(1:maxVarFreq) = 1
+  ! initialize flags to reset/finalize statistics
+  hru_data%resetStats%dat(:)    = .true.   ! start by resetting statistics
+  hru_data%finalizeStats%dat(:) = .false.  ! do not finalize stats on the first time step
+  ! set stats flag for the timestep-level output
+  hru_data%finalizeStats%dat(iLookFreq%timestep)=.true.
+
   ! identify the end of the initialization
   call date_and_time(values=endInit)
 
diff --git a/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90 b/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90
index 76b02b2..37254cf 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90
@@ -84,11 +84,6 @@ subroutine runPhysics(&
               indxHRU,             &
               modelTimeStep,       &
               handle_hru_data,     &
-              fracJulDay,          & ! fraction of the current Julian day
-              tmZoneOffsetFracDay, & ! time zone offset (fraction of the day)
-              yearLength,          & ! number of days in the current year
-              ! run time variables
-              computeVegFlux,      & ! flag to indicate if we are computing fluxes over vegetation
               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
               err) bind(C, name='RunPhysics')
@@ -114,20 +109,16 @@ subroutine runPhysics(&
   ! ---------------------------------------------------------------------------------------
   ! Dummy Variables
   ! ---------------------------------------------------------------------------------------
-  integer(c_long),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)             :: fracJulDay                    ! fractional julian days since the start of year
-  real(c_double),intent(inout)             :: tmZoneOffsetFracDay 
-  integer(c_int),intent(inout)             :: yearLength                    ! number of days in the current year
-  integer(c_int),intent(inout)             :: computeVegFlux         ! flag to indicate if we are computing fluxes over vegetation
-  real(c_double),intent(inout)             :: dt_init                ! used to initialize the length of the sub-step for each HRU
-  integer(c_int),intent(in)                :: dt_init_factor         ! used to adjust the length of the timestep in the event of a failure
-  integer(c_int),intent(inout)             :: err                    ! error code
+  integer(c_long),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
+  integer(c_int), intent(in)                :: dt_init_factor         ! used to adjust the length of the timestep in the event of a failure
+  integer(c_int), intent(inout)             :: err                    ! error code
   ! ---------------------------------------------------------------------------------------
   ! FORTRAN POINTERS
   ! ---------------------------------------------------------------------------------------
-  type(hru_type),pointer                   :: hru_data                ! hru data
+  type(hru_type),pointer                    :: hru_data               ! hru data
 
   ! ---------------------------------------------------------------------------------------
   ! local variables: general
@@ -160,8 +151,8 @@ subroutine runPhysics(&
       ! (compute the exposed LAI and SAI and whether veg is buried by snow)
       call vegPhenlgy(&
                       ! model control
-                      fracJulDay,             & ! intent(in):    fractional julian days since the start of year
-                      yearLength,             & ! intent(in):    number of days in the current year
+                      hru_data%fracJulDay,             & ! intent(in):    fractional julian days since the start of year
+                      hru_data%yearLength,             & ! intent(in):    number of days in the current year
                       ! input/output: data structures
                       model_decisions,        & ! intent(in):    model decisions
                       hru_data%typeStruct,    & ! intent(in):    type of vegetation and soil
@@ -178,8 +169,8 @@ subroutine runPhysics(&
 
     
       ! save the flag for computing the vegetation fluxes
-      if(computeVegFluxFlag)      computeVegFlux = yes
-      if(.not.computeVegFluxFlag) computeVegFlux = no
+      if(computeVegFluxFlag)      hru_data%computeVegFlux = yes
+      if(.not.computeVegFluxFlag) hru_data%computeVegFlux = no
       
       ! define the green vegetation fraction of the grid box (used to compute LAI)
       hru_data%diagStruct%var(iLookDIAG%scalarGreenVegFraction)%dat(1) = greenVegFrac_monthly(hru_data%timeStruct%var(iLookTIME%im))
@@ -216,7 +207,7 @@ subroutine runPhysics(&
   nSoil   = hru_data%indxStruct%var(iLookINDEX%nSoil)%dat(1)    ! number of soil layers
   nLayers = hru_data%indxStruct%var(iLookINDEX%nLayers)%dat(1)  ! total number of layers
   
-  computeVegFluxFlag = (ComputeVegFlux == yes)
+  computeVegFluxFlag = (hru_data%ComputeVegFlux == yes)
 
   !******************************************************************************
   !****************************** From run_oneHRU *******************************
@@ -267,7 +258,7 @@ subroutine runPhysics(&
         hru_data%progStruct,         & ! data structure of model prognostic variables
         hru_data%diagStruct,         & ! data structure of model diagnostic variables
         hru_data%fluxStruct,         & ! data structure of model fluxes
-        tmZoneOffsetFracDay,         & ! time zone offset in fractional days
+        hru_data%tmZoneOffsetFracDay,         & ! time zone offset in fractional days
         err,cmessage)                  ! error control
   if(err/=0)then;err=20; message=trim(message)//cmessage; print*, message; return; endif
  
@@ -281,8 +272,8 @@ subroutine runPhysics(&
                   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
-                  fracJulDay,                  & ! intent(in):    fractional julian days since the start of year
-                  yearLength,                  & ! intent(in):    number of days in the current year
+                  hru_data%fracJulDay,                  & ! intent(in):    fractional julian days since the start of year
+                  hru_data%yearLength,                  & ! intent(in):    number of days in the current year
                   ! data structures (input)
                   hru_data%typeStruct,         & ! intent(in):    local classification of soil veg etc. for each HRU
                   hru_data%attrStruct,         & ! intent(in):    local attributes for each HRU
@@ -304,8 +295,8 @@ subroutine runPhysics(&
 
   !************************************* End of run_oneHRU *****************************************
   ! save the flag for computing the vegetation fluxes
-  if(computeVegFluxFlag)      ComputeVegFlux = yes
-  if(.not.computeVegFluxFlag) ComputeVegFlux = no
+  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)
 
diff --git a/build/source/actors/hru_actor/fortran_code/hru_read.f90 b/build/source/actors/hru_actor/fortran_code/hru_read.f90
new file mode 100644
index 0000000..8d75099
--- /dev/null
+++ b/build/source/actors/hru_actor/fortran_code/hru_read.f90
@@ -0,0 +1,277 @@
+module hru_read
+
+
+USE,intrinsic :: iso_c_binding
+USE nrtype
+USE data_types,only:&
+                    var_i,          &  
+                    var_i8,         &
+                    var_d,          &
+                    var_ilength,    &
+                    var_dlength,    &
+                    flagVec,        &
+                    hru_type
+implicit none
+public::setTimeZoneOffset
+public::HRU_readForcing
+private::getFirstTimeStep
+
+
+real(dp),parameter  :: verySmall=1e-3_rkind      ! tiny number
+real(dp),parameter  :: smallOffset=1.e-8_rkind   ! small offset (units=days) to force ih=0 at the start of the day
+
+contains
+
+! set the refTimeString and extract the time to set the tmZonOffsetFracDay
+subroutine setTimeZoneOffset(iFile, handle_hru_data, err) bind(C, name="setTimeZoneOffset")
+  USE globalData,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
+  USE summafilemanager,only:NC_TIME_ZONE
+  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
+  integer(c_int),intent(out)            :: err
+
+  ! 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
+  call extractTime(forcingDataStruct(iFile)%refTimeString, & ! input  = units string for time data
+                   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
+  
+  ! set the timezone offset
+  select case(trim(NC_TIME_ZONE))
+    case('ncTime');     hru_data%tmZoneOffsetFracDay = sign(1, ih_tz) * fracDay(ih_tz,   & ! time zone hour
+                                                                                imin_tz, & ! time zone minute
+                                                                                dsec_tz)   ! time zone second
+    case('utcTime');   hru_data%tmZoneOffsetFracDay = 0._dp
+    case('localTime'); hru_data%tmZoneOffsetFracDay = 0._dp
+    case default; err=20; message=trim(message)//'unable to identify time zone info option'; return
+  end select ! (option time zone option)
+
+end subroutine setTimeZoneOffset
+
+subroutine HRU_readForcing(indxGRU, iStep, iRead, iFile, handle_hru_data, err) bind(C, name="HRU_readForcing")
+  USE multiconst,only:secprday                  ! number of seconds in a day
+  USE time_utils_module,only:compJulDay         ! convert calendar date to julian day
+  ! global Data
+  USE globalData,only:data_step                 ! length of the data step (s)
+  USE globalData,only:dJulianStart              ! julian day of start time of simulation
+  USE globalData,only:refJulDay_data            ! reference time for data files (fractional julian days)
+  USE globalData,only:integerMissing            ! integer missing value
+  USE globalData,only:vecTime
+  USE globalData,only:forcingDataStruct
+  USE globalData,only:time_meta,forc_meta
+  USE var_lookup,only:iLookTIME,iLookFORCE
+  USE data_types,only:var_i,var_d
+  USE netcdf,only:nf90_max_name                                   ! used for nf90_max_name
+  USE time_utils_module,only:compcalday                 ! convert julian day to calendar date
+  USE globalData,only:refJulDay                 ! reference time (fractional julian days)
+
+  implicit none
+
+  integer(c_int),intent(in)               :: indxGRU          ! Index of the GRU in gru_struc
+  integer(c_int),intent(in)               :: istep            ! Model Timestep
+  integer(c_int),intent(inout)            :: iRead            ! Model Timestep 
+  integer(c_int),intent(in)               :: iFile            ! index of current forcing file from forcing file list 
+  type(c_ptr),intent(in),value            :: handle_hru_data  ! vector of time data for a given time step
+  integer(c_int),intent(out)              :: err              ! Model Timestep
+  ! 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
+  
+  ! Counters
+  integer(i4b)                            :: iline            ! loop through lines in the file
+  integer(i4b)                            :: iVar
+  integer(i4b)                            :: iNC
+  ! other
+  logical(lgt),dimension(size(forc_meta)) :: checkForce       ! flags to check forcing data variables exist
+  logical(lgt),parameter                  :: checkTime=.false.! flag to check the time
+  
+  real(dp)                                :: dsec             ! double precision seconds (not used)
+  real(dp),parameter                      :: dataMin=-1._dp   ! minimum allowable data value (all forcing variables should be positive)
+  character(len = nf90_max_name)          :: varName          ! dimenison name
+
+  character(len=256)                      :: message          ! error message
+  character(len=256)                      :: cmessage         ! error message
+
+  call c_f_pointer(handle_hru_data, hru_data)
+  err=0;message="hru_actor.f90 - readForcingHRU";
+
+  if(istep == 1) then
+    call getFirstTimestep(iFile, iRead, err)
+    if(err/=0)then; message=trim(message)//"getFirstTimestep"; print*,message;return; end if
+  endif
+  
+  ! determine the julDay of current model step (istep) we need to read
+  currentJulDay = dJulianStart + (data_step*real(iStep-1,dp))/secprday
+
+  hru_data%timeStruct%var(:) = integerMissing
+  dataJulDay = vecTime(iFile)%dat(iRead)/forcingDataStruct(iFile)%convTime2Days + refJulDay_data
+  if(abs(currentJulDay - dataJulDay) > verySmall)then
+    write(message,'(a,f18.8,a,f18.8)') trim(message)//'date for time step: ',dataJulDay,' differs from the expected date: ',currentJulDay
+    print*, message
+    err=40
+    return
+  end if
+  
+  ! convert julian day to time vector
+  ! NOTE: use small offset to force ih=0 at the start of the day
+  call compcalday(dataJulDay+smallOffset,         & ! input  = julian day
+                  hru_data%timeStruct%var(iLookTIME%iyyy),      & ! output = year
+                  hru_data%timeStruct%var(iLookTIME%im),        & ! output = month
+                  hru_data%timeStruct%var(iLookTIME%id),        & ! output = day
+                  hru_data%timeStruct%var(iLookTIME%ih),        & ! output = hour
+                  hru_data%timeStruct%var(iLookTIME%imin),dsec, & ! output = minute/second
+                  err,cmessage)                     ! output = error control
+  if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
+    
+  ! check to see if any of the time data is missing -- note that it is OK if ih_tz or imin_tz are missing
+  if((hru_data%timeStruct%var(iLookTIME%iyyy)==integerMissing) .or. (hru_data%timeStruct%var(iLookTIME%im)==integerMissing) .or. (hru_data%timeStruct%var(iLookTIME%id)==integerMissing) .or. (hru_data%timeStruct%var(iLookTIME%ih)==integerMissing) .or. (hru_data%timeStruct%var(iLookTIME%imin)==integerMissing))then
+      do iline=1,size(hru_data%timeStruct%var)
+          if(hru_data%timeStruct%var(iline)==integerMissing)then; err=40; message=trim(message)//"variableMissing[var='"//trim(time_meta(iline)%varname)//"']"; return; end if
+      end do
+  end if
+
+  ! initialize flags for forcing data
+  checkForce(:) = .false.
+  checkForce(iLookFORCE%time) = .true.  ! time is handled separately
+
+  do iNC=1,forcingDataStruct(iFile)%nVars
+    ! check variable is desired
+    if(forcingDataStruct(iFile)%var_ix(iNC)==integerMissing) cycle
+
+    ! get index in forcing structure
+    iVar = forcingDataStruct(iFile)%var_ix(iNC)
+    checkForce(iVar) = .true.
+
+    ! check individual data value
+    if(forcingDataStruct(iFile)%var(ivar)%dataFromFile(indxGRU,iRead)<dataMin)then
+      write(message,'(a,f13.5)') trim(message)//'forcing data for variable '//trim(varname)//' is less than minimum allowable value ', dataMin
+      err=20; return
+    endif
+    ! put the data into structures
+    hru_data%forcStruct%var(ivar) = forcingDataStruct(iFile)%var(ivar)%dataFromFile(indxGRU,iRead)
+  end do  ! loop through forcing variables
+    
+  ! check if any forcing data is missing
+  if(count(checkForce)<size(forc_meta))then
+    do iline=1,size(forc_meta)
+    if(.not.checkForce(iline))then
+      message=trim(message)//"checkForce_variableMissing[var='"//trim(forc_meta(iline)%varname)//"']"
+      err=20; return
+    endif    ! if variable is missing
+    end do   ! looping through variables
+  end if   ! if any variables are missing
+
+
+ ! **********************************************************************************************
+ ! ***** part 2: compute time
+ ! **********************************************************************************************
+  ! compute the julian day at the start of the year
+  call compjulday(hru_data%timeStruct%var(iLookTIME%iyyy),          & ! input  = year
+                  1, 1, 1, 1, 0._dp,                  & ! input  = month, day, hour, minute, second
+                  startJulDay,err,cmessage)             ! output = julian day (fraction of day) + error control
+  if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
+
+  ! compute the fractional julian day for the current time step
+  call compjulday(hru_data%timeStruct%var(iLookTIME%iyyy),           & ! input  = year
+                  hru_data%timeStruct%var(iLookTIME%im),             & ! input  = month
+                  hru_data%timeStruct%var(iLookTIME%id),             & ! input  = day
+                  hru_data%timeStruct%var(iLookTIME%ih),             & ! input  = hour
+                  hru_data%timeStruct%var(iLookTIME%imin),0._dp,     & ! input  = minute/second
+                  currentJulDay,err,cmessage)            ! output = julian day (fraction of day) + error control
+  if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
+  ! compute the time since the start of the year (in fractional days)
+  hru_data%fracJulDay = currentJulDay - startJulDay
+  ! set timing of current forcing vector (in seconds since reference day)
+  ! NOTE: It is a bit silly to have time information for each HRU and GRU
+  hru_data%forcStruct%var(iLookFORCE%time) = (currentJulDay-refJulDay)*secprday
+
+  ! compute the number of days in the current year
+  hru_data%yearLength = 365
+  if(mod(hru_data%timeStruct%var(iLookTIME%iyyy),4) == 0)then
+    hru_data%yearLength = 366
+    if(mod(hru_data%timeStruct%var(iLookTIME%iyyy),100) == 0)then
+    hru_data%yearLength = 365
+    if(mod(hru_data%timeStruct%var(iLookTIME%iyyy),400) == 0)then
+      hru_data%yearLength = 366
+    end if
+    end if
+  end if
+
+  ! test
+  if(checkTime)then
+    write(*,'(i4,1x,4(i2,1x),f9.3,1x,i4)')  hru_data%timeStruct%var(iLookTIME%iyyy),           & ! year
+                                            hru_data%timeStruct%var(iLookTIME%im),             & ! month
+                                            hru_data%timeStruct%var(iLookTIME%id),             & ! day
+                                            hru_data%timeStruct%var(iLookTIME%ih),             & ! hour
+                                            hru_data%timeStruct%var(iLookTIME%imin),           & ! minute
+                                            hru_data%fracJulDay,                          & ! fractional julian day for the current time step
+                                            hru_data%yearLength                             ! number of days in the current year
+    !pause ' checking time'
+  end if
+
+
+
+
+end subroutine HRU_readForcing 
+
+ ! Find the first timestep within the forcing file
+subroutine getFirstTimestep(iFile, iRead, err)
+  USE globalData,only:forcingDataStruct         ! forcing structure
+  USE globalData,only:vecTime                   ! time structure for forcing 
+  USE globalData,only:dJulianStart              ! julian day of start time of simulation
+  USE globalData,only:data_step                 ! length of the data step (s)
+  USE globalData,only:refJulDay_data            ! reference time for data files (fractional julian days)
+    
+  USE multiconst,only:secprday                  ! number of seconds in a day
+  
+  USE nr_utility_module,only:arth               ! get a sequence of numbers
+
+  implicit none
+
+  integer(i4b),intent(in)                                 :: iFile
+  integer(i4b),intent(out)                                :: iRead
+  integer(i4b),intent(out)                                :: err
+  ! local variables
+  character(len=256)                                      :: message
+  real(dp)                                                :: timeVal(1)    ! single time value (restrict time read)
+  real(dp),dimension(forcingDataStruct(iFile)%nTimeSteps) :: fileTime      ! array of time from netcdf file
+  real(dp),dimension(forcingDataStruct(iFile)%nTimeSteps) :: diffTime      ! array of time differences
+
+  err=0; message="hru_actor.f90 - getFirstTimeStep"
+
+  ! get time vector & convert units based on offset and data step
+  timeVal(1) = vecTime(iFile)%dat(1)
+  fileTime = arth(0,1,forcingDataStruct(iFile)%nTimeSteps) * data_step/secprday + refJulDay_data &
+              + timeVal(1)/forcingDataStruct(iFile)%convTime2Days
+
+  ! find difference of fileTime from currentJulDay
+  diffTime=abs(fileTime-dJulianStart)
+  
+  if(any(diffTime < verySmall))then
+    iRead=minloc(diffTime,1)
+  else
+    iRead=-1 ! set to -1 to designinate this forcing file is not the start
+  endif
+  
+end subroutine getFirstTimestep
+
+end module hru_read
+
diff --git a/build/source/actors/hru_actor/fortran_code/hru_setup.f90 b/build/source/actors/hru_actor/fortran_code/hru_setup.f90
index e014839..50256ba 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_setup.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_setup.f90
@@ -83,7 +83,7 @@ subroutine setupHRUParam(&
   USE paramCheck_module,only:paramCheck                       ! module to check consistency of model parameters
   USE pOverwrite_module,only:pOverwrite                       ! module to overwrite default parameter values with info from the Noah tables
   USE ConvE2Temp_module,only:E2T_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
-#ifdef SUNDIALS_ACTIVE  
+#ifdef V4_ACTIVE  
   USE t2enthalpy_module,only:T2E_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
 #endif
   USE var_derive_module,only:fracFuture                       ! module to calculate the fraction of runoff in future time steps (time delay histogram)
@@ -112,7 +112,7 @@ subroutine setupHRUParam(&
   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
+  real(c_double),intent(inout)             :: upArea
   integer(c_int),intent(inout)             :: err
 
   ! local variables
diff --git a/build/source/actors/hru_actor/fortran_code/hru_writeOutput.f90 b/build/source/actors/hru_actor/fortran_code/hru_writeOutput.f90
index 0000875..abf05d8 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_writeOutput.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_writeOutput.f90
@@ -55,61 +55,14 @@ USE get_ixname_module,only:get_freqName       ! get name of frequency from frequ
 
 implicit none
 private
-public::initStatisticsFlags
 public::writeHRUToOutputStructure
 
 contains
-subroutine initStatisticsFlags(handle_statCounter, handle_outputTimeStep, &
-    handle_resetStats, handle_finalizeStats, err) bind(C, name="initStatisticsFlags")
-  USE var_lookup,only:maxvarFreq                ! maximum number of output files
-  USE var_lookup,only:iLookFreq                 ! named variables for the frequency structure
-  implicit none
-  ! dummy variables
-  type(c_ptr), intent(in), value       :: handle_statCounter
-  type(c_ptr), intent(in), value       :: handle_outputTimeStep
-  type(c_ptr), intent(in), value       :: handle_resetStats
-  type(c_ptr), intent(in), value       :: handle_finalizeStats
-  integer(c_int), intent(out)          :: err
-  ! local variables
-  type(var_i), pointer                 :: statCounter
-  type(var_i), pointer                 :: outputTimeStep
-  type(flagVec), pointer               :: resetStats
-  type(flagVec), pointer               :: finalizeStats
-  ! Convert C pointers to Fortran pointers
-  call c_f_pointer(handle_statCounter, statCounter)
-  call c_f_pointer(handle_outputTimeStep, outputTimeStep)
-  call c_f_pointer(handle_resetStats, resetStats)
-  call c_f_pointer(handle_finalizeStats, finalizeStats)
-  ! Start of Subroutine
-
-  ! initialize the statistics flags
-  allocate(statCounter%var(maxVarFreq), stat=err)
-  allocate(outputTimeStep%var(maxVarFreq), stat=err)
-  statCounter%var(1:maxVarFreq) = 1
-  outputTimeStep%var(1:maxVarFreq) = 1
-
-  allocate(resetStats%dat(maxVarFreq), stat=err)
-  allocate(finalizeStats%dat(maxVarFreq), stat=err)
-  ! initialize flags to reset/finalize statistics
-  resetStats%dat(:)    = .true.   ! start by resetting statistics
-  finalizeStats%dat(:) = .false.  ! do not finalize stats on the first time step
-
-  ! set stats flag for the timestep-level output
-  finalizeStats%dat(iLookFreq%timestep)=.true.
-
-
-end subroutine initStatisticsFlags
-
 subroutine writeHRUToOutputStructure(&
                             indxHRU,                   &
                             indxGRU,                   &
                             outputStep,                & ! index into the output Struc
-                            handle_hru_data,           & ! local HRU data
-                            handle_statCounter,        & ! x%var(:)
-                            handle_outputTimeStep,     & ! x%var(:)
-                            handle_resetStats,         & ! x%var(:)
-                            handle_finalizeStats,      & ! x%var(:)
-                            ! run time variables
+                            handle_hru_data,           & ! local HRU data  
                             err) bind(C, name="writeHRUToOutputStructure") 
   USE nrtype
   USE globalData,only:structInfo
@@ -142,18 +95,10 @@ subroutine writeHRUToOutputStructure(&
   integer(c_int),intent(in)             :: indxGRU               ! index of the GRU
   integer(c_int),intent(in)             :: outputStep            ! index into the output Struc
   type(c_ptr),intent(in),value          :: handle_hru_data       ! local HRU data
-  type(c_ptr),intent(in),value          :: handle_statCounter    ! x%var(:)
-  type(c_ptr),intent(in),value          :: handle_outputTimeStep ! x%var(:)
-  type(c_ptr),intent(in),value          :: handle_resetStats     ! x%var(:)
-  type(c_ptr),intent(in),value          :: handle_finalizeStats  ! x%var(:)
   integer(c_int),intent(out)            :: err
 
   ! local pointers
   type(hru_type), pointer               :: hru_data              ! local HRU data
-  type(var_i),pointer                   :: statCounter     ! time counter for stats
-  type(var_i),pointer                   :: outputTimeStep  ! timestep in output files
-  type(flagVec),pointer                 :: resetStats      ! flags to reset statistics
-  type(flagVec),pointer                 :: finalizeStats   ! flags to finalize statistics
   ! local variables
   character(len=256)                    :: cmessage
   character(len=256)                    :: message 
@@ -166,10 +111,6 @@ subroutine writeHRUToOutputStructure(&
   integer(i4b)                          :: iFreq             ! index of the output frequency
   ! convert the C pointers to Fortran pointers
   call c_f_pointer(handle_hru_data, hru_data)
-  call c_f_pointer(handle_statCounter, statCounter)
-  call c_f_pointer(handle_outputTimeStep, outputTimeStep)
-  call c_f_pointer(handle_resetStats, resetStats)
-  call c_f_pointer(handle_finalizeStats, finalizeStats)
   err=0; message='summa_manageOutputFiles/'
   ! identify the start of the writing
 
@@ -178,41 +119,41 @@ subroutine writeHRUToOutputStructure(&
                             newOutputFile,  defNewOutputFile,            &
                             ixRestart,      printRestart,                &   ! flag to print the restart file
                             ixProgress,     printProgress,               &   ! flag to print simulation progress
-                            resetStats%dat, finalizeStats%dat,           &   ! flags to reset and finalize stats
-                            statCounter%var,                             &   ! statistics counter
+                            hru_data%resetStats%dat, hru_data%finalizeStats%dat,           &   ! flags to reset and finalize stats
+                            hru_data%statCounter%var,                             &   ! statistics counter
                             err, cmessage)                                  ! error control
   if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
 
  ! If we do not do this looping we segfault - I am not sure why
-  outputStructure(1)%finalizeStats%gru(indxGRU)%hru(indxHRU)%tim(outputStep)%dat(:) = finalizeStats%dat(:)
+  outputStructure(1)%finalizeStats%gru(indxGRU)%hru(indxHRU)%tim(outputStep)%dat(:) = hru_data%finalizeStats%dat(:)
 
  ! ****************************************************************************
  ! *** calculate output statistics
  ! ****************************************************************************
   do iStruct=1,size(structInfo)
     select case(trim(structInfo(iStruct)%structName))
-      case('forc'); call calcStats(hru_data%forcStat%var, hru_data%forcStruct%var, statForc_meta, resetStats%dat, finalizeStats%dat, statCounter%var, err, cmessage)
-      case('prog'); call calcStats(hru_data%progStat%var, hru_data%progStruct%var, statProg_meta, resetStats%dat, finalizeStats%dat, statCounter%var, err, cmessage)
-      case('diag'); call calcStats(hru_data%diagStat%var, hru_data%diagStruct%var, statDiag_meta, resetStats%dat, finalizeStats%dat, statCounter%var, err, cmessage)
-      case('flux'); call calcStats(hru_data%fluxStat%var, hru_data%fluxStruct%var, statFlux_meta, resetStats%dat, finalizeStats%dat, statCounter%var, err, cmessage)
-      case('indx'); call calcStats(hru_data%indxStat%var, hru_data%indxStruct%var, statIndx_meta, resetStats%dat, finalizeStats%dat, statCounter%var, err, cmessage)     
+      case('forc'); call calcStats(hru_data%forcStat%var, hru_data%forcStruct%var, statForc_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
+      case('prog'); call calcStats(hru_data%progStat%var, hru_data%progStruct%var, statProg_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
+      case('diag'); call calcStats(hru_data%diagStat%var, hru_data%diagStruct%var, statDiag_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
+      case('flux'); call calcStats(hru_data%fluxStat%var, hru_data%fluxStruct%var, statFlux_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
+      case('indx'); call calcStats(hru_data%indxStat%var, hru_data%indxStruct%var, statIndx_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)     
     end select
     if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
   end do  ! (looping through structures)
     
   ! calc basin stats
-  call calcStats(hru_data%bvarStat%var(:), hru_data%bvarStruct%var(:), statBvar_meta, resetStats%dat, finalizeStats%dat, statCounter%var, err, cmessage)
+  call calcStats(hru_data%bvarStat%var(:), hru_data%bvarStruct%var(:), statBvar_meta, hru_data%resetStats%dat, hru_data%finalizeStats%dat, hru_data%statCounter%var, err, cmessage)
   if(err/=0)then; message=trim(message)//trim(cmessage)//'[bvar stats]'; return; endif
   
   ! write basin-average variables
-  call writeBasin(indxGRU,indxHRU,outputStep,finalizeStats%dat, &
-                  outputTimeStep%var,bvar_meta,hru_data%bvarStat%var,hru_data%bvarStruct%var,bvarChild_map,err,cmessage)
+  call writeBasin(indxGRU,indxHRU,outputStep,hru_data%finalizeStats%dat, &
+                  hru_data%outputTimeStep%var,bvar_meta,hru_data%bvarStat%var,hru_data%bvarStruct%var,bvarChild_map,err,cmessage)
   if(err/=0)then; message=trim(message)//trim(cmessage)//'[bvar]'; return; endif
 
   ! ****************************************************************************
   ! *** write data
   ! ****************************************************************************
-  call writeTime(indxGRU,indxHRU,outputStep,finalizeStats%dat, &
+  call writeTime(indxGRU,indxHRU,outputStep,hru_data%finalizeStats%dat, &
                 time_meta,hru_data%timeStruct%var,err,message)
 
   ! write the model output to the OutputStructure
@@ -221,15 +162,15 @@ subroutine writeHRUToOutputStructure(&
   ! Thus, we must also pass the stats parent->child maps from childStruct.
   do iStruct=1,size(structInfo)
     select case(trim(structInfo(iStruct)%structName))
-      case('forc'); call writeData(indxGRU,indxHRU,outputStep,"forc",finalizeStats%dat,&
+      case('forc'); call writeData(indxGRU,indxHRU,outputStep,"forc",hru_data%finalizeStats%dat,&
                     maxLayers,forc_meta,hru_data%forcStat,hru_data%forcStruct,forcChild_map,hru_data%indxStruct,err,cmessage)
-      case('prog'); call writeData(indxGRU,indxHRU,outputStep,"prog",finalizeStats%dat,&
+      case('prog'); call writeData(indxGRU,indxHRU,outputStep,"prog",hru_data%finalizeStats%dat,&
                     maxLayers,prog_meta,hru_data%progStat,hru_data%progStruct,progChild_map,hru_data%indxStruct,err,cmessage)
-      case('diag'); call writeData(indxGRU,indxHRU,outputStep,"diag",finalizeStats%dat,&
+      case('diag'); call writeData(indxGRU,indxHRU,outputStep,"diag",hru_data%finalizeStats%dat,&
                     maxLayers,diag_meta,hru_data%diagStat,hru_data%diagStruct,diagChild_map,hru_data%indxStruct,err,cmessage)
-      case('flux'); call writeData(indxGRU,indxHRU,outputStep,"flux",finalizeStats%dat,&
+      case('flux'); call writeData(indxGRU,indxHRU,outputStep,"flux",hru_data%finalizeStats%dat,&
                     maxLayers,flux_meta,hru_data%fluxStat,hru_data%fluxStruct,fluxChild_map,hru_data%indxStruct,err,cmessage)
-      case('indx'); call writeData(indxGRU,indxHRU,outputStep,"indx",finalizeStats%dat,&
+      case('indx'); call writeData(indxGRU,indxHRU,outputStep,"indx",hru_data%finalizeStats%dat,&
                     maxLayers,indx_meta,hru_data%indxStat,hru_data%indxStruct,indxChild_map,hru_data%indxStruct,err,cmessage)
     end select
     if(err/=0)then 
@@ -244,12 +185,12 @@ subroutine writeHRUToOutputStructure(&
 
   ! increment output file timestep
   do iFreq = 1,maxvarFreq
-    statCounter%var(iFreq) = statCounter%var(iFreq)+1
-    if(finalizeStats%dat(iFreq)) outputTimeStep%var(iFreq) = outputTimeStep%var(iFreq) + 1
+    hru_data%statCounter%var(iFreq) = hru_data%statCounter%var(iFreq)+1
+    if(hru_data%finalizeStats%dat(iFreq)) hru_data%outputTimeStep%var(iFreq) = hru_data%outputTimeStep%var(iFreq) + 1
   end do
 
   ! if finalized stats, then reset stats on the next time step
-  resetStats%dat(:) = finalizeStats%dat(:)
+  hru_data%resetStats%dat(:) = hru_data%finalizeStats%dat(:)
 
   ! save time vector
   hru_data%oldTime_hru%var(:) = hru_data%timeStruct%var(:)
diff --git a/build/source/dshare/data_types.f90 b/build/source/dshare/data_types.f90
index 82d092c..fa351c6 100755
--- a/build/source/dshare/data_types.f90
+++ b/build/source/dshare/data_types.f90
@@ -714,6 +714,21 @@ MODULE data_types
     type(var_i),pointer                        :: finishTime_hru             ! end time for the model simulation
     type(var_i),pointer                        :: refTime_hru                ! reference time for the model simulation
     type(var_i),pointer                        :: oldTime_hru                ! time from previous step
+
+    ! Statistic flags 
+    type(var_i), pointer                       :: statCounter                ! time counter for stats
+    type(var_i), pointer                       :: outputTimeStep             ! timestep in output files
+    type(flagVec), pointer                     :: resetStats                 ! flags to finalize statistics
+    type(flagVec), pointer                     :: finalizeStats              ! flags to finalize statistics
+
+    ! Julian Day Variables
+    real(c_double)                             :: fracJulDay                 ! fractional julian days since the start of year
+    real(c_double)                             :: tmZoneOffsetFracDay        ! time zone offset in fractional days
+    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
+    
+
   end type hru_type
 
 
-- 
GitLab