diff --git a/build/includes/file_access_actor/file_access_actor.hpp b/build/includes/file_access_actor/file_access_actor.hpp
index 87af7b566778e0d99ecc7a825928f7507e546dd0..7a4e21cfc46599cf9c5a8df87ab3c22dee24a377 100644
--- a/build/includes/file_access_actor/file_access_actor.hpp
+++ b/build/includes/file_access_actor/file_access_actor.hpp
@@ -94,14 +94,6 @@ struct file_access_state {
 behavior file_access_actor(stateful_actor<file_access_state>* self, int startGRU, int numGRU, 
    File_Access_Actor_Settings file_access_actor_settings, actor parent);
 
-// Read in the attributes for all HRUs that are in the run-domain
-void readAttributes(stateful_actor<file_access_state>* self); 
-
-// read in the parameters for all HRUs that are in the run-domain
-void readParameters(stateful_actor<file_access_state>* self);
-
-// Read in the inital conditions for all the HRUs that are in the run-domain
-void readInitConditions(stateful_actor<file_access_state>* self);
 
 void initalizeOutputHandles(stateful_actor<file_access_state>* self);
 
diff --git a/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp b/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp
index 0d97c855add3a1630fa8fb0bd134409bd41ca12e..002e4c97cb8de6ea6f787fd6d0ea2364866c29ff 100644
--- a/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp
+++ b/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp
@@ -15,38 +15,8 @@ extern "C" {
                                int* err);
 
   // OutputStructure and Output functions
-  // void initOutputStructure(void* handle_forcFileInfo, int* max_steps, int* numGRU, int* err);
   void deallocateOutputStructure(int* err);
   void writeOutput_fortran(void* handle_ncid, int* num_steps, int* start_gru, int* max_gru, int* err);
-  // void initOutputTimeStep(int* num_gru, int* err);
-
-  // Attributes Files- called inside initalizeFileAccessActor
-  void allocateAttributeStructures(int* index_gru, int* index_hru, void* handle_attr_struct, 
-      void* handle_type_struct, void* handle_id_struct, int* err);
-  void openAttributeFile(int* att_ncid, int* err);
-  void getNumVarAttr(int* attr_ncid, int* num_var_attr, int* err);
-  void closeAttributeFile(int* attr_ncid, int* err);
-  void readAttributeFromNetCDF(int* attr_ncid, int* index_gru, int* index_hru, int* num_var,
-    void* attr_array, void* type_array, void* id_array, int* err);
-  
-  // Parameters File - called inside initalizeFileAccessActor
-  void allocateParamStructures(int* index_gru, int* index_hru, void* handle_dpar_struct,
-      void* handle_mpar_struct, void* handle_bpar_struct, int* err);
-  void openParamFile(int* param_ncid, bool* param_file_exists, int* err);
-  void getNumVarParam(int* param_ncid, int* num_var_param, int* err);
-  void closeParamFile(int* param_ncid, int* err);
-  void getParamSizes(int* dpar_array_size, int* bpar_array_size, int* type_array_size);
-  void overwriteParam(int* index_gru, int* index_hru,
-    void* handle_type_struct, void* handle_dpar_struct, void* handle_mpar_struct, 
-    void* handle_bpar_struct, int* err);
-  void readParamFromNetCDF(int* param_ncid, int* index_gru, int* index_hru, int* start_index_gru,
-    int* num_var_param, void* handle_mpar_struct, void* _handle_bpar_struct, int* err);
-
-  // Set up global initial conditions
-  void openInitCondFile(int* init_cond_ncid, int* err);
-  void closeInitCondFile(int* init_cond_ncid, int* err);
-  void readInitCond_prog(int* init_cond_ncid, int* start_gru, int* num_gru, int* err);
-  void readInitCond_bvar(int* init_cond_ncid, int* start_gru, int* num_gru, int* err);
 
   void updateFailed(int* indxHRU);
 
diff --git a/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp b/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp
index cee5ba4cc3810f0a373421db1a163dbdf9fb0288..5a1656efb72b766a563b0cb62ce663691eb231d0 100644
--- a/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp
+++ b/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp
@@ -7,11 +7,13 @@ extern "C" {
             // Statistics Structures
             void* forcStat, void* progStat, void* diagStat, void* fluxStat, void* indxStat, void* bvarStat,
             // Primary Data Structures (scalars) 
-            void* timeStruct, void* forcStruct,
+            void* timeStruct, void* forcStruct, void* attrStruct, void* typeStruct, void* idStruct,
             // primary data structures (variable length vectors)
-            void* indxStruct, void* progStruct, void* diagStruct, void* fluxStruct,
+            void* indxStruct, void* mparStruct, void* progStruct, void* diagStruct, void* fluxStruct,
             // basin-average structures
-            void* bvarStruct,
+            void* bvarStruct, void* bparStruct,
+            // ancillary data structures
+            void* dparStruct,
             // local HRU data 
             void* startTime, void* finshTime, void* refTime, void* oldTime, int* err);
 
@@ -23,7 +25,7 @@ extern "C" {
             // primary data structures (scalars)
             void* attrStruct, void* typeStruct, void* idStruct,
             // primary data structures (variable length vectors)
-            void* mparStruct, void* bparStruct, void* bvarStruct, void* dparStruct,
+            void* indexStruct, void* mparStruct, void* progStruct, void* bparStruct, void* bvarStruct, void* dparStruct,
             // lookup tables
             void* lookupStruct,
             // local HRU data
diff --git a/build/source/actors/file_access_actor/cpp_code/file_access_actor.cpp b/build/source/actors/file_access_actor/cpp_code/file_access_actor.cpp
index 2284a8a48459beb114324b136eee2757f191fe3b..3b0b55a860446ff410f66207934866c4c2c92330 100644
--- a/build/source/actors/file_access_actor/cpp_code/file_access_actor.cpp
+++ b/build/source/actors/file_access_actor/cpp_code/file_access_actor.cpp
@@ -44,12 +44,6 @@ behavior file_access_actor(stateful_actor<file_access_state>* self, int start_gr
 
     aout(self) << "Simluations Steps: " << self->state.num_steps << "\n";
 
-    // Read in the attribute and parameter information for the HRUs to request
-    readAttributes(self);
-    readParameters(self);
-
-    // read in the inital conditions for the grus/hrus
-    readInitConditions(self);
     
     // Inital Files Have Been Loaded - Send Message to Job_Actor to Start Simulation
     self->send(self->state.parent, init_gru_v);
@@ -66,11 +60,10 @@ behavior file_access_actor(stateful_actor<file_access_state>* self, int start_gr
     }
 
     // Set up the output container
-    self->state.output_container = new Output_Container(
-        self->state.file_access_actor_settings.num_partitions_in_output_buffer,
-        self->state.num_gru,
-        self->state.file_access_actor_settings.num_timesteps_in_output_buffer,
-        self->state.num_steps); 
+    self->state.output_container = new Output_Container(self->state.file_access_actor_settings.num_partitions_in_output_buffer,
+                                                        self->state.num_gru,
+                                                        self->state.file_access_actor_settings.num_timesteps_in_output_buffer,
+                                                        self->state.num_steps); 
 
     return {
         [=](write_param, int index_gru, int index_hru, std::vector<double> attr_struct, 
@@ -248,113 +241,6 @@ behavior file_access_actor(stateful_actor<file_access_state>* self, int start_gr
 }
 
 
-void readAttributes(stateful_actor<file_access_state>* self) {
-
-    int err = 0;
-    openAttributeFile(&self->state.attribute_ncid, &err);
-    
-    getNumVarAttr(&self->state.attribute_ncid, &self->state.num_var_in_attributes_file, &err);
-    
-    for (int index_gru = 1; index_gru < self->state.num_gru + 1; index_gru++) {
-
-        void* handle_attr_struct = new_handle_var_d();
-        void* handle_type_struct = new_handle_var_i();
-        void* handle_id_struct   = new_handle_var_i8();
-        int index_hru = 1;
-
-        allocateAttributeStructures(&index_gru, &index_hru, handle_attr_struct, handle_type_struct,
-            handle_id_struct, &err);
-
-        readAttributeFromNetCDF(&self->state.attribute_ncid, &index_gru, &index_hru,
-            &self->state.num_var_in_attributes_file, handle_attr_struct, handle_type_struct,
-            handle_id_struct, &err);
-        
-        // attr struct
-        std::vector<double> attr_struct_to_push = get_var_d(handle_attr_struct);
-        self->state.attr_structs_for_hrus.push_back(attr_struct_to_push);
-        delete_handle_var_d(handle_attr_struct);
-        // type struct
-        std::vector<int> type_struct_to_push = get_var_i(handle_type_struct);
-        self->state.type_structs_for_hrus.push_back(type_struct_to_push);
-        delete_handle_var_i(handle_type_struct);
-        // id struct
-        std::vector<long int> id_struct_to_push = get_var_i8(handle_id_struct);
-        self->state.id_structs_for_hrus.push_back(id_struct_to_push);
-        delete_handle_var_i8(handle_id_struct);
-    }
-
-    closeAttributeFile(&self->state.attribute_ncid, &err);
-}
-
-void readParameters(stateful_actor<file_access_state>* self) {
-
-    int err = 0;
-    int index_hru = 1;
-
-    openParamFile(&self->state.param_ncid, &self->state.param_file_exists, 
-        &err);
-
-    getParamSizes(&self->state.dpar_array_size, &self->state.bpar_array_size,
-        &self->state.type_array_size);
-    
-
-    if (self->state.param_file_exists) {
-        getNumVarParam(&self->state.param_ncid, &self->state.num_var_in_param_file,
-            &err);
-    } else {
-        self->state.num_var_in_param_file = self->state.type_array_size;
-    }
-
-    for (int index_gru = 1; index_gru < self->state.num_gru + 1; index_gru++) {
-
-        std::vector<double> dpar_array(self->state.dpar_array_size);
-        void* handle_type_struct = new_handle_var_i();
-        void* handle_dpar_struct = new_handle_var_d();
-        void* handle_mpar_struct = new_handle_var_dlength();      
-        void* handle_bpar_struct = new_handle_var_d();  
-        std::vector<double> bpar_array(self->state.dpar_array_size);        
-
-        allocateParamStructures(&index_gru, &index_hru, handle_dpar_struct, 
-            handle_mpar_struct, handle_bpar_struct, &err);
-
-        // need to convert attr_struct to FORTRAN format   
-        set_var_i(self->state.type_structs_for_hrus[index_gru-1], handle_type_struct); 
-    
-        overwriteParam(&index_gru, &index_hru, 
-            handle_type_struct,
-            handle_dpar_struct, 
-            handle_mpar_struct, 
-            handle_bpar_struct, 
-            &err);
-
-        if (self->state.param_file_exists) {
-            readParamFromNetCDF(&self->state.param_ncid, &index_gru, &index_hru,
-                &self->state.start_gru, 
-                &self->state.num_var_in_param_file, 
-                handle_mpar_struct, 
-                handle_bpar_struct,
-                &err);
-        }
-
-        // type_struct
-        delete_handle_var_i(handle_type_struct);
-        
-        // dpar_struct
-        std::vector<double> dpar_struct_to_push = get_var_d(handle_dpar_struct);
-        self->state.dpar_structs_for_hrus.push_back(dpar_struct_to_push);
-        delete_handle_var_d(handle_dpar_struct);
-        // mpar_struct
-        std::vector<std::vector<double>> mpar_struct_to_push = get_var_dlength(handle_mpar_struct);
-        self->state.mpar_structs_for_hrus.push_back(mpar_struct_to_push);
-        delete_handle_var_dlength(handle_mpar_struct);
-        // bpar_struct
-        std::vector<double> bpar_struct_to_push = get_var_d(handle_bpar_struct);
-        self->state.bpar_structs_for_hrus.push_back(bpar_struct_to_push);
-        delete_handle_var_d(handle_bpar_struct);
-    }
-    closeParamFile(&self->state.param_ncid, &err);
-}
-
 void writeOutput(stateful_actor<file_access_state>* self, Output_Partition* partition) {
                 
     int num_timesteps_to_write = partition->getNumStoredTimesteps();
@@ -378,12 +264,4 @@ void writeOutput(stateful_actor<file_access_state>* self, Output_Partition* part
     partition->resetReadyToWriteList();
 }
 
-void readInitConditions(stateful_actor<file_access_state>* self) {
-    int err;
-    openInitCondFile(&self->state.init_cond_ncid, &err);
-    readInitCond_prog(&self->state.init_cond_ncid, &self->state.start_gru, &self->state.num_gru, &err);
-    readInitCond_bvar(&self->state.init_cond_ncid, &self->state.start_gru, &self->state.num_gru, &err);
-    closeInitCondFile(&self->state.init_cond_ncid, &err); 
-}
-
 } // end namespace
\ 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 ac0b842092de90aa9c27f75f1b8dd433a96413d0..f4b881f01c25e8d51dc9aee20628405f5f207324 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
@@ -6,6 +6,9 @@ module cppwrap_fileAccess
   USE nrtype
   USE data_types
   USE globalData
+  USE globalData,only:integerMissing      ! missing integer value
+  USE globalData,only:realMissing         ! missing double precision value
+
   USE var_lookup,only:maxvarFreq                ! maximum number of output files
 
 
@@ -36,18 +39,52 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   USE def_output_actors_module,only:def_output                ! module to define output variables
   USE output_structure_module,only:initOutputStructure        ! module to initialize output structure
   USE output_structure_module,only:initOutputTimeStep         ! module to initialize output timestep structure (tracks GRUs timestep for output)
-
+  USE read_attrb_module,only:read_attrb                       ! module to read local attributes
+  USE read_param_module,only:read_param                       ! module to read model parameter sets
+  USE pOverwrite_module,only:pOverwrite                       ! module to overwrite default parameter values with info from the Noah tables
+  USE var_derive_module,only:fracFuture                       ! module to calculate the fraction of runoff in future time steps (time delay histogram)
+  USE paramCheck_module,only:paramCheck                       ! module to check consistency of model parameters
+  USE read_icond_module,only:read_icond                       ! module to read initial conditions
+  USE check_icond_module,only:check_icond                     ! module to check initial conditions
+
+  USE mDecisions_module,only:&
+                              sameRulesAllLayers, & ! SNTHERM option: same combination/sub-dividion rules applied to all layers
+                              rulesDependLayerIndex ! CLM option: combination/sub-dividion rules depend on layer index
   USE globalData,only:localParFallback                        ! local column default parameters
   USE globalData,only:basinParFallback                        ! basin-average default parameters
   USE summaFileManager,only:LOCALPARAM_INFO,BASINPARAM_INFO   ! files defining the default values and constraints for model parameters
   USE globalData,only:mpar_meta,bpar_meta                     ! parameter metadata structures
   USE summaFileManager,only:SETTINGS_PATH                     ! define path to settings files (e.g., parameters, soil and veg. tables)
+  USE summaFileManager,only:LOCAL_ATTRIBUTES                  ! name of model initial attributes file
   USE summaFileManager,only:GENPARM,VEGPARM,SOILPARM,MPTABLE  ! files defining the noah tables
+  USE summaFileManager,only:MODEL_INITCOND                    ! name of model initial conditions file
+  USE summaFileManager,only:STATE_PATH                        ! optional path to state/init. condition files (defaults to SETTINGS_PATH)
   USE globalData,only:model_decisions                         ! model decision structure
   USE var_lookup,only:iLookDECISIONS                          ! look-up values for model decisions
-
+  USE var_lookup,only:iLookTYPE                               ! look-up values for model types
+  USE var_lookup,only:iLookID                                 ! look-up values for model IDs
+  USE var_lookup,only:iLookPARAM
+  USE var_lookup,only:iLookATTR                               ! look-up values for model attributes
+  USE var_lookup,only:iLookBVAR                               ! look-up values for basin-average variables
+  USE output_structure_module,only:outputStructure            ! output structure
   USE globalData,only:failedHRUs                              ! Flag for file access actor to know which GRUs have failed
+  
+  USE globalData,only:iRunModeFull,iRunModeGRU,iRunModeHRU
+  USE globalData,only:iRunMode                                ! define the current running mode
+  USE globalData,only:checkHRU                                ! index of the HRU for a single HRU run
+  
+  ! look-up values for the choice of heat capacity computation
+  USE mDecisions_module,only:  &
+  enthalpyFD                             ! heat capacity using enthalpy
+  USE mDecisions_module,only:&
+ 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
+
+  USE ConvE2Temp_module,only:E2T_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
+  USE t2enthalpy_module,only:T2E_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
 
+  USE NOAHMP_VEG_PARAMETERS,only:SAIM,LAIM                    ! 2-d tables for stem area index and leaf area index (vegType,month)
+  USE NOAHMP_VEG_PARAMETERS,only:HVT,HVB                      ! height at the top and bottom of vegetation (vegType)
 
   USE globalData,only:numtim                 ! number of time steps in the simulation
 
@@ -67,9 +104,16 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
 
   ! local Variables
   type(file_info_array),pointer          :: forcFileInfo
-  type(var_i),pointer                    :: output_ncid                        ! id of output file
+  type(var_i),pointer                    :: output_ncid        ! id of output file
+  integer(i4b)                           :: iGRU               ! counter for GRUs
+  integer(i4b)                           :: iHRU               ! counter for HRUs
+  integer(i4b)                           :: jHRU,kHRU          ! HRU indices
+  integer(i4b)                           :: ivar               ! counter for variables
+  character(len=256)                     :: attrFile           ! attributes file name
+  character(LEN=256)                    :: restartFile        ! restart file name
   integer(i4b)                           :: indxGRU=1
-  character(len=256)                     :: message         ! error message for downwind routine
+  character(len=256)                     :: message            ! error message for downwind routine
+
 
   err=0; message="fileAccessActor_init_fortran/"
 
@@ -85,6 +129,16 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   if(err/=0)then; print*,trim(message); return; endif
   num_timesteps = numtim
 
+    ! get the maximum number of snow layers
+  select case(model_decisions(iLookDECISIONS%snowLayers)%iDecision)
+    case(sameRulesAllLayers);    maxSnowLayers = 100
+    case(rulesDependLayerIndex); maxSnowLayers = 5
+    case default; err=20; message=trim(message)//'unable to identify option to combine/sub-divide snow layers';print*,message;return
+  end select ! (option to combine/sub-divide snow layers)
+
+
+  maxLayers = gru_struc(1)%hruInfo(1)%nSoil + maxSnowLayers
+
   ! *****************************************************************************
   ! *** read default model parameters
   ! *****************************************************************************
@@ -96,20 +150,37 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   call read_pinit(BASINPARAM_INFO,.FALSE.,bpar_meta,basinParFallback,err,message)
   if(err/=0)then; print*,trim(message); return; endif
   
+  
+  ! *****************************************************************************
+  ! *** read Noah vegetation and soil tables
+  ! *****************************************************************************
 
- ! read Noah soil and vegetation tables
+  greenVegFrac_monthly = (/0.01_dp, 0.02_dp, 0.03_dp, 0.07_dp, 0.50_dp, 0.90_dp, 0.95_dp, 0.96_dp, 0.65_dp, 0.24_dp, 0.11_dp, 0.02_dp/)
+
+
+  ! read Noah soil and vegetation tables
   call soil_veg_gen_parm(trim(SETTINGS_PATH)//trim(VEGPARM),      & ! filename for vegetation table
-      trim(SETTINGS_PATH)//trim(SOILPARM),                        & ! filename for soils table
-      trim(SETTINGS_PATH)//trim(GENPARM),                         & ! filename for general table
-      trim(model_decisions(iLookDECISIONS%vegeParTbl)%cDecision), & ! classification system used for vegetation
-      trim(model_decisions(iLookDECISIONS%soilCatTbl)%cDecision))   ! classification system used for soils
+                         trim(SETTINGS_PATH)//trim(SOILPARM),                        & ! filename for soils table
+                         trim(SETTINGS_PATH)//trim(GENPARM),                         & ! filename for general table
+                         trim(model_decisions(iLookDECISIONS%vegeParTbl)%cDecision), & ! classification system used for vegetation
+                         trim(model_decisions(iLookDECISIONS%soilCatTbl)%cDecision))   ! classification system used for soils
   if(err/=0)then; print*,trim(message); return; endif
 
   ! read Noah-MP vegetation tables
   call read_mp_veg_parameters(trim(SETTINGS_PATH)//trim(MPTABLE),                       & ! filename for Noah-MP table
-       trim(model_decisions(iLookDECISIONS%vegeParTbl)%cDecision)) ! classification system used for vegetation
+                              trim(model_decisions(iLookDECISIONS%vegeParTbl)%cDecision)) ! classification system used for vegetation
   if(err/=0)then; print*,trim(message); return; endif
 
+    ! define urban vegetation category
+  select case(trim(model_decisions(iLookDECISIONS%vegeParTbl)%cDecision))
+    case('USGS');                     urbanVegCategory =    1
+    case('MODIFIED_IGBP_MODIS_NOAH'); urbanVegCategory =   13
+    case('plumberCABLE');             urbanVegCategory = -999
+    case('plumberCHTESSEL');          urbanVegCategory = -999
+    case('plumberSUMMA');             urbanVegCategory = -999
+    case default; message=trim(message)//'unable to identify vegetation category';print*,message;return
+  end select
+
   ! *****************************************************************************
   ! *** Initalize failed HRU tracker
   ! *****************************************************************************
@@ -136,6 +207,172 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   call initOutputTimeStep(num_gru, err)
   if(err/=0)then; print*,trim(message); return; endif
 
+
+  ! *****************************************************************************
+  ! *** Read Attributes
+  ! *****************************************************************************
+
+  attrFile = trim(SETTINGS_PATH)//trim(LOCAL_ATTRIBUTES)
+  call read_attrb(trim(attrFile),num_gru,outputStructure(1)%attrStruct,&
+                  outputStructure(1)%typeStruct,outputStructure(1)%idStruct,err,message)
+  if(err/=0)then; print*,trim(message); return; endif
+
+
+  ! set default model parameters
+  do iGRU=1, num_gru
+    do iHRU=1, gru_struc(iGRU)%hruCount
+      ! set parmameters to their default value
+      outputStructure(1)%dparStruct%gru(iGRU)%hru(iHRU)%var(:) = localParFallback(:)%default_val         ! x%hru(:)%var(:)
+
+      ! overwrite default model parameters with information from the Noah-MP tables
+      call pOverwrite(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex),  &  ! vegetation category
+                      outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%soilTypeIndex), &  ! soil category
+                      outputStructure(1)%dparStruct%gru(iGRU)%hru(iHRU)%var,                          &  ! default model parameters
+                      err,message)                                                   ! error control
+      if(err/=0)then; print*, trim(message); return; endif
+
+
+   ! copy over to the parameter structure
+   ! NOTE: constant for the dat(:) dimension (normally depth)
+      do ivar=1,size(localParFallback)
+        outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(ivar)%dat(:) = outputStructure(1)%dparStruct%gru(iGRU)%hru(iHRU)%var(ivar)
+      end do  ! looping through variables
+    
+    end do  ! looping through HRUs
+    
+    ! set default for basin-average parameters
+    outputStructure(1)%bparStruct%gru(iGRU)%var(:) = basinParFallback(:)%default_val
+    
+  end do  ! looping through GRUs
+
+
+  ! *****************************************************************************
+  ! *** Read Parameters
+  ! *****************************************************************************
+  checkHRU = integerMissing
+  call read_param(iRunMode,checkHRU,start_gru,num_hru,num_gru,outputStructure(1)%idStruct,&
+                  outputStructure(1)%mparStruct,outputStructure(1)%bparStruct,err,message)
+  if(err/=0)then; print*,trim(message); return; endif
+
+  ! *****************************************************************************
+  ! *** compute derived model variables that are pretty much constant for the basin as a whole
+  ! *****************************************************************************
+  ! ! loop through GRUs
+  do iGRU=1,num_gru
+    ! calculate the fraction of runoff in future time steps
+    call fracFuture(outputStructure(1)%bparStruct%gru(iGRU)%var,    &  ! vector of basin-average model parameters
+                    outputStructure(1)%bvarStruct_init%gru(iGRU),    &  ! data structure of basin-average variables
+                    err,message)                   ! error control
+    if(err/=0)then; print*, trim(message); return; endif
+
+    ! loop through local HRUs
+    do iHRU=1,gru_struc(iGRU)%hruCount
+
+    
+      kHRU=0
+      ! check the network topology (only expect there to be one downslope HRU)
+      do jHRU=1,gru_struc(iGRU)%hruCount
+      if(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%downHRUindex) == outputStructure(1)%idStruct%gru(iGRU)%hru(jHRU)%var(iLookID%hruId))then
+        if(kHRU==0)then  ! check there is a unique match
+        kHRU=jHRU
+        else
+        message=trim(message)//'only expect there to be one downslope HRU'; print*, message; return
+        end if  ! (check there is a unique match)
+      end if  ! (if identified a downslope HRU)
+      end do
+
+
+      ! check that the parameters are consistent
+      call paramCheck(outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU),err,message)
+      if(err/=0)then; print*, message; return; endif
+
+
+      ! calculate a look-up table for the temperature-enthalpy conversion: snow
+      ! NOTE1: this should eventually be replaced by the more general routine below
+      ! NOTE2: this does not actually need to be called for each HRU and GRU
+      call E2T_lookup(outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU),err,message)
+      if(err/=0)then; print*,message; return; endif
+
+      ! calculate a lookup table to compute enthalpy from temperature, only for enthalpyFD
+      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
+                        outputStructure(1)%lookupStruct%gru(iGRU)%hru(iHRU),      &   ! intent(inout): lookup table data structure
+                        err,message)                              ! intent(out):   error control
+        if(err/=0)then; print*, message; return; endif
+      endif
+
+      ! overwrite the vegetation height
+      HVT(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex)) = outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%heightCanopyTop)%dat(1)
+      HVB(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex)) = outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%heightCanopyBottom)%dat(1)
+         
+      ! overwrite the tables for LAI and SAI
+      if(model_decisions(iLookDECISIONS%LAI_method)%iDecision == specified)then
+        SAIM(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex),:) = outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%winterSAI)%dat(1)
+        LAIM(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex),:) = outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%summerLAI)%dat(1)*greenVegFrac_monthly
+      endif
+
+    end do ! HRU
+    
+    ! compute total area of the upstream HRUS that flow into each HRU
+    do iHRU=1,gru_struc(iGRU)%hruCount
+      outputStructure(1)%upArea%gru(iGRU)%hru(iHRU) = 0._rkind
+      do jHRU=1,gru_struc(iGRU)%hruCount
+       ! check if jHRU flows into iHRU; assume no exchange between GRUs
+       if(outputStructure(1)%typeStruct%gru(iGRU)%hru(jHRU)%var(iLookTYPE%downHRUindex)==outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookID%hruId))then
+        outputStructure(1)%upArea%gru(iGRU)%hru(iHRU) = outputStructure(1)%upArea%gru(iGRU)%hru(iHRU) + outputStructure(1)%attrStruct%gru(iGRU)%hru(jHRU)%var(iLookATTR%HRUarea)
+       endif   ! (if jHRU is an upstream HRU)
+      end do  ! jHRU
+    end do  ! iHRU
+  
+    ! identify the total basin area for a GRU (m2)  
+    outputStructure(1)%bvarStruct_init%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) = 0._rkind
+    do iHRU=1,gru_struc(iGRU)%hruCount
+      outputStructure(1)%bvarStruct_init%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) = &
+      outputStructure(1)%bvarStruct_init%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) + outputStructure(1)%attrStruct%gru(iGRU)%hru(iHRU)%var(iLookATTR%HRUarea)
+    end do
+  
+  end do ! GRU
+
+
+
+
+
+
+  ! *****************************************************************************
+  ! Restart File
+  ! *****************************************************************************
+  ! define restart file path/name
+  if(STATE_PATH == '') then
+    restartFile = trim(SETTINGS_PATH)//trim(MODEL_INITCOND)
+  else
+    restartFile = trim(STATE_PATH)//trim(MODEL_INITCOND)
+  endif
+
+ ! read initial conditions
+  call read_icond(restartFile,                        & ! intent(in):    name of initial conditions file
+                  num_gru,                            & ! intent(in):    number of response units
+                  outputStructure(1)%mparStruct,      & ! intent(in):    model parameters
+                  outputStructure(1)%progStruct_init, & ! intent(inout): model prognostic variables
+                  outputStructure(1)%bvarStruct_init, & ! intent(inout): model basin (GRU) variables
+                  outputStructure(1)%indxStruct_init, & ! intent(inout): model indices
+                  err,message)                          ! intent(out):   error control
+  if(err/=0)then; print*, message; return; endif
+
+  call check_icond(num_gru,                            &
+                   outputStructure(1)%progStruct_init, &  ! intent(inout): model prognostic variables
+                   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
+  
+
+
+
+
+
+
+  
 end subroutine fileAccessActor_init_fortran
 
 
@@ -164,8 +401,6 @@ subroutine FileAccessActor_DeallocateStructures(handle_forcFileInfo, handle_ncid
   USE globalData,only:forcingDataStruct
   USE globalData,only:vectime
   USE globalData,only:outputTimeStep
-  USE globalData,only:init_cond_prog
-  USE globalData,only:init_cond_bvar
   implicit none
   type(c_ptr),intent(in), value        :: handle_forcFileInfo
   type(c_ptr),intent(in), value        :: handle_ncid
@@ -193,8 +428,6 @@ subroutine FileAccessActor_DeallocateStructures(handle_forcFileInfo, handle_ncid
   deallocate(ncid)
   deallocate(failedHRUs)
   deallocate(outputTimeStep)
-  deallocate(init_cond_prog)
-  if (allocated(init_cond_bvar))then; deallocate(init_cond_bvar); endif;
 end subroutine FileAccessActor_DeallocateStructures
 
 
diff --git a/build/source/actors/file_access_actor/fortran_code/output_structure.f90 b/build/source/actors/file_access_actor/fortran_code/output_structure.f90
index b497bb59717ae80a6e3f542c11c72b41eced4af4..c1e85b115115a8455257afd42c3ff28e6ed9283f 100644
--- a/build/source/actors/file_access_actor/fortran_code/output_structure.f90
+++ b/build/source/actors/file_access_actor/fortran_code/output_structure.f90
@@ -13,6 +13,7 @@ module output_structure_module
                       var_ilength,         & ! x%var(:)%dat        (i4b)
                       var_dlength,         & ! x%var(:)%dat        (rkind)
                       ! gru dimension
+                      gru_d,               & ! x%gru(:)%var(:)     (rkind)  
                       gru_int,             & ! x%gru(:)%var(:)     (i4b)
                       gru_int8,            & ! x%gru(:)%var(:)     integer(8)
                       gru_double,          & ! x%gru(:)%var(:)     (rkind)
@@ -59,43 +60,50 @@ module output_structure_module
   private::is_var_desired
 
   type, public :: summa_output_type
+    type(gru_hru_z_vLookup)                          :: lookupStruct                   ! x%gru(:)%hru(:)%z(:)%var(:)%lookup(:) -- lookup tables
+
     ! define the statistics structures
-    type(gru_hru_time_doubleVec),allocatable          :: forcStat(:)                   ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model forcing data
-    type(gru_hru_time_doubleVec),allocatable          :: progStat(:)                   ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model prognostic (state) variables
-    type(gru_hru_time_doubleVec),allocatable          :: diagStat(:)                   ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model diagnostic variables
-    type(gru_hru_time_doubleVec),allocatable          :: fluxStat(:)                   ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model fluxes
-    type(gru_hru_time_doubleVec),allocatable          :: indxStat(:)                   ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model indices
-    type(gru_hru_time_doubleVec),allocatable          :: bvarStat(:)                   ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- basin-average variabl
+    type(gru_hru_time_doubleVec)                      :: forcStat                      ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model forcing data
+    type(gru_hru_time_doubleVec)                      :: progStat                      ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model prognostic (state) variables
+    type(gru_hru_time_doubleVec)                      :: diagStat                      ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model diagnostic variables
+    type(gru_hru_time_doubleVec)                      :: fluxStat                      ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model fluxes
+    type(gru_hru_time_doubleVec)                      :: indxStat                      ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model indices
+    type(gru_hru_time_doubleVec)                      :: bvarStat                      ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- basin-average variabl
 
     ! define the primary data structures (scalars)
-    type(gru_hru_time_int),allocatable                :: timeStruct(:)                 ! x%gru(:)%hru(:)%var(:)%tim(:)     -- model time data
-    type(gru_hru_time_double),allocatable             :: forcStruct(:)                 ! x%gru(:)%hru(:)%var(:)%tim(:)     -- model forcing data
-    type(gru_hru_double),allocatable                  :: attrStruct(:)                 ! x%gru(:)%hru(:)%var(:)            -- local attributes for each HRU, DOES NOT CHANGE OVER TIMESTEPS
-    type(gru_hru_int),allocatable                     :: typeStruct(:)                 ! x%gru(:)%hru(:)%var(:)%tim(:)     -- local classification of soil veg etc. for each HRU, DOES NOT CHANGE OVER TIMESTEPS
-    type(gru_hru_int8),allocatable                    :: idStruct(:)                   ! x%gru(:)%hru(:)%var(:)
+    type(gru_hru_time_int)                            :: timeStruct                    ! x%gru(:)%hru(:)%var(:)%tim(:)     -- model time data
+    type(gru_hru_time_double)                         :: forcStruct                    ! x%gru(:)%hru(:)%var(:)%tim(:)     -- model forcing data
+    type(gru_hru_double)                              :: attrStruct                    ! x%gru(:)%hru(:)%var(:)            -- local attributes for each HRU, DOES NOT CHANGE OVER TIMESTEPS
+    type(gru_hru_int)                                 :: typeStruct                    ! x%gru(:)%hru(:)%var(:)            -- local classification of soil veg etc. for each HRU, DOES NOT CHANGE OVER TIMESTEPS
+    type(gru_hru_int8)                                :: idStruct                      ! x%gru(:)%hru(:)%var(:)
 
     ! define the primary data structures (variable length vectors)
-    type(gru_hru_time_intVec),allocatable             :: indxStruct(:)                 ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model indices
-    type(gru_hru_doubleVec),allocatable               :: mparStruct(:)                 ! x%gru(:)%hru(:)%var(:)%dat        -- model parameters, DOES NOT CHANGE OVER TIMESTEPS TODO: MAYBE
-    type(gru_hru_time_doubleVec),allocatable          :: progStruct(:)                 ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model prognostic (state) variables
-    type(gru_hru_time_doubleVec),allocatable          :: diagStruct(:)                 ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model diagnostic variables
-    type(gru_hru_time_doubleVec),allocatable          :: fluxStruct(:)                 ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model fluxes
+    type(gru_hru_time_intVec)                         :: indxStruct                    ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model indices
+    type(gru_hru_intVec)                              :: indxStruct_init               ! x%gru(:)%hru(:)%var(:)%dat        -- model indices
+    type(gru_hru_doubleVec)                           :: mparStruct                    ! x%gru(:)%hru(:)%var(:)%dat        -- model parameters, DOES NOT CHANGE OVER TIMESTEPS TODO: MAYBE
+    type(gru_hru_time_doubleVec)                      :: progStruct                    ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model prognostic (state) variables
+    type(gru_hru_doubleVec)                           :: progStruct_init               ! x%gru(:)%hru(:)%var(:)%dat        -- model prognostic (state) variables
+    type(gru_hru_time_doubleVec)                      :: diagStruct                    ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model diagnostic variables
+    type(gru_hru_time_doubleVec)                      :: fluxStruct                    ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- model fluxes
 
     ! define the basin-average structures
-    type(gru_double),allocatable                      :: bparStruct(:)                 ! x%gru(:)%var(:)                   -- basin-average parameters, DOES NOT CHANGE OVER TIMESTEPS
-    type(gru_hru_time_doubleVec),allocatable          :: bvarStruct(:)                 ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- basin-average variables
-
+    type(gru_double)                                  :: bparStruct                    ! x%gru(:)%var(:)                   -- basin-average parameters, DOES NOT CHANGE OVER TIMESTEPS
+    type(gru_hru_time_doubleVec)                      :: bvarStruct                    ! x%gru(:)%hru(:)%var(:)%tim(:)%dat -- basin-average variables
+    type(gru_doubleVec)                               :: bvarStruct_init               ! x%gru(:)%hru(:)%var(:)%dat        -- basin-average variables
     ! define the ancillary data structures
-    type(gru_hru_double),allocatable                  :: dparStruct(:)                 ! x%gru(:)%hru(:)%var(:)
+    type(gru_hru_double)                              :: dparStruct                    ! x%gru(:)%hru(:)%var(:)
 
     ! finalize stats structure
-    type(gru_hru_time_flagVec),allocatable            :: finalizeStats(:)              ! x%gru(:)%hru(:)%tim(:)%dat -- flags on when to write to file
+    type(gru_hru_time_flagVec)                        :: finalizeStats                 ! x%gru(:)%hru(:)%tim(:)%dat -- flags on when to write to file
+
+
+    type(gru_d)                                       :: upArea
 
     integer(i4b)                                      :: nTimeSteps
   end type summa_output_type  
   
   
-  type(summa_output_type),allocatable,save,public :: outputStructure(:) ! summa_OutputStructure(iFile)%struc%var(:)%dat(nTimeSteps) 
+  type(summa_output_type),allocatable,save,public :: outputStructure(:) ! summa_OutputStructure(1)%struc%var(:)%dat(nTimeSteps) 
   
   contains
 
@@ -124,6 +132,7 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
   USE globalData,only:prog_meta,diag_meta,flux_meta,id_meta   ! metadata structures
   USE globalData,only:mpar_meta,indx_meta                     ! metadata structures
   USE globalData,only:bpar_meta,bvar_meta                     ! metadata structures
+  USE globalData,only:lookup_meta
   USE globalData,only:statForc_meta                           ! child metadata for stats
   USE globalData,only:statProg_meta                           ! child metadata for stats
   USE globalData,only:statDiag_meta                           ! child metadata for stats
@@ -132,9 +141,11 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
   USE globalData,only:statBvar_meta                           ! child metadata for stats
   USE globalData,only:gru_struc
   USE globalData,only:structInfo                              ! information on the data structures
-  USE multiconst,only:secprday               ! number of seconds in a day
+  USE multiconst,only:secprday                                ! number of seconds in a day
   USE data_types,only:file_info_array
-  USE var_lookup,only:maxvarFreq                ! maximum number of output files
+  USE var_lookup,only:maxvarFreq                              ! maximum number of output files
+
+  USE allocspace_module,only:allocGlobal                      ! module to allocate space for global data structures
   
   implicit none
   type(file_info_array),intent(in)      :: forcFileInfo
@@ -159,97 +170,92 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
   if (.not.allocated(outputStructure))then
     allocate(outputStructure(1))
   else
-    print*, "Already Allocated"
-    return;
+    print*, "Already Allocated"; return;
   end if
 
-  ! Statistics Structures
-  allocate(outputStructure(1)%forcStat(1))
-  allocate(outputStructure(1)%progStat(1))
-  allocate(outputStructure(1)%diagStat(1))
-  allocate(outputStructure(1)%fluxStat(1))
-  allocate(outputStructure(1)%indxStat(1))
-  allocate(outputStructure(1)%bvarStat(1))
-  allocate(outputStructure(1)%forcStat(1)%gru(num_gru))
-  allocate(outputStructure(1)%progStat(1)%gru(num_gru))
-  allocate(outputStructure(1)%diagStat(1)%gru(num_gru))
-  allocate(outputStructure(1)%fluxStat(1)%gru(num_gru))
-  allocate(outputStructure(1)%indxStat(1)%gru(num_gru))
-  allocate(outputStructure(1)%bvarStat(1)%gru(num_gru))
+  ! LookupStructure
 
+  ! Statistics Structures
+  allocate(outputStructure(1)%forcStat%gru(num_gru))
+  allocate(outputStructure(1)%progStat%gru(num_gru))
+  allocate(outputStructure(1)%diagStat%gru(num_gru))
+  allocate(outputStructure(1)%fluxStat%gru(num_gru))
+  allocate(outputStructure(1)%indxStat%gru(num_gru))
+  allocate(outputStructure(1)%bvarStat%gru(num_gru))
   ! Primary Data Structures (scalars)
-  allocate(outputStructure(1)%timeStruct(1))
-  allocate(outputStructure(1)%forcStruct(1))
-  allocate(outputStructure(1)%attrStruct(1))
-  allocate(outputStructure(1)%typeStruct(1))
-  allocate(outputStructure(1)%idStruct(1))
-  allocate(outputStructure(1)%timeStruct(1)%gru(num_gru))
-  allocate(outputStructure(1)%forcStruct(1)%gru(num_gru))
-  allocate(outputStructure(1)%attrStruct(1)%gru(num_gru))
-  allocate(outputStructure(1)%typeStruct(1)%gru(num_gru))
-  allocate(outputStructure(1)%idStruct(1)%gru(num_gru))
-  
+  allocate(outputStructure(1)%timeStruct%gru(num_gru))
+  allocate(outputStructure(1)%forcStruct%gru(num_gru))
   ! Primary Data Structures (variable length vectors)
-  allocate(outputStructure(1)%indxStruct(1))
-  allocate(outputStructure(1)%mparStruct(1))
-  allocate(outputStructure(1)%progStruct(1))
-  allocate(outputStructure(1)%diagStruct(1))
-  allocate(outputStructure(1)%fluxStruct(1))
-  allocate(outputStructure(1)%indxStruct(1)%gru(num_gru))
-  allocate(outputStructure(1)%mparStruct(1)%gru(num_gru))
-  allocate(outputStructure(1)%progStruct(1)%gru(num_gru))
-  allocate(outputStructure(1)%diagStruct(1)%gru(num_gru))
-  allocate(outputStructure(1)%fluxStruct(1)%gru(num_gru))
-
+  allocate(outputStructure(1)%indxStruct%gru(num_gru))
+  allocate(outputStructure(1)%progStruct%gru(num_gru))
+  allocate(outputStructure(1)%diagStruct%gru(num_gru))
+  allocate(outputStructure(1)%fluxStruct%gru(num_gru))
   ! Basin-Average structures
-  allocate(outputStructure(1)%bparStruct(1))
-  allocate(outputStructure(1)%bvarStruct(1))
-  allocate(outputStructure(1)%bparStruct(1)%gru(num_gru))
-  allocate(outputStructure(1)%bvarStruct(1)%gru(num_gru))
-
-  ! define the ancillary data structures
-  allocate(outputStructure(1)%dparStruct(1))
-  allocate(outputStructure(1)%dparStruct(1)%gru(num_gru))
-
+  allocate(outputStructure(1)%bvarStruct%gru(num_gru))
   ! Finalize Stats for writing
-  allocate(outputStructure(1)%finalizeStats(1))
-  allocate(outputStructure(1)%finalizeStats(1)%gru(num_gru))
+  allocate(outputStructure(1)%finalizeStats%gru(num_gru))
+  ! Extras
+  allocate(outputStructure(1)%upArea%gru(num_gru))
   
   
   do iGRU = 1, num_gru
     num_hru = gru_struc(iGRU)%hruCount
     ! Statistics Structures
-    allocate(outputStructure(1)%forcStat(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%progStat(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%diagStat(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%fluxStat(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%indxStat(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%bvarStat(1)%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%forcStat%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%progStat%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%diagStat%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%fluxStat%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%indxStat%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%bvarStat%gru(iGRU)%hru(num_hru))
 
     ! Primary Data Structures (scalars)
-    allocate(outputStructure(1)%timeStruct(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%forcStruct(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%attrStruct(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%typeStruct(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%idStruct(1)%gru(iGRU)%hru(num_hru))
-  
+    allocate(outputStructure(1)%timeStruct%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%forcStruct%gru(iGRU)%hru(num_hru))
+
     ! Primary Data Structures (variable length vectors)
-    allocate(outputStructure(1)%indxStruct(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%mparStruct(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%progStruct(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%diagStruct(1)%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%fluxStruct(1)%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%indxStruct%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%progStruct%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%diagStruct%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%fluxStruct%gru(iGRU)%hru(num_hru))
   
     ! Basin-Average structures
-    allocate(outputStructure(1)%bvarStruct(1)%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%bvarStruct%gru(iGRU)%hru(num_hru))
+
 
    ! define the ancillary data structures
-    allocate(outputStructure(1)%dparStruct(1)%gru(iGRU)%hru(num_hru))
 
     ! Finalize Stats for writing
-    allocate(outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(num_hru))
+    allocate(outputStructure(1)%finalizeStats%gru(iGRU)%hru(num_hru))
+
+    allocate(outputStructure(1)%upArea%gru(iGRU)%hru(num_hru))
+
+  end do
+
 
+  ! Allocate variables that do not require time
+  do iStruct=1,size(structInfo)
+    select case(trim(structInfo(iStruct)%structName))
+    case('time'); cycle;
+    case('forc'); cycle;
+    case('attr'); call allocGlobal(attr_meta,  outputStructure(1)%attrStruct,  err, message)
+    case('type'); call allocGlobal(type_meta,  outputStructure(1)%typeStruct,  err, message)
+    case('id'  ); call allocGlobal(id_meta,    outputStructure(1)%idStruct,    err, message)
+    case('mpar'); call allocGlobal(mpar_meta,  outputStructure(1)%mparStruct,  err, message); 
+    case('indx'); call allocGlobal(indx_meta,  outputStructure(1)%indxStruct_init, err, message);
+    case('prog'); call allocGlobal(prog_meta,  outputStructure(1)%progStruct_init, err, message);
+    case('diag'); cycle;
+    case('flux'); cycle;
+    case('bpar'); call allocGlobal(bpar_meta,outputStructure(1)%bparStruct    ,err, message);  ! basin-average params 
+    case('bvar'); call allocGlobal(bvar_meta,outputStructure(1)%bvarStruct_init,err,message);  ! basin-average variables
+    case('deriv'); cycle;
+    case('lookup'); call allocGlobal(lookup_meta,outputStructure(1)%lookupStruct,err, message);
+    end select
   end do
+  
+
+  call allocGlobal(mpar_meta,outputStructure(1)%dparStruct,err,message);
+
+
 
   do iGRU=1,num_gru
     do iHRU=1,gru_struc(iGRU)%hruCount
@@ -263,70 +269,58 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
         ! allocate space structures
           select case(trim(structInfo(iStruct)%structName))    
             case('time')
-              call alloc_outputStruc(time_meta,outputStructure(1)%timeStruct(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(time_meta,outputStructure(1)%timeStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,err=err,message=message)     ! model forcing data
             case('forc')
               ! Structure
-              call alloc_outputStruc(forc_meta,outputStructure(1)%forcStruct(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(forc_meta,outputStructure(1)%forcStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model forcing data
               ! Statistics
-              call alloc_outputStruc(statForc_meta(:)%var_info,outputStructure(1)%forcStat(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(statForc_meta(:)%var_info,outputStructure(1)%forcStat%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model forcing data
-            case('attr')
-              call alloc_outputStruc(attr_meta,outputStructure(1)%attrStruct(1)%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! local attributes for each HRU
-            case('type')
-              call alloc_outputStruc(type_meta,outputStructure(1)%typeStruct(1)%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! classification of soil veg etc.
-            case('id'  )
-              call alloc_outputStruc(id_meta,outputStructure(1)%idStruct(1)%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);        ! local values of hru gru IDs
-            case('mpar') ! model parameters
-              call alloc_outputStruc(mpar_meta,outputStructure(1)%mparStruct(1)%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message); 
-
-              call alloc_outputStruc(mpar_meta, outputStructure(1)%dparStruct(1)%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,err=err,message=message)
+            case('attr'); cycle;
+              ! call allocGlobal(attr_meta, outputStructure(1)%attrStruct(1)%gru(iGRU)%hru(iHRU), err, message)
+            case('type'); cycle;
+            case('id'  ); cycle;
+            case('mpar'); cycle;
             case('indx')
               ! Structure
-              call alloc_outputStruc(indx_meta,outputStructure(1)%indxStruct(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(indx_meta,outputStructure(1)%indxStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,str_name='indx',message=message);    ! model variables
               ! Statistics
-              call alloc_outputStruc(statIndx_meta(:)%var_info,outputStructure(1)%indxStat(1)%gru(iGRU)%hru(1), &
+              call alloc_outputStruc(statIndx_meta(:)%var_info,outputStructure(1)%indxStat%gru(iGRU)%hru(1), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! index vars
             case('prog')
               ! Structure
-              call alloc_outputStruc(prog_meta,outputStructure(1)%progStruct(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(prog_meta,outputStructure(1)%progStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model prognostic (state) variables
               ! Statistics
-              call alloc_outputStruc(statProg_meta(:)%var_info,outputStructure(1)%progStat(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(statProg_meta(:)%var_info,outputStructure(1)%progStat%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model prognostic 
             case('diag')
               ! Structure
-              call alloc_outputStruc(diag_meta,outputStructure(1)%diagStruct(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(diag_meta,outputStructure(1)%diagStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model diagnostic variables
               ! Statistics
-              call alloc_outputStruc(statDiag_meta(:)%var_info,outputStructure(1)%diagStat(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(statDiag_meta(:)%var_info,outputStructure(1)%diagStat%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model diagnostic
             case('flux')
               ! Structure
-              call alloc_outputStruc(flux_meta,outputStructure(1)%fluxStruct(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(flux_meta,outputStructure(1)%fluxStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model fluxes
               ! Statistics
-              call alloc_outputStruc(statFlux_meta(:)%var_info,outputStructure(1)%fluxStat(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(statFlux_meta(:)%var_info,outputStructure(1)%fluxStat%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model fluxes
-            case('bpar')
-              call alloc_outputStruc(bpar_meta,outputStructure(1)%bparStruct(1)%gru(iGRU), &
-                                      nSteps=maxSteps,nSnow=0,nSoil=0,err=err,message=message);  ! basin-average params 
+            case('bpar'); cycle;
             case('bvar')
               ! Structure
-              call alloc_outputStruc(bvar_meta,outputStructure(1)%bvarStruct(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(bvar_meta,outputStructure(1)%bvarStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=0,nSoil=0,err=err,message=message);  ! basin-average variables
               ! Statistics
-              call alloc_outputStruc(statBvar_meta(:)%var_info,outputStructure(1)%bvarStat(1)%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(statBvar_meta(:)%var_info,outputStructure(1)%bvarStat%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=0,nSoil=0,err=err,message=message);  ! basin-average variables
             case('deriv');  cycle
-            case('lookup'); cycle
+            case('lookup'); cycle  ! lookup tables
             case default; err=20; message='unable to find structure name: '//trim(structInfo(iStruct)%structName)
         end select
 
@@ -339,14 +333,15 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
       end do  ! looping through data structures
     
       ! Finalize stats structure for writing to output file
-      allocate(outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(iHRU)%tim(maxSteps))
+      allocate(outputStructure(1)%finalizeStats%gru(iGRU)%hru(iHRU)%tim(maxSteps))
       do iStep = 1, maxSteps
-        allocate(outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(iHRU)%tim(iStep)%dat(1:maxVarFreq))
+        allocate(outputStructure(1)%finalizeStats%gru(iGRU)%hru(iHRU)%tim(iStep)%dat(1:maxVarFreq))
       end do ! timeSteps
     end do ! Looping through GRUs
   end do
 
 
+
 end subroutine initOutputStructure
 
 subroutine deallocateOutputStructure(err) bind(C, name="deallocateOutputStructure")
diff --git a/build/source/actors/file_access_actor/fortran_code/read_attrb.f90 b/build/source/actors/file_access_actor/fortran_code/read_attrb.f90
deleted file mode 100644
index 7c58fabe9f92c8bfe986aac057a2f6d6a6dd25c3..0000000000000000000000000000000000000000
--- a/build/source/actors/file_access_actor/fortran_code/read_attrb.f90
+++ /dev/null
@@ -1,356 +0,0 @@
-module read_attrb_module
-USE, intrinsic :: iso_c_binding
-USE nrtype
-
-implicit none
-private
-public::allocateAttributeStructures
-public::openAttributeFile
-public::getNumVarAttr
-public::closeAttributeFile
-public::readAttributeFromNetCDF
-
-contains
-
-subroutine allocateAttributeStructures(index_gru, index_hru, & ! indexes into gru_struc
-    handle_attr_struct, handle_type_struct, handle_id_struct, err) bind(C, name="allocateAttributeStructures")
-  USE data_types,only:var_d, var_i, var_i8
-  USE globalData,only:gru_struc
-  USE globalData,only:attr_meta,type_meta,id_meta
-  USE allocspace_module,only:allocLocal
-  implicit none
-  integer(c_int),intent(in)           :: index_gru
-  integer(c_int),intent(in)           :: index_hru
-  type(c_ptr), intent(in), value      :: handle_attr_struct
-  type(c_ptr), intent(in), value      :: handle_type_struct
-  type(c_ptr), intent(in), value      :: handle_id_struct
-  integer(c_int), intent(out)         :: err
-  type(var_d), pointer                :: attr_struct
-  type(var_i), pointer                :: type_struct
-  type(var_i8), pointer               :: id_struct
-  integer(i4b)                        :: nSoil
-  integer(i4b)                        :: nSnow
-  character(len=256)                  :: message
-  ! ---------------------------------------------------------------------------------------
-  ! * Convert From C++ to Fortran
-  ! ---------------------------------------------------------------------------------------
-  call c_f_pointer(handle_attr_struct, attr_struct)
-  call c_f_pointer(handle_type_struct, type_struct)
-  call c_f_pointer(handle_id_struct,   id_struct)
-  ! Start subroutine
-  err=0; message="read_attrb.f90 - allocateAttributeStructures"
-
-  nSnow = gru_struc(index_gru)%hruInfo(index_hru)%nSnow
-  nSoil = gru_struc(index_gru)%hruInfo(index_hru)%nSoil
-
-  call allocLocal(attr_meta,attr_struct,nSnow,nSoil,err,message);
-  call allocLocal(type_meta,type_struct,nSnow,nSoil,err,message);
-  call allocLocal(id_meta,id_struct,nSnow,nSoil,err,message);
-  if(err/=0)then; message=trim(message); print*, message; return; endif;
-  
-end subroutine allocateAttributeStructures
-
-
-subroutine openAttributeFile(attr_ncid, err) bind(C, name="openAttributeFile")
-  USE netcdf
-  USE netcdf_util_module,only:nc_file_open                   ! open netcdf file
-  ! Attribute File
-  USE summaFileManager,only:SETTINGS_PATH                     ! define path to settings files (e.g., parameters, soil and veg. tables)
-  USE summaFileManager,only:LOCAL_ATTRIBUTES                  ! name of model initial attributes file
-  implicit none
-  integer(c_int),intent(out)                :: attr_ncid
-  integer(c_int),intent(out)                :: err
-
-  ! local variables
-  character(len=256)                        :: message       ! error message
-  character(len=256)                        :: attrFile           ! attributes file name
-
-  err=0; message="read_attrb.f90 - openAttributesFile"
-  attrFile = trim(SETTINGS_PATH)//trim(LOCAL_ATTRIBUTES)
-
-  call nc_file_open(trim(attrFile),nf90_noWrite,attr_ncid,err,message)
-  if(err/=0)then
-      message=trim(message)
-      print*, message
-      return
-  endif
-
-end subroutine
-
-subroutine getNumVarAttr(attr_ncid, num_var, err) bind(C, name="getNumVarAttr")
-  USE netcdf
-  USE netcdf_util_module,only:netcdf_err                     ! netcdf error handling function
-  implicit none
-  integer(c_int),intent(in)                 :: attr_ncid
-  integer(c_int),intent(out)                :: num_var
-  integer(c_int),intent(out)                :: err
-
-  ! local variables
-  character(len=256)                        :: message       ! error message
-  err=0; message="read_attrb.f90 - getNumVar"
-  ! get number of variables total in netcdf file
-  err = nf90_inquire(attr_ncid,nvariables=num_var)
-  call netcdf_err(err,message)
-  if (err/=0) then
-      message=trim(message)//'problem with nf90_inquire'
-      return
-  endif
-
-end subroutine getNumVarAttr
-
-subroutine closeAttributeFile(attr_ncid, err) bind(C, name="closeAttributeFile")
-  USE netcdf_util_module,only:nc_file_close
-  implicit none
-  integer(c_int),intent(in)         :: attr_ncid
-  integer(c_int),intent(out)        :: err
-  ! local variables
-  character(len=256)                :: message
-  err=0; message="read_attrb.f90 - closeAttributeFile"
-
-  call nc_file_close(attr_ncid,err,message)
-  if (err/=0)then
-    message=trim(message)
-    return
-  end if
-
-end subroutine closeAttributeFile
-
-
-
-! Read in the local attributes for an HRU
-subroutine readAttributeFromNetCDF(ncid, index_gru, index_hru, num_var, &
-  handle_attr_struct, handle_type_struct, handle_id_struct, err) bind(C, name="readAttributeFromNetCDF")
-  ! netcdf utilities
-  USE netcdf
-  USE netcdf_util_module,only:nc_file_open                   ! open netcdf file
-  USE netcdf_util_module,only:nc_file_close                  ! close netcdf file
-  USE nr_utility_module ,only:arth
-  ! needed global data structures and metadata
-  USE globalData,only:gru_struc                              ! gru-hru mapping structure
-  USE globalData,only:attr_meta,type_meta,id_meta            ! metadata structures
-  USE get_ixname_module,only:get_ixAttr,get_ixType,get_ixId  ! access function to find index of elements in structure
-  ! Information to make up the attributes file
-  USE summaFileManager,only:SETTINGS_PATH                     ! define path to settings files (e.g., parameters, soil and veg. tables)
-  USE summaFileManager,only:LOCAL_ATTRIBUTES                  ! name of model initial attributes file
-  ! Fortran Data Type Structures
-  USE data_types,only:var_d, var_i, var_i8
-  implicit none
-  ! indexes into gru_struc
-  integer(c_int), intent(in)            :: ncid
-  integer(c_int), intent(in)            :: index_gru
-  integer(c_int), intent(in)            :: index_hru
-  ! number of variables from the netCDF file
-  integer(c_int), intent(in)            :: num_var
-  ! data structures to populate
-  type(c_ptr), intent(in), value        :: handle_attr_struct
-  type(c_ptr), intent(in), value        :: handle_type_struct
-  type(c_ptr), intent(in), value        :: handle_id_struct
-  ! error control
-  integer(c_int), intent(out)           :: err
-  ! local variables
-  integer(i4b)                          :: iVar               ! loop through varibles in the netcdf file
-  integer(i4b)                          :: varType            ! type of variable (categorica, numerical, idrelated)
-  integer(i4b)                          :: varIndx            ! index of variable within its data structure
-  ! Fortran structures
-  type(var_d), pointer                  :: attr_struct
-  type(var_i), pointer                  :: type_struct
-  type(var_i8), pointer                 :: id_struct
-  ! check structures - to verify input
-  integer(i4b)                          :: iCheck             ! index of an attribute name
-  logical(lgt),allocatable              :: checkType(:)       ! vector to check if we have all desired categorical values
-  logical(lgt),allocatable              :: checkId(:)         ! vector to check if we have all desired IDs
-  logical(lgt),allocatable              :: checkAttr(:)       ! vector to check if we have all desired local attributes
-  ! netcdf variables
-  character(LEN=nf90_max_name)          :: varName            ! character array of netcdf variable name
-  integer(i4b),parameter                :: categorical=101    ! named variable to denote categorical data
-  integer(i4b),parameter                :: numerical=102      ! named variable to denote numerical data
-  integer(i4b),parameter                :: idrelated=103      ! named variable to denote ID related data
-  integer(i4b)                          :: categorical_var(1) ! temporary categorical variable from local attributes netcdf file
-  real(rkind)                           :: numeric_var(1)     ! temporary numeric variable from local attributes netcdf file
-  integer(8)                            :: idrelated_var(1)   ! temporary ID related variable from local attributes netcdf file
-  character(len=256)                    :: attr_file          ! attributes file name
-  character(len=256)                    :: message           
-  ! ---------------------------------------------------------------------------------------
-  ! * Convert From C++ to Fortran
-  ! ---------------------------------------------------------------------------------------
-  call c_f_pointer(handle_attr_struct, attr_struct)
-  call c_f_pointer(handle_type_struct, type_struct)
-  call c_f_pointer(handle_id_struct,   id_struct)
-
-  err=0; message="read_attrb_file_access_actor - read_attrb.f90"
-
-  attr_file= trim(SETTINGS_PATH)//trim(LOCAL_ATTRIBUTES)
-
-  ! **********************************************************************************************
-  ! (1) prepare check vectors
-  ! **********************************************************************************************
-  allocate(checkType(size(type_meta)),checkAttr(size(attr_meta)),checkId(size(id_meta)),stat=err)
-  if(err/=0)then
-      err=20
-      message=trim(message)//'problem allocating space for variable check vectors'
-      print*, message
-      return
-  endif
-  checkType(:) = .false.
-  checkAttr(:) = .false.
-  checkId(:)   = .false.
-
-  iCheck = 1
-  do iVar = 1,num_var
-    ! inqure about current variable name, type, number of dimensions
-    err = nf90_inquire_variable(ncID,iVar,name=varName)
-    if(err/=nf90_noerr)then; 
-        message=trim(message)//'problem inquiring variable: '//trim(varName)//'/'//trim(nf90_strerror(err)); 
-        print*, message
-        return 
-    endif
-    
-    ! find attribute name
-    select case(trim(varName))
-   
-      ! ** categorical data
-      case('vegTypeIndex','soilTypeIndex','slopeTypeIndex','downHRUindex')
-
-          ! get the index of the variable
-          varType = categorical
-          varIndx = get_ixType(varName)
-          checkType(varIndx) = .true.
-
-          ! check that the variable could be identified in the data structure
-          if(varIndx < 1)then
-              err=20; 
-              message=trim(message)//'unable to find variable ['//trim(varName)//'] in data structure'; 
-              print*, message
-              return; 
-          endif
-  
-
-          err = nf90_get_var(ncID,iVar,categorical_var,start=(/gru_struc(index_gru)%hruInfo(index_hru)%hru_nc/),count=(/1/))
-          if(err/=nf90_noerr)then
-              message=trim(message)//'problem reading: '//trim(varName)
-              print*, message
-              return
-          end if
-          type_struct%var(varIndx) = categorical_var(1)
-
-          ! ** ID related data
-      case('hruId')
-          ! get the index of the variable
-          varType = idrelated
-          varIndx = get_ixId(varName)
-          checkId(varIndx) = .true.
-  
-          ! check that the variable could be identified in the data structure
-          if(varIndx < 1)then
-              err=20
-              message=trim(message)//'unable to find variable ['//trim(varName)//'] in data structure'
-              print*, message
-              return
-          endif
-  
-          ! get data from netcdf file and store in vector
-
-          err = nf90_get_var(ncID,iVar,idrelated_var,start=(/gru_struc(index_gru)%hruInfo(index_hru)%hru_nc/),count=(/1/))
-          if(err/=nf90_noerr)then
-              message=trim(message)//'problem reading: '//trim(varName)
-              print*, message
-              return
-          end if
-          id_struct%var(varIndx) = idrelated_var(1)
-
-
-      ! ** numerical data
-      case('latitude','longitude','elevation','tan_slope','contourLength','HRUarea','mHeight')
-
-          ! get the index of the variable
-          varType = numerical
-          varIndx = get_ixAttr(varName)
-          checkAttr(varIndx) = .true.
-  
-          ! check that the variable could be identified in the data structure
-          if(varIndx < 1)then
-              err=20; message=trim(message)//'unable to find variable ['//trim(varName)//'] in data structure'
-              print*, message
-              return 
-          endif
-          ! get data from netcdf file and store in vector
-                 
-          err = nf90_get_var(ncID,iVar,numeric_var,start=(/gru_struc(index_gru)%hruInfo(index_hru)%hru_nc/),count=(/1/))
-          if(err/=nf90_noerr)then
-              message=trim(message)//'problem reading: '//trim(varName)
-              print*, message
-              return
-          end if
-          attr_struct%var(varIndx) = numeric_var(1)
-    
-      
-      ! for mapping varibles, do nothing (information read above)
-      case('hru2gruId','gruId'); cycle
-  
-      ! check that variables are what we expect
-      case default
-          message=trim(message)//'unknown variable ['//trim(varName)//'] in local attributes file'
-          print*,message
-          err=20
-          return
-
-    end select ! select variable
-
-  end do ! (looping through netcdf local attribute file)
-
-  ! ** now handle the optional aspect variable if it's missing
-  varIndx = get_ixAttr('aspect')
-  ! check that the variable was not found in the attribute file
-  if(.not. checkAttr(varIndx)) then
-      if (index_gru == 1) then
-        write(*,*) NEW_LINE('A')//'INFO: aspect not found in the input attribute file, continuing ...'//NEW_LINE('A')
-      endif
-      attr_struct%var(varIndx) = nr_realMissing      ! populate variable with out-of-range value, used later
-      checkAttr(varIndx) = .true.
-  endif
-
-  ! **********************************************************************************************
-  ! (4) check that we have all the desired varaibles
-  ! **********************************************************************************************
-  ! check that we have all desired categorical variables
-  if(any(.not.checkType))then
-    do iCheck = 1,size(type_meta)
-        if(.not.checkType(iCheck))then 
-            err=20
-            message=trim(message)//'missing variable ['//trim(type_meta(iCheck)%varname)//'] in local attributes file'
-            print*, message
-            return
-        endif
-    end do
-  endif
-
-  ! check that we have all desired ID variables
-  if(any(.not.checkId))then
-    do iCheck = 1,size(id_meta)
-        if(.not.checkId(iCheck))then
-            err=20
-            message=trim(message)//'missing variable ['//trim(id_meta(iCheck)%varname)//'] in local attributes file'
-            print*, message
-            return 
-        endif
-    end do
-  endif
-
-  ! check that we have all desired local attributes
-  if(any(.not.checkAttr))then
-    do iCheck = 1,size(attr_meta)
-        if(.not.checkAttr(iCheck))then; 
-            err=20
-            message=trim(message)//'missing variable ['//trim(attr_meta(iCheck)%varname)//'] in local attributes file'
-            print*, message
-            return 
-        endif
-    end do
-  endif
-
-  ! free memory
-  deallocate(checkType)
-  deallocate(checkId)
-  deallocate(checkAttr)
-end subroutine readAttributeFromNetCDF
-end module read_attrb_module
diff --git a/build/source/actors/file_access_actor/fortran_code/read_icondFromStructure.f90 b/build/source/actors/file_access_actor/fortran_code/read_icondFromStructure.f90
deleted file mode 100644
index 6eabadd6706558e8b8c79cdee7978d90cdb1edc5..0000000000000000000000000000000000000000
--- a/build/source/actors/file_access_actor/fortran_code/read_icondFromStructure.f90
+++ /dev/null
@@ -1,280 +0,0 @@
-module read_icondFromStructure_module
-USE, intrinsic :: iso_c_binding
-USE nrtype
-
-
-implicit none
-private
-public :: openInitCondFile
-public :: closeInitCondFile
-public :: readInitCond_prog
-public :: readInitCond_bvar
-
-! define single HRU restart file
-integer(i4b), parameter :: singleHRU=1001
-integer(i4b), parameter :: multiHRU=1002
-integer(i4b), parameter :: restartFileType=multiHRU
-contains
-
-subroutine openInitCondFile(init_cond_ncid, err) bind(C, name="openInitCondFile")
-  USE netcdf_util_module,only:nc_file_open               ! open netcdf file
-  USE netcdf
-  ! file paths
-  USE summaFileManager,only:SETTINGS_PATH           ! path to settings files (e.g., Noah vegetation tables)
-  USE summaFileManager,only:STATE_PATH              ! optional path to state/init. condition files (defaults to SETTINGS_PATH)
-  USE summaFileManager,only:MODEL_INITCOND          ! name of model initial conditions file
-  implicit none
-
-  integer(c_int),intent(out)          :: init_cond_ncid
-  integer(c_int),intent(out)          :: err
-
-  character(len=256)                  :: init_cond_fileName
-  character(len=256)                  :: message          ! error message
-  character(len=1024)                 :: cmessage         ! error message for downwind routine
-  ! --------------------------------------------------------------------------------------------------------
-
-  ! define restart file path/name
-  if(STATE_PATH == '') then
-    init_cond_fileName = trim(SETTINGS_PATH)//trim(MODEL_INITCOND)
-  else
-    init_cond_fileName = trim(STATE_PATH)//trim(MODEL_INITCOND)
-  endif
-
-  call nc_file_open(init_cond_fileName,nf90_nowrite,init_cond_ncid,err,cmessage)
-  if (err/=0) then
-    message=trim(message)//trim(cmessage)
-    print(message)
-    return
-  end if
-
-end subroutine openInitCondFile
-
-subroutine closeInitCondFile(init_cond_ncid,err) bind(C, name="closeInitCondFile")
-  USE netcdf_util_module,only:nc_file_close
-  implicit none
-  
-  integer(c_int),intent(out)          :: init_cond_ncid
-  integer(c_int),intent(out)          :: err
-  
-  ! local variables
-  character(len=256)                  :: message
-  ! --------------------------------------------------------------------------------------------------------
-  err=0; message="read_icondFromStructure.f90 - closeInitCondFile/"
-
-  call nc_file_close(init_cond_ncid,err,message)
-  if (err/=0) then
-    message=trim(message)
-    print(message)
-    return
-  end if
-
-end subroutine closeInitCondFile
-
-subroutine readInitCond_prog(init_cond_ncid, start_gru, num_gru, err) bind(C, name="readInitCond_prog")
-  ! Structures to populate
-  USE globalData,only:init_cond_prog
-  ! Netcdf
-  USE netcdf
-  USE netcdf_util_module,only:nc_file_close              ! close netcdf file
-  USE netcdf_util_module,only:netcdf_err                 ! netcdf error handling
-
-  ! metadata
-  USE globalData,only:prog_meta                     ! metadata for prognostic variables
-  USE globalData,only:bvar_meta                     ! metadata for basin (GRU) variables
-  ! var_lookup
-  USE var_lookup,only:iLookVarType                  ! variable lookup structure
-  USE var_lookup,only:iLookPROG                     ! variable lookup structure
-  USE var_lookup,only:iLookPARAM                    ! variable lookup structure
-  USE var_lookup,only:iLookINDEX                    ! variable lookup structure
-  ! access type string
-  USE get_ixName_module,only:get_varTypeName        ! to access type strings for error messages
-  implicit none
-
-  integer(c_int), intent(in)             :: init_cond_ncid
-  integer(c_int), intent(in)             :: start_gru
-  integer(c_int), intent(in)             :: num_gru
-  integer(c_int), intent(out)            :: err
- 
-  ! local variables
-  integer(i4b)                           :: iVar
-  integer(i4b)                           :: dimID        ! varible dimension ids
-  integer(i4b)                           :: ncVarID      ! variable ID in netcdf file
-  integer(i4b)                           :: dimLen       ! data dimensions
-  character(256)                         :: dimName      ! not used except as a placeholder in call to inq_dim function
-  integer(i4b)                           :: fileHRU      ! number of HRUs in file
-
-  character(LEN=256)                     :: icond_file  ! restart/icond_file file name
-  character(len=256)                     :: message
-  character(len=256)                     :: cmessage
-
-  character(len=32),parameter            :: scalDimName   ='scalarv'  ! dimension name for scalar data
-  character(len=32),parameter            :: midSoilDimName='midSoil'  ! dimension name for soil-only layers
-  character(len=32),parameter            :: midTotoDimName='midToto'  ! dimension name for layered varaiables
-  character(len=32),parameter            :: ifcTotoDimName='ifcToto'  ! dimension name for layered varaiables
-  ! --------------------------------------------------------------------------------------------------------
-
-  err=0; message="read_icondFromStructure.f90 - readInitCond_prog"
-  
-  ! get number of HRUs in file
-  err = nf90_inq_dimid(init_cond_ncid,"hru",dimID);
-  if(err/=nf90_noerr)then; message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err)); return; end if
-  err = nf90_inquire_dimension(init_cond_ncid,dimID,len=fileHRU);
-  if(err/=nf90_noerr)then; message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err)); return; end if
-
-  allocate(init_cond_prog(size(prog_meta)))
-
-  ! loop through prognostic variables
-  do iVar = 1,size(prog_meta)
-    ! skip variables that are computed later
-    if(prog_meta(iVar)%varName=='scalarCanopyWat'           .or. &
-       prog_meta(iVar)%varName=='spectralSnowAlbedoDiffuse' .or. &
-       prog_meta(iVar)%varName=='scalarSurfaceTemp'         .or. &
-       prog_meta(iVar)%varName=='mLayerVolFracWat'          .or. &
-       prog_meta(iVar)%varName=='mLayerHeight'                   )then 
-       cycle
-    endif
-
-    ! get variable id
-    err = nf90_inq_varid(init_cond_ncid,trim(prog_meta(iVar)%varName),ncVarID); call netcdf_err(err,message)
-    if(err/=0)then
-      message=trim(message)//': problem with getting variable id, var='//trim(prog_meta(iVar)%varName)
-      print*,message
-      return
-    endif
-
-    ! get variable dimension IDs
-    select case (prog_meta(iVar)%varType)
-      case (iLookVarType%scalarv); err = nf90_inq_dimid(init_cond_ncid,trim(scalDimName)   ,dimID)
-        call netcdf_err(err,message)
-      case (iLookVarType%midSoil); err = nf90_inq_dimid(init_cond_ncid,trim(midSoilDimName),dimID)
-        call netcdf_err(err,message)
-      case (iLookVarType%midToto); err = nf90_inq_dimid(init_cond_ncid,trim(midTotoDimName),dimID)
-        call netcdf_err(err,message)
-      case (iLookVarType%ifcToto); err = nf90_inq_dimid(init_cond_ncid,trim(ifcTotoDimName),dimID)
-        call netcdf_err(err,message)
-    case default
-      message=trim(message)//"unexpectedVariableType[name='"//trim(prog_meta(iVar)%varName)//"';type='"//trim(get_varTypeName(prog_meta(iVar)%varType))//"']"
-      print*, message
-      err=20; return
-    end select
-
-    ! check errors
-    if(err/=0)then
-      message=trim(message)//': problem with dimension ids, var='//trim(prog_meta(iVar)%varName)
-      print*, message
-      return
-    endif
-
-    ! get the dimension length
-    err = nf90_inquire_dimension(init_cond_ncid,dimID,dimName,dimLen); call netcdf_err(err,message)
-    if(err/=0)then;message=trim(message)//': problem getting the dimension length';print*, message;return;endif
-
-    ! initialize the variable data
-    allocate(init_cond_prog(iVar)%var_data(num_gru,dimLen),stat=err)
-    if(err/=0)then;message=trim(message)//'problem allocating HRU variable data';print*, message;return;endif
-
-     ! get data
-    err = nf90_get_var(init_cond_ncid,ncVarID,init_cond_prog(iVar)%var_data, start=(/start_gru,1/),count=(/num_gru,dimLen/)) 
-    call netcdf_err(err,message)
-    if(err/=0)then; message=trim(message)//': problem getting the data for variable '//trim(prog_meta(iVar)%varName); return; endif
-
-  end do
-
-end subroutine readInitCond_prog
-
-subroutine readInitCond_bvar(init_cond_ncid, start_gru, num_gru, err) bind(C, name="readInitCond_bvar")
-  USE globalData,only:init_cond_bvar
-  USE globalData,only: nTimeDelay   ! number of hours in the time delay histogram
-  ! var_lookup
-  USE var_lookup,only:iLookBVAR                     ! variable lookup structure
-  ! metadata structures
-  USE globalData,only:bvar_meta                     ! metadata for basin (GRU) variables
-  ! netcdf
-  USE netcdf
-  USE netcdf_util_module,only:netcdf_err                 ! netcdf error handling
-  implicit none
-  
-  integer(c_int), intent(in)             :: init_cond_ncid
-  integer(c_int), intent(in)             :: start_gru
-  integer(c_int), intent(in)             :: num_gru
-  integer(c_int), intent(out)            :: err
-
-  ! local variables
-  integer(i4b)                           :: nTDH          ! number of points in time-delay histogram
-  integer(i4b)                           :: dimID        ! varible dimension ids
-  integer(i4b)                           :: fileGRU      ! number of GRUs in file
-  integer(i4b)                           :: i          
-  integer(i4b)                           :: iVar     
-  integer(i4b)                           :: dimLen       ! data dimensions
-  character(256)                         :: dimName      ! not used except as a placeholder in call to inq_dim function
-  integer(i4b)                           :: ncVarID      ! variable ID in netcdf file
-  integer(i4b),dimension(1)              :: ndx          ! intermediate array of loop indices
-   
-
-  character(len=256)                     :: message
-  
-  character(len=32),parameter            :: tdhDimName    ='tdh'      ! dimension name for time-delay basin variables
-
-  ! --------------------------------------------------------------------------------------------------------
-  err = 0; message="read_icondFromStructure.f90 - readInitCond_bvar/"
-  if(restartFileType/=singleHRU)then
-    ! get dimension of time delay histogram (TDH) from initial conditions file
-    err = nf90_inq_dimid(init_cond_ncid,"tdh",dimID);
-    
-    if(err/=nf90_noerr)then
-      write(*,*) 'WARNING: routingRunoffFuture is not in the initial conditions file ... using zeros'  ! previously created in var_derive.f90
-      err=nf90_noerr    ! reset this err
-    
-    else 
-      ! the state file *does* have the basin variable(s), so process them
-      err = nf90_inquire_dimension(init_cond_ncid,dimID,len=nTDH);
-      if(err/=nf90_noerr)then
-        message=trim(message)//'problem reading tdh dimension from initial condition file/'//trim(nf90_strerror(err))
-        print*, message
-        return
-      end if
-      
-      ! get number of GRUs in file
-      err = nf90_inq_dimid(init_cond_ncid,"gru",dimID);               if(err/=nf90_noerr)then; message=trim(message)//'problem finding gru dimension/'//trim(nf90_strerror(err)); return; end if
-      err = nf90_inquire_dimension(init_cond_ncid,dimID,len=fileGRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading gru dimension/'//trim(nf90_strerror(err)); return; end if
-
-      ! check vs hardwired value set in globalData.f90
-      if(nTDH /= nTimeDelay)then
-        write(*,*) 'tdh=',nTDH,' nTimeDelay=',nTimeDelay
-        message=trim(message)//': state file time delay dimension tdh does not match summa expectation of nTimeDelay set in globalData()'
-        return
-      endif
-
-      ndx = (/iLookBVAR%routingRunoffFuture/)   ! array of desired variable indices
-      allocate(init_cond_bvar(size(ndx)))
-      do i = 1, size(ndx)
-        iVar = ndx(i)
-        ! get tdh dimension Id in file (should be 'tdh')
-        err = nf90_inq_dimid(init_cond_ncid,trim(tdhDimName), dimID);
-        if(err/=0)then;message=trim(message)//': problem with dimension ids for tdh vars';print*,message;return;endif
-
-        ! get the tdh dimension length (dimName and dimLen are outputs of this call)
-        err = nf90_inquire_dimension(init_cond_ncid,dimID,dimName,dimLen); call netcdf_err(err,message)
-        if(err/=0)then;message=trim(message)//': problem getting the dimension length for tdh vars';print*,message;return;endif
-
-        ! get tdh-based variable id
-        err = nf90_inq_varid(init_cond_ncid,trim(bvar_meta(iVar)%varName),ncVarID); call netcdf_err(err,message)
-        if(err/=0)then; message=trim(message)//': problem with getting basin variable id, var='//trim(bvar_meta(iVar)%varName); return; endif
-
-        allocate(init_cond_bvar(i)%var_data(num_gru,dimLen),stat=err)
-        if(err/=0)then; print*, 'err= ',err; message=trim(message)//'problem allocating GRU variable data'; return; endif
-
-        ! get data
-        err = nf90_get_var(init_cond_ncid,ncVarID,init_cond_bvar(i)%var_data, start=(/start_gru,1/),count=(/num_gru,dimLen/)); call netcdf_err(err,message)
-        if(err/=0)then; message=trim(message)//': problem getting the data'; return; endif
-      end do
-    endif
-  endif
-
-
-end subroutine readInitCond_bvar
-
-
-
-
-end module read_icondFromStructure_module
\ No newline at end of file
diff --git a/build/source/actors/file_access_actor/fortran_code/read_param.f90 b/build/source/actors/file_access_actor/fortran_code/read_param.f90
deleted file mode 100644
index e8d8c6784f1caff7919c2d2f1021747431b74713..0000000000000000000000000000000000000000
--- a/build/source/actors/file_access_actor/fortran_code/read_param.f90
+++ /dev/null
@@ -1,401 +0,0 @@
-module read_param_module
-  USE, intrinsic :: iso_c_binding
-  USE nrtype
-  implicit none
-  private
-  public::allocateParamStructures
-  public::openParamFile
-  public::getNumVarParam
-  public::closeParamFile
-  public::getParamSizes
-  public::overwriteParam
-  public::readParamFromNetCDF
-  contains
-subroutine allocateParamStructures(index_gru, index_hru, handle_dpar_struct, &
-    handle_mpar_struct, handle_bpar_struct, err) bind(C, name="allocateParamStructures")
-  
-  USE globalData,only:mpar_meta,bpar_meta
-  USE globalData,only:gru_struc
-  USE data_types,only:var_dlength,var_d
-  USE allocspace_module,only:allocLocal
-
-  implicit none
-  integer(c_int),intent(in)       :: index_gru
-  integer(c_int),intent(in)       :: index_hru
-  type(c_ptr),intent(in),value    :: handle_dpar_struct
-  type(c_ptr),intent(in),value    :: handle_mpar_struct
-  type(c_ptr),intent(in),value    :: handle_bpar_struct
-  integer(c_int),intent(out)      :: err
-
-  type(var_d), pointer            :: dpar_struct
-  type(var_dlength), pointer      :: mpar_struct
-  type(var_d), pointer            :: bpar_struct
-
-  integer(i4b)                    :: nSnow  
-  integer(i4b)                    :: nSoil
-
-  character(len=256)              :: message
-
-  ! ---------------------------------------------------------------------------------------
-  ! * Convert From C++ to Fortran
-  ! ---------------------------------------------------------------------------------------
-  call c_f_pointer(handle_dpar_struct, dpar_struct)
-  call c_f_pointer(handle_mpar_struct, mpar_struct)
-  call c_f_pointer(handle_bpar_struct, bpar_struct)
-  ! start of subroutine
-  err=0; message="read_attrb.f90 - allocateAttributeStructures"
-
-  nSnow = gru_struc(index_gru)%hruInfo(index_hru)%nSnow
-  nSoil = gru_struc(index_gru)%hruInfo(index_hru)%nSoil
-
-  ! initalize the structure with allocatable components
-  call allocLocal(mpar_meta,dpar_struct,nSnow,nSoil,err,message); 
-  call allocLocal(mpar_meta,mpar_struct,nSnow,nSoil,err,message); 
-  call allocLocal(bpar_meta,bpar_struct,nSnow=0,nSoil=0,err=err,message=message); 
-  if(err/=0)then; message=trim(message); print*, message; return; endif;
-end subroutine
-
-subroutine openParamFile(param_ncid, param_file_exists, err) bind(C, name="openParamFile")
-  USE netcdf
-  USE netcdf_util_module,only:nc_file_open
-  ! Parameters File
-  USE summaFileManager,only:SETTINGS_PATH             ! path for metadata files
-  USE summaFileManager,only:PARAMETER_TRIAL           ! file with parameter trial values
-  
-  implicit none
-  integer(c_int),intent(out)          :: param_ncid
-  logical(c_bool),intent(out)         :: param_file_exists
-  integer(c_int),intent(out)          :: err
-  
-  ! local variables
-  character(LEN=1024)                 :: infile           ! input filename
-  character(len=256)                  :: message          ! error message
-  character(len=1024)                 :: cmessage         ! error message for downwind routine
-  err=0; message="read_param.f90 - openParamFile"    
-  ! **********************************************************************************************
-  ! * open files, etc.
-  ! **********************************************************************************************
-  
-  infile = trim(SETTINGS_PATH)//trim(PARAMETER_TRIAL) ! build filename
-
-  ! check whether the user-specified file exists and warn if it does not
-  inquire(file=trim(infile),exist=param_file_exists)
-  if (.not.param_file_exists) then
-      write(*,'(A)') NEW_LINE('A')//'!! WARNING:  trial parameter file not found; proceeding instead with other default parameters; check path in file manager input if this was not the desired behavior'//NEW_LINE('A')
-      return
-  endif
-
-  ! open trial parameters file if it exists
-  call nc_file_open(trim(infile),nf90_nowrite,param_ncid,err,cmessage)
-  if(err/=0)then
-    message=trim(message)//trim(cmessage)
-    print*, message
-    return
-  end if
-
-end subroutine openParamFile
-
-subroutine getNumVarParam(param_ncid, num_var, err) bind(C, name="getNumVarParam")
-  USE netcdf
-  USE netcdf_util_module,only:netcdf_err                     ! netcdf error handling function
-  implicit none
-  integer(c_int),intent(in)            :: param_ncid
-  integer(c_int),intent(out)           :: num_var
-  integer(c_int),intent(out)           :: err
-
-  ! local variables
-  character(len=256)                   :: message       ! error message
-  
-  err=0; message="read_param.f90 - getNumVarAndDims/"
-  ! get the number of variables in the parameter file
-  err=nf90_inquire(param_ncid, nVariables=num_var)
-  call netcdf_err(err,message)
-  if (err/=0) then
-    err=20
-    print*, message
-    return
-  end if
-
-end subroutine getNumVarParam
-
-subroutine closeParamFile(param_ncid, err) bind(C, name="closeParamFile")
-  USE netcdf_util_module,only:nc_file_close
-  implicit none
-  integer(c_int),intent(in)            :: param_ncid
-  integer(c_int),intent(out)           :: err
-  ! local variables
-  character(len=256)                   :: message
-  ! --------------------------------------------------------------------------------------------------------
- 
-  err=0; message="read_param.f90 - closeParamFile/"
-
-  call nc_file_close(param_ncid,err,message)
-  if(err/=0)then
-    message=trim(message)
-    print*, message
-    return
-  end if
-
-end subroutine closeParamFile
-
-! get the sizes of the arrays for dpar_array bpar_array
-subroutine getParamSizes(dpar_array_size, bpar_array_size, type_array_size) bind(C, name="getParamSizes")
-  USE var_lookup,only:maxvarMpar      ! model parameters: maximum number variables
-  USE var_lookup,only:maxvarBpar      ! model parameters: maximum number variables
-  USE var_lookup,only:maxvarType
-
-  implicit none
-  integer(c_int),intent(out)    :: dpar_array_size
-  integer(c_int),intent(out)    :: bpar_array_size
-  integer(c_int),intent(out)    :: type_array_size
-
-
-  dpar_array_size = maxvarMpar
-  bpar_array_size = maxvarBpar
-  type_array_size = maxvarType
-
-
-end subroutine getParamSizes
-
-subroutine overwriteParam(index_gru, index_hru, handle_type_struct, &
-  handle_dpar_struct, handle_mpar_struct, handle_bpar_struct, err) bind(C, name="overwriteParam")
-  USE var_lookup,only:maxvarMpar      ! model parameters: maximum number variables
-  USE var_lookup,only:maxvarBpar      ! model parameters: maximum number variables
-  USE var_lookup,only:iLookTYPE      ! named variables to index elements of the data vectors
-  ! global data
-  USE globalData,only:gru_struc
-  USE globalData,only:localParFallback                        ! local column default parameters
-  USE globalData,only:basinParFallback                        ! basin-average default parameter
-  USE data_types,only:var_dlength,var_i,var_d
-
-  USE pOverwrite_module,only:pOverwrite                       ! module to overwrite default parameter values with info from the Noah tables
-  USE allocspace_module,only:allocLocal
-  implicit none
-  integer(c_int),intent(in)     :: index_gru
-  integer(c_int),intent(in)     :: index_hru
-  ! structures
-  type(c_ptr),intent(in),value  :: handle_type_struct
-  type(c_ptr),intent(in),value  :: handle_dpar_struct
-  type(c_ptr),intent(in),value  :: handle_mpar_struct
-  type(c_ptr),intent(in),value  :: handle_bpar_struct
-
-  ! error control
-  integer(c_int), intent(out)   :: err
-
-  ! local variables
-  type(var_i),pointer           :: type_struct                 !  model parameters
-  type(var_d),pointer           :: dpar_struct                 !  model parameters
-  type(var_dlength),pointer     :: mpar_struct                 !  model parameters
-  type(var_d),pointer           :: bpar_struct                 !  model parameters
-
-  integer(i4b)                  :: iVar
-  integer(i4b)                  :: iDat
-
-  character(len=256)            :: message
-  ! ---------------------------------------------------------------------------------------
-  ! * Convert From C++ to Fortran
-  ! ---------------------------------------------------------------------------------------
-  call c_f_pointer(handle_type_struct, type_struct)
-  call c_f_pointer(handle_dpar_struct, dpar_struct)
-  call c_f_pointer(handle_mpar_struct, mpar_struct)
-  call c_f_pointer(handle_bpar_struct, bpar_struct)
-  ! Start subroutine
-  err=0; message="read_param.f90 - overwriteParam"
-
-  ! Set the basin parameters with the default values
-  do ivar=1, size(localParFallback)
-    dpar_struct%var(iVar) = localParFallback(iVar)%default_val
-  end do
-
-  call pOverwrite(type_struct%var(iLookTYPE%vegTypeIndex), &  ! vegetation category
-                  type_struct%var(iLookTYPE%soilTypeIndex),&  ! soil category
-                  dpar_struct%var(:),err,message)             ! default model parameters
-  
-  do ivar=1, size(localParFallback)
-    do iDat=1, size(mpar_struct%var(iVar)%dat)
-      mpar_struct%var(iVar)%dat(iDat) = dpar_struct%var(iVar)
-    end do
-  end do
-
-  do iVar=1, size(basinParFallback)
-    bpar_struct%var(iVar) = basinParFallback(iVar)%default_val
-  end do
-
-end subroutine overwriteParam
-
-
-subroutine readParamFromNetCDF(param_ncid, index_gru, index_hru, start_index_gru, &
-    num_vars, handle_mpar_struct, handle_bpar_struct, err) bind(C, name="readParamFromNetCDF")
-  USE netcdf
-  USE netcdf_util_module,only:netcdf_err     ! netcdf error handling function
-
-  USE data_types,only:var_dlength,var_d
-  USE get_ixname_module,only:get_ixparam,get_ixbpar   ! access function to find index of elements in structure
-  
-  USE globalData,only:index_map,gru_struc             ! mapping from global HRUs to the elements in the data structures
-  USE globalData,only:integerMissing  ! missing integer
-  
-  implicit none
-  ! dummy variables
-  integer(c_int),intent(in)     :: param_ncid
-  integer(c_int),intent(in)     :: index_gru
-  integer(c_int),intent(in)     :: index_hru
-  integer(c_int),intent(in)     :: start_index_gru
-  integer(c_int),intent(in)     :: num_vars
-  type(c_ptr),intent(in),value  :: handle_mpar_struct
-  type(c_ptr),intent(in),value  :: handle_bpar_struct
-  integer(c_int), intent(out)   :: err
-  ! define local variables
-  type(var_dlength),pointer     :: mpar_struct                 !  model parameters
-  type(var_d),pointer           :: bpar_struct                 !  model parameters
-
-  character(len=256)            :: message          ! error message
-  character(len=1024)           :: cmessage         ! error message for downwind routine
-  integer(i4b)                  :: localHRU_ix      ! index of HRU within data structure
-  integer(i4b)                  :: ixParam          ! index of the model parameter in the data structure
-  
-  ! indices/metadata in the NetCDF file
-  integer(i4b)                  :: num_dims            ! number of dimensions
-
-  integer(i4b)                  :: ivarid           ! variable index
-  character(LEN=64)             :: dimName          ! dimension name
-
-  character(LEN=64)             :: parName          ! parameter name
-  integer(i4b)                  :: nSoil_file       ! number of soil layers in the file
-  integer(i4b)                  :: idim_list(2)     ! list of dimension ids
-  ! data in the netcdf file
-  integer(i4b)                  :: parLength        ! length of the parameter data
-  real(dp),allocatable          :: parVector(:)     ! model parameter vector
-  integer(i4b)                  :: fHRU             ! index of HRU in input file
-  integer(i4b)                  :: netcdf_index
-
-  ! ---------------------------------------------------------------------------------------
-  ! * Convert From C++ to Fortran
-  ! ---------------------------------------------------------------------------------------
-  call c_f_pointer(handle_mpar_struct, mpar_struct)
-  call c_f_pointer(handle_bpar_struct, bpar_struct)
-  err=0; message="read_param.f90 - readParamFromNetCDF/"
-
-
-  ! **********************************************************************************************
-  ! * read the local parameters and the basin parameters
-  ! ********************************************************************************************** 
-  do ivarid=1,num_vars
-    ! get the parameter name
-    err=nf90_inquire_variable(param_ncid, ivarid, name=parName)
-    call netcdf_err(err,message)
-    if (err/=0) then
-      err=20
-      print*,message
-      return
-    end if
-
-    ! get the local parameters
-    ixParam = get_ixparam( trim(parName) )
-    if(ixParam/=integerMissing)then
-      ! **********************************************************************************************
-      ! * read the local parameters
-      ! **********************************************************************************************
-      ! get the variable shape
-      err=nf90_inquire_variable(param_ncid, ivarid, nDims=num_dims, dimids=idim_list)
-      if(err/=0)then
-          message=trim(message)//trim(cmessage)
-          print*, message
-          return
-      end if
-        
-        ! get the length of the depth dimension (if it exists)
-        if(num_dims==2)then
-            ! get the information on the 2nd dimension for 2-d variables
-            err=nf90_inquire_dimension(param_ncid, idim_list(2), dimName, nSoil_file)
-            if(err/=0)then
-                message=trim(message)//trim(cmessage)
-                print*, message
-                return
-            end if
-            
-            ! check that it is the depth dimension
-            if(trim(dimName)/='depth')then
-                message=trim(message)//'expect 2nd dimension of 2-d variable to be depth (dimension name = '//trim(dimName)//')'
-                err=20; return
-            endif
-
-            ! ! check that the dimension length is correct
-            if(size(mpar_struct%var(ixParam)%dat) /= nSoil_file)then
-                message=trim(message)//'unexpected number of soil layers in parameter file'
-                err=20; return
-            endif
-
-            ! define parameter length
-            parLength = nSoil_file
-        else
-            parLength = 1
-        endif ! if two dimensions
-
-        ! allocate space for model parameters
-        allocate(parVector(parLength),stat=err)
-        if(err/=0)then
-            message=trim(message)//'problem allocating space for parameter vector'
-            err=20; return
-        endif
-        
-
-        localHRU_ix=index_map(index_hru)%localHRU_ix
-        fHRU = gru_struc(index_gru)%hruInfo(localHRU_ix)%hru_nc
-        ! read parameter data
-        select case(num_dims)
-            case(1); err=nf90_get_var(param_ncid, ivarid, parVector, start=(/fHRU/), count=(/1/) )
-            case(2); err=nf90_get_var(param_ncid, ivarid, parVector, start=(/fHRU,1/), count=(/1,nSoil_file/) )
-            case default; err=20; message=trim(message)//'unexpected number of dimensions for parameter '//trim(parName)
-        end select
-          
-        ! error check for the parameter read
-        if(err/=0)then
-            message=trim(message)//trim(cmessage)
-            print*, message
-            return
-        end if
-          
-          ! populate parameter structures
-        select case(num_dims)
-            case(1); mpar_struct%var(ixParam)%dat(:) = parVector(1)  ! also distributes scalar across depth dimension
-            case(2); mpar_struct%var(ixParam)%dat(:) = parVector(:)
-            case default; err=20; message=trim(message)//'unexpected number of dimensions for parameter '//trim(parName)
-        end select
-
-        ! deallocate space for model parameters
-        deallocate(parVector,stat=err)
-        if(err/=0)then
-            message=trim(message)//'problem deallocating space for parameter vector'
-            print*, message
-            err=20; return
-        endif
-    
-    ! **********************************************************************************************
-    ! * read the basin parameters
-    ! **********************************************************************************************
-    
-    ! get the basin parameters
-    else
-        ! get the parameter index
-        ixParam = get_ixbpar( trim(parName) )
-
-        ! allow extra variables in the file that are not used
-        if(ixParam==integerMissing) cycle
-
-        ! read parameter data
-        netcdf_index = start_index_gru + index_gru - 1
-        err=nf90_get_var(param_ncid, ivarid, bpar_struct%var(ixParam), start=(/netcdf_index/))
-        if(err/=0)then
-          message=trim(message)//trim(cmessage)
-          print*, message
-          return
-        end if
-    endif
-
-  end do ! (looping through the parameters in the NetCDF file)
-
-end subroutine readParamFromNetCDF
-
-
-end module read_param_module
\ No newline at end of file
diff --git a/build/source/actors/file_access_actor/fortran_code/writeOutputFromOutputStructure.f90 b/build/source/actors/file_access_actor/fortran_code/writeOutputFromOutputStructure.f90
index 3e159d536b9c98bbd3cc8cd19deec82b6e82a657..364d14537077dde693c61e89a9237a351175c6ee 100644
--- a/build/source/actors/file_access_actor/fortran_code/writeOutputFromOutputStructure.f90
+++ b/build/source/actors/file_access_actor/fortran_code/writeOutputFromOutputStructure.f90
@@ -126,7 +126,7 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, err)
     stepCounter(:) = outputTimeStep(iGRU)%dat(:) ! We want to avoid updating outputTimeStep
     do iStep=1, num_steps
       call writeTime(ncid,outputTimeStep(iGRU)%dat(:),iStep,time_meta, &
-              outputStructure(1)%timeStruct(1)%gru(iGRU)%hru(indxHRU)%var,err,cmessage)
+              outputStructure(1)%timeStruct%gru(iGRU)%hru(indxHRU)%var,err,cmessage)
     end do ! istep
   end do ! iGRU
 
@@ -137,7 +137,7 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, err)
   ! ****************************************************************************
   call writeBasin(ncid,outputTimeStep(start_gru)%dat(:),outputTimeStepUpdate,num_steps,&
                   start_gru, max_gru, numGRU, &
-                  bvar_meta,outputStructure(1)%bvarStat(1),outputStructure(1)%bvarStruct(1), &
+                  bvar_meta,outputStructure(1)%bvarStat,outputStructure(1)%bvarStruct, &
                   bvarChild_map,err,cmessage)
 
   ! ****************************************************************************
@@ -148,28 +148,28 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, err)
       case('forc')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, & 
-                        forc_meta,outputStructure(1)%forcStat(1),outputStructure(1)%forcStruct(1),'forc', &
-                        forcChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
+                        forc_meta,outputStructure(1)%forcStat,outputStructure(1)%forcStruct,'forc', &
+                        forcChild_map,outputStructure(1)%indxStruct,err,cmessage)
       case('prog')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, &
-                        prog_meta,outputStructure(1)%progStat(1),outputStructure(1)%progStruct(1),'prog', &
-                        progChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
+                        prog_meta,outputStructure(1)%progStat,outputStructure(1)%progStruct,'prog', &
+                        progChild_map,outputStructure(1)%indxStruct,err,cmessage)
       case('diag')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, &
-                        diag_meta,outputStructure(1)%diagStat(1),outputStructure(1)%diagStruct(1),'diag', &
-                        diagChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
+                        diag_meta,outputStructure(1)%diagStat,outputStructure(1)%diagStruct,'diag', &
+                        diagChild_map,outputStructure(1)%indxStruct,err,cmessage)
       case('flux')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, &
-                        flux_meta,outputStructure(1)%fluxStat(1),outputStructure(1)%fluxStruct(1),'flux', &
-                        fluxChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
+                        flux_meta,outputStructure(1)%fluxStat,outputStructure(1)%fluxStruct,'flux', &
+                        fluxChild_map,outputStructure(1)%indxStruct,err,cmessage)
       case('indx')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, &
-                        indx_meta,outputStructure(1)%indxStat(1),outputStructure(1)%indxStruct(1),'indx', &
-                        indxChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
+                        indx_meta,outputStructure(1)%indxStat,outputStructure(1)%indxStruct,'indx', &
+                        indxChild_map,outputStructure(1)%indxStruct,err,cmessage)
     end select
     if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
   end do  ! (looping through structures)
@@ -325,9 +325,9 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,nSteps,
 
         do iStep = 1, nSteps
           ! check if we want this timestep
-          if(.not.outputStructure(1)%finalizeStats(1)%gru(verifiedGRUIndex)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+          if(.not.outputStructure(1)%finalizeStats%gru(verifiedGRUIndex)%hru(1)%tim(iStep)%dat(iFreq)) cycle
           stepCounter = stepCounter+1
-          timeVec(stepCounter) = outputStructure(1)%forcStruct(1)%gru(verifiedGRUIndex)%hru(1)%var(iVar)%tim(iStep)
+          timeVec(stepCounter) = outputStructure(1)%forcStruct%gru(verifiedGRUIndex)%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
@@ -396,7 +396,7 @@ subroutine writeScalar(ncid, outputTimestep, outputTimestepUpdate, nSteps, minGR
         stepCounter = 0
         gruCounter = gruCounter + 1
         do iStep = 1, nSteps
-          if(.not.outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+          if(.not.outputStructure(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
@@ -412,9 +412,8 @@ subroutine writeScalar(ncid, outputTimestep, outputTimestepUpdate, nSteps, minGR
         print*, "   maxGRU = ", maxGRU
         err = 20
         return
-
-      err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realVec(1:gruCounter, 1:stepCounter),start=(/minGRU,outputTimestep(iFreq)/),count=(/numGRU,stepCounter/))
       endif
+      err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realVec(1:gruCounter, 1:stepCounter),start=(/minGRU,outputTimestep(iFreq)/),count=(/numGRU,stepCounter/))
     class default; err=20; message=trim(message)//'stats must be scalarv and of type gru_hru_doubleVec'; return
   end select  ! stat
 
@@ -492,11 +491,11 @@ subroutine writeVector(ncid, outputTimestep, maxLayers, nSteps, minGRU, maxGRU,
       ! get the data vectors
       select type (dat)
           class is (gru_hru_time_doubleVec)
-              if(.not.outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+              if(.not.outputStructure(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(:)
 
           class is (gru_hru_time_intVec)
-              if(.not.outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+              if(.not.outputStructure(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(:)
           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
@@ -567,6 +566,7 @@ subroutine writeBasin(ncid,outputTimestep,outputTimestepUpdate,nSteps,&
   ! initialize error control
   err=0;message="f-writeBasin/"
 
+
   ! loop through output frequencies
   do iFreq=1,maxvarFreq
 
@@ -632,7 +632,7 @@ subroutine writeTime(ncid,outputTimestep,iStep,meta,dat,err,message)
   do iFreq=1,maxvarFreq
 
     ! check that we have finalized statistics for a given frequency
-    if(.not.outputStructure(1)%finalizeStats(1)%gru(1)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+    if(.not.outputStructure(1)%finalizeStats%gru(1)%hru(1)%tim(iStep)%dat(iFreq)) cycle
 
     ! loop through model variables
     do iVar = 1,size(meta)
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 959fd3e4118870ae21da8fc4f08dc02278596abd..4dc4668b296fbb897b3ddffcb270583e2f16e481 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
@@ -8,9 +8,6 @@ USE data_types
 
 implicit none
 public::writeParamToNetCDF
-public::writeDataToNetCDF
-public::writeBasinToNetCDF
-public::writeTimeToNetCDF
 public::writeGRUStatistics
 
 contains
@@ -199,83 +196,6 @@ subroutine writeDataToNetCDF(handle_ncid,           &
   end do  ! (looping through structures)
 end subroutine writeDataToNetCDF
 
-subroutine writeBasinToNetCDF(handle_ncid, index_gru, handle_finalize_stats, &
-  handle_output_timestep, handle_bvar_stat, handle_bvar_struct, err) bind(C, name="writeBasinToNetCDF")
-  USE modelwrite_module,only:writeBasin
-  USE globalData,only:bvar_meta                 ! metadata on basin-average variables
-  USE globalData,only:bvarChild_map             ! index of the child data structure: stats bvar
- 
-  implicit none
-  ! dummy variables
-  type(c_ptr),    intent(in), value  :: handle_ncid
-  integer(c_int), intent(in)         :: index_gru
-  type(c_ptr),    intent(in), value  :: handle_finalize_stats
-  type(c_ptr),    intent(in), value  :: handle_output_timestep
-  type(c_ptr),    intent(in), value  :: handle_bvar_stat
-  type(c_ptr),    intent(in), value  :: handle_bvar_struct
-  integer(c_int), intent(out)        :: err
-  ! local pointers for dummy variables
-  type(var_i), pointer               :: ncid
-  type(flagVec), pointer             :: finalize_stats
-  type(var_i),pointer                :: output_timestep
-  type(var_dlength),pointer          :: bvar_stat
-  type(var_dlength),pointer          :: bvar_struct
-  ! local Variables
-  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_output_timestep, output_timestep)
-  call c_f_pointer(handle_bvar_stat, bvar_stat)
-  call c_f_pointer(handle_bvar_struct, bvar_struct)
-  message="file_access_actor.f90 - writeBasinToNetCDF"
-
-
-  call writeBasin(ncid,index_gru,finalize_stats%dat(:),output_timestep%var(:),&
-      bvar_meta, bvar_stat%var, bvar_struct%var, bvarChild_map, err, message)
-  if(err/=0)then 
-    message=trim(message)//'[bvar]'
-    print*, message
-    return
-  endif 
-
-end subroutine writeBasinToNetCDF
-
-subroutine writeTimeToNetCDF(handle_ncid, handle_finalize_stats, handle_output_timestep, &
-    handle_time_struct, err) bind(C, name="writeTimeToNetCDF")
-  USE modelwrite_module,only:writeTime
-  USE globalData,only:time_meta
-
-  implicit none
-  type(c_ptr), intent(in), value :: handle_ncid
-  type(c_ptr), intent(in), value :: handle_finalize_stats
-  type(c_ptr), intent(in), value :: handle_output_timestep
-  type(c_ptr), intent(in), value :: handle_time_struct
-  integer(c_int), intent(out)    :: err
-
-  type(var_i), pointer           :: ncid
-  type(flagVec), pointer         :: finalize_stats
-  type(var_i),pointer            :: output_timestep
-  type(var_i),pointer            :: time_struct
-  character(len=256)             :: message
-
-  call c_f_pointer(handle_ncid, ncid)
-  call c_f_pointer(handle_finalize_stats, finalize_stats)
-  call c_f_pointer(handle_output_timestep, output_timestep)
-  call c_f_pointer(handle_time_struct, time_struct)
-  message="file_access_actor.f90 - writeTimeToNetCDF"
-
-  call writeTime(ncid, finalize_stats%dat(:),output_timestep%var(:),&
-      time_meta, time_struct%var,err,message)
-  if(err/=0)then 
-    message=trim(message)//'writeTime'
-    print*, message
-    return
-  endif 
-
-end subroutine writeTimeToNetCDF
 
 subroutine writeGRUStatistics(handle_ncid,      &
                               gru_var_ids,      &
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 28e9feaf11d440d13ff4f1bbadc856ccd96bec34..38a39d2d811126119a244e1a153cd4ddc5b21820 100644
--- a/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
+++ b/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
@@ -41,12 +41,18 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
             self->state.handle_indxStat, 
             self->state.handle_bvarStat, 
             self->state.handle_timeStruct, 
-            self->state.handle_forcStruct, 
+            self->state.handle_forcStruct,
+            self->state.handle_attrStruct,
+            self->state.handle_typeStruct,
+            self->state.handle_idStruct,
             self->state.handle_indxStruct,
+            self->state.handle_mparStruct,
             self->state.handle_progStruct, 
             self->state.handle_diagStruct, 
             self->state.handle_fluxStruct,
+            self->state.handle_bparStruct,
             self->state.handle_bvarStruct, 
+            self->state.handle_dparStruct,
             self->state.handle_startTime, 
             self->state.handle_finshTime, 
             self->state.handle_refTime,
@@ -76,30 +82,34 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
                   .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);
                   });
 
-    // Get the attributes and parameters for the HRU
-    self->request(self->state.file_access_actor,
-                  caf::infinite,
-                  get_attributes_params_v,
-                  self->state.indxGRU)
-                  .await([=](std::tuple<std::vector<double>,
-                                        std::vector<int>,
-                                        std::vector<long int>, 
-                                        std::vector<double>, 
-                                        std::vector<double>, 
-                                        std::vector<std::vector<double>>> attr_and_params) {
-                                int err = 0;
-                                set_var_d(std::get<0>(attr_and_params), self->state.handle_attrStruct);
-                                set_var_i(std::get<1>(attr_and_params), self->state.handle_typeStruct);
-                                set_var_i8(std::get<2>(attr_and_params), self->state.handle_idStruct);
-                                set_var_d(std::get<3>(attr_and_params), self->state.handle_bparStruct);
-                                set_var_d(std::get<4>(attr_and_params), self->state.handle_dparStruct);
-                                set_var_dlength(std::get<5>(attr_and_params), self->state.handle_mparStruct);
-
-                                Initialize_HRU(self);
-
-                                self->send(self, start_hru_v); });
+
+
+    // // Get the attributes and parameters for the HRU
+    // self->request(self->state.file_access_actor,
+    //               caf::infinite,
+    //               get_attributes_params_v,
+    //               self->state.indxGRU)
+    //               .await([=](std::tuple<std::vector<double>,
+    //                                     std::vector<int>,
+    //                                     std::vector<long int>, 
+    //                                     std::vector<double>, 
+    //                                     std::vector<double>, 
+    //                                     std::vector<std::vector<double>>> attr_and_params) {
+    //                             int err = 0;
+    //                             set_var_d(std::get<0>(attr_and_params), self->state.handle_attrStruct);
+    //                             set_var_i(std::get<1>(attr_and_params), self->state.handle_typeStruct);
+    //                             set_var_i8(std::get<2>(attr_and_params), self->state.handle_idStruct);
+    //                             set_var_d(std::get<3>(attr_and_params), self->state.handle_bparStruct);
+    //                             set_var_d(std::get<4>(attr_and_params), self->state.handle_dparStruct);
+    //                             set_var_dlength(std::get<5>(attr_and_params), self->state.handle_mparStruct);
+
+    //                             Initialize_HRU(self);
+
+    //                             self->send(self, start_hru_v); });
 
 
     return {
@@ -242,12 +252,14 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
 
 void Initialize_HRU(stateful_actor<hru_state>* self) {
 
-    setupHRUParam(&self->state.indxHRU, 
-                  &self->state.indxGRU,
+    setupHRUParam(&self->state.indxGRU,
+                  &self->state.indxHRU, 
                   self->state.handle_attrStruct, 
                   self->state.handle_typeStruct, 
                   self->state.handle_idStruct,
-                  self->state.handle_mparStruct, 
+                  self->state.handle_indxStruct,
+                  self->state.handle_mparStruct,
+                  self->state.handle_progStruct, 
                   self->state.handle_bparStruct, 
                   self->state.handle_bvarStruct,
                   self->state.handle_dparStruct, 
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 09a63f8e1788f85c25dce3632a954cbe44276333..f8ac820e1d66f229bd2b2ecd5d3ccad8e4a2f00e 100755
--- a/build/source/actors/hru_actor/fortran_code/hru_init.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_init.f90
@@ -55,13 +55,20 @@ contains
                         ! primary data structures (scalars)
                         handle_timeStruct,  & !  model time data
                         handle_forcStruct,  & !  model forcing data
+                        handle_attrStruct,  & !  model attribute data
+                        handle_typeStruct,  & !  model type data
+                        handle_idStruct,    & !  model id data
                         ! primary data structures (variable length vectors)
                         handle_indxStruct,  & !  model indices
+                        handle_mparStruct,  & !  model parameters
                         handle_progStruct,  & !  model prognostic (state) variables
                         handle_diagStruct,  & !  model diagnostic variables
                         handle_fluxStruct,  & !  model fluxes
                         ! basin-average structures
+                        handle_bparStruct,  & !  basin-average variables
                         handle_bvarStruct,  & !  basin-average variables
+                        ! ancillary data structures
+                        handle_dparStruct,  & !  default model parameters
                         ! local HRU data structures
                         handle_startTime,   & ! start time for the model simulation
                         handle_finshTime,   & ! end time for the model simulation
@@ -108,14 +115,20 @@ contains
   ! primary data structures (scalars)
   type(c_ptr), intent(in), value             :: handle_timeStruct !  model time data
   type(c_ptr), intent(in), value             :: handle_forcStruct !  model forcing data
+  type(c_ptr), intent(in), value             :: handle_attrStruct !  model attribute data
+  type(c_ptr), intent(in), value             :: handle_typeStruct !  model type data
+  type(c_ptr), intent(in), value             :: handle_idStruct !  model id data
   ! primary data structures (variable length vectors)
   type(c_ptr), intent(in), value             :: handle_indxStruct !  model indices
+  type(c_ptr), intent(in), value             :: handle_mparStruct !  model parameters
   type(c_ptr), intent(in), value             :: handle_progStruct !  model prognostic (state) variables
   type(c_ptr), intent(in), value             :: handle_diagStruct !  model diagnostic variables
   type(c_ptr), intent(in), value             :: handle_fluxStruct !  model fluxes
   ! basin-average structures
+  type(c_ptr), intent(in), value             :: handle_bparStruct !  basin-average variables
   type(c_ptr), intent(in), value             :: handle_bvarStruct !  basin-average variables
   ! ancillary data structures
+  type(c_ptr), intent(in), value             :: handle_dparStruct !  ancillary data structures 
   ! local hru data structures
   type(c_ptr), intent(in), value             :: handle_startTime  ! start time for the model simulation
   type(c_ptr), intent(in), value             :: handle_finshTime ! end time for the model simulation
@@ -125,23 +138,29 @@ contains
   ! ---------------------------------------------------------------------------------------
   ! * Fortran Variables For Conversion
   ! ---------------------------------------------------------------------------------------
-  type(zLookup),pointer                      :: lookupStruct               !  z(:)%var(:)%lookup(:) -- lookup tables
-  type(var_dlength),pointer                  :: forcStat                   !  model forcing data
-  type(var_dlength),pointer                  :: progStat                   !  model prognostic (state) variables
-  type(var_dlength),pointer                  :: diagStat                   !  model diagnostic variables
-  type(var_dlength),pointer                  :: fluxStat                   !  model fluxes
-  type(var_dlength),pointer                  :: indxStat                   !  model indices
-  type(var_dlength),pointer                  :: bvarStat                   !  basin-average variabl
+  type(zLookup),pointer                      :: lookupStruct               ! z(:)%var(:)%lookup(:) -- lookup tables
+  type(var_dlength),pointer                  :: forcStat                   ! model forcing data
+  type(var_dlength),pointer                  :: progStat                   ! model prognostic (state) variables
+  type(var_dlength),pointer                  :: diagStat                   ! model diagnostic variables
+  type(var_dlength),pointer                  :: fluxStat                   ! model fluxes
+  type(var_dlength),pointer                  :: indxStat                   ! model indices
+  type(var_dlength),pointer                  :: bvarStat                   ! basin-average variabl
   ! primary data structures (scalars)
-  type(var_i),pointer                        :: timeStruct                 !  model time data
-  type(var_d),pointer                        :: forcStruct                 !  model forcing data
+  type(var_i),pointer                        :: timeStruct                 ! model time data
+  type(var_d),pointer                        :: forcStruct                 ! model forcing data
+  type(var_d),pointer                        :: attrStruct                 ! model attribute data
+  type(var_i),pointer                        :: typeStruct                 ! model type data
+  type(var_i8),pointer                       :: idStruct                   ! model id data
   ! primary data structures (variable length vectors)
-  type(var_ilength),pointer                  :: indxStruct                 !  model indices
-  type(var_dlength),pointer                  :: progStruct                 !  model prognostic (state) variables
-  type(var_dlength),pointer                  :: diagStruct                 !  model diagnostic variables
-  type(var_dlength),pointer                  :: fluxStruct                 !  model fluxes
+  type(var_ilength),pointer                  :: indxStruct                 ! model indices
+  type(var_dlength),pointer                  :: mparStruct                 ! model parameters
+  type(var_dlength),pointer                  :: progStruct                 ! model prognostic (state) variables
+  type(var_dlength),pointer                  :: diagStruct                 ! model diagnostic variables
+  type(var_dlength),pointer                  :: fluxStruct                 ! model fluxes
   ! basin-average structures
-  type(var_dlength),pointer                  :: bvarStruct                 !  basin-average variables
+  type(var_d),pointer                        :: bparStruct                 ! basin-average variables
+  type(var_dlength),pointer                  :: bvarStruct                 ! basin-average variables
+  type(var_d),pointer                        :: dparStruct                 ! default model parameters
   ! local HRU data structures
   type(var_i),pointer                        :: startTime_hru              ! start time for the model simulation
   type(var_i),pointer                        :: finishTime_hru             ! end time for the model simulation
@@ -165,11 +184,17 @@ contains
   call c_f_pointer(handle_bvarStat,   bvarStat)
   call c_f_pointer(handle_timeStruct, timeStruct)
   call c_f_pointer(handle_forcStruct, forcStruct)
+  call c_f_pointer(handle_attrStruct, attrStruct)
+  call c_f_pointer(handle_typeStruct, typeStruct)
+  call c_f_pointer(handle_idStruct,   idStruct)
   call c_f_pointer(handle_indxStruct, indxStruct)
+  call c_f_pointer(handle_mparStruct, mparStruct)
   call c_f_pointer(handle_progStruct, progStruct)
   call c_f_pointer(handle_diagStruct, diagStruct)
   call c_f_pointer(handle_fluxStruct, fluxStruct)
+  call c_f_pointer(handle_bparStruct, bparStruct)
   call c_f_pointer(handle_bvarStruct, bvarStruct)
+  call c_f_pointer(handle_dparStruct, dparStruct)
   call c_f_pointer(handle_startTime,  startTime_hru)
   call c_f_pointer(handle_finshTime,  finishTime_hru)
   call c_f_pointer(handle_refTime,    refTime_hru)
@@ -221,17 +246,17 @@ contains
   do iStruct=1,size(structInfo)
   ! allocate space  
   select case(trim(structInfo(iStruct)%structName))    
-    case('time'); call allocLocal(time_meta,timeStruct,err=err,message=cmessage)     ! model forcing data
+    case('time'); call allocLocal(time_meta,timeStruct,err=err,message=cmessage)     ! model time data
     case('forc'); call allocLocal(forc_meta,forcStruct,nSnow,nSoil,err,cmessage);    ! model forcing data
-    case('attr'); cycle ! set by file_access_actor  
-    case('type'); cycle ! set by file_access_actor
-    case('id'  ); cycle ! set by file_access_actor   
-    case('mpar'); cycle ! set by file_access_actor  
+    case('attr'); call allocLocal(attr_meta,attrStruct,nSnow,nSoil,err,cmessage);    ! model attribute data
+    case('type'); call allocLocal(type_meta,typeStruct,nSnow,nSoil,err,cmessage);    ! model type data
+    case('id'  ); call allocLocal(id_meta,idStruct,nSnow,nSoil,err,cmessage);        ! model id data
+    case('mpar'); call allocLocal(mpar_meta,mparStruct,nSnow,nSoil,err,cmessage);    ! model parameters  
     case('indx'); call allocLocal(indx_meta,indxStruct,nSnow,nSoil,err,cmessage);    ! model variables
     case('prog'); call allocLocal(prog_meta,progStruct,nSnow,nSoil,err,cmessage);    ! model prognostic (state) variables
     case('diag'); call allocLocal(diag_meta,diagStruct,nSnow,nSoil,err,cmessage);    ! model diagnostic variables
     case('flux'); call allocLocal(flux_meta,fluxStruct,nSnow,nSoil,err,cmessage);    ! model fluxes
-    case('bpar'); cycle ! set by file_access_actor
+    case('bpar'); call allocLocal(bpar_meta,bparStruct,nSnow=0,nSoil=0,err=err,message=cmessage);  ! basin-average variables
     case('bvar'); call allocLocal(bvar_meta,bvarStruct,nSnow=0,nSoil=0,err=err,message=cmessage);  ! basin-average variables
     case('lookup'); cycle ! allocated in t2enthalpy.f90
     case('deriv'); cycle
@@ -245,6 +270,12 @@ contains
   endif
   end do  ! looping through data structures
 
+  ! allocate space for default model parameters
+	! NOTE: This is done here, rather than in the loop above, because dpar is not one of the "standard" data structures
+	call allocLocal(mpar_meta,dparStruct,nSnow,nSoil,err,cmessage);    ! default model parameters
+	if(err/=0)then; message=trim(message)//trim(cmessage)//' [problem allocating dparStruct]'; print*,message;return;endif
+	 
+
 
   ! *****************************************************************************
   ! *** allocate space for output statistics data structures
diff --git a/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90 b/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90
index a71f284c1aa12e48119083f928956e95ba743c38..96c95688dc159e3527bf898de4f190044ca7278d 100755
--- a/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90
@@ -116,17 +116,17 @@ subroutine writeParm(indxGRU,indxHRU,ispatial,struct,meta,structName,err,message
       select type (struct)
         class is (var_i)
         if (structName == "type")then
-          outputStructure(1)%typeStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar) = struct%var(iVar)
+          outputStructure(1)%typeStruct%gru(indxGRU)%hru(indxHRU)%var(iVar) = struct%var(iVar)
         end if
         class is (var_i8)
         
         class is (var_d)
         if (structName == "attr")then
-          outputStructure(1)%attrStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar) = struct%var(iVar)
+          outputStructure(1)%attrStruct%gru(indxGRU)%hru(indxHRU)%var(iVar) = struct%var(iVar)
         end if
         class is (var_dlength)
         if (structName == "mpar")then
-          outputStructure(1)%mparStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar) = struct%var(iVar)
+          outputStructure(1)%mparStruct%gru(indxGRU)%hru(indxHRU)%var(iVar) = struct%var(iVar)
         end if
         
         class default; err=20; message=trim(message)//'unknown variable type (with HRU)'; return
@@ -138,7 +138,7 @@ subroutine writeParm(indxGRU,indxHRU,ispatial,struct,meta,structName,err,message
       select type (struct)
         class is (var_d)
         if (structName == "bpar")then
-          outputStructure(1)%bparStruct(1)%gru(indxGRU)%var(iVar) = struct%var(iVar) ! this will overwrite data
+          outputStructure(1)%bparStruct%gru(indxGRU)%var(iVar) = struct%var(iVar) ! this will overwrite data
           print*, "bpar"
         end if
         class is (var_i8)
@@ -213,7 +213,7 @@ subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
         ! Write the time step values
         select type(dat)      ! forcStruc
           class is (var_d)    ! x%var(:)
-            outputStructure(1)%forcStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep) = dat%var(iVar)
+            outputStructure(1)%forcStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep) = dat%var(iVar)
           class default; err=20; message=trim(message)//'time variable must be of type var_d (forcing data structure)'; return
         end select
       end if  ! id time
@@ -230,15 +230,15 @@ subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
           class is (var_dlength)
             select case(trim(structName))
             case('forc')
-              outputStructure(1)%forcStat(1)%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
+              outputStructure(1)%forcStat%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
             case('prog')
-              outputStructure(1)%progStat(1)%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
+              outputStructure(1)%progStat%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
             case('diag')
-              outputStructure(1)%diagStat(1)%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
+              outputStructure(1)%diagStat%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
             case('flux')
-              outputStructure(1)%fluxStat(1)%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
+              outputStructure(1)%fluxStat%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
             case('indx')
-              outputStructure(1)%indxStat(1)%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
+              outputStructure(1)%indxStat%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
             case default
               err=21; message=trim(message)//"Stats structure not found"; return
             end select
@@ -250,11 +250,11 @@ subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
 
         ! get the model layers
         nSoil   = indx%var(iLookIndex%nSoil)%dat(1)
-        outputStructure(1)%indxStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iLookIndex%nSoil)%tim(iStep)%dat(1) = nSoil
+        outputStructure(1)%indxStruct%gru(indxGRU)%hru(indxHRU)%var(iLookIndex%nSoil)%tim(iStep)%dat(1) = nSoil
         nSnow   = indx%var(iLookIndex%nSnow)%dat(1)
-        outputStructure(1)%indxStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iLookIndex%nSnow)%tim(iStep)%dat(1) = nSnow
+        outputStructure(1)%indxStruct%gru(indxGRU)%hru(indxHRU)%var(iLookIndex%nSnow)%tim(iStep)%dat(1) = nSnow
         nLayers = indx%var(iLookIndex%nLayers)%dat(1)
-        outputStructure(1)%indxStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iLookIndex%nLayers)%tim(iStep)%dat(1) = nLayers
+        outputStructure(1)%indxStruct%gru(indxGRU)%hru(indxHRU)%var(iLookIndex%nLayers)%tim(iStep)%dat(1) = nLayers
 
         ! get the length of each data vector
         select case (meta(iVar)%varType)
@@ -273,16 +273,16 @@ subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
           class is (var_dlength)
             select case(trim(structName))
               case('prog')
-                outputStructure(1)%progStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
+                outputStructure(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
               case('diag')
-                outputStructure(1)%diagStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
+                outputStructure(1)%diagStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
               case('flux')
-                outputStructure(1)%fluxStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
+                outputStructure(1)%fluxStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
               case default
                 err=21; message=trim(message)//'data structure not found for output'
             end select
           class is (var_ilength) 
-            outputStructure(1)%indxStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
+            outputStructure(1)%indxStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
           class default; err=20; message=trim(message)//'data must not be scalarv and either of type var_dlength or var_ilength'; return
         end select
 
@@ -362,10 +362,10 @@ subroutine writeBasin(indxGRU,indxHRU,iStep,finalizeStats,&
    select case (meta(iVar)%varType)
 
     case (iLookVarType%scalarv)
-      outputStructure(1)%bvarStat(1)%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat(map(iVar))%dat(iFreq)
+      outputStructure(1)%bvarStat%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat(map(iVar))%dat(iFreq)
     case (iLookVarType%routing)
      if (iFreq==1 .and. outputTimestep(iFreq)==1) then
-      outputStructure(1)%bvarStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(iFreq) = dat(iVar)%dat(iFreq)
+      outputStructure(1)%bvarStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(iFreq) = dat(iVar)%dat(iFreq)
      end if
 
     case default
@@ -417,7 +417,7 @@ subroutine writeTime(indxGRU,indxHRU,iStep,finalizeStats,meta,dat,err,message)
    if (meta(iVar)%statIndex(iFreq)/=iLookStat%inst) cycle
 
    ! add to outputStructure
-   outputStructure(1)%timeStruct(1)%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep) = dat(iVar)
+   outputStructure(1)%timeStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep) = dat(iVar)
    if (err/=0) message=trim(message)//trim(meta(iVar)%varName)
    if (err/=0) then; err=20; return; end if
 
diff --git a/build/source/actors/hru_actor/fortran_code/hru_restart.f90 b/build/source/actors/hru_actor/fortran_code/hru_restart.f90
index ab487caaa673df5bfad133cc72204a2c96f4afaa..cc604d0448875cd52fcb9b4edc162d5bff8cc912 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_restart.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_restart.f90
@@ -69,8 +69,8 @@ subroutine summa_readRestart(&
   USE nrtype                                                  ! variable types, etc.
   ! functions and subroutines
   USE time_utils_module,only:elapsedSec                       ! calculate the elapsed time
-  USE read_icond_actors_module,only:read_icond               ! module to read initial conditions
-  USE check_icond4chm_module,only:check_icond4chm             ! module to check initial conditions
+  ! USE read_icond_module,only:read_icond               ! module to read initial conditions
+  ! USE check_icond4chm_module,only:check_icond4chm             ! module to check initial conditions
   USE var_derive_module,only:calcHeight                       ! module to calculate height at layer interfaces and layer mid-point
   USE var_derive_module,only:v_shortcut                       ! module to calculate "short-cut" variables
   USE var_derive_module,only:rootDensty                       ! module to calculate the vertical distribution of roots
@@ -131,31 +131,6 @@ subroutine summa_readRestart(&
 
   ! initialize error control
   err=0; message='hru_actor_readRestart/'
-  nGRU = 1
-
-  ! *****************************************************************************
-  ! *** read/check initial conditions
-  ! *****************************************************************************
-  ! read initial conditions
-  call read_icond(&
-                  indxGRU,                       & ! intent(in):    index of GRU in gru_struc
-                  indxHRU,                       & ! intent(in):    index of HRU in gru_struc
-                  mparStruct,                    & ! intent(in):    model parameters
-                  progStruct,                    & ! intent(inout): model prognostic variables
-                  bvarStruct,                    & ! intent(inout): model basin (GRU) variables
-                  indxStruct,                    & ! intent(inout): model indices
-                  err,cmessage)                    ! intent(out):   error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
-
-  ! check initial conditions
-  call check_icond4chm(&
-                    indxGRU,                        & ! intent(in):   index of GRU in gru_struc
-                    indxHRU,                        & ! intent(in):   index of HRU in gru_struc           
-                    progStruct,                     & ! intent(in):   model prognostic (state) variables
-                    mparStruct,                     & ! intent(in):   model parameters
-                    indxStruct,                     & ! intent(in):   layer indexes
-                    err,cmessage)                   ! intent(out):   error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
 
 
   ! *****************************************************************************
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 323c4991b1ce6ef9ad6897fcc853e579a401dc96..c9edd4c1c34105cf5b1cec374711b15c9743d902 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_setup.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_setup.f90
@@ -49,9 +49,7 @@ USE globalData,only:urbanVegCategory                        ! vegetation categor
 USE globalData,only:mpar_meta,bpar_meta                     ! parameter metadata structures
 
 ! named variables to define the decisions for snow layers
-USE mDecisions_module,only:&
-  sameRulesAllLayers, & ! SNTHERM option: same combination/sub-dividion rules applied to all layers
-  rulesDependLayerIndex ! CLM option: combination/sub-dividion rules depend on layer index
+
 
 ! named variables to define LAI decisions
 USE mDecisions_module,only:&
@@ -74,7 +72,9 @@ subroutine setupHRUParam(&
                   handle_typeStruct,              & ! local classification of soil veg etc. for each HRU
                   handle_idStruct,                & ! local classification of soil veg etc. for each HRU
                   ! primary data structures (variable length vectors)
+                  handle_indxStruct,              & ! model indices
                   handle_mparStruct,              & ! model parameters
+                  handle_progStruct,              & ! model prognostic (state) variables
                   handle_bparStruct,              & ! basin-average parameters
                   handle_bvarStruct,              & ! basin-average variables
                   handle_dparStruct,              & ! default model parameters
@@ -89,10 +89,10 @@ subroutine setupHRUParam(&
   ! * desired modules
   ! ---------------------------------------------------------------------------------------
   USE nrtype                                                  ! variable types, etc.
+  USE output_structure_module,only:outputStructure
   ! subroutines and functions
   use time_utils_module,only:elapsedSec                       ! calculate the elapsed time
   USE mDecisions_module,only:mDecisions                       ! module to read model decisions
-  ! USE read_attrb_module,only:read_attrb               ! module to read local attributes
   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
@@ -125,7 +125,9 @@ subroutine setupHRUParam(&
   type(c_ptr), intent(in), value           :: handle_attrStruct    ! local attributes for each HRU
   type(c_ptr), intent(in), value           :: handle_typeStruct    ! local classification of soil veg etc. for each HRU
   type(c_ptr), intent(in), value           :: handle_idStruct      !  
+  type(c_ptr), intent(in), value           :: handle_indxStruct !  model indices
   type(c_ptr), intent(in), value           :: handle_mparStruct    ! model parameters
+  type(c_ptr), intent(in), value           :: handle_progStruct    !  model prognostic (state) variables
   type(c_ptr), intent(in), value           :: handle_bparStruct    ! basin-average parameters
   type(c_ptr), intent(in), value           :: handle_bvarStruct    ! basin-average variables
   type(c_ptr), intent(in), value           :: handle_dparStruct    ! default model parameters
@@ -139,13 +141,18 @@ subroutine setupHRUParam(&
   type(var_d),pointer                      :: attrStruct           ! local attributes for each HRU
   type(var_i),pointer                      :: typeStruct           ! local classification of soil veg etc. for each HRU
   type(var_i8),pointer                     :: idStruct             !
+  type(var_ilength),pointer                :: indxStruct           ! model indices
   type(var_dlength),pointer                :: mparStruct           ! model parameters
+  type(var_dlength),pointer                :: progStruct           ! model prognostic (state) variables
   type(var_d),pointer                      :: bparStruct           ! basin-average parameters
   type(var_dlength),pointer                :: bvarStruct           ! basin-average variables
   type(var_d),pointer                      :: dparStruct           ! default model parameters
   type(zLookup),pointer                    :: lookupStruct         ! default model parameters
   type(var_i),pointer                      :: startTime            ! start time for the model simulation
   type(var_i),pointer                      :: oldTime              ! time for the previous model time step
+
+  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
    
@@ -157,7 +164,9 @@ subroutine setupHRUParam(&
   call c_f_pointer(handle_attrStruct, attrStruct)
   call c_f_pointer(handle_typeStruct, typeStruct)
   call c_f_pointer(handle_idStruct, idStruct)
+  call c_f_pointer(handle_indxStruct, indxStruct)
   call c_f_pointer(handle_mparStruct, mparStruct)
+  call c_f_pointer(handle_progStruct, progStruct)
   call c_f_pointer(handle_bparStruct, bparStruct)
   call c_f_pointer(handle_bvarStruct, bvarStruct)
   call c_f_pointer(handle_dparStruct, dparStruct)
@@ -169,81 +178,39 @@ subroutine setupHRUParam(&
   
   oldTime%var(:) = startTime%var(:)
 
-  ! get the maximum number of snow layers
-  select case(model_decisions(iLookDECISIONS%snowLayers)%iDecision)
-    case(sameRulesAllLayers);    maxSnowLayers = 100
-    case(rulesDependLayerIndex); maxSnowLayers = 5
-    case default; err=20; 
-        message=trim(message)//'unable to identify option to combine/sub-divide snow layers'
-        print*, message
-        return
-  end select ! (option to combine/sub-divide snow layers)
-
-  ! get the maximum number of layers
-  maxLayers = gru_struc(1)%hruInfo(1)%nSoil + maxSnowLayers
-
-  ! define monthly fraction of green vegetation
-  greenVegFrac_monthly = (/0.01_dp, 0.02_dp, 0.03_dp, 0.07_dp, 0.50_dp, 0.90_dp, 0.95_dp, 0.96_dp, 0.65_dp, 0.24_dp, 0.11_dp, 0.02_dp/)
-
-  ! define urban vegetation category
-  select case(trim(model_decisions(iLookDECISIONS%vegeParTbl)%cDecision))
-    case('USGS');                     urbanVegCategory =    1
-    case('MODIFIED_IGBP_MODIS_NOAH'); urbanVegCategory =   13
-    case('plumberCABLE');             urbanVegCategory = -999
-    case('plumberCHTESSEL');          urbanVegCategory = -999
-    case('plumberSUMMA');             urbanVegCategory = -999
-    case default
-      message=trim(message)//'unable to identify vegetation category'
-      print*, message
-      return
-  end select
-
-  ! *****************************************************************************
-  ! *** compute derived model variables that are pretty much constant for the basin as a whole
-  ! *****************************************************************************
-  ! calculate the fraction of runoff in future time steps
-  call fracFuture(bparStruct%var,    &  ! vector of basin-average model parameters
-                  bvarStruct,        &  ! data structure of basin-average variables
-                  err,cmessage)                   ! error control
-  if(err/=0)then;message=trim(message)//trim(cmessage);print*,message;return;endif
-
-  ! check that the parameters are consistent
-  call paramCheck(mparStruct,err,cmessage)
-  if(err/=0)then;message=trim(message)//trim(cmessage);print*,message;return;endif
-
-  ! calculate a look-up table for the temperature-enthalpy conversion
-  call E2T_lookup(mparStruct,err,cmessage)
-  if(err/=0)then;message=trim(message)//trim(cmessage);print*, message;return;endif
-
-  ! calculate a lookup table to compute enthalpy from temperature
-  call T2E_lookup(gru_struc(indxGRU)%hruInfo(1)%nSoil,   &   ! intent(in):    number of soil layers
-                  mparStruct,        &   ! intent(in):    parameter data structure
-                  lookupStruct,      &   ! intent(inout): lookup table data structure
-                  err,cmessage)                              ! intent(out):   error control
-  if(err/=0)then; message=trim(message)//trim(cmessage);print*,message;return;endif
-
-  ! overwrite the vegetation height
-  HVT(typeStruct%var(iLookTYPE%vegTypeIndex)) = mparStruct%var(iLookPARAM%heightCanopyTop)%dat(1)
-  HVB(typeStruct%var(iLookTYPE%vegTypeIndex)) = mparStruct%var(iLookPARAM%heightCanopyBottom)%dat(1)
-
-  ! overwrite the tables for LAI and SAI
-  if(model_decisions(iLookDECISIONS%LAI_method)%iDecision == specified)then
-    SAIM(typeStruct%var(iLookTYPE%vegTypeIndex),:) = mparStruct%var(iLookPARAM%winterSAI)%dat(1)
-    LAIM(typeStruct%var(iLookTYPE%vegTypeIndex),:) = mparStruct%var(iLookPARAM%summerLAI)%dat(1)*greenVegFrac_monthly
-  endif
-
-  ! compute total area of the upstream HRUS that flow into each HRU
-  upArea = 0._dp
-  ! Check if lateral flows exists within the HRU
-  if(typeStruct%var(iLookTYPE%downHRUindex)==typeStruct%var(iLookID%hruId))then
-    upArea = upArea + attrStruct%var(iLookATTR%HRUarea)
+  ! Copy the attrStruct
+  attrStruct%var(:) = outputStructure(1)%attrStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  ! Copy the typeStruct
+  typeStruct%var(:) = outputStructure(1)%typeStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  ! Copy the idStruct
+  idStruct%var(:) = outputStructure(1)%idStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+
+  ! Copy the mparStruct
+  mparStruct%var(:) = outputStructure(1)%mparStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  ! Copy the bparStruct
+  bparStruct%var(:) = outputStructure(1)%bparStruct%gru(indxGRU)%var(:)
+  ! Copy the dparStruct
+  dparStruct%var(:) = outputStructure(1)%dparStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  ! Copy the bvarStruct
+  do ivar=1, size(outputStructure(1)%bvarStruct_init%gru(indxGRU)%var(:))
+    bvarStruct%var(ivar)%dat(:) = outputStructure(1)%bvarStruct_init%gru(indxGRU)%var(ivar)%dat(:)
+  enddo
+  ! Copy the lookup Struct if its allocated
+  if (allocated(outputStructure(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z)) then
+    do i_z=1, size(outputStructure(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z(:))
+      do iVar=1, size(outputStructure(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z(i_z)%var(:))
+        lookupStruct%z(i_z)%var(ivar)%lookup(:) = outputStructure(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z(i_z)%var(iVar)%lookup(:)
+      end do
+    end do
   endif
-
-  ! identify the total basin area for a GRU (m2)
-  associate(totalArea => bvarStruct%var(iLookBVAR%basin__totalArea)%dat(1) )
-    totalArea = 0._dp
-    totalArea = totalArea + attrStruct%var(iLookATTR%HRUarea)
-  end associate
+  ! Copy the progStruct_init
+  do ivar=1, size(outputStructure(1)%progStruct_init%gru(indxGRU)%hru(indxHRU)%var(:))
+    progStruct%var(ivar)%dat(:) = outputStructure(1)%progStruct_init%gru(indxGRU)%hru(indxHRU)%var(ivar)%dat(:)
+  enddo
+  ! copy the indexStruct_init
+  do ivar=1, size(outputStructure(1)%indxStruct_init%gru(indxGRU)%hru(indxHRU)%var(:))
+    indxStruct%var(ivar)%dat(:) = outputStructure(1)%indxStruct_init%gru(indxGRU)%hru(indxHRU)%var(ivar)%dat(:)
+  enddo
 
 end subroutine setupHRUParam
 
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 3fd07dae1aecf20350d8131a72b70e24dbadd087..91f9d6794903a301c98ee7227bb07277d28d35e5 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_writeOutput.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_writeOutput.f90
@@ -265,7 +265,7 @@ subroutine writeHRUToOutputStructure(&
   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(1)%gru(indxGRU)%hru(indxHRU)%tim(outputStep)%dat(:) = finalizeStats%dat(:)
+  outputStructure(1)%finalizeStats%gru(indxGRU)%hru(indxHRU)%tim(outputStep)%dat(:) = finalizeStats%dat(:)
 
  ! ****************************************************************************
  ! *** calculate output statistics
diff --git a/build/source/actors/job_actor/job_actor.f90 b/build/source/actors/job_actor/job_actor.f90
index e335c7c9e7cf0bfe3316f741723f89944533be2a..bee10b60c47b3d8ed26dd5327ed6887802a4f7df 100644
--- a/build/source/actors/job_actor/job_actor.f90
+++ b/build/source/actors/job_actor/job_actor.f90
@@ -28,8 +28,8 @@ subroutine job_init_fortran(file_manager, start_gru, num_gru,&
   USE summaFileManager,only:LOCAL_ATTRIBUTES                  ! name of model initial attributes file
   
   ! subroutines and functions: read dimensions (NOTE: NetCDF)
-  USE read_attrb_actors_module,only:read_dimension              ! module to read dimensions of GRU and HRU
-  USE read_icond_actors_module,only:read_icond_nlayers               ! module to read initial condition dimensions
+  USE read_attrb_module,only:read_dimension              ! module to read dimensions of GRU and HRU
+  USE read_icond_module,only:read_icond_nlayers               ! module to read initial condition dimensions
 
   USE globalData,only:indx_meta                     ! metadata structures
   USE globalData,only:startTime,finshTime,refTime,oldTime
@@ -46,7 +46,7 @@ subroutine job_init_fortran(file_manager, start_gru, num_gru,&
 
   ! dummy variables
   character(kind=c_char,len=1),intent(in)   :: file_manager
-  integer(c_int),intent(in)                 :: start_gru
+  integer(c_int),intent(inout)              :: start_gru
   integer(c_int),intent(inout)              :: num_gru
   integer(c_int),intent(inout)              :: num_hru
   integer(c_int),intent(out)                :: err
@@ -68,7 +68,7 @@ subroutine job_init_fortran(file_manager, start_gru, num_gru,&
 
   ! Set variables that were previosuly set by getCommandArguments()
   startGRU=start_gru
-  iRunMode=iRunModeFull
+  iRunMode=iRunModeGRU
   checkHRU=integerMissing
 
   call summa_SetTimesDirsAndFiles(summaFileManagerIn,err,message)
@@ -83,8 +83,8 @@ subroutine job_init_fortran(file_manager, start_gru, num_gru,&
   ! obtain the HRU and GRU dimensions in the LocalAttribute file
   attrFile = trim(SETTINGS_PATH)//trim(LOCAL_ATTRIBUTES)
   select case (iRunMode)
-    case(iRunModeFull); call read_dimension(attrFile,num_gru,num_hru,start_gru,err)
-    case(iRunModeGRU ); err=20; message='iRunModeGRU not implemented for Actors Code'
+    case(iRunModeFull); err=20; message='iRunModeFull not implemented for Actors Code'
+    case(iRunModeGRU ); call read_dimension(trim(attrFile),fileGRU,fileHRU,num_gru,num_hru,err,message,startGRU=start_gru)
     case(iRunModeHRU ); err=20; message='iRunModeHRU not implemented for Actors Code'
   end select
   if(err/=0)then; print*, trim(message); return; endif
diff --git a/build/source/engine/read_attrb.f90 b/build/source/engine/read_attrb.f90
deleted file mode 100644
index 39b1c95ee3cb51a5649f24568b8d989561ede0d2..0000000000000000000000000000000000000000
--- a/build/source/engine/read_attrb.f90
+++ /dev/null
@@ -1,220 +0,0 @@
-! SUMMA - Structure for Unifying Multiple Modeling Alternatives
-! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
-!
-! This file is part of SUMMA
-!
-! For more information see: http://www.ral.ucar.edu/projects/summa
-!
-! This program is free software: you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License
-! along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-module read_attrb_actors_module
-USE nrtype
-implicit none
-private
-public::read_dimension
-contains
-
-! ************************************************************************************************
-! public subroutine read_dimension: read HRU and GRU dimension information on local attributes
-! ************************************************************************************************
-subroutine read_dimension(attrFile,numGRUs,numHRUs,startGRU,err)
-
-  USE netcdf
-  USE netcdf_util_module,only:nc_file_open                   ! open netcdf file
-  USE netcdf_util_module,only:nc_file_close                  ! close netcdf file
-  USE nr_utility_module ,only:arth
-  ! provide access to global data
-  USE globalData,only:gru_struc                              ! gru->hru mapping structure
-  USE globalData,only:index_map                              ! hru->gru mapping structure
-
-
-  implicit none
-
-  ! Dummy Variables
-  character(*),intent(in)                :: attrFile           ! name of attributed file
-  integer(i4b),intent(in)                :: numGRUs            ! number of GRUs for the run domain
-  integer(i4b),intent(out)               :: numHRUs            ! number of HRUs for the run domain (value filled in this subroutine)
-  integer(i4b),intent(in)                :: startGRU           ! Index of the starting GRU
-  integer(i4b),intent(out)               :: err                ! error code
-  
-  ! Local Variables
-  integer(i4b)                           :: fileGRU            ! number of GRUs in the input file
-  integer(i4b)                           :: fileHRU            ! number of HRUs in the input file
-  integer(i4b)                           :: iHRU               ! HRU couinting index
-  integer(i4b)                           :: iGRU               ! GRU loop index
-  integer(8),allocatable                 :: gru_id(:),hru_id(:)! read gru/hru IDs in from attributes file
-  integer(8),allocatable                 :: hru2gru_id(:)      ! read hru->gru mapping in from attributes file
-  integer(i4b),allocatable               :: hru_ix(:)          ! hru index for search
-  character(len=256)                     :: message            ! error message
-
-
-  ! define variables for NetCDF file operation
-  integer(i4b)                           :: ncID               ! NetCDF file ID
-  integer(i4b)                           :: varID              ! NetCDF variable ID
-  integer(i4b)                           :: gruDimId           ! variable id of GRU dimension from netcdf file
-  integer(i4b)                           :: hruDimId           ! variable id of HRU dimension from netcdf file
-  character(len=256)                     :: cmessage           ! error message for downwind routine
-
-  err=0; message="read_dimension/"
-
-  ! open nc file
-  call nc_file_open(trim(attrFile),nf90_noWrite,ncID,err,cmessage)
-  if(err/=0)then
-    message=trim(message)//trim(cmessage) 
-    print*, message
-    print*, attrFile
-    return
-  end if
-  
-  ! *********************************************************************************************
-  ! read and set GRU dimensions
-  ! **********************************************************************************************
-  ! get gru dimension of whole file
-  err = nf90_inq_dimid(ncID,"gru",gruDimId)
-  if(err/=nf90_noerr)then
-    message=trim(message)//'problem finding gru dimension/'//trim(nf90_strerror(err))
-    print*, message
-    return
-  end if
-
-  err = nf90_inquire_dimension(ncID, gruDimId, len = fileGRU)
-  if(err/=nf90_noerr)then; 
-    message=trim(message)//'problem reading gru dimension/'//trim(nf90_strerror(err))
-    print*, message
-    return
-  end if
-
-  ! get hru dimension of whole file
-  err = nf90_inq_dimid(ncID,"hru",hruDimId)
-  if(err/=nf90_noerr)then
-    message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err))
-    print*, message
-    return
-  end if
-
-  err = nf90_inquire_dimension(ncID, hruDimId, len = fileHRU)
-  if(err/=nf90_noerr)then
-    message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err))
-    print*, message
-    return
-  end if
-
-  ! check dimensions
-  if(numGRUs > fileGRU .or. numGRUs < 1) then; 
-    err=20
-    message=trim(message)//"numGRUs is out of range"
-    print*, message
-    return
-  end if
-
-  ! *********************************************************************************************
-  ! read mapping vectors and populate mapping structures
-  ! **********************************************************************************************
-  ! allocate space for GRU indices and HRU indices
-  allocate(gru_id(fileGRU))
-  allocate(hru_ix(fileHRU),hru_id(fileHRU),hru2gru_id(fileHRU))
-
-  ! read gru_id from netcdf file
-  err = nf90_inq_varid(ncID,"gruId",varID)
-  if (err/=0) then
-    message=trim(message)//'problem finding gruId'
-    print*, message
-    return
-  end if
-
-  err = nf90_get_var(ncID,varID,gru_id)
-  if (err/=0) then
-    message=trim(message)//'problem reading gruId'
-    print*, message
-    return
-  end if
-
-  ! read hru_id from netcdf file
-  err = nf90_inq_varid(ncID,"hruId",varID)
-  if (err/=0) then
-    message=trim(message)//'problem finding hruId'
-    print*, message
-    return
-  end if
-
-  err = nf90_get_var(ncID,varID,hru_id)
-  if (err/=0) then
-    message=trim(message)//'problem reading hruId'
-    print*, message
-    return
-  end if
-
-  ! read hru2gru_id from netcdf file
-  err = nf90_inq_varid(ncID,"hru2gruId",varID)
-  if (err/=0) then
-    message=trim(message)//'problem finding hru2gruId'
-    print*, message
-    return
-  end if
-
-  err = nf90_get_var(ncID,varID,hru2gru_id)
-  if (err/=0) then
-    message=trim(message)//'problem reading hru2gruId'
-    print*, message
-    return
-  end if
-  
-  ! array from 1 to total # of HRUs in attributes file
-  hru_ix=arth(1,1,fileHRU)
-
-  ! check that the mappings are not alreaday allocated
-  if (allocated(gru_struc)) then
-    deallocate(gru_struc)
-  endif
-
-  if (allocated(index_map)) then
-    deallocate(index_map)
-  endif
-
-  ! allocate first level of gru to hru mapping
-  allocate(gru_struc(numGRUs))
-
-  ! allocate space for the run
-  iHRU = 1
-  do iGRU = 1,numGRUs
-    if (count(hru2gru_Id == gru_id(iGRU+startGRU-1)) < 1) then; err=20; message=trim(message)//'problem finding HRUs belonging to GRU'; return; end if
-    gru_struc(iGRU)%hruCount          = count(hru2gru_Id == gru_id(iGRU+startGRU-1))                 ! number of HRUs in each GRU
-    gru_struc(iGRU)%gru_id            = gru_id(iGRU+startGRU-1)                                  ! set gru id
-    gru_struc(iGRU)%gru_nc            = iGRU+startGRU-1                                          ! set gru index in the netcdf file
-
-    allocate(gru_struc(iGRU)%hruInfo(gru_struc(iGRU)%hruCount))                                  ! allocate second level of gru to hru map
-    gru_struc(iGRU)%hruInfo(:)%hru_nc = pack(hru_ix,hru2gru_id == gru_struc(iGRU)%gru_id)        ! set hru id in attributes netcdf file
-    gru_struc(iGRU)%hruInfo(:)%hru_ix = arth(iHRU,1,gru_struc(iGRU)%hruCount)                    ! set index of hru in run domain
-    gru_struc(iGRU)%hruInfo(:)%hru_id = hru_id(gru_struc(iGRU)%hruInfo(:)%hru_nc)                ! set id of hru
-    iHRU = iHRU + gru_struc(iGRU)%hruCount
-  end do ! iGRU = 1,nGRU
-
-  ! set hru to gru mapping
-  numHRUs = sum(gru_struc%hruCount)   ! Total number of HRUs
-  allocate(index_map(numHRUs))        ! allocate first level of hru to gru mapping
-
-  do iGRU = 1, numGRUs
-    index_map(gru_struc(iGRU)%hruInfo(:)%hru_ix)%gru_ix   = iGRU                                 ! index of gru in run domain to which the hru belongs
-    index_map(gru_struc(iGRU)%hruInfo(:)%hru_ix)%localHRU_ix = hru_ix(1:gru_struc(iGRU)%hruCount)! index of hru within the gru
-  end do ! iGRU =1, numGRUs
-
-
-  deallocate(gru_id, hru_ix, hru_id, hru2gru_id)
-  ! close netcdf file
-  call nc_file_close(ncID,err,cmessage)
-  if (err/=0) then; message=trim(message)//trim(cmessage); return; end if
-
-end subroutine read_dimension
-
-end module read_attrb_actors_module
diff --git a/build/source/netcdf/modelwrite.f90 b/build/source/netcdf/modelwrite.f90
index 36a5900ecd0dbca7e2c43a5999361c4590b94f84..a4f45e6de0844d8132e6005d0cbc7048f42ac331 100644
--- a/build/source/netcdf/modelwrite.f90
+++ b/build/source/netcdf/modelwrite.f90
@@ -205,32 +205,21 @@ subroutine writeData(ncid, finalize_stats, output_timestep, max_layers, index_gr
         ! get variable index
         err = nf90_inq_varid(ncid%var(iFreq),trim(meta(iVar)%varName),ncVarID)
         call netcdf_err(err,message)
-        if (err/=0) then
-          print*, message
-          return
+        if (err/=0) then; print*, message; return
         endif
 
         select type(dat)
           class is(var_d)
             err = nf90_put_var(ncid%var(iFreq),ncVarID,dat%var(iVar),start=(/output_timestep(iFreq)/))
             call netcdf_err(err,message)
-            if (err/=0) then
-              print*, message
-              return
-            endif
+            if (err/=0) then; print*, message; return; endif
             cycle
-            class default
-              err=20
-              message=trim(message)//'time variable must be of type var_dlength (forcing data structure)'
-              print*, message
-              return
+          class default
+            err=20;message=trim(message)//'time variable must be of type var_dlength (forcing data structure)';print*, message;return
         end select
 
         call netcdf_err(err,message)
-        if (err/=0) then
-          print*, message
-          return
-        endif
+        if (err/=0) then; print*, message;return;endif
       endif
 
       ! define the statistics index
@@ -245,17 +234,10 @@ subroutine writeData(ncid, finalize_stats, output_timestep, max_layers, index_gr
           class is (var_dlength)
             realVec(1) = stat%var(map(iVar))%dat(iFreq)
             err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realVec,start=(/index_gru,output_timestep(iFreq)/),count=(/num_gru,1/))
-            if (err/=0) then
-              print*, message
-              return
-            endif
+            if (err/=0) then; print*, message; return; endif
           class default
-            err=20
-            message=trim(message)//'stats must be scalarv and of type var_dlength'
-            print*, message
-            return
+            err=20; message=trim(message)//'stats must be scalarv and of type var_dlength'; print*, message; return
         end select ! stat
-
       else
 
          ! Write the data
@@ -366,15 +348,19 @@ subroutine writeBasin(ncid,iGRU,finalizeStats,outputTimestep,meta,stat,dat,map,e
   ! initialize error control
   err=0;message="f-writeBasin/"
 
+  print*, "WE should see this"
+
   ! loop through output frequencies
   do iFreq=1,maxvarFreq
 
     ! skip frequencies that are not needed
+    print*, "Before outputFreq"
     if(.not.outFreq(iFreq)) cycle
-
+    print*, "After outputFreq"
     ! check that we have finalized statistics for a given frequency
+    print*, "Before stats"
     if(.not.finalizeStats(iFreq)) cycle
-
+    print*, "After stats"
     ! loop through model variables
     do iVar = 1,size(meta)
 
@@ -388,6 +374,7 @@ subroutine writeBasin(ncid,iGRU,finalizeStats,outputTimestep,meta,stat,dat,map,e
       select case (meta(iVar)%varType)
 
         case (iLookVarType%scalarv)
+          print*, "output", stat(map(iVar))%dat(iFreq)
           err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),(/stat(map(iVar))%dat(iFreq)/),start=(/iGRU,outputTimestep(iFreq)/),count=(/1,1/))
 
         case (iLookVarType%routing)
diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90
deleted file mode 100644
index c95af1cc5f1ff2069a31eb8ee8eb9b28284e7ba4..0000000000000000000000000000000000000000
--- a/build/source/netcdf/read_icond.f90
+++ /dev/null
@@ -1,344 +0,0 @@
-! SUMMA - Structure for Unifying Multiple Modeling Alternatives
-! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
-!
-! This file is part of SUMMA
-!
-! For more information see: http://www.ral.ucar.edu/projects/summa
-!
-! This program is free software: you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License
-! along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-module read_icond_actors_module
-USE, intrinsic :: iso_c_binding
-USE nrtype
-USE netcdf
-USE globalData,only: ixHRUfile_min,ixHRUfile_max
-USE globalData,only: nTimeDelay   ! number of hours in the time delay histogram
-implicit none
-private
-public::read_icond
-public::read_icond_nlayers
-! define single HRU restart file
-integer(i4b), parameter :: singleHRU=1001
-integer(i4b), parameter :: multiHRU=1002
-integer(i4b), parameter :: restartFileType=multiHRU
-contains
-
- ! ************************************************************************************************
- ! public subroutine read_icond_nlayers: read model initial conditions file for number of snow/soil layers
- ! ************************************************************************************************
-subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message)
-  ! --------------------------------------------------------------------------------------------------------
-  ! modules
-  USE nrtype
-  USE var_lookup,only:iLookIndex                        ! variable lookup structure
-  USE globalData,only:gru_struc                         ! gru-hru mapping structures
-  USE netcdf_util_module,only:nc_file_close             ! close netcdf file
-  USE netcdf_util_module,only:nc_file_open              ! close netcdf file
-  USE netcdf_util_module,only:netcdf_err                ! netcdf error handling
-  USE data_types,only:gru_hru_intVec                    ! actual data
-  USE data_types,only:var_info                          ! metadata
-  implicit none
-
-  ! --------------------------------------------------------------------------------------------------------
-  ! variable declarations
-  ! dummies
-  character(*)        ,intent(in)     :: iconFile       ! name of input (restart) file
-  integer(i4b)        ,intent(in)     :: nGRU           ! total # of GRUs in run domain
-  type(var_info)      ,intent(in)     :: indx_meta(:)   ! metadata
-  integer(i4b)        ,intent(out)    :: err            ! error code
-  character(*)        ,intent(out)    :: message        ! returned error message
-
-  ! locals
-  integer(i4b)                        :: ncID                       ! netcdf file id
-  integer(i4b)                        :: dimID                      ! netcdf file dimension id
-  integer(i4b)                        :: fileHRU                    ! number of HRUs in netcdf file
-  integer(i4b)                        :: snowID, soilID             ! netcdf variable ids
-  integer(i4b)                        :: iGRU, iHRU                 ! loop indexes
-  integer(i4b)                        :: iHRU_global                ! index of HRU in the netcdf file
-  integer(i4b),allocatable            :: snowData(:)                ! number of snow layers in all HRUs
-  integer(i4b),allocatable            :: soilData(:)                ! number of soil layers in all HRUs  
-  character(len=256)                  :: cmessage                   ! downstream error message
-
-  ! --------------------------------------------------------------------------------------------------------
-  ! initialize error message
-  err=0
-  message = 'read_icond_nlayers/'
-
-  ! open netcdf file
-  call nc_file_open(iconFile,nf90_nowrite,ncid,err,cmessage);
-  if (err/=0) then; message=trim(message)//trim(cmessage); return; end if
-  
-
-  ! get number of HRUs in file (the GRU variable(s), if present, are processed at the end)
-  err = nf90_inq_dimid(ncID,"hru",dimId);               if(err/=nf90_noerr)then; message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err)); return; end if
-  err = nf90_inquire_dimension(ncID,dimId,len=fileHRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err)); return; end if
-
-  ! allocate storage for reading from file (allocate entire file size, even when doing subdomain run)
-  allocate(snowData(fileHRU))
-  allocate(soilData(fileHRU))
-  snowData = 0
-  soilData = 0
-
-  ! get netcdf ids for the variables holding number of snow and soil layers in each hru
-  err = nf90_inq_varid(ncid,trim(indx_meta(iLookIndex%nSnow)%varName),snowid); call netcdf_err(err,message)
-  err = nf90_inq_varid(ncid,trim(indx_meta(iLookIndex%nSoil)%varName),soilid); call netcdf_err(err,message)
-
-  ! get nSnow and nSoil data (reads entire state file)
-  err = nf90_get_var(ncid,snowid,snowData); call netcdf_err(err,message)
-  err = nf90_get_var(ncid,soilid,soilData); call netcdf_err(err,message)
-
-  ixHRUfile_min=huge(1)
-  ixHRUfile_max=0
-  ! find the min and max hru indices in the state file
-  do iGRU = 1,nGRU
-  do iHRU = 1,gru_struc(iGRU)%hruCount
-  if(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc < ixHRUfile_min) ixHRUfile_min = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc
-  if(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc > ixHRUfile_max) ixHRUfile_max = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc
-  end do
-  end do
-  
-
-  ! loop over grus in current run to update snow/soil layer information
-  do iGRU = 1,nGRU
-  do iHRU = 1,gru_struc(iGRU)%hruCount
-  iHRU_global = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc
-
-  ! single HRU (Note: 'restartFileType' is hardwired above to multiHRU)
-  if(restartFileType==singleHRU) then
-      gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(1)
-      gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(1)
-
-  ! multi HRU
-  else
-      gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(iHRU_global)
-      gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(iHRU_global)
-  endif
-
-  end do
-  end do
-
-
-  ! close file
-  call nc_file_close(ncid,err,cmessage)
-  if(err/=0)then;message=trim(message)//trim(cmessage);return;end if
-
-  ! cleanup
-  deallocate(snowData,soilData)
-
-end subroutine read_icond_nlayers
-
-
-! ************************************************************************************************
-! public subroutine read_icond: read model initial conditions
-! ************************************************************************************************
-subroutine read_icond(&
-                    indxGRU,                       & ! intent(in):    Index of GRU in gru_struc
-                    indxHRU,                       & ! intent(in):    Index of HRU in gru_struc
-                    mparData,                      & ! intent(in):    model parameters
-                    progData,                      & ! intent(inout): model prognostic variables
-                    bvarData,                      & ! intent(inout): model basin (GRU) variables
-                    indxData,                      & ! intent(inout): model indices
-                    err,message)                     ! intent(out):   error control
-  ! --------------------------------------------------------------------------------------------------------
-  ! modules
-  USE nrtype
-  USE var_lookup,only:iLookVarType                       ! variable lookup structure
-  USE var_lookup,only:iLookPROG                          ! variable lookup structure
-  USE var_lookup,only:iLookPARAM                         ! variable lookup structure
-  USE var_lookup,only:iLookBVAR                          ! variable lookup structure
-  USE var_lookup,only:iLookINDEX                         ! variable lookup structure
-  USE globalData,only:prog_meta                          ! metadata for prognostic variables
-  USE globalData,only:bvar_meta                          ! metadata for basin (GRU) variables
-  USE globalData,only:gru_struc                          ! gru-hru mapping structures
-  USE globaldata,only:iname_soil,iname_snow              ! named variables to describe the type of layer
-  USE data_types,only:var_ilength                        ! full integer structure
-  USE data_types,only:var_dlength                        ! double precision structure for a single HRU
-  USE data_types,only:var_info                           ! metadata
-  USE get_ixName_module,only:get_varTypeName             ! to access type strings for error messages
-  USE updatState_module,only:updateSoil                  ! update soil states
-
-  USE netcdf
-
-  USE globalData,only:init_cond_prog
-  USE globalData,only:init_cond_bvar
-
-
-  implicit none
-
-  ! --------------------------------------------------------------------------------------------------------
-  ! variable declarations
-  ! dummies
-  integer(i4b)     ,intent(in)             :: indxGRU      ! index of GRU in gru_struc
-  integer(i4b)     ,intent(in)             :: indxHRU      ! index of HRU in hru_struc
-  type(var_dlength),intent(in)             :: mparData     ! model parameters
-  type(var_dlength),intent(inout)          :: progData     ! model prognostic variables
-  type(var_dlength),intent(inout)          :: bvarData     ! model basin (GRU) variables
-  type(var_ilength),intent(inout)          :: indxData     ! model indices
-  integer(i4b)     ,intent(out)            :: err          ! error code
-  character(*)     ,intent(out)            :: message      ! returned error message
-
-  ! locals
-  character(len=256)                     :: cmessage     ! downstream error message
-  integer(i4b)                           :: fileHRU      ! number of HRUs in file
-  integer(i4b)                           :: fileGRU      ! number of GRUs in file
-  integer(i4b)                           :: iVar, i      ! loop indices
-  integer(i4b),dimension(1)              :: ndx          ! intermediate array of loop indices
-  integer(i4b)                           :: dimID        ! varible dimension ids
-  integer(i4b)                           :: ncVarID      ! variable ID in netcdf file
-  character(256)                         :: dimName      ! not used except as a placeholder in call to inq_dim function
-  integer(i4b)                           :: dimLen       ! data dimensions
-  integer(i4b)                           :: ncID         ! netcdf file ID
-  integer(i4b)                           :: ixFile       ! index in file
-  integer(i4b)                           :: nSoil, nSnow, nToto ! # layers
-  integer(i4b)                           :: nTDH          ! number of points in time-delay histogram
-  integer(i4b)                           :: iLayer,jLayer ! layer indices
-  integer(i4b),parameter                 :: nBand=2       ! number of spectral bands
-
-  character(len=32),parameter            :: scalDimName   ='scalarv'  ! dimension name for scalar data
-  character(len=32),parameter            :: midSoilDimName='midSoil'  ! dimension name for soil-only layers
-  character(len=32),parameter            :: midTotoDimName='midToto'  ! dimension name for layered varaiables
-  character(len=32),parameter            :: ifcTotoDimName='ifcToto'  ! dimension name for layered varaiables
-  character(len=32),parameter            :: tdhDimName    ='tdh'      ! dimension name for time-delay basin variables
-
-  ! --------------------------------------------------------------------------------------------------------
-
-  ! Start procedure here
-  err=0; message="read_icondActors.f90 - read_icond/"
-
-  ! loop through prognostic variables
-  do iVar = 1,size(prog_meta)
-
-    ! skip variables that are computed later
-    if(prog_meta(iVar)%varName=='scalarCanopyWat'          .or. &
-      prog_meta(iVar)%varName=='spectralSnowAlbedoDiffuse' .or. &
-      prog_meta(iVar)%varName=='scalarSurfaceTemp'         .or. &
-      prog_meta(iVar)%varName=='mLayerVolFracWat'          .or. &
-      prog_meta(iVar)%varName=='mLayerHeight'                   ) cycle
-
-
-    ! get the number of layers
-    nSnow = gru_struc(indxGRU)%hruInfo(indxHRU)%nSnow
-    nSoil = gru_struc(indxGRU)%hruInfo(indxHRU)%nSoil
-    nToto = nSnow + nSoil
-
-    ! put the data into data structures and check that none of the values are set to nf90_fill_double
-    select case (prog_meta(iVar)%varType)
-      case (iLookVarType%scalarv)
-        progData%var(iVar)%dat(1)       = init_cond_prog(iVar)%var_data(indxGRU,1)
-        if(abs(progData%var(iVar)%dat(1) - nf90_fill_double) < epsilon(init_cond_prog(iVar)%var_data))then; err=20; endif
-      case (iLookVarType%midSoil)
-        progData%var(iVar)%dat(1:nSoil) = init_cond_prog(iVar)%var_data(indxGRU,1:nSoil)
-        if(any(abs(progData%var(iVar)%dat(1:nSoil) - nf90_fill_double) < epsilon(init_cond_prog(iVar)%var_data)))then; err=20; endif
-      case (iLookVarType%midToto)
-        progData%var(iVar)%dat(1:nToto) = init_cond_prog(iVar)%var_data(indxGRU,1:nToto)
-        if(any(abs(progData%var(iVar)%dat(1:nToto) - nf90_fill_double) < epsilon(init_cond_prog(iVar)%var_data)))then; err=20; endif
-      case (iLookVarType%ifcToto)
-        progData%var(iVar)%dat(0:nToto) = init_cond_prog(iVar)%var_data(indxGRU,1:nToto+1)
-        if(any(abs(progData%var(iVar)%dat(0:nToto) - nf90_fill_double) < epsilon(init_cond_prog(iVar)%var_data)))then; err=20; endif
-      case default
-        message=trim(message)//"unexpectedVariableType[name='"//trim(prog_meta(iVar)%varName)//"';type='"//trim(get_varTypeName(prog_meta(iVar)%varType))//"']"
-        print*,message
-        err=20; return
-    end select
-
-    if(err==20)then; 
-      message=trim(message)//"data set to the fill value (name='"//trim(prog_meta(iVar)%varName)//"')"; 
-      print*, message
-      return;
-    endif
-
-    ! fix the snow albedo
-    if(progData%var(iLookPROG%scalarSnowAlbedo)%dat(1) < 0._dp)then
-      progData%var(iLookPROG%scalarSnowAlbedo)%dat(1) = mparData%var(iLookPARAM%albedoMax)%dat(1)
-    endif
-
-    ! make sure canopy water is positive
-    if(progData%var(iLookPROG%scalarCanopyWat)%dat(1) < 0.0001_rkind)then
-      progData%var(iLookPROG%scalarCanopyliq)%dat(1) = 0.0001_rkind
-    endif
-
-    ! initialize the spectral albedo
-    progData%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1:nBand) = progData%var(iLookPROG%scalarSnowAlbedo)%dat(1)
-
-  end do ! end looping through prognostic variables (iVar)
-
-  ! --------------------------------------------------------------------------------------------------------
-  ! (2) set number of layers
-  ! --------------------------------------------------------------------------------------------------------
-
-  ! save the number of layers
-  indxData%var(iLookINDEX%nSnow)%dat(1)   = gru_struc(indxGRU)%hruInfo(indxHRU)%nSnow
-  indxData%var(iLookINDEX%nSoil)%dat(1)   = gru_struc(indxGRU)%hruInfo(indxHRU)%nSoil
-  indxData%var(iLookINDEX%nLayers)%dat(1) = gru_struc(indxGRU)%hruInfo(indxHRU)%nSnow + gru_struc(indxGRU)%hruInfo(indxHRU)%nSoil
-
-  ! set layer type
-  indxData%var(iLookINDEX%layerType)%dat(1:gru_struc(indxGRU)%hruInfo(indxHRU)%nSnow) = iname_snow
-  indxData%var(iLookINDEX%layerType)%dat((gru_struc(indxGRU)%hruInfo(indxHRU)%nSnow+1):(gru_struc(indxGRU)%hruInfo(indxHRU)%nSnow+gru_struc(indxGRU)%hruInfo(indxHRU)%nSoil)) = iname_soil
-
-
-  ! --------------------------------------------------------------------------------------------------------
-  ! (3) update soil layers (diagnostic variables)
-  ! --------------------------------------------------------------------------------------------------------
-  ! loop through soil layers
-  do iLayer = 1,indxData%var(iLookINDEX%nSoil)%dat(1)
-
-    ! get layer in the total vector
-    jLayer = iLayer+indxData%var(iLookINDEX%nSnow)%dat(1)
-
-    ! update soil layers
-    call updateSoil(&
-                    ! input
-                    progData%var(iLookPROG%mLayerTemp          )%dat(jLayer),& ! intent(in): temperature vector (K)
-                    progData%var(iLookPROG%mLayerMatricHead    )%dat(iLayer),& ! intent(in): matric head (m)
-                    mparData%var(iLookPARAM%vGn_alpha          )%dat(iLayer),& ! intent(in): van Genutchen "alpha" parameter
-                    mparData%var(iLookPARAM%vGn_n              )%dat(iLayer),& ! intent(in): van Genutchen "n" parameter
-                    mparData%var(iLookPARAM%theta_sat          )%dat(iLayer),& ! intent(in): soil porosity (-)
-                    mparData%var(iLookPARAM%theta_res          )%dat(iLayer),& ! intent(in): soil residual volumetric water content (-)
-                    1._dp - 1._dp/mparData%var(iLookPARAM%vGn_n)%dat(iLayer),& ! intent(in): van Genutchen "m" parameter (-)
-                    ! output
-                    progData%var(iLookPROG%mLayerVolFracWat    )%dat(jLayer),& ! intent(out): volumetric fraction of total water (-)
-                    progData%var(iLookPROG%mLayerVolFracLiq    )%dat(jLayer),& ! intent(out): volumetric fraction of liquid water (-)
-                    progData%var(iLookPROG%mLayerVolFracIce    )%dat(jLayer),& ! intent(out): volumetric fraction of ice (-)
-                    err,message)                                                                   ! intent(out): error control
-    if (err/=0) then; message=trim(message)//trim(cmessage); return; end if
-
-  end do  ! looping through soil layers
-
-  ! --------------------------------------------------------------------------------------------------------
-  ! (2) now get the basin variable(s)
-  ! --------------------------------------------------------------------------------------------------------
-
-  ! get the index in the file: single HRU
-  if(allocated(init_cond_bvar))then
-
-    ! loop through specific basin variables (currently 1 but loop provided to enable inclusion of others)
-    ndx = (/iLookBVAR%routingRunoffFuture/)   ! array of desired variable indices
-    do i = 1,size(ndx)
-      iVar = ndx(i)
-
-      ! store data in basin var (bvar) structure
-
-      ! put the data into data structures
-      bvarData%var(iVar)%dat(1:nTDH) = init_cond_bvar(i)%var_data(indxGRU,1:nTDH)
-      ! check whether the first values is set to nf90_fill_double
-      if(any(abs(bvarData%var(iVar)%dat(1:nTDH) - nf90_fill_double) < epsilon(init_cond_bvar(i)%var_data)))then; err=20; endif
-      if(err==20)then; message=trim(message)//"data set to the fill value (name='"//trim(bvar_meta(iVar)%varName)//"')"; return; endif
-
-    end do ! end looping through basin variables
-  endif  ! end if case for not being a singleHRU run - gaurded by the structure not being allocated
-
-
-end subroutine read_icond
-
-end module read_icond_actors_module