diff --git a/build/CMakeLists.txt b/build/CMakeLists.txt
index 1ca4bff4573ecf072bc26e2c14491031843483c3..fce75751284dfe875a6bcb26076f7ea9b687d593 100644
--- a/build/CMakeLists.txt
+++ b/build/CMakeLists.txt
@@ -224,13 +224,9 @@ set(FILE_ACCESS_INTERFACE
 set(JOB_INTERFACE
     ${JOB_ACTOR_DIR}/job_actor.f90)
 set(HRU_INTERFACE
-    ${HRU_ACTOR_DIR}/fortran_code/hru_actor.f90
     ${HRU_ACTOR_DIR}/fortran_code/hru_init.f90
     ${HRU_ACTOR_DIR}/fortran_code/hru_read.f90
     ${HRU_ACTOR_DIR}/fortran_code/hru_modelRun.f90
-    ${HRU_ACTOR_DIR}/fortran_code/hru_modelwrite.f90
-    ${HRU_ACTOR_DIR}/fortran_code/hru_restart.f90
-    ${HRU_ACTOR_DIR}/fortran_code/hru_setup.f90
     ${HRU_ACTOR_DIR}/fortran_code/hru_writeOutput.f90)
 
 # Actors actual actor modules
diff --git a/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp b/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp
index 396fcbf64ac9f42f4281f74799aef76e6666f557..fea12f36ca86ca5abef28803a3a4b6bccc94ccd7 100644
--- a/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp
+++ b/build/includes/hru_actor/hru_actor_subroutine_wrappers.hpp
@@ -2,41 +2,35 @@
 
 extern "C" {
   // Initialize HRU data_structures
-	void initHRU(int* indxGRU, int* num_steps, void* hru_data, int* err);
+    void initHRU(int* indxGRU, int* num_steps, void* hru_data, int* err);
     
-      void setupHRUParam( int* indxGRU, int* indxHRU, void* hru_data, double* upArea, int* err);
-
-  
-      // Setup summa_readRestart File if this option has been chosen 
-	void summa_readRestart(int* indxGRU, int* indxHRU, void* hru_data, double* dtInit, int* err);
-
-
-      // Run the model for one timestep
-      void RunPhysics(int* id, int* stepIndex, void* hru_data, double* dt, int* dt_int_factor, int* err);
-
-      void writeHRUToOutputStructure(int* index_hru, int* index_gru, int* output_step, void* hru_data, int* err);
-      
-      
-      void setTimeZoneOffset(int* iFile, void* hru_data, int* err);
-
-
-      void HRU_readForcing(int* index_gru, int* iStep, int* iRead, int* iFile, void* hru_data,  int* err);
-
-
-      void setIDATolerances(void* hru_data,
-                            double* relTolTempCas,
-                            double* absTolTempCas,
-                            double* relTolTempVeg,
-                            double* absTolTempVeg,
-                            double* relTolWatVeg,
-                            double* absTolWatVeg,
-                            double* relTolTempSoilSnow,
-                            double* absTolTempSoilSnow,
-                            double* relTolWatSnow,
-                            double* absTolWatSnow,
-                            double* relTolMatric,
-                            double* absTolMatric,
-                            double* relTolAquifr,
-                            double* absTolAquifr);
-
+    void setupHRUParam( int* indxGRU, int* indxHRU, void* hru_data, double* upArea, int* err);
+    
+    // Setup summa_readRestart File if this option has been chosen 
+    void summa_readRestart(int* indxGRU, int* indxHRU, void* hru_data, double* dtInit, int* err);
+    
+    // Run the model for one timestep
+    void RunPhysics(int* id, int* stepIndex, void* hru_data, double* dt, int* dt_int_factor, int* err);
+    
+    void hru_writeOutput(int* index_hru, int* index_gru, int* output_step, void* hru_data, int* err);
+    
+    void setTimeZoneOffset(int* iFile, void* hru_data, int* err);
+
+    void HRU_readForcing(int* index_gru, int* iStep, int* iRead, int* iFile, void* hru_data,  int* err);
+
+    void setIDATolerances(void* hru_data,
+                        double* relTolTempCas,
+                        double* absTolTempCas,
+                        double* relTolTempVeg,
+                        double* absTolTempVeg,
+                        double* relTolWatVeg,
+                        double* absTolWatVeg,
+                        double* relTolTempSoilSnow,
+                        double* absTolTempSoilSnow,
+                        double* relTolWatSnow,
+                        double* absTolWatSnow,
+                        double* relTolMatric,
+                        double* absTolMatric,
+                        double* relTolAquifr,
+                        double* absTolAquifr);
 }
\ 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 535722293c77b61893111f0e184967c067c32aab..d37946d44dd709022c1e478ab45b7a1f8cd06267 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
@@ -34,7 +34,7 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   USE ffile_info_actors_module,only:ffile_info
   USE mDecisions_module,only:mDecisions                       ! module to read model decisions
   USE read_pinit_module,only:read_pinit                       ! module to read initial model parameter values
-  USE SummaActors_setup,only:SOIL_VEG_GEN_PARM
+  USE INIT_HRU_ACTOR,only:SOIL_VEG_GEN_PARM
   USE module_sf_noahmplsm,only:read_mp_veg_parameters         ! module to read NOAH vegetation tables
   USE def_output_actors_module,only:def_output                ! module to define output variables
   USE output_structure_module,only:initOutputStructure        ! module to initialize output structure
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 8da70f15035d4799e699f69edb09d91e30076536..61b493db9fee4a33edf50d95eaa60b189c416e1b 100644
--- a/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
+++ b/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
@@ -86,11 +86,6 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
         [=](run_hru) {
             int err = 0;
 
-            // if (self->state.timestep == 1) {
-            //     getFirstTimestep(&self->state.iFile, &self->state.forcingStep, &err);
-            //     if (self->state.forcingStep == -1) { aout(self) << "HRU - Wrong starting forcing file\n";} 
-            // }
-
             while(self->state.num_steps_until_write > 0) {
                 if (self->state.forcingStep > self->state.stepsInCurrentFFile) {
                     self->send(self->state.file_access_actor, access_forcing_v, self->state.iFile+1, self);
@@ -100,24 +95,8 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
                 self->state.num_steps_until_write--;
 
                 err = Run_HRU(self); // Simulate a Timestep
-
                 if (err != 0) {
                     // We should have already printed the error to the screen if we get here
-
-                    self->send(self->state.parent, hru_error::run_physics_unhandleable, self);
-                    self->quit();
-                    return;
-                }
-
-                writeHRUToOutputStructure(&self->state.indxHRU, 
-                                          &self->state.indxGRU, 
-                                          &self->state.output_structure_step_index,
-                                          self->state.hru_data,
-                                          &err);
-                if (err != 0) {
-                    aout(self) << "Error: HRU_Actor - writeHRUToOutputStructure - HRU = " << self->state.indxHRU
-                               << " - indxGRU = " << self->state.indxGRU << " - refGRU = " << self->state.refGRU
-                               << "\nError = " << err << "\n";
                     self->send(self->state.parent, hru_error::run_physics_unhandleable, self);
                     self->quit();
                     return;
@@ -126,7 +105,7 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
                 self->state.timestep++;
                 self->state.forcingStep++;
                 self->state.output_structure_step_index++;
- 
+                
                 if (self->state.timestep > self->state.num_steps) {
                     self->send(self, done_hru_v);
                     break;
@@ -260,6 +239,19 @@ int Run_HRU(stateful_actor<hru_state>* self) {
         return 20;
     }
 
+    hru_writeOutput(&self->state.indxHRU, 
+                &self->state.indxGRU, 
+                &self->state.output_structure_step_index,
+                self->state.hru_data,
+                &self->state.err);
+    if (self->state.err != 0) {
+        aout(self) << "Error: HRU_Actor - writeHRUToOutputStructure - HRU = " << self->state.indxHRU
+                << " - indxGRU = " << self->state.indxGRU << " - refGRU = " << self->state.refGRU
+                << "\nError = " << self->state.err  << "\n";
+        self->quit();
+        return 21;
+    }
+
     return 0;      
 }
 
diff --git a/build/source/actors/hru_actor/fortran_code/hru_actor.f90 b/build/source/actors/hru_actor/fortran_code/hru_actor.f90
deleted file mode 100644
index a022349a90f026fd452192cffc10946e096576f7..0000000000000000000000000000000000000000
--- a/build/source/actors/hru_actor/fortran_code/hru_actor.f90
+++ /dev/null
@@ -1,85 +0,0 @@
-module hru_actor
-USE,intrinsic :: iso_c_binding
-USE nrtype
-USE data_types,only:&
-                    var_i,          &  
-                    var_i8,         &
-                    var_d,          &
-                    var_ilength,    &
-                    var_dlength,    &
-                    flagVec,        &
-                    hru_type
-implicit none
-
-
-public::setIDATolerances
-
-real(dp),parameter  :: verySmall=1e-3_rkind      ! tiny number
-real(dp),parameter  :: smallOffset=1.e-8_rkind   ! small offset (units=days) to force ih=0 at the start of the day
-
-contains
-
-
-
-
-
-! Set the HRU's relative and absolute tolerances
-subroutine setIDATolerances(handle_hru_data,    &
-                            relTolTempCas,      &
-                            absTolTempCas,      &
-                            relTolTempVeg,      &
-                            absTolTempVeg,      &
-                            relTolWatVeg,       &
-                            absTolWatVeg,       &
-                            relTolTempSoilSnow, &
-                            absTolTempSoilSnow, &
-                            relTolWatSnow,      &
-                            absTolWatSnow,      &
-                            relTolMatric,       &
-                            absTolMatric,       &
-                            relTolAquifr,       &
-                            absTolAquifr) bind(C, name="setIDATolerances")
-  USE data_types,only:var_dlength
-  USE var_lookup,only:iLookPARAM
-
-  implicit none
-
-  type(c_ptr), intent(in), value          :: handle_hru_data    !  model time data
-  real(c_double),intent(in)               :: relTolTempCas
-  real(c_double),intent(in)               :: absTolTempCas
-  real(c_double),intent(in)               :: relTolTempVeg
-  real(c_double),intent(in)               :: absTolTempVeg
-  real(c_double),intent(in)               :: relTolWatVeg
-  real(c_double),intent(in)               :: absTolWatVeg
-  real(c_double),intent(in)               :: relTolTempSoilSnow
-  real(c_double),intent(in)               :: absTolTempSoilSnow
-  real(c_double),intent(in)               :: relTolWatSnow
-  real(c_double),intent(in)               :: absTolWatSnow
-  real(c_double),intent(in)               :: relTolMatric
-  real(c_double),intent(in)               :: absTolMatric
-  real(c_double),intent(in)               :: relTolAquifr
-  real(c_double),intent(in)               :: absTolAquifr
-  ! local variables
-  type(hru_type),pointer                  :: hru_data          !  model time data
-
-  call c_f_pointer(handle_hru_data, hru_data)
-  
-#ifdef SUNDIALS_ACTIVE
-  hru_data%mparStruct%var(iLookPARAM%relTolTempCas)%dat(1)       = relTolTempCas 
-  hru_data%mparStruct%var(iLookPARAM%absTolTempCas)%dat(1)       = absTolTempCas
-  hru_data%mparStruct%var(iLookPARAM%relTolTempVeg)%dat(1)       = relTolTempVeg
-  hru_data%mparStruct%var(iLookPARAM%absTolTempVeg)%dat(1)       = absTolTempVeg
-  hru_data%mparStruct%var(iLookPARAM%relTolWatVeg)%dat(1)        = relTolWatVeg
-  hru_data%mparStruct%var(iLookPARAM%absTolWatVeg)%dat(1)        = absTolWatVeg
-  hru_data%mparStruct%var(iLookPARAM%relTolTempSoilSnow)%dat(1)  = relTolTempSoilSnow
-  hru_data%mparStruct%var(iLookPARAM%absTolTempSoilSnow)%dat(1)  = absTolTempSoilSnow
-  hru_data%mparStruct%var(iLookPARAM%relTolWatSnow)%dat(1)       = relTolWatSnow
-  hru_data%mparStruct%var(iLookPARAM%absTolWatSnow)%dat(1)       = absTolWatSnow
-  hru_data%mparStruct%var(iLookPARAM%relTolMatric)%dat(1)        = relTolMatric
-  hru_data%mparStruct%var(iLookPARAM%absTolMatric)%dat(1)        = absTolMatric
-  hru_data%mparStruct%var(iLookPARAM%relTolAquifr)%dat(1)        = relTolAquifr
-  hru_data%mparStruct%var(iLookPARAM%absTolAquifr)%dat(1)        = absTolAquifr
-#endif
-  end subroutine setIDATolerances
-
-end module hru_actor
\ No newline at end of file
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 6718e021fa1edfaeb7c405e9c605b07009caecea..c1050a8b4896dffd8fca7134d7ff4b4ac88c3add 100755
--- a/build/source/actors/hru_actor/fortran_code/hru_init.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_init.f90
@@ -23,8 +23,9 @@ USE globalData,only:prog_meta,diag_meta,flux_meta,id_meta   ! metadata structure
 USE globalData,only:mpar_meta,indx_meta                     ! metadata structures
 USE globalData,only:bpar_meta,bvar_meta                     ! metadata structures
 USE globalData,only:averageFlux_meta                        ! metadata for time-step average fluxes
-! USE globalData,only:lookup_meta 
-
+#ifdef V4_ACTIVE
+USE globalData,only:lookup_meta 
+#endif
 ! statistics metadata structures
 USE globalData,only:statForc_meta                           ! child metadata for stats
 USE globalData,only:statProg_meta                           ! child metadata for stats
@@ -34,18 +35,40 @@ USE globalData,only:statIndx_meta                           ! child metadata for
 USE globalData,only:statBvar_meta                           ! child metadata for stats
 ! maxvarFreq 
 USE var_lookup,only:maxVarFreq                               ! # of available output frequencies
+! named variables
+USE var_lookup,only:iLookATTR                               ! look-up values for local attributes
+USE var_lookup,only:iLookTYPE                               ! look-up values for classification of veg, soils etc.
+USE var_lookup,only:iLookPARAM                              ! look-up values for local column model parameters
+USE var_lookup,only:iLookID                                 ! look-up values for local column model parameters
+
+USE var_lookup,only:iLookPROG                               ! look-up values for local column model prognostic (state) variables
+USE var_lookup,only:iLookDIAG                               ! look-up values for local column model diagnostic variables
+USE var_lookup,only:iLookFLUX                               ! look-up values for local column model fluxes
+USE var_lookup,only:iLookBVAR                               ! look-up values for basin-average model variables
+USE var_lookup,only:iLookDECISIONS                          ! look-up values for model decisions
+USE globalData,only:urbanVegCategory                        ! vegetation category for urban areas
+
+! named variables to define LAI decisions
+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
 
 ! safety: set private unless specified otherwise
 implicit none
 private
 public::initHRU
+public::setupHRUParam
+public::SOIL_VEG_GEN_PARM
+public::summa_readRestart
+public::setIDATolerances
 contains
-
- ! used to declare and allocate summa data structures and initialize model state to known values
- subroutine initHRU(indxGRU,            & !  Index of HRU's GRU parent
-                    num_steps,          &
-                    handle_hru_data,    &
-                    err) bind(C,name='initHRU')
+! **************************************************************************************************
+! public subroutine initHRU: ! used to declare and allocate summa data structures and initialize model state to known values
+! **************************************************************************************************
+subroutine initHRU(indxGRU,            & !  Index of HRU's GRU parent
+                   num_steps,          &
+                   handle_hru_data,    &
+                   err) bind(C,name='initHRU')
   ! ---------------------------------------------------------------------------------------
   ! * desired modules
   ! ---------------------------------------------------------------------------------------
@@ -210,6 +233,584 @@ contains
   ! end association to info in data structures
   end associate
 
- end subroutine initHRU
+end subroutine initHRU
+
+
+
+! **************************************************************************************************
+! public subroutine setupHRUParam: initializes parameter data structures (e.g. vegetation and soil parameters).
+! **************************************************************************************************
+subroutine setupHRUParam(indxGRU,                 & ! ID of hru
+                         indxHRU,                 & ! Index of the parent GRU of the HRU 
+                         handle_hru_data,         &
+                         upArea,                  & ! area upslope of each HRU,
+                         err) bind(C, name='setupHRUParam')
+  ! ---------------------------------------------------------------------------------------
+  ! * 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 paramCheck_module,only:paramCheck                       ! module to check consistency of model parameters
+  USE pOverwrite_module,only:pOverwrite                       ! module to overwrite default parameter values with info from the Noah tables
+  USE ConvE2Temp_module,only:E2T_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
+#ifdef V4_ACTIVE  
+  USE t2enthalpy_module,only:T2E_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
+#endif
+  USE var_derive_module,only:fracFuture                       ! module to calculate the fraction of runoff in future time steps (time delay histogram)
+  USE module_sf_noahmplsm,only:read_mp_veg_parameters         ! module to read NOAH vegetation tables
+  ! global data structures
+  USE globalData,only:gru_struc                               ! gru-hru mapping structures
+  USE globalData,only:localParFallback                        ! local column default parameters
+  USE globalData,only:model_decisions                         ! model decision structure
+  USE globalData,only:greenVegFrac_monthly                    ! fraction of green vegetation in each month (0-1)
+  ! output constraints
+  USE globalData,only:maxLayers                               ! maximum number of layers
+  USE globalData,only:maxSnowLayers                           ! maximum number of snow layers
+  ! timing variables
+  USE globalData,only:startSetup,endSetup                     ! date/time for the start and end of the parameter setup
+  USE globalData,only:elapsedSetup                            ! elapsed time for the parameter setup
+  ! Noah-MP parameters
+  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)
+
+  ! ---------------------------------------------------------------------------------------
+  ! * variables
+  ! ---------------------------------------------------------------------------------------
+  implicit none
+  ! dummy variables
+  ! calling variables
+  integer(c_int),intent(in)                :: indxGRU              ! Index of the parent GRU of the HRU
+  integer(c_int),intent(in)                :: indxHRU              ! ID to locate correct HRU from netcdf file 
+  type(c_ptr), intent(in), value           :: handle_hru_data      ! pointer to the hru data structure (for error messages
+  real(c_double),intent(inout)             :: upArea
+  integer(c_int),intent(inout)             :: err
+
+  ! local variables
+  type(hru_type),pointer                   :: hru_data             ! local hru data structure
+
+  integer(i4b)                             :: ivar                 ! loop counter
+  integer(i4b)                             :: i_z                  ! loop counter
+  character(len=256)                       :: message              ! error message
+  character(len=256)                       :: cmessage             ! error message of downwind routine
+
+  ! ---------------------------------------------------------------------------------------
+  ! initialize error control
+  err=0; message='setupHRUParam/'
+
+  ! convert to fortran pointer from C++ pointer
+  call c_f_pointer(handle_hru_data, hru_data)
+
+  ! ffile_info and mDecisions moved to their own seperate subroutine call
+
+  hru_data%oldTime_hru%var(:) = hru_data%startTime_hru%var(:)
+
+  ! Copy the attrStruct
+  hru_data%attrStruct%var(:) = outputStructure(1)%attrStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  ! Copy the typeStruct
+  hru_data%typeStruct%var(:) = outputStructure(1)%typeStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  ! Copy the idStruct
+  hru_data%idStruct%var(:) = outputStructure(1)%idStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+
+  ! Copy the mparStruct
+  hru_data%mparStruct%var(:) = outputStructure(1)%mparStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  ! Copy the bparStruct
+  hru_data%bparStruct%var(:) = outputStructure(1)%bparStruct%gru(indxGRU)%var(:)
+  ! Copy the dparStruct
+  hru_data%dparStruct%var(:) = outputStructure(1)%dparStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  ! Copy the bvarStruct
+  do ivar=1, size(outputStructure(1)%bvarStruct_init%gru(indxGRU)%var(:))
+    hru_data%bvarStruct%var(ivar)%dat(:) = outputStructure(1)%bvarStruct_init%gru(indxGRU)%var(ivar)%dat(:)
+  enddo
+  ! Copy the lookup Struct if its allocated
+#ifdef V4_ACTIVE
+  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(:))
+        hru_data%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
+#endif
+  ! Copy the progStruct_init
+  do ivar=1, size(outputStructure(1)%progStruct_init%gru(indxGRU)%hru(indxHRU)%var(:))
+    hru_data%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(:))
+    hru_data%indxStruct%var(ivar)%dat(:) = outputStructure(1)%indxStruct_init%gru(indxGRU)%hru(indxHRU)%var(ivar)%dat(:)
+  enddo
+end subroutine setupHRUParam
+
+
+
+! **************************************************************************************************
+! private subroutine SOIL_VEG_GEN_PARM: Read soil, vegetation and other model parameters (from NOAH)
+! **************************************************************************************************
+SUBROUTINE SOIL_VEG_GEN_PARM(FILENAME_VEGTABLE, FILENAME_SOILTABLE, FILENAME_GENERAL, MMINLU, MMINSL)
+  !-----------------------------------------------------------------
+  use module_sf_noahlsm, only : shdtbl, nrotbl, rstbl, rgltbl, &
+  &                        hstbl, snuptbl, maxalb, laimintbl, &
+  &                        bb, drysmc, f11, maxsmc, laimaxtbl, &
+  &                        emissmintbl, emissmaxtbl, albedomintbl, &
+  &                        albedomaxtbl, wltsmc, qtz, refsmc, &
+  &                        z0mintbl, z0maxtbl, &
+  &                        satpsi, satdk, satdw, &
+  &                        theta_res, theta_sat, vGn_alpha, vGn_n, k_soil, &  ! MPC add van Genutchen parameters
+  &                        fxexp_data, lvcoef_data, &
+  &                        lutype, maxalb, &
+  &                        slope_data, frzk_data, bare, cmcmax_data, &
+  &                        cfactr_data, csoil_data, czil_data, &
+  &                        refkdt_data, natural, refdk_data, &
+  &                        rsmax_data, salp_data, sbeta_data, &
+  &                        zbot_data, smhigh_data, smlow_data, &
+  &                        lucats, topt_data, slcats, slpcats, sltype
+
+  IMPLICIT NONE
+
+  CHARACTER(LEN=*), INTENT(IN) :: FILENAME_VEGTABLE, FILENAME_SOILTABLE, FILENAME_GENERAL
+  CHARACTER(LEN=*), INTENT(IN) :: MMINLU, MMINSL
+  integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
+  integer :: ierr
+  INTEGER , PARAMETER :: OPEN_OK = 0
+
+  character*128 :: mess , message
+
+  !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
+  !             ALBBCK: SFC albedo (in percentage)
+  !                 Z0: Roughness length (m)
+  !             SHDFAC: Green vegetation fraction (in percentage)
+  !  Note: The ALBEDO, Z0, and SHDFAC values read from the following table
+  !          ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
+  !          the monthly green vegetation data
+  !             CMXTBL: MAX CNPY Capacity (m)
+  !             NROTBL: Rooting depth (layer)
+  !              RSMIN: Mimimum stomatal resistance (s m-1)
+  !              RSMAX: Max. stomatal resistance (s m-1)
+  !                RGL: Parameters used in radiation stress function
+  !                 HS: Parameter used in vapor pressure deficit functio
+  !               TOPT: Optimum transpiration air temperature. (K)
+  !             CMCMAX: Maximum canopy water capacity
+  !             CFACTR: Parameter used in the canopy inteception calculati
+  !               SNUP: Threshold snow depth (in water equivalent m) that
+  !                     implies 100% snow cover
+  !                LAI: Leaf area index (dimensionless)
+  !             MAXALB: Upper bound on maximum albedo over deep snow
+  !
+  !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
+  !
+
+  OPEN(19, FILE=trim(FILENAME_VEGTABLE),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
+  IF(ierr .NE. OPEN_OK ) THEN
+  WRITE(message,FMT='(A)') &
+  'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
+  CALL wrf_error_fatal ( message )
+  END IF
+
+  LUMATCH=0
+
+  FIND_LUTYPE : DO WHILE (LUMATCH == 0)
+  READ (19,*,END=2002)
+  READ (19,*,END=2002)LUTYPE
+  READ (19,*)LUCATS,IINDEX
+
+  IF(LUTYPE.EQ.MMINLU)THEN
+  WRITE( mess , * ) 'LANDUSE TYPE = ' // TRIM ( LUTYPE ) // ' FOUND', LUCATS,' CATEGORIES'
+  ! CALL wrf_message( mess )
+  LUMATCH=1
+  ELSE
+  ! call wrf_message ( "Skipping over LUTYPE = " // TRIM ( LUTYPE ) )
+  DO LC = 1, LUCATS+12
+  read(19,*)
+  ENDDO
+  ENDIF
+  ENDDO FIND_LUTYPE
+  ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
+  IF ( SIZE(SHDTBL)       < LUCATS .OR. &
+  SIZE(NROTBL)       < LUCATS .OR. &
+  SIZE(RSTBL)        < LUCATS .OR. &
+  SIZE(RGLTBL)       < LUCATS .OR. &
+  SIZE(HSTBL)        < LUCATS .OR. &
+  SIZE(SNUPTBL)      < LUCATS .OR. &
+  SIZE(MAXALB)       < LUCATS .OR. &
+  SIZE(LAIMINTBL)    < LUCATS .OR. &
+  SIZE(LAIMAXTBL)    < LUCATS .OR. &
+  SIZE(Z0MINTBL)     < LUCATS .OR. &
+  SIZE(Z0MAXTBL)     < LUCATS .OR. &
+  SIZE(ALBEDOMINTBL) < LUCATS .OR. &
+  SIZE(ALBEDOMAXTBL) < LUCATS .OR. &
+  SIZE(EMISSMINTBL ) < LUCATS .OR. &
+  SIZE(EMISSMAXTBL ) < LUCATS ) THEN
+  CALL wrf_error_fatal('Table sizes too small for value of LUCATS in module_sf_noahdrv.F')
+  ENDIF
+
+  IF(LUTYPE.EQ.MMINLU)THEN
+  DO LC=1,LUCATS
+  READ (19,*)IINDEX,SHDTBL(LC),                        &
+  NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
+  SNUPTBL(LC),MAXALB(LC), LAIMINTBL(LC),     &
+  LAIMAXTBL(LC),EMISSMINTBL(LC),             &
+  EMISSMAXTBL(LC), ALBEDOMINTBL(LC),         &
+  ALBEDOMAXTBL(LC), Z0MINTBL(LC), Z0MAXTBL(LC)
+  ENDDO
+
+  READ (19,*)
+  READ (19,*)TOPT_DATA
+  READ (19,*)
+  READ (19,*)CMCMAX_DATA
+  READ (19,*)
+  READ (19,*)CFACTR_DATA
+  READ (19,*)
+  READ (19,*)RSMAX_DATA
+  READ (19,*)
+  READ (19,*)BARE
+  READ (19,*)
+  READ (19,*)NATURAL
+  ENDIF
+
+  2002 CONTINUE
+
+  CLOSE (19)
+  IF (LUMATCH == 0) then
+  CALL wrf_error_fatal ("Land Use Dataset '"//MMINLU//"' not found in VEGPARM.TBL.")
+  ENDIF
+
+  !
+  !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
+  !
+  OPEN(19, FILE=trim(FILENAME_SOILTABLE),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
+  IF(ierr .NE. OPEN_OK ) THEN
+  WRITE(message,FMT='(A)') &
+  'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
+  CALL wrf_error_fatal ( message )
+  END IF
+
+  WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ', TRIM ( MMINSL )
+  ! CALL wrf_message( mess )
+
+  LUMATCH=0
+
+  ! MPC add a new soil table
+  FIND_soilTYPE : DO WHILE (LUMATCH == 0)
+  READ (19,*)
+  READ (19,*,END=2003)SLTYPE
+  READ (19,*)SLCATS,IINDEX
+  IF(SLTYPE.EQ.MMINSL)THEN
+  WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ', TRIM ( SLTYPE ) , ' FOUND', &
+  SLCATS,' CATEGORIES'
+  ! CALL wrf_message ( mess )
+  LUMATCH=1
+  ELSE
+  ! call wrf_message ( "Skipping over SLTYPE = " // TRIM ( SLTYPE ) )
+  DO LC = 1, SLCATS
+  read(19,*)
+  ENDDO
+  ENDIF
+  ENDDO FIND_soilTYPE
+  ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
+  IF ( SIZE(BB    ) < SLCATS .OR. &
+  SIZE(DRYSMC) < SLCATS .OR. &
+  SIZE(F11   ) < SLCATS .OR. &
+  SIZE(MAXSMC) < SLCATS .OR. &
+  SIZE(REFSMC) < SLCATS .OR. &
+  SIZE(SATPSI) < SLCATS .OR. &
+  SIZE(SATDK ) < SLCATS .OR. &
+  SIZE(SATDW ) < SLCATS .OR. &
+  SIZE(WLTSMC) < SLCATS .OR. &
+  SIZE(QTZ   ) < SLCATS  ) THEN
+  CALL wrf_error_fatal('Table sizes too small for value of SLCATS in module_sf_noahdrv.F')
+  ENDIF
+
+  ! MPC add new soil table
+  select case(trim(SLTYPE))
+  case('STAS','STAS-RUC')  ! original soil tables
+  DO LC=1,SLCATS
+  READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
+  REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC),   &
+  WLTSMC(LC), QTZ(LC)
+  ENDDO
+  case('ROSETTA')          ! new soil table
+  DO LC=1,SLCATS
+  READ (19,*) IINDEX,&
+  ! new soil parameters (from Rosetta)
+  theta_res(LC), theta_sat(LC),        &
+  vGn_alpha(LC), vGn_n(LC), k_soil(LC), &
+  ! original soil parameters
+  BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
+  REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC),   &
+  WLTSMC(LC), QTZ(LC)
+  ENDDO
+  case default
+  CALL wrf_message( 'SOIL TEXTURE IN INPUT FILE DOES NOT ' )
+  CALL wrf_message( 'MATCH SOILPARM TABLE'                 )
+  CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
+  end select
+
+  2003 CONTINUE
+
+  CLOSE (19)
+
+  IF(LUMATCH.EQ.0)THEN
+  CALL wrf_message( 'SOIL TEXTURE IN INPUT FILE DOES NOT ' )
+  CALL wrf_message( 'MATCH SOILPARM TABLE'                 )
+  CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
+  ENDIF
+
+  !
+  !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
+  !
+  OPEN(19, FILE=trim(FILENAME_GENERAL),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
+  IF(ierr .NE. OPEN_OK ) THEN
+  WRITE(message,FMT='(A)') &
+  'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
+  CALL wrf_error_fatal ( message )
+  END IF
+
+  READ (19,*)
+  READ (19,*)
+  READ (19,*) NUM_SLOPE
+
+  SLPCATS=NUM_SLOPE
+  ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
+  IF ( SIZE(slope_data) < NUM_SLOPE ) THEN
+  CALL wrf_error_fatal('NUM_SLOPE too large for slope_data array in module_sf_noahdrv')
+  ENDIF
+
+  DO LC=1,SLPCATS
+  READ (19,*)SLOPE_DATA(LC)
+  ENDDO
+
+  READ (19,*)
+  READ (19,*)SBETA_DATA
+  READ (19,*)
+  READ (19,*)FXEXP_DATA
+  READ (19,*)
+  READ (19,*)CSOIL_DATA
+  READ (19,*)
+  READ (19,*)SALP_DATA
+  READ (19,*)
+  READ (19,*)REFDK_DATA
+  READ (19,*)
+  READ (19,*)REFKDT_DATA
+  READ (19,*)
+  READ (19,*)FRZK_DATA
+  READ (19,*)
+  READ (19,*)ZBOT_DATA
+  READ (19,*)
+  READ (19,*)CZIL_DATA
+  READ (19,*)
+  READ (19,*)SMLOW_DATA
+  READ (19,*)
+  READ (19,*)SMHIGH_DATA
+  READ (19,*)
+  READ (19,*)LVCOEF_DATA
+  CLOSE (19)
+
+END SUBROUTINE SOIL_VEG_GEN_PARM
+
+
+
+! **************************************************************************************************
+! public subroutine summa_readRestart: read restart data and reset the model state
+! **************************************************************************************************
+subroutine summa_readRestart(indxGRU,         & ! index of GRU in gru_struc
+                            indxHRU,          & ! index of HRU in gru_struc
+                            handle_hru_data,  & ! data structure for the HRU
+                            dt_init,          & ! used to initialize the length of the sub-step for each HRU
+                            err) bind(C,name='summa_readRestart')
+  ! ---------------------------------------------------------------------------------------
+  ! * desired modules
+  ! ---------------------------------------------------------------------------------------
+  ! data types
+  USE nrtype                                                  ! variable types, etc.
+  ! functions and subroutines
+  USE time_utils_module,only:elapsedSec                       ! calculate the elapsed time
+  ! USE read_icond_module,only:read_icond               ! module to read initial conditions
+  ! USE check_icond4chm_module,only:check_icond4chm             ! module to check initial conditions
+  USE var_derive_module,only:calcHeight                       ! module to calculate height at layer interfaces and layer mid-point
+  USE var_derive_module,only:v_shortcut                       ! module to calculate "short-cut" variables
+  USE var_derive_module,only:rootDensty                       ! module to calculate the vertical distribution of roots
+  USE var_derive_module,only:satHydCond                       ! module to calculate the saturated hydraulic conductivity in each soil layer
+  ! global data structures
+  USE globalData,only:model_decisions                         ! model decision structure
+  ! timing variables
+  USE globalData,only:startRestart,endRestart                 ! date/time for the start and end of reading model restart files
+  USE globalData,only:elapsedRestart                          ! elapsed time to read model restart files
+  ! model decisions
+  USE mDecisions_module,only:&                                ! look-up values for the choice of method for the spatial representation of groundwater
+  localColumn, & ! separate groundwater representation in each local soil column
+  singleBasin    ! single groundwater store over the entire basin
+  implicit none
+  ! ---------------------------------------------------------------------------------------
+  ! Dummy variables
+  ! ---------------------------------------------------------------------------------------
+  integer(c_int),intent(in)               :: indxGRU            !  index of GRU in gru_struc
+  integer(c_int),intent(in)               :: indxHRU            !  index of HRU in gru_struc
+  type(c_ptr), intent(in), value          :: handle_hru_data   !  data structure for the HRU
+
+  real(c_double), intent(inout)           :: dt_init
+  integer(c_int), intent(inout)           :: err
+  ! ---------------------------------------------------------------------------------------
+  ! Fortran Pointers
+  ! ---------------------------------------------------------------------------------------
+  type(hru_type),pointer                 :: hru_data
+  ! ---------------------------------------------------------------------------------------
+  ! local variables
+  ! ---------------------------------------------------------------------------------------
+  character(len=256)                     :: message            ! error message
+  character(LEN=256)                     :: cmessage           ! error message of downwind routine
+  character(LEN=256)                     :: restartFile        ! restart file name
+  integer(i4b)                           :: nGRU
+  ! ---------------------------------------------------------------------------------------
+
+  call c_f_pointer(handle_hru_data, hru_data)
+
+  ! initialize error control
+  err=0; message='hru_actor_readRestart/'
+
+
+  ! *****************************************************************************
+  ! *** compute ancillary variables
+  ! *****************************************************************************
+
+  ! re-calculate height of each layer
+  call calcHeight(hru_data%indxStruct,   & ! layer type
+      hru_data%progStruct,   & ! model prognostic (state) variables for a local HRU
+      err,cmessage)                       ! error control
+  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+
+  ! calculate vertical distribution of root density
+  call rootDensty(hru_data%mparStruct,   & ! vector of model parameters
+      hru_data%indxStruct,   & ! data structure of model indices
+      hru_data%progStruct,   & ! data structure of model prognostic (state) variables
+      hru_data%diagStruct,   & ! data structure of model diagnostic variables
+      err,cmessage)                       ! error control
+  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+
+  ! calculate saturated hydraulic conductivity in each soil layer
+  call satHydCond(hru_data%mparStruct,   & ! vector of model parameters
+      hru_data%indxStruct,   & ! data structure of model indices
+      hru_data%progStruct,   & ! data structure of model prognostic (state) variables
+      hru_data%fluxStruct,   & ! data structure of model fluxes
+      err,cmessage)                       ! error control
+  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+
+  ! calculate "short-cut" variables such as volumetric heat capacity
+  call v_shortcut(hru_data%mparStruct,   & ! vector of model parameters
+      hru_data%diagStruct,   & ! data structure of model diagnostic variables
+      err,cmessage)                       ! error control
+  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+
+  ! initialize canopy drip
+  ! NOTE: canopy drip from the previous time step is used to compute throughfall for the current time step
+  hru_data%fluxStruct%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) = 0._dp  ! not used
+
+  ! *****************************************************************************
+  ! *** initialize aquifer storage
+  ! *****************************************************************************
+
+  ! initialize aquifer storage
+  ! NOTE: this is ugly: need to add capabilities to initialize basin-wide state variables
+
+  ! There are two options for groundwater:
+  !  (1) where groundwater is included in the local column (i.e., the HRUs); and
+  !  (2) where groundwater is included for the single basin (i.e., the GRUS, where multiple HRUS drain into a GRU).
+
+  ! For water balance calculations it is important to ensure that the local aquifer storage is zero if groundwater is treated as a basin-average state variable (singleBasin);
+  !  and ensure that basin-average aquifer storage is zero when groundwater is included in the local columns (localColumn).
+
+  ! select groundwater option
+  select case(model_decisions(iLookDECISIONS%spatial_gw)%iDecision)
+
+  ! the basin-average aquifer storage is not used if the groundwater is included in the local column
+  case(localColumn)
+  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 0._dp ! set to zero to be clear that there is no basin-average aquifer storage in this configuration
+
+  ! the local column aquifer storage is not used if the groundwater is basin-average
+  ! (i.e., where multiple HRUs drain to a basin-average aquifer)
+  case(singleBasin)
+  hru_data%bvarStruct%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 1._dp
+  hru_data%progStruct%var(iLookPROG%scalarAquiferStorage)%dat(1) = 0._dp  ! set to zero to be clear that there is no local aquifer storage in this configuration
+
+  ! error check
+  case default
+  message=trim(message)//'unable to identify decision for regional representation of groundwater'
+  return
+
+  end select  ! groundwater option
+
+  ! *****************************************************************************
+  ! *** initialize time step
+  ! *****************************************************************************
+
+  ! initialize time step length
+  dt_init = hru_data%progStruct%var(iLookPROG%dt_init)%dat(1) ! seconds
+
+
+  ! *****************************************************************************
+  ! *** finalize
+  ! *****************************************************************************
+
+end subroutine summa_readRestart
+
+! Set the HRU's relative and absolute tolerances
+subroutine setIDATolerances(handle_hru_data,    &
+                            relTolTempCas,      &
+                            absTolTempCas,      &
+                            relTolTempVeg,      &
+                            absTolTempVeg,      &
+                            relTolWatVeg,       &
+                            absTolWatVeg,       &
+                            relTolTempSoilSnow, &
+                            absTolTempSoilSnow, &
+                            relTolWatSnow,      &
+                            absTolWatSnow,      &
+                            relTolMatric,       &
+                            absTolMatric,       &
+                            relTolAquifr,       &
+                            absTolAquifr) bind(C, name="setIDATolerances")
+  USE data_types,only:var_dlength
+  USE var_lookup,only:iLookPARAM
+
+  implicit none
+
+  type(c_ptr), intent(in), value          :: handle_hru_data    !  model time data
+  real(c_double),intent(in)               :: relTolTempCas
+  real(c_double),intent(in)               :: absTolTempCas
+  real(c_double),intent(in)               :: relTolTempVeg
+  real(c_double),intent(in)               :: absTolTempVeg
+  real(c_double),intent(in)               :: relTolWatVeg
+  real(c_double),intent(in)               :: absTolWatVeg
+  real(c_double),intent(in)               :: relTolTempSoilSnow
+  real(c_double),intent(in)               :: absTolTempSoilSnow
+  real(c_double),intent(in)               :: relTolWatSnow
+  real(c_double),intent(in)               :: absTolWatSnow
+  real(c_double),intent(in)               :: relTolMatric
+  real(c_double),intent(in)               :: absTolMatric
+  real(c_double),intent(in)               :: relTolAquifr
+  real(c_double),intent(in)               :: absTolAquifr
+  ! local variables
+  type(hru_type),pointer                  :: hru_data          !  model time data
+
+  call c_f_pointer(handle_hru_data, hru_data)
 
+#ifdef SUNDIALS_ACTIVE
+  hru_data%mparStruct%var(iLookPARAM%relTolTempCas)%dat(1)       = relTolTempCas 
+  hru_data%mparStruct%var(iLookPARAM%absTolTempCas)%dat(1)       = absTolTempCas
+  hru_data%mparStruct%var(iLookPARAM%relTolTempVeg)%dat(1)       = relTolTempVeg
+  hru_data%mparStruct%var(iLookPARAM%absTolTempVeg)%dat(1)       = absTolTempVeg
+  hru_data%mparStruct%var(iLookPARAM%relTolWatVeg)%dat(1)        = relTolWatVeg
+  hru_data%mparStruct%var(iLookPARAM%absTolWatVeg)%dat(1)        = absTolWatVeg
+  hru_data%mparStruct%var(iLookPARAM%relTolTempSoilSnow)%dat(1)  = relTolTempSoilSnow
+  hru_data%mparStruct%var(iLookPARAM%absTolTempSoilSnow)%dat(1)  = absTolTempSoilSnow
+  hru_data%mparStruct%var(iLookPARAM%relTolWatSnow)%dat(1)       = relTolWatSnow
+  hru_data%mparStruct%var(iLookPARAM%absTolWatSnow)%dat(1)       = absTolWatSnow
+  hru_data%mparStruct%var(iLookPARAM%relTolMatric)%dat(1)        = relTolMatric
+  hru_data%mparStruct%var(iLookPARAM%absTolMatric)%dat(1)        = absTolMatric
+  hru_data%mparStruct%var(iLookPARAM%relTolAquifr)%dat(1)        = relTolAquifr
+  hru_data%mparStruct%var(iLookPARAM%absTolAquifr)%dat(1)        = absTolAquifr
+#endif
+end subroutine setIDATolerances
 end module INIT_HRU_ACTOR
diff --git a/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90 b/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90
deleted file mode 100755
index f562fd11e99abe10f19aa2cfd9ec76f6061c584e..0000000000000000000000000000000000000000
--- a/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90
+++ /dev/null
@@ -1,702 +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 contains subroutines for writing to the global output strucutre
-module hru_modelwrite_module
-
-! NetCDF types
-USE netcdf
-USE netcdf_util_module,only:netcdf_err                    ! netcdf error handling function
-
-! top-level data types
-USE nrtype
-
-! missing values
-USE globalData,only: integerMissing, realMissing
-
-! provide access to global data
-USE globalData,only:gru_struc                             ! gru->hru mapping structure
-USE output_structure_module,only:outputStructure
-
-
-! provide access to the derived types to define the data structures
-USE data_types,only:&
-                    ! final data vectors
-                    dlength,             & ! var%dat
-                    ilength,             & ! var%dat
-                    ! no spatial dimension
-                    var_i,               & ! x%var(:)            (i4b)
-                    var_i8,              & ! x%var(:)            integer(8)
-                    var_d,               & ! x%var(:)            (dp)
-                    var_ilength,         & ! x%var(:)%dat        (i4b)
-                    var_dlength,         & ! x%var(:)%dat        (dp)
-                    ! no variable dimension
-                    hru_i,               & ! x%hru(:)            (i4b)
-                    hru_d,               & ! x%hru(:)            (dp)
-                    ! gru dimension
-                    gru_int,             & ! x%gru(:)%var(:)     (i4b)
-                    gru_double,          & ! x%gru(:)%var(:)     (dp)
-                    gru_intVec,          & ! x%gru(:)%var(:)%dat (i4b)
-                    gru_doubleVec,       & ! x%gru(:)%var(:)%dat (dp)
-                    ! gru+hru dimension
-                    gru_hru_int,         & ! x%gru(:)%hru(:)%var(:)     (i4b)
-                    gru_hru_int8,        & ! x%gru(:)%hru(:)%var(:)     integer(8)
-                    gru_hru_double,      & ! x%gru(:)%hru(:)%var(:)     (dp)
-                    gru_hru_intVec,      & ! x%gru(:)%hru(:)%var(:)%dat (i4b)
-                    gru_hru_doubleVec      ! x%gru(:)%hru(:)%var(:)%dat (dp)
-
-! vector lengths
-USE var_lookup, only: maxvarFreq ! number of output frequencies
-USE var_lookup, only: maxvarStat ! number of statistics
-
-implicit none
-private
-public::writeParm
-public::writeData
-public::writeBasin
-public::writeTime
-public::writeRestart
-! define dimension lengths
-integer(i4b),parameter      :: maxSpectral=2              ! maximum number of spectral bands
-contains
-
-! **********************************************************************************************************
-! public subroutine writeParm: write model parameters
-! **********************************************************************************************************
-subroutine writeParm(indxGRU,indxHRU,ispatial,struct,meta,structName,err,message)
-  ! USE globalData,only:ncid                      ! netcdf file ids
-  USE data_types,only:var_info                    ! metadata info
-  USE var_lookup,only:iLookStat                   ! index in statistics vector
-  USE var_lookup,only:iLookFreq                   ! index in vector of model output frequencies
-  implicit none
-
-  ! declare input variables
-  integer(i4b)  ,intent(in)   :: indxGRU          ! Index into output Structure
-  integer(i4b)  ,intent(in)   :: indxHRU          ! Index into output Structure
-  integer(i4b)  ,intent(in)   :: iSpatial         ! hydrologic response unit
-  class(*)      ,intent(in)   :: struct           ! data structure
-  type(var_info),intent(in)   :: meta(:)          ! metadata structure
-  character(*)  ,intent(in)   :: structName       ! Name to know which global struct to write to
-  integer(i4b)  ,intent(out)  :: err              ! error code
-  character(*)  ,intent(out)  :: message          ! error message
-  ! local variables
-  integer(i4b)                :: iVar             ! loop through variables
-
-  ! initialize error control
-  err=0;message="hru_modelwrite.f90-writeParm/"
-
-  ! loop through local column model parameters
-  do iVar = 1,size(meta)
-
-    ! check that the variable is desired
-    if (meta(iVar)%statIndex(iLookFREQ%timestep)==integerMissing) cycle
-
-    ! initialize message
-    message=trim(message)//trim(meta(iVar)%varName)//'/'
-
-    ! HRU data
-    if (iSpatial/=integerMissing) then
-      select type (struct)
-        class is (var_i)
-        if (structName == "type")then
-          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%gru(indxGRU)%hru(indxHRU)%var(iVar) = struct%var(iVar)
-        end if
-        class is (var_dlength)
-        if (structName == "mpar")then
-          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
-      end select
-      call netcdf_err(err,message); if (err/=0) return
-
-    ! GRU data
-    else
-      select type (struct)
-        class is (var_d)
-        if (structName == "bpar")then
-          outputStructure(1)%bparStruct%gru(indxGRU)%var(iVar) = struct%var(iVar) ! this will overwrite data
-          print*, "bpar"
-        end if
-        class is (var_i8)
-        class default; err=20; message=trim(message)//'unknown variable type (no HRU)'; return
-      end select
-    end if
-    call netcdf_err(err,message); if (err/=0) return
-
-    ! re-initialize message
-    message="writeParm/"
-  end do  ! looping through local column model parameters
-
-end subroutine writeParm
-
-  ! **************************************************************************************
-  ! public subroutine writeData: write model time-dependent data
-  ! **************************************************************************************
-subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
-                      maxLayers,meta,stat,dat,map,indx,err,message)
-  USE data_types,only:var_info                       ! metadata type
-  USE var_lookup,only:maxVarStat                     ! index into stats structure
-  USE var_lookup,only:iLookVarType                   ! index into type structure
-  USE var_lookup,only:iLookIndex                     ! index into index structure
-  USE var_lookup,only:iLookStat                      ! index into stat structure
-  USE globalData,only:outFreq                        ! output file information
-  USE get_ixName_module,only:get_varTypeName         ! to access type strings for error messages
-  USE get_ixName_module,only:get_statName            ! to access type strings for error messages
-  implicit none
-  ! declare dummy variables
-  integer(i4b)  ,intent(in)        :: indxGRU
-  integer(i4b)  ,intent(in)        :: indxHRU
-  integer(i4b)  ,intent(in)        :: iStep
-  character(*)  ,intent(in)        :: structName
-  logical(lgt)  ,intent(in)        :: finalizeStats(:)  ! flags to finalize statistics
-  integer(i4b)  ,intent(in)        :: maxLayers         ! maximum number of layers
-  type(var_info),intent(in)        :: meta(:)           ! meta data
-  class(*)      ,intent(in)        :: stat              ! stats data
-  class(*)      ,intent(in)        :: dat               ! timestep data
-  integer(i4b)  ,intent(in)        :: map(:)            ! map into stats child struct
-  type(var_ilength) ,intent(in)    :: indx              ! index data
-  integer(i4b)  ,intent(out)       :: err               ! error code
-  character(*)  ,intent(out)       :: message           ! error message
-  ! local variables
-  integer(i4b)                     :: iVar              ! variable index
-  integer(i4b)                     :: iStat             ! statistics index
-  integer(i4b)                     :: iFreq             ! frequency index
-  integer(i4b)                     :: nSnow             ! number of snow layers
-  integer(i4b)                     :: nSoil             ! number of soil layers
-  integer(i4b)                     :: nLayers           ! total number of layers
-  ! output arrays
-  integer(i4b)                     :: datLength         ! length of each data vector
-  integer(i4b)                     :: maxLength         ! maximum length of each data vector
-  integer(i4b),parameter           :: ixInteger=1001    ! named variable for integer
-  integer(i4b),parameter           :: ixReal=1002       ! named variable for real
-  ! initialize error control
-  err=0;message="hru_modelwrite.f90-writeData/"
-
-  ! loop through output frequencies
-  do iFreq=1,maxvarFreq
-
-    ! skip frequencies that are not needed
-    if(.not.outFreq(iFreq)) cycle
-
-    ! check that we have finalized statistics for a given frequency
-    if(.not.finalizeStats(iFreq)) cycle
-
-    ! loop through model variables
-    do iVar = 1,size(meta)
-
-      ! handle time first
-      if (meta(iVar)%varName=='time')then
-        ! Write the time step values
-        select type(dat)      ! forcStruc
-          class is (var_d)    ! x%var(:)
-            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
-
-      ! define the statistics index
-      iStat = meta(iVar)%statIndex(iFreq)
-
-      ! check that the variable is desired
-      if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle
-
-      ! stats output: only scalar variable type
-      if(meta(iVar)%varType==iLookVarType%scalarv) then
-        select type(stat)
-          class is (var_dlength)
-            select case(trim(structName))
-            case('forc')
-              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%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
-            case('diag')
-              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%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
-            case('indx')
-              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
-          class default; err=20; message=trim(message)//'stats must be scalarv and of type gru_hru_doubleVec'; return
-        end select  ! stat
-
-        ! non-scalar variables: regular data structures
-      else
-
-        ! get the model layers
-        nSoil   = indx%var(iLookIndex%nSoil)%dat(1)
-        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%gru(indxGRU)%hru(indxHRU)%var(iLookIndex%nSnow)%tim(iStep)%dat(1) = nSnow
-        nLayers = indx%var(iLookIndex%nLayers)%dat(1)
-        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)
-          case(iLookVarType%wLength); datLength = maxSpectral
-          case(iLookVarType%midToto); datLength = nLayers
-          case(iLookVarType%midSnow); datLength = nSnow
-          case(iLookVarType%midSoil); datLength = nSoil
-          case(iLookVarType%ifcToto); datLength = nLayers+1
-          case(iLookVarType%ifcSnow); datLength = nSnow+1
-          case(iLookVarType%ifcSoil); datLength = nSoil+1
-          case default; cycle
-        end select ! vartype
-      
-        ! get the data vectors
-        select type (dat)
-          class is (var_dlength)
-            select case(trim(structName))
-              case('prog')
-                outputStructure(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = dat%var(iVar)%dat(:)
-              case('diag')
-                outputStructure(1)%diagStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = dat%var(iVar)%dat(:)
-              case('flux')
-                outputStructure(1)%fluxStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = 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%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = 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
-
-        ! get the maximum length of each data vector
-        select case (meta(iVar)%varType)
-          case(iLookVarType%wLength); maxLength = maxSpectral
-          case(iLookVarType%midToto); maxLength = maxLayers
-          case(iLookVarType%midSnow); maxLength = maxLayers-nSoil
-          case(iLookVarType%midSoil); maxLength = nSoil
-          case(iLookVarType%ifcToto); maxLength = maxLayers+1
-          case(iLookVarType%ifcSnow); maxLength = (maxLayers-nSoil)+1
-          case(iLookVarType%ifcSoil); maxLength = nSoil+1
-          case default; cycle
-        end select ! vartype
-
-      end if ! not scalarv
-
-        ! process error code
-      if (err/=0) message=trim(message)//trim(meta(iVar)%varName)//'_'//trim(get_statName(iStat))
-        call netcdf_err(err,message); if (err/=0) return
-
-    end do ! iVar
-  end do ! iFreq
-end subroutine writeData
-
-! **************************************************************************************
-! public subroutine writeBasin: write basin-average variables
-! **************************************************************************************
-subroutine writeBasin(indxGRU,indxHRU,iStep,finalizeStats,&
-                      outputTimestep,meta,stat,dat,map,err,message)
- USE data_types,only:var_info                       ! metadata type
- USE var_lookup,only:maxVarStat                     ! index into stats structure
- USE var_lookup,only:iLookVarType                   ! index into type structure
- USE globalData,only:outFreq                   ! output file information
- USE get_ixName_module,only:get_varTypeName         ! to access type strings for error messages
- USE get_ixName_module,only:get_statName            ! to access type strings for error messages
- implicit none
-
- ! declare dummy variables
- integer(i4b)  ,intent(in)     :: indxGRU
- integer(i4b)  ,intent(in)     :: indxHRU
- integer(i4b)  ,intent(in)     :: iStep
- logical(lgt)  ,intent(in)     :: finalizeStats(:)  ! flags to finalize statistics
- integer(i4b)  ,intent(in)     :: outputTimestep(:) ! output time step
- type(var_info),intent(in)     :: meta(:)           ! meta data
- type(dlength) ,intent(in)     :: stat(:)           ! stats data
- type(dlength) ,intent(in)     :: dat(:)            ! timestep data
- integer(i4b)  ,intent(in)     :: map(:)            ! map into stats child struct
- integer(i4b)  ,intent(out)    :: err               ! error code
- character(*)  ,intent(out)    :: message           ! error message
- ! local variables
- integer(i4b)                  :: iVar              ! variable index
- integer(i4b)                  :: iStat             ! statistics index
- integer(i4b)                  :: iFreq             ! frequency index
- ! initialize error control
- err=0;message="hru_modelwrite.f90-writeBasin/"
-
- ! loop through output frequencies
- do iFreq=1,maxvarFreq
-
-  ! skip frequencies that are not needed
-  if(.not.outFreq(iFreq)) cycle
-
-  ! check that we have finalized statistics for a given frequency
-  if(.not.finalizeStats(iFreq)) cycle
-
-  ! loop through model variables
-  do iVar = 1,size(meta)
-
-   ! define the statistics index
-   iStat = meta(iVar)%statIndex(iFreq)
-
-   ! check that the variable is desired
-   if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle
-
-   ! stats/data output - select data type
-   select case (meta(iVar)%varType)
-
-    case (iLookVarType%scalarv)
-      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%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(iFreq) = dat(iVar)%dat(iFreq)
-     end if
-
-    case default
-     err=40; message=trim(message)//"unknownVariableType[name='"//trim(meta(iVar)%varName)//"';type='"//trim(get_varTypeName(meta(iVar)%varType))//    "']"; return
-   end select ! variable type
-
-   ! process error code
-   if (err.ne.0) message=trim(message)//trim(meta(iVar)%varName)//'_'//trim(get_statName(iStat))
-   if (err/=0) return
-
-  end do ! iVar
- end do ! iFreq
-
-end subroutine writeBasin
-
- ! **************************************************************************************
- ! public subroutine writeTime: write current time to all files
- ! **************************************************************************************
-subroutine writeTime(indxGRU,indxHRU,iStep,finalizeStats,meta,dat,err,message)
- USE data_types,only:var_info                       ! metadata type
- USE var_lookup,only:iLookStat                      ! index into stat structure
- implicit none
-
- ! declare dummy variables
- integer(i4b)  ,intent(in)     :: indxGRU
- integer(i4b)  ,intent(in)     :: indxHRU
- integer(i4b)  ,intent(in)     :: iStep
- logical(lgt)  ,intent(in)     :: finalizeStats(:)  ! flags to finalize statistics
- type(var_info),intent(in)     :: meta(:)           ! meta data
- integer       ,intent(in)     :: dat(:)            ! timestep data
- integer(i4b)  ,intent(out)    :: err               ! error code
- character(*)  ,intent(out)    :: message           ! error message
- ! local variables
- integer(i4b)                  :: iVar              ! variable index
- integer(i4b)                  :: iFreq             ! frequency index
- ! initialize error control
- err=0;message="hru_modelwrite.f90-writeTime/"
-
- ! loop through output frequencies
- do iFreq=1,maxvarFreq
-
-  ! check that we have finalized statistics for a given frequency
-  if(.not.finalizeStats(iFreq)) cycle
-
-  ! loop through model variables
-  do iVar = 1,size(meta)
-
-   ! check instantaneous
-   if (meta(iVar)%statIndex(iFreq)/=iLookStat%inst) cycle
-
-   ! add to outputStructure
-   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
-
-  end do ! iVar
- end do ! iFreq
-
-end subroutine writeTime
-
-! *********************************************************************************************************
-! public subroutine printRestartFile: print a re-start file
-! *********************************************************************************************************
-subroutine writeRestart(filename,         & ! intent(in): name of restart file
-                         nGRU,             & ! intent(in): number of GRUs
-                         nHRU,             & ! intent(in): number of HRUs
-                         prog_meta,        & ! intent(in): prognostics metadata
-                         prog_data,        & ! intent(in): prognostics data
-                         bvar_meta,        & ! intent(in): basin (gru) variable metadata
-                         bvar_data,        & ! intent(in): basin (gru) variable data
-                         maxLayers,        & ! intent(in): maximum number of layers
-                         maxSnowLayers,    & ! intent(in): maximum number of snow layers
-                         indx_meta,        & ! intent(in): index metadata
-                         indx_data,        & ! intent(in): index data
-                         err,message)        ! intent(out): error control
- ! --------------------------------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------------------------------
- ! access the derived types to define the data structures
- USE data_types,only:var_info               ! metadata
- ! access named variables defining elements in the data structures
- USE var_lookup,only:iLookINDEX             ! named variables for structure elements
- USE var_lookup,only:iLookVarType           ! named variables for structure elements
- USE var_lookup,only:iLookBVAR              ! named variables for structure elements
- ! constants
- USE globalData,only:gru_struc              ! gru-hru mapping structures
- ! external routines
- USE netcdf_util_module,only:nc_file_close  ! close netcdf file
- USE netcdf_util_module,only:nc_file_open   ! open netcdf file
- USE globalData,only:nTimeDelay             ! number of timesteps in the time delay histogram
- 
- implicit none
- ! --------------------------------------------------------------------------------------------------------
- ! input
- character(len=256),intent(in)      :: filename      ! name of the restart file
- integer(i4b),intent(in)            :: nGRU          ! number of GRUs
- integer(i4b),intent(in)            :: nHRU          ! number of HRUs
- type(var_info),intent(in)          :: prog_meta(:)  ! prognostic variable metadata
- type(var_dlength),intent(in)       :: prog_data     ! prognostic vars
- type(var_info),intent(in)          :: bvar_meta(:)  ! basin variable metadata
- type(var_dlength),intent(in)       :: bvar_data     ! basin variables
- type(var_info),intent(in)          :: indx_meta(:)  ! metadata
- type(var_ilength),intent(in)       :: indx_data     ! indexing vars
- ! output: error control
- integer(i4b),intent(out)           :: err           ! error code
- character(*),intent(out)           :: message       ! error message
- ! --------------------------------------------------------------------------------------------------------
- ! dummy variables
- integer(i4b), intent(in)           :: maxLayers     ! maximum number of total layers
- integer(i4b), intent(in)           :: maxSnowLayers ! maximum number of snow layers
-
- ! local variables
- integer(i4b)                       :: ncid          ! netcdf file id
- integer(i4b),allocatable           :: ncVarID(:)    ! netcdf variable id
- integer(i4b)                       :: ncSnowID      ! index variable id
- integer(i4b)                       :: ncSoilID      ! index variable id
-
- integer(i4b)                       :: nSoil         ! number of soil layers
- integer(i4b)                       :: nSnow         ! number of snow layers
- integer(i4b)                       :: maxSnow       ! maximum number of snow layers
- integer(i4b)                       :: maxSoil       ! maximum number of soil layers
- integer(i4b)                       :: nLayers       ! number of total layers
- integer(i4b),parameter             :: nSpectral=2   ! number of spectal bands
- integer(i4b),parameter             :: nScalar=1     ! size of a scalar
- integer(i4b)                       :: nProgVars     ! number of prognostic variables written to state file
-
- integer(i4b)                       :: hruDimID      ! variable dimension ID
- integer(i4b)                       :: gruDimID      ! variable dimension ID
- integer(i4b)                       :: tdhDimID      ! variable dimension ID
- integer(i4b)                       :: scalDimID     ! variable dimension ID
- integer(i4b)                       :: specDimID     ! variable dimension ID
- integer(i4b)                       :: midSnowDimID  ! variable dimension ID
- integer(i4b)                       :: midSoilDimID  ! variable dimension ID
- integer(i4b)                       :: midTotoDimID  ! variable dimension ID
- integer(i4b)                       :: ifcSnowDimID  ! variable dimension ID
- integer(i4b)                       :: ifcSoilDimID  ! variable dimension ID
- integer(i4b)                       :: ifcTotoDimID  ! variable dimension ID
-
- character(len=32),parameter        :: hruDimName    ='hru'      ! dimension name for HRUs
- character(len=32),parameter        :: gruDimName    ='gru'      ! dimension name for GRUs
- character(len=32),parameter        :: tdhDimName    ='tdh'      ! dimension name for time-delay basin variables
- character(len=32),parameter        :: scalDimName   ='scalarv'  ! dimension name for scalar data
- character(len=32),parameter        :: specDimName   ='spectral' ! dimension name for spectral bands
- character(len=32),parameter        :: midSnowDimName='midSnow'  ! dimension name for snow-only layers
- 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        :: ifcSnowDimName='ifcSnow'  ! dimension name for snow-only layers
- character(len=32),parameter        :: ifcSoilDimName='ifcSoil'  ! dimension name for soil-only layers
- character(len=32),parameter        :: ifcTotoDimName='ifcToto'  ! dimension name for layered variables
-
- integer(i4b)                       :: cHRU          ! count of HRUs
- integer(i4b)                       :: iHRU          ! index of HRUs
- integer(i4b)                       :: iGRU          ! index of GRUs
- integer(i4b)                       :: iVar          ! variable index
- logical(lgt)                       :: okLength      ! flag to check if the vector length is OK
- character(len=256)                 :: cmessage      ! downstream error message
- ! --------------------------------------------------------------------------------------------------------
-
- ! initialize error control
- err=0; message='writeRestart/'
-
- ! size of prognostic variable vector
- nProgVars = size(prog_meta)
- allocate(ncVarID(nProgVars+1))     ! include 1 additional basin variable in ID array (possibly more later)
-
- ! maximum number of soil layers
- maxSoil = gru_struc(1)%hruInfo(1)%nSoil
-
- ! maximum number of snow layers
- maxSnow = maxSnowLayers
- 
- ! create file
- !TODO: Verify if this is needed
- err = nf90_create(trim(filename),nf90_classic_model,ncid)
- message='iCreate[create]'; call netcdf_err(err,message); if(err/=0)return
-
- ! define dimensions
- !TODO: Verify if this is needed
-                err = nf90_def_dim(ncid,trim(hruDimName)    ,nHRU       ,    hruDimID); message='iCreate[hru]'     ; call netcdf_err(err,message); if(err/=0)return
-                err = nf90_def_dim(ncid,trim(gruDimName)    ,nGRU       ,    gruDimID); message='iCreate[gru]'     ; call netcdf_err(err,message); if(err/=0)return
-                err = nf90_def_dim(ncid,trim(tdhDimName)    ,nTimeDelay ,    tdhDimID); message='iCreate[tdh]'     ; call netcdf_err(err,message); if(err/=0)return
-                err = nf90_def_dim(ncid,trim(scalDimName)   ,nScalar    ,   scalDimID); message='iCreate[scalar]'  ; call netcdf_err(err,message); if(err/=0)return
-                err = nf90_def_dim(ncid,trim(specDimName)   ,nSpectral  ,   specDimID); message='iCreate[spectral]'; call netcdf_err(err,message); if(err/=0)return
-                err = nf90_def_dim(ncid,trim(midSoilDimName),maxSoil    ,midSoilDimID); message='iCreate[ifcSoil]' ; call netcdf_err(err,message); if(err/=0)return
-                err = nf90_def_dim(ncid,trim(midTotoDimName),maxLayers  ,midTotoDimID); message='iCreate[midToto]' ; call netcdf_err(err,message); if(err/=0)return
-                err = nf90_def_dim(ncid,trim(ifcSoilDimName),maxSoil+1  ,ifcSoilDimID); message='iCreate[ifcSoil]' ; call netcdf_err(err,message); if(err/=0)return
-                err = nf90_def_dim(ncid,trim(ifcTotoDimName),maxLayers+1,ifcTotoDimID); message='iCreate[ifcToto]' ; call netcdf_err(err,message); if(err/=0)return
- if (maxSnow>0) err = nf90_def_dim(ncid,trim(midSnowDimName),maxSnow    ,midSnowDimID); message='iCreate[ifcSnow]' ; call netcdf_err(err,message); if(err/=0)return
- if (maxSnow>0) err = nf90_def_dim(ncid,trim(ifcSnowDimName),maxSnow+1  ,ifcSnowDimID); message='iCreate[ifcSnow]' ; call netcdf_err(err,message); if(err/=0)return
- ! re-initialize error control
- err=0; message='writeRestart/'
-
- ! define prognostic variables
- do iVar = 1,nProgVars
-  if (prog_meta(iVar)%varType==iLookvarType%unknown) cycle
-
-  ! define variable
-  select case(prog_meta(iVar)%varType)
-    !TODO: Verify if this is needed
-   case(iLookvarType%scalarv);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,  scalDimID /),ncVarID(iVar))
-   case(iLookvarType%wLength);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,  specDimID /),ncVarID(iVar))
-   case(iLookvarType%midSoil);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,midSoilDimID/),ncVarID(iVar))
-   case(iLookvarType%midToto);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,midTotoDimID/),ncVarID(iVar))
-   case(iLookvarType%ifcSoil);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,ifcSoilDimID/),ncVarID(iVar))
-   case(iLookvarType%ifcToto);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,ifcTotoDimID/),ncVarID(iVar))
-   case(iLookvarType%midSnow); if (maxSnow>0) err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,midSnowDimID/),ncVarID(iVar))
-   case(iLookvarType%ifcSnow); if (maxSnow>0) err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,ifcSnowDimID/),ncVarID(iVar))
-  end select
-
-  ! check errors
-  if(err/=0)then
-   message=trim(message)//trim(cmessage)//' [variable '//trim(prog_meta(iVar)%varName)//']'
-   return
-  end if
-
-  ! add parameter description
-  !TODO: Verify if this is needed
-  err = nf90_put_att(ncid,ncVarID(iVar),'long_name',trim(prog_meta(iVar)%vardesc))
-  call netcdf_err(err,message)
-
-  ! add parameter units
-  !TODO: Verify if this is needed
-  err = nf90_put_att(ncid,ncVarID(iVar),'units',trim(prog_meta(iVar)%varunit))
-  call netcdf_err(err,message)
-
- end do ! iVar
- 
- ! define selected basin variables (derived) -- e.g., hillslope routing
- !TODO: Verify if this is needed
- err = nf90_def_var(ncid, trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varName), nf90_double, (/gruDimID, tdhDimID /), ncVarID(nProgVars+1))
- err = nf90_put_att(ncid,ncVarID(nProgVars+1),'long_name',trim(bvar_meta(iLookBVAR%routingRunoffFuture)%vardesc));   call netcdf_err(err,message)
- err = nf90_put_att(ncid,ncVarID(nProgVars+1),'units'    ,trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varunit));   call netcdf_err(err,message)
-
- ! define index variables - snow
- !TODO: Verify if this is needed
- err = nf90_def_var(ncid,trim(indx_meta(iLookIndex%nSnow)%varName),nf90_int,(/hruDimID/),ncSnowID); call netcdf_err(err,message)
- err = nf90_put_att(ncid,ncSnowID,'long_name',trim(indx_meta(iLookIndex%nSnow)%vardesc));           call netcdf_err(err,message)
- err = nf90_put_att(ncid,ncSnowID,'units'    ,trim(indx_meta(iLookIndex%nSnow)%varunit));           call netcdf_err(err,message)
-
- ! define index variables - soil
- !TODO: Verify if this is needed
- err = nf90_def_var(ncid,trim(indx_meta(iLookIndex%nSoil)%varName),nf90_int,(/hruDimID/),ncSoilID); call netcdf_err(err,message)
- err = nf90_put_att(ncid,ncSoilID,'long_name',trim(indx_meta(iLookIndex%nSoil)%vardesc));           call netcdf_err(err,message)
- err = nf90_put_att(ncid,ncSoilID,'units'    ,trim(indx_meta(iLookIndex%nSoil)%varunit));           call netcdf_err(err,message)
-
- ! end definition phase
- !TODO: Verify if this is needed
- err = nf90_enddef(ncid); call netcdf_err(err,message); if (err/=0) return
-
- ! write variables
- do iGRU = 1,nGRU
-  do iHRU = 1,gru_struc(iGRU)%hruCount
-   cHRU = gru_struc(iGRU)%hruInfo(iHRU)%hru_ix
-   do iVar = 1,size(prog_meta)
-
-    ! excape if this variable is not used
-    if (prog_meta(iVar)%varType==iLookvarType%unknown) cycle
-
-    ! actual number of layers
-    nSnow = gru_struc(iGRU)%hruInfo(iHRU)%nSnow
-    nSoil = gru_struc(iGRU)%hruInfo(iHRU)%nSoil
-    nLayers = nSoil + nSnow
-
-    ! check size
-    ! NOTE: this may take time that we do not wish to use
-    okLength=.true.
-    select case (prog_meta(iVar)%varType)
-    ! case(iLookVarType%scalarv);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nScalar  )
-    ! case(iLookVarType%wlength);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSpectral)
-    ! case(iLookVarType%midSoil);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSoil    )
-    ! case(iLookVarType%midToto);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nLayers  )
-    ! case(iLookVarType%ifcSoil);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSoil+1  )
-    ! case(iLookVarType%ifcToto);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nLayers+1)
-    ! case(iLookVarType%midSnow); if (nSnow>0) okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSnow    )
-    ! case(iLookVarType%ifcSnow); if (nSnow>0) okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSnow+1  )
-     case(iLookVarType%scalarv);              okLength = (size(prog_data%var(iVar)%dat) == nScalar  )
-     case(iLookVarType%wlength);              okLength = (size(prog_data%var(iVar)%dat) == nSpectral)
-     case(iLookVarType%midSoil);              okLength = (size(prog_data%var(iVar)%dat) == nSoil    )
-     case(iLookVarType%midToto);              okLength = (size(prog_data%var(iVar)%dat) == nLayers  )
-     case(iLookVarType%ifcSoil);              okLength = (size(prog_data%var(iVar)%dat) == nSoil+1  )
-     case(iLookVarType%ifcToto);              okLength = (size(prog_data%var(iVar)%dat) == nLayers+1)
-     case(iLookVarType%midSnow); if (nSnow>0) okLength = (size(prog_data%var(iVar)%dat) == nSnow    )
-     case(iLookVarType%ifcSnow); if (nSnow>0) okLength = (size(prog_data%var(iVar)%dat) == nSnow+1  )
-     case default; err=20; message=trim(message)//'unknown var type'; return
-    end select
-
-    ! error check
-    if(.not.okLength)then
-     message=trim(message)//'bad vector length for variable '//trim(prog_meta(iVar)%varname)
-     err=20; return
-    endif
-
-    ! write data
-    select case (prog_meta(iVar)%varType)
-     case(iLookVarType%scalarv);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nScalar  /))
-     case(iLookVarType%wlength);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSpectral/))
-     case(iLookVarType%midSoil);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSoil    /))
-     case(iLookVarType%midToto);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nLayers  /))
-     case(iLookVarType%ifcSoil);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSoil+1  /))
-     case(iLookVarType%ifcToto);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nLayers+1/))
-     case(iLookVarType%midSnow); if (nSnow>0) err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSnow    /))
-     case(iLookVarType%ifcSnow); if (nSnow>0) err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSnow+1  /))
-     case default; err=20; message=trim(message)//'unknown var type'; return
-    end select
-
-    ! error check
-    if (err.ne.0) message=trim(message)//'writing variable:'//trim(prog_meta(iVar)%varName)
-    call netcdf_err(err,message); if (err/=0) return
-    err=0; message='writeRestart/'
-
-   end do ! iVar loop
-
-   ! write index variables
-   !TODO: Verify if this is needed
-   err=nf90_put_var(ncid,ncSnowID,(/indx_data%var(iLookIndex%nSnow)%dat/),start=(/cHRU/),count=(/1/))
-   err=nf90_put_var(ncid,ncSoilID,(/indx_data%var(iLookIndex%nSoil)%dat/),start=(/cHRU/),count=(/1/))
-
-  end do ! iHRU loop
-  
-  ! write selected basin variables
- !TODO: Verify if this is needed
-  err=nf90_put_var(ncid,ncVarID(nProgVars+1),(/bvar_data%var(iLookBVAR%routingRunoffFuture)%dat/),  start=(/iGRU/),count=(/1,nTimeDelay/))
-
- end do  ! iGRU loop
-
- ! close file
- call nc_file_close(ncid,err,cmessage)
- if(err/=0)then;message=trim(message)//trim(cmessage);return;end if
-
- ! cleanup
- deallocate(ncVarID)
-
-end subroutine writeRestart
-
-end module hru_modelwrite_module
diff --git a/build/source/actors/hru_actor/fortran_code/hru_restart.f90 b/build/source/actors/hru_actor/fortran_code/hru_restart.f90
deleted file mode 100644
index 523e501b8e406204401a1f18ab14c56b4fb02795..0000000000000000000000000000000000000000
--- a/build/source/actors/hru_actor/fortran_code/hru_restart.f90
+++ /dev/null
@@ -1,198 +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 summa_restart
-! read restart data and reset the model state
-USE,intrinsic :: iso_c_binding
-
-
-USE data_types,only:&
-                    ! no spatial dimension
-                    var_i,               & ! x%var(:)            (i4b)
-                    var_i8,              & ! x%var(:)            (i8b)
-                    var_d,               & ! x%var(:)            (dp)
-                    var_ilength,         & ! x%var(:)%dat        (i4b)
-                    var_dlength,&            ! x%var(:)%dat        (dp)
-                    hru_type
-! access missing values
-USE globalData,only:integerMissing   ! missing integer
-USE globalData,only:realMissing      ! missing double precision number
-
-! named variables
-USE var_lookup,only:iLookPROG                               ! look-up values for local column model prognostic (state) variables
-USE var_lookup,only:iLookDIAG                               ! look-up values for local column model diagnostic variables
-USE var_lookup,only:iLookFLUX                               ! look-up values for local column model fluxes
-USE var_lookup,only:iLookBVAR                               ! look-up values for basin-average model variables
-USE var_lookup,only:iLookDECISIONS                          ! look-up values for model decisions
-
-! safety: set private unless specified otherwise
-implicit none
-private
-public::summa_readRestart
-contains
-
-! read restart data and reset the model state
-subroutine summa_readRestart(&
-                indxGRU,    & ! index of GRU in gru_struc
-                indxHRU,    & ! index of HRU in gru_struc
-                handle_hru_data,   & ! data structure for the HRU
-                ! primary data structures (variable length vectors)
-                dt_init,    & ! used to initialize the length of the sub-step for each HRU
-                err) bind(C,name='summa_readRestart')
-  ! ---------------------------------------------------------------------------------------
-  ! * desired modules
-  ! ---------------------------------------------------------------------------------------
-  ! data types
-  USE nrtype                                                  ! variable types, etc.
-  ! functions and subroutines
-  USE time_utils_module,only:elapsedSec                       ! calculate the elapsed time
-  ! USE read_icond_module,only:read_icond               ! module to read initial conditions
-  ! USE check_icond4chm_module,only:check_icond4chm             ! module to check initial conditions
-  USE var_derive_module,only:calcHeight                       ! module to calculate height at layer interfaces and layer mid-point
-  USE var_derive_module,only:v_shortcut                       ! module to calculate "short-cut" variables
-  USE var_derive_module,only:rootDensty                       ! module to calculate the vertical distribution of roots
-  USE var_derive_module,only:satHydCond                       ! module to calculate the saturated hydraulic conductivity in each soil layer
-  ! global data structures
-  USE globalData,only:model_decisions                         ! model decision structure
-  ! timing variables
-  USE globalData,only:startRestart,endRestart                 ! date/time for the start and end of reading model restart files
-  USE globalData,only:elapsedRestart                          ! elapsed time to read model restart files
-  ! model decisions
-  USE mDecisions_module,only:&                                ! look-up values for the choice of method for the spatial representation of groundwater
-    localColumn, & ! separate groundwater representation in each local soil column
-    singleBasin    ! single groundwater store over the entire basin
-  implicit none
-  ! ---------------------------------------------------------------------------------------
-  ! Dummy variables
-  ! ---------------------------------------------------------------------------------------
-  integer(c_int),intent(in)               :: indxGRU            !  index of GRU in gru_struc
-  integer(c_int),intent(in)               :: indxHRU            !  index of HRU in gru_struc
-  type(c_ptr), intent(in), value          :: handle_hru_data   !  data structure for the HRU
-
-  real(c_double), intent(inout)           :: dt_init
-  integer(c_int), intent(inout)           :: err
-  ! ---------------------------------------------------------------------------------------
-  ! Fortran Pointers
-  ! ---------------------------------------------------------------------------------------
-  type(hru_type),pointer                 :: hru_data
-  ! ---------------------------------------------------------------------------------------
-  ! local variables
-  ! ---------------------------------------------------------------------------------------
-  character(len=256)                     :: message            ! error message
-  character(LEN=256)                     :: cmessage           ! error message of downwind routine
-  character(LEN=256)                     :: restartFile        ! restart file name
-  integer(i4b)                           :: nGRU
-  ! ---------------------------------------------------------------------------------------
-
-  call c_f_pointer(handle_hru_data, hru_data)
-
-  ! initialize error control
-  err=0; message='hru_actor_readRestart/'
-
-
-  ! *****************************************************************************
-  ! *** compute ancillary variables
-  ! *****************************************************************************
-
-  ! re-calculate height of each layer
-  call calcHeight(hru_data%indxStruct,   & ! layer type
-                  hru_data%progStruct,   & ! model prognostic (state) variables for a local HRU
-                  err,cmessage)                       ! error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
-
-  ! calculate vertical distribution of root density
-  call rootDensty(hru_data%mparStruct,   & ! vector of model parameters
-                  hru_data%indxStruct,   & ! data structure of model indices
-                  hru_data%progStruct,   & ! data structure of model prognostic (state) variables
-                  hru_data%diagStruct,   & ! data structure of model diagnostic variables
-                  err,cmessage)                       ! error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
-
-  ! calculate saturated hydraulic conductivity in each soil layer
-  call satHydCond(hru_data%mparStruct,   & ! vector of model parameters
-                  hru_data%indxStruct,   & ! data structure of model indices
-                  hru_data%progStruct,   & ! data structure of model prognostic (state) variables
-                  hru_data%fluxStruct,   & ! data structure of model fluxes
-                  err,cmessage)                       ! error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
-
-  ! calculate "short-cut" variables such as volumetric heat capacity
-  call v_shortcut(hru_data%mparStruct,   & ! vector of model parameters
-                  hru_data%diagStruct,   & ! data structure of model diagnostic variables
-                  err,cmessage)                       ! error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
-
-  ! initialize canopy drip
-  ! NOTE: canopy drip from the previous time step is used to compute throughfall for the current time step
-  hru_data%fluxStruct%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) = 0._dp  ! not used
-
-  ! *****************************************************************************
-  ! *** initialize aquifer storage
-  ! *****************************************************************************
-
-  ! initialize aquifer storage
-  ! NOTE: this is ugly: need to add capabilities to initialize basin-wide state variables
-
-  ! There are two options for groundwater:
-  !  (1) where groundwater is included in the local column (i.e., the HRUs); and
-  !  (2) where groundwater is included for the single basin (i.e., the GRUS, where multiple HRUS drain into a GRU).
-
-  ! For water balance calculations it is important to ensure that the local aquifer storage is zero if groundwater is treated as a basin-average state variable (singleBasin);
-  !  and ensure that basin-average aquifer storage is zero when groundwater is included in the local columns (localColumn).
-
-  ! select groundwater option
-  select case(model_decisions(iLookDECISIONS%spatial_gw)%iDecision)
-
-   ! the basin-average aquifer storage is not used if the groundwater is included in the local column
-   case(localColumn)
-    hru_data%bvarStruct%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 0._dp ! set to zero to be clear that there is no basin-average aquifer storage in this configuration
-
-   ! the local column aquifer storage is not used if the groundwater is basin-average
-   ! (i.e., where multiple HRUs drain to a basin-average aquifer)
-   case(singleBasin)
-    hru_data%bvarStruct%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 1._dp
-    hru_data%progStruct%var(iLookPROG%scalarAquiferStorage)%dat(1) = 0._dp  ! set to zero to be clear that there is no local aquifer storage in this configuration
-
-   ! error check
-   case default
-    message=trim(message)//'unable to identify decision for regional representation of groundwater'
-    return
-
-  end select  ! groundwater option
-
-  ! *****************************************************************************
-  ! *** initialize time step
-  ! *****************************************************************************
-
-  ! initialize time step length
-  dt_init = hru_data%progStruct%var(iLookPROG%dt_init)%dat(1) ! seconds
-   
-
-  ! *****************************************************************************
-  ! *** finalize
-  ! *****************************************************************************
-
-end subroutine summa_readRestart
-
-end module summa_restart
-
-
-
-
diff --git a/build/source/actors/hru_actor/fortran_code/hru_setup.f90 b/build/source/actors/hru_actor/fortran_code/hru_setup.f90
deleted file mode 100644
index 50256ba3b23f6051d6fc59a3d3da3e5ebf596e1d..0000000000000000000000000000000000000000
--- a/build/source/actors/hru_actor/fortran_code/hru_setup.f90
+++ /dev/null
@@ -1,442 +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 SummaActors_setup
-USE,intrinsic :: iso_c_binding
-
-! initializes parameter data structures (e.g. vegetation and soil parameters).
-
-USE data_types,only:&
-                    ! no spatial dimension
-                    var_i,               & ! x%var(:)            (i4b)
-                    var_i8,              & ! x%var(:)            (i8b)
-                    var_d,               & ! x%var(:)            (dp)
-                    var_ilength,         & ! x%var(:)%dat        (i4b)
-                    var_dlength,         & ! x%var(:)%dat        (dp)
-                    zLookup,            & ! x%z(:)%var(:)%lookup(dp)
-                    hru_type
-
-! access missing values
-USE globalData,only:integerMissing   ! missing integer
-USE globalData,only:realMissing      ! missing double precision number
-
-! named variables
-USE var_lookup,only:iLookATTR                               ! look-up values for local attributes
-USE var_lookup,only:iLookTYPE                               ! look-up values for classification of veg, soils etc.
-USE var_lookup,only:iLookPARAM                              ! look-up values for local column model parameters
-USE var_lookup,only:iLookID                                 ! look-up values for local column model parameters
-USE var_lookup,only:iLookBVAR                               ! look-up values for basin-average model variables
-USE var_lookup,only:iLookDECISIONS                          ! look-up values for model decisions
-USE globalData,only:urbanVegCategory                        ! vegetation category for urban areas
-
-! metadata structures
-USE globalData,only:mpar_meta,bpar_meta                     ! parameter metadata structures
-
-! named variables to define the decisions for snow layers
-
-
-! named variables to define LAI decisions
-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
-
-! safety: set private unless specified otherwise
-implicit none
-private
-public::setupHRUParam
-public::SOIL_VEG_GEN_PARM
-contains
-
-! initializes parameter data structures (e.g. vegetation and soil parameters).
-subroutine setupHRUParam(&
-                  indxGRU,                 & ! ID of hru
-                  indxHRU,                 & ! Index of the parent GRU of the HRU 
-                  handle_hru_data,         &
-                  ! miscellaneous variables
-                  upArea,                  & ! area upslope of each HRU,
-                  err) bind(C, name='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 paramCheck_module,only:paramCheck                       ! module to check consistency of model parameters
-  USE pOverwrite_module,only:pOverwrite                       ! module to overwrite default parameter values with info from the Noah tables
-  USE ConvE2Temp_module,only:E2T_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
-#ifdef V4_ACTIVE  
-  USE t2enthalpy_module,only:T2E_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
-#endif
-  USE var_derive_module,only:fracFuture                       ! module to calculate the fraction of runoff in future time steps (time delay histogram)
-  USE module_sf_noahmplsm,only:read_mp_veg_parameters         ! module to read NOAH vegetation tables
-  ! global data structures
-  USE globalData,only:gru_struc                               ! gru-hru mapping structures
-  USE globalData,only:localParFallback                        ! local column default parameters
-  USE globalData,only:model_decisions                         ! model decision structure
-  USE globalData,only:greenVegFrac_monthly                    ! fraction of green vegetation in each month (0-1)
-  ! output constraints
-  USE globalData,only:maxLayers                               ! maximum number of layers
-  USE globalData,only:maxSnowLayers                           ! maximum number of snow layers
-  ! timing variables
-  USE globalData,only:startSetup,endSetup                     ! date/time for the start and end of the parameter setup
-  USE globalData,only:elapsedSetup                            ! elapsed time for the parameter setup
-  ! Noah-MP parameters
-  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)
-
-  ! ---------------------------------------------------------------------------------------
-  ! * variables
-  ! ---------------------------------------------------------------------------------------
-  implicit none
-  ! dummy variables
-  ! calling variables
-  integer(c_int),intent(in)                :: indxGRU              ! Index of the parent GRU of the HRU
-  integer(c_int),intent(in)                :: indxHRU              ! ID to locate correct HRU from netcdf file 
-  type(c_ptr), intent(in), value           :: handle_hru_data      ! pointer to the hru data structure (for error messages
-  real(c_double),intent(inout)             :: upArea
-  integer(c_int),intent(inout)             :: err
-
-  ! local variables
-  type(hru_type),pointer                   :: hru_data             ! local hru data structure
-
-  integer(i4b)                             :: ivar                 ! loop counter
-  integer(i4b)                             :: i_z                  ! loop counter
-  character(len=256)                       :: message              ! error message
-  character(len=256)                       :: cmessage             ! error message of downwind routine
-   
-  ! ---------------------------------------------------------------------------------------
-  ! initialize error control
-  err=0; message='setupHRUParam/'
-
-  ! convert to fortran pointer from C++ pointer
-  call c_f_pointer(handle_hru_data, hru_data)
-
-  ! ffile_info and mDecisions moved to their own seperate subroutine call
-  
-  hru_data%oldTime_hru%var(:) = hru_data%startTime_hru%var(:)
-
-  ! Copy the attrStruct
-  hru_data%attrStruct%var(:) = outputStructure(1)%attrStruct%gru(indxGRU)%hru(indxHRU)%var(:)
-  ! Copy the typeStruct
-  hru_data%typeStruct%var(:) = outputStructure(1)%typeStruct%gru(indxGRU)%hru(indxHRU)%var(:)
-  ! Copy the idStruct
-  hru_data%idStruct%var(:) = outputStructure(1)%idStruct%gru(indxGRU)%hru(indxHRU)%var(:)
-
-  ! Copy the mparStruct
-  hru_data%mparStruct%var(:) = outputStructure(1)%mparStruct%gru(indxGRU)%hru(indxHRU)%var(:)
-  ! Copy the bparStruct
-  hru_data%bparStruct%var(:) = outputStructure(1)%bparStruct%gru(indxGRU)%var(:)
-  ! Copy the dparStruct
-  hru_data%dparStruct%var(:) = outputStructure(1)%dparStruct%gru(indxGRU)%hru(indxHRU)%var(:)
-  ! Copy the bvarStruct
-  do ivar=1, size(outputStructure(1)%bvarStruct_init%gru(indxGRU)%var(:))
-    hru_data%bvarStruct%var(ivar)%dat(:) = outputStructure(1)%bvarStruct_init%gru(indxGRU)%var(ivar)%dat(:)
-  enddo
-  ! Copy the lookup Struct if its allocated
-#ifdef SUNDIALS_ACTIVE
-  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(:))
-        hru_data%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
-#endif
-  ! Copy the progStruct_init
-  do ivar=1, size(outputStructure(1)%progStruct_init%gru(indxGRU)%hru(indxHRU)%var(:))
-    hru_data%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(:))
-    hru_data%indxStruct%var(ivar)%dat(:) = outputStructure(1)%indxStruct_init%gru(indxGRU)%hru(indxHRU)%var(ivar)%dat(:)
-  enddo
-
-end subroutine setupHRUParam
-
-! **************************************************************************************************
-! private subroutine SOIL_VEG_GEN_PARM: Read soil, vegetation and other model parameters (from NOAH)
-! **************************************************************************************************
-!-----------------------------------------------------------------
-SUBROUTINE SOIL_VEG_GEN_PARM(FILENAME_VEGTABLE, FILENAME_SOILTABLE, FILENAME_GENERAL, MMINLU, MMINSL)
- !-----------------------------------------------------------------
-  use module_sf_noahlsm, only : shdtbl, nrotbl, rstbl, rgltbl, &
-       &                        hstbl, snuptbl, maxalb, laimintbl, &
-       &                        bb, drysmc, f11, maxsmc, laimaxtbl, &
-       &                        emissmintbl, emissmaxtbl, albedomintbl, &
-       &                        albedomaxtbl, wltsmc, qtz, refsmc, &
-       &                        z0mintbl, z0maxtbl, &
-       &                        satpsi, satdk, satdw, &
-       &                        theta_res, theta_sat, vGn_alpha, vGn_n, k_soil, &  ! MPC add van Genutchen parameters
-       &                        fxexp_data, lvcoef_data, &
-       &                        lutype, maxalb, &
-       &                        slope_data, frzk_data, bare, cmcmax_data, &
-       &                        cfactr_data, csoil_data, czil_data, &
-       &                        refkdt_data, natural, refdk_data, &
-       &                        rsmax_data, salp_data, sbeta_data, &
-       &                        zbot_data, smhigh_data, smlow_data, &
-       &                        lucats, topt_data, slcats, slpcats, sltype
-
-  IMPLICIT NONE
-
-  CHARACTER(LEN=*), INTENT(IN) :: FILENAME_VEGTABLE, FILENAME_SOILTABLE, FILENAME_GENERAL
-  CHARACTER(LEN=*), INTENT(IN) :: MMINLU, MMINSL
-  integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
-  integer :: ierr
-  INTEGER , PARAMETER :: OPEN_OK = 0
-
-  character*128 :: mess , message
-
-  !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
-  !             ALBBCK: SFC albedo (in percentage)
-  !                 Z0: Roughness length (m)
-  !             SHDFAC: Green vegetation fraction (in percentage)
-  !  Note: The ALBEDO, Z0, and SHDFAC values read from the following table
-  !          ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
-  !          the monthly green vegetation data
-  !             CMXTBL: MAX CNPY Capacity (m)
-  !             NROTBL: Rooting depth (layer)
-  !              RSMIN: Mimimum stomatal resistance (s m-1)
-  !              RSMAX: Max. stomatal resistance (s m-1)
-  !                RGL: Parameters used in radiation stress function
-  !                 HS: Parameter used in vapor pressure deficit functio
-  !               TOPT: Optimum transpiration air temperature. (K)
-  !             CMCMAX: Maximum canopy water capacity
-  !             CFACTR: Parameter used in the canopy inteception calculati
-  !               SNUP: Threshold snow depth (in water equivalent m) that
-  !                     implies 100% snow cover
-  !                LAI: Leaf area index (dimensionless)
-  !             MAXALB: Upper bound on maximum albedo over deep snow
-  !
-  !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
-  !
-
-  OPEN(19, FILE=trim(FILENAME_VEGTABLE),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
-  IF(ierr .NE. OPEN_OK ) THEN
-     WRITE(message,FMT='(A)') &
-          'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
-     CALL wrf_error_fatal ( message )
-  END IF
-
-  LUMATCH=0
-
-  FIND_LUTYPE : DO WHILE (LUMATCH == 0)
-     READ (19,*,END=2002)
-     READ (19,*,END=2002)LUTYPE
-     READ (19,*)LUCATS,IINDEX
-
-     IF(LUTYPE.EQ.MMINLU)THEN
-        WRITE( mess , * ) 'LANDUSE TYPE = ' // TRIM ( LUTYPE ) // ' FOUND', LUCATS,' CATEGORIES'
-        ! CALL wrf_message( mess )
-        LUMATCH=1
-     ELSE
-       ! call wrf_message ( "Skipping over LUTYPE = " // TRIM ( LUTYPE ) )
-        DO LC = 1, LUCATS+12
-           read(19,*)
-        ENDDO
-     ENDIF
-  ENDDO FIND_LUTYPE
-  ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
-  IF ( SIZE(SHDTBL)       < LUCATS .OR. &
-       SIZE(NROTBL)       < LUCATS .OR. &
-       SIZE(RSTBL)        < LUCATS .OR. &
-       SIZE(RGLTBL)       < LUCATS .OR. &
-       SIZE(HSTBL)        < LUCATS .OR. &
-       SIZE(SNUPTBL)      < LUCATS .OR. &
-       SIZE(MAXALB)       < LUCATS .OR. &
-       SIZE(LAIMINTBL)    < LUCATS .OR. &
-       SIZE(LAIMAXTBL)    < LUCATS .OR. &
-       SIZE(Z0MINTBL)     < LUCATS .OR. &
-       SIZE(Z0MAXTBL)     < LUCATS .OR. &
-       SIZE(ALBEDOMINTBL) < LUCATS .OR. &
-       SIZE(ALBEDOMAXTBL) < LUCATS .OR. &
-       SIZE(EMISSMINTBL ) < LUCATS .OR. &
-       SIZE(EMISSMAXTBL ) < LUCATS ) THEN
-     CALL wrf_error_fatal('Table sizes too small for value of LUCATS in module_sf_noahdrv.F')
-  ENDIF
-
-  IF(LUTYPE.EQ.MMINLU)THEN
-     DO LC=1,LUCATS
-        READ (19,*)IINDEX,SHDTBL(LC),                        &
-             NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
-             SNUPTBL(LC),MAXALB(LC), LAIMINTBL(LC),     &
-             LAIMAXTBL(LC),EMISSMINTBL(LC),             &
-             EMISSMAXTBL(LC), ALBEDOMINTBL(LC),         &
-             ALBEDOMAXTBL(LC), Z0MINTBL(LC), Z0MAXTBL(LC)
-     ENDDO
-
-     READ (19,*)
-     READ (19,*)TOPT_DATA
-     READ (19,*)
-     READ (19,*)CMCMAX_DATA
-     READ (19,*)
-     READ (19,*)CFACTR_DATA
-     READ (19,*)
-     READ (19,*)RSMAX_DATA
-     READ (19,*)
-     READ (19,*)BARE
-     READ (19,*)
-     READ (19,*)NATURAL
-  ENDIF
-
-2002 CONTINUE
-
-  CLOSE (19)
-  IF (LUMATCH == 0) then
-     CALL wrf_error_fatal ("Land Use Dataset '"//MMINLU//"' not found in VEGPARM.TBL.")
-  ENDIF
-
-!
-!-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
-!
-  OPEN(19, FILE=trim(FILENAME_SOILTABLE),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
-  IF(ierr .NE. OPEN_OK ) THEN
-     WRITE(message,FMT='(A)') &
-          'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
-     CALL wrf_error_fatal ( message )
-  END IF
-
-  WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ', TRIM ( MMINSL )
-  ! CALL wrf_message( mess )
-
-  LUMATCH=0
-
-  ! MPC add a new soil table
-  FIND_soilTYPE : DO WHILE (LUMATCH == 0)
-   READ (19,*)
-   READ (19,*,END=2003)SLTYPE
-   READ (19,*)SLCATS,IINDEX
-   IF(SLTYPE.EQ.MMINSL)THEN
-     WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ', TRIM ( SLTYPE ) , ' FOUND', &
-          SLCATS,' CATEGORIES'
-     ! CALL wrf_message ( mess )
-     LUMATCH=1
-   ELSE
-    ! call wrf_message ( "Skipping over SLTYPE = " // TRIM ( SLTYPE ) )
-    DO LC = 1, SLCATS
-     read(19,*)
-    ENDDO
-   ENDIF
-  ENDDO FIND_soilTYPE
-  ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
-  IF ( SIZE(BB    ) < SLCATS .OR. &
-       SIZE(DRYSMC) < SLCATS .OR. &
-       SIZE(F11   ) < SLCATS .OR. &
-       SIZE(MAXSMC) < SLCATS .OR. &
-       SIZE(REFSMC) < SLCATS .OR. &
-       SIZE(SATPSI) < SLCATS .OR. &
-       SIZE(SATDK ) < SLCATS .OR. &
-       SIZE(SATDW ) < SLCATS .OR. &
-       SIZE(WLTSMC) < SLCATS .OR. &
-       SIZE(QTZ   ) < SLCATS  ) THEN
-     CALL wrf_error_fatal('Table sizes too small for value of SLCATS in module_sf_noahdrv.F')
-  ENDIF
-
-  ! MPC add new soil table
-  select case(trim(SLTYPE))
-   case('STAS','STAS-RUC')  ! original soil tables
-     DO LC=1,SLCATS
-        READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
-             REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC),   &
-             WLTSMC(LC), QTZ(LC)
-     ENDDO
-   case('ROSETTA')          ! new soil table
-     DO LC=1,SLCATS
-        READ (19,*) IINDEX,&
-             ! new soil parameters (from Rosetta)
-             theta_res(LC), theta_sat(LC),        &
-             vGn_alpha(LC), vGn_n(LC), k_soil(LC), &
-             ! original soil parameters
-             BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
-             REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC),   &
-             WLTSMC(LC), QTZ(LC)
-     ENDDO
-   case default
-     CALL wrf_message( 'SOIL TEXTURE IN INPUT FILE DOES NOT ' )
-     CALL wrf_message( 'MATCH SOILPARM TABLE'                 )
-     CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
-  end select
-
-2003 CONTINUE
-
-  CLOSE (19)
-
-  IF(LUMATCH.EQ.0)THEN
-     CALL wrf_message( 'SOIL TEXTURE IN INPUT FILE DOES NOT ' )
-     CALL wrf_message( 'MATCH SOILPARM TABLE'                 )
-     CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
-  ENDIF
-
-!
-!-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
-!
-  OPEN(19, FILE=trim(FILENAME_GENERAL),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
-  IF(ierr .NE. OPEN_OK ) THEN
-     WRITE(message,FMT='(A)') &
-          'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
-     CALL wrf_error_fatal ( message )
-  END IF
-
-  READ (19,*)
-  READ (19,*)
-  READ (19,*) NUM_SLOPE
-
-  SLPCATS=NUM_SLOPE
-! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
-  IF ( SIZE(slope_data) < NUM_SLOPE ) THEN
-     CALL wrf_error_fatal('NUM_SLOPE too large for slope_data array in module_sf_noahdrv')
-  ENDIF
-
-  DO LC=1,SLPCATS
-     READ (19,*)SLOPE_DATA(LC)
-  ENDDO
-
-  READ (19,*)
-  READ (19,*)SBETA_DATA
-  READ (19,*)
-  READ (19,*)FXEXP_DATA
-  READ (19,*)
-  READ (19,*)CSOIL_DATA
-  READ (19,*)
-  READ (19,*)SALP_DATA
-  READ (19,*)
-  READ (19,*)REFDK_DATA
-  READ (19,*)
-  READ (19,*)REFKDT_DATA
-  READ (19,*)
-  READ (19,*)FRZK_DATA
-  READ (19,*)
-  READ (19,*)ZBOT_DATA
-  READ (19,*)
-  READ (19,*)CZIL_DATA
-  READ (19,*)
-  READ (19,*)SMLOW_DATA
-  READ (19,*)
-  READ (19,*)SMHIGH_DATA
-  READ (19,*)
-  READ (19,*)LVCOEF_DATA
-  CLOSE (19)
-
-!-----------------------------------------------------------------
-END SUBROUTINE SOIL_VEG_GEN_PARM
-!-----------------------------------------------------------------
-
-end module SummaActors_setup
\ No newline at end of file
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 abf05d8b980ee7e3c4c6c457ee032896786537dd..8d0e33effd1762d5df9ee0646153fc8638d75939 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_writeOutput.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_writeOutput.f90
@@ -2,9 +2,12 @@ module HRUwriteoOutput_module
 USE,intrinsic :: iso_c_binding
 
 USE data_types,only:&
+                    dlength,        & ! var%dat
+                    ilength,        & ! var%dat
                     var_i,          &  
                     var_i8,         &
                     var_d,          &
+                    
                     var_ilength,    &
                     var_dlength,    &
                     flagVec,        &
@@ -50,20 +53,29 @@ USE var_lookup,only:iLookPROG                 ! named variables for local column
 USE var_lookup,only:iLookINDEX                ! named variables for local column index variables
 USE var_lookup,only:iLookFreq                 ! named variables for the frequency structure
 USE var_lookup,only:iLookBVAR                 ! named variables for basin parameters
+! vector lengths
+USE var_lookup, only: maxvarFreq ! number of output frequencies
+USE var_lookup, only: maxvarStat ! number of statistics
 USE get_ixname_module,only:get_freqName       ! get name of frequency from frequency index
+USE output_structure_module,only:outputStructure
 
 
 implicit none
 private
-public::writeHRUToOutputStructure
-
+public::hru_writeOutput
+public::writeParm
+public::writeData
+public::writeBasin
+public::writeTime
+public::writeRestart
+integer(i4b),parameter      :: maxSpectral=2              ! maximum number of spectral bands
 contains
-subroutine writeHRUToOutputStructure(&
+subroutine hru_writeOutput(&
                             indxHRU,                   &
                             indxGRU,                   &
                             outputStep,                & ! index into the output Struc
                             handle_hru_data,           & ! local HRU data  
-                            err) bind(C, name="writeHRUToOutputStructure") 
+                            err) bind(C, name="hru_writeOutput") 
   USE nrtype
   USE globalData,only:structInfo
   USE globalData,only:startWrite,endWrite
@@ -79,10 +91,6 @@ subroutine writeHRUToOutputStructure(&
 
   USE globalData,only:forc_meta,attr_meta,type_meta           ! metaData structures
   USE output_stats,only:calcStats                             ! module for compiling output statistics
-  USE hru_modelwrite_module,only:writeData,writeBasin       ! module to write model output
-  USE hru_modelwrite_module,only:writeTime                  ! module to write model time
-  USE hru_modelwrite_module,only:writeRestart               ! module to write model summa_readRestart
-  USE hru_modelwrite_module,only:writeParm                  ! module to write model parameters
   USE time_utils_module,only:elapsedSec                       ! calculate the elapsed time
   USE globalData,only:elapsedWrite                            ! elapsed time to write data
   USE output_structure_module,only:outputStructure
@@ -195,6 +203,628 @@ subroutine writeHRUToOutputStructure(&
   ! save time vector
   hru_data%oldTime_hru%var(:) = hru_data%timeStruct%var(:)
 
-end subroutine writeHRUToOutputStructure
+end subroutine hru_writeOutput
+
+! **********************************************************************************************************
+! public subroutine writeParm: write model parameters
+! **********************************************************************************************************
+subroutine writeParm(indxGRU,indxHRU,ispatial,struct,meta,structName,err,message)
+  ! USE globalData,only:ncid                      ! netcdf file ids
+  USE data_types,only:var_info                    ! metadata info
+  USE var_lookup,only:iLookStat                   ! index in statistics vector
+  USE var_lookup,only:iLookFreq                   ! index in vector of model output frequencies
+  implicit none
+
+  ! declare input variables
+  integer(i4b)  ,intent(in)   :: indxGRU          ! Index into output Structure
+  integer(i4b)  ,intent(in)   :: indxHRU          ! Index into output Structure
+  integer(i4b)  ,intent(in)   :: iSpatial         ! hydrologic response unit
+  class(*)      ,intent(in)   :: struct           ! data structure
+  type(var_info),intent(in)   :: meta(:)          ! metadata structure
+  character(*)  ,intent(in)   :: structName       ! Name to know which global struct to write to
+  integer(i4b)  ,intent(out)  :: err              ! error code
+  character(*)  ,intent(out)  :: message          ! error message
+  ! local variables
+  integer(i4b)                :: iVar             ! loop through variables
+
+  ! initialize error control
+  err=0;message="hru_modelwrite.f90-writeParm/"
+
+  ! loop through local column model parameters
+  do iVar = 1,size(meta)
+
+    ! check that the variable is desired
+    if (meta(iVar)%statIndex(iLookFREQ%timestep)==integerMissing) cycle
+
+    ! initialize message
+    message=trim(message)//trim(meta(iVar)%varName)//'/'
+
+    ! HRU data
+    if (iSpatial/=integerMissing) then
+      select type (struct)
+        class is (var_i)
+        if (structName == "type")then
+          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%gru(indxGRU)%hru(indxHRU)%var(iVar) = struct%var(iVar)
+        end if
+        class is (var_dlength)
+        if (structName == "mpar")then
+          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
+      end select
+      call netcdf_err(err,message); if (err/=0) return
+
+    ! GRU data
+    else
+      select type (struct)
+        class is (var_d)
+        if (structName == "bpar")then
+          outputStructure(1)%bparStruct%gru(indxGRU)%var(iVar) = struct%var(iVar) ! this will overwrite data
+          print*, "bpar"
+        end if
+        class is (var_i8)
+        class default; err=20; message=trim(message)//'unknown variable type (no HRU)'; return
+      end select
+    end if
+    call netcdf_err(err,message); if (err/=0) return
+
+    ! re-initialize message
+    message="writeParm/"
+  end do  ! looping through local column model parameters
+
+end subroutine writeParm
+
+  ! **************************************************************************************
+  ! public subroutine writeData: write model time-dependent data
+  ! **************************************************************************************
+subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
+                      maxLayers,meta,stat,dat,map,indx,err,message)
+  USE data_types,only:var_info                       ! metadata type
+  USE var_lookup,only:maxVarStat                     ! index into stats structure
+  USE var_lookup,only:iLookVarType                   ! index into type structure
+  USE var_lookup,only:iLookIndex                     ! index into index structure
+  USE var_lookup,only:iLookStat                      ! index into stat structure
+  USE globalData,only:outFreq                        ! output file information
+  USE get_ixName_module,only:get_varTypeName         ! to access type strings for error messages
+  USE get_ixName_module,only:get_statName            ! to access type strings for error messages
+  implicit none
+  ! declare dummy variables
+  integer(i4b)  ,intent(in)        :: indxGRU
+  integer(i4b)  ,intent(in)        :: indxHRU
+  integer(i4b)  ,intent(in)        :: iStep
+  character(*)  ,intent(in)        :: structName
+  logical(lgt)  ,intent(in)        :: finalizeStats(:)  ! flags to finalize statistics
+  integer(i4b)  ,intent(in)        :: maxLayers         ! maximum number of layers
+  type(var_info),intent(in)        :: meta(:)           ! meta data
+  class(*)      ,intent(in)        :: stat              ! stats data
+  class(*)      ,intent(in)        :: dat               ! timestep data
+  integer(i4b)  ,intent(in)        :: map(:)            ! map into stats child struct
+  type(var_ilength) ,intent(in)    :: indx              ! index data
+  integer(i4b)  ,intent(out)       :: err               ! error code
+  character(*)  ,intent(out)       :: message           ! error message
+  ! local variables
+  integer(i4b)                     :: iVar              ! variable index
+  integer(i4b)                     :: iStat             ! statistics index
+  integer(i4b)                     :: iFreq             ! frequency index
+  integer(i4b)                     :: nSnow             ! number of snow layers
+  integer(i4b)                     :: nSoil             ! number of soil layers
+  integer(i4b)                     :: nLayers           ! total number of layers
+  ! output arrays
+  integer(i4b)                     :: datLength         ! length of each data vector
+  integer(i4b)                     :: maxLength         ! maximum length of each data vector
+  integer(i4b),parameter           :: ixInteger=1001    ! named variable for integer
+  integer(i4b),parameter           :: ixReal=1002       ! named variable for real
+  ! initialize error control
+  err=0;message="hru_modelwrite.f90-writeData/"
+
+  ! loop through output frequencies
+  do iFreq=1,maxvarFreq
+
+    ! skip frequencies that are not needed
+    if(.not.outFreq(iFreq)) cycle
+
+    ! check that we have finalized statistics for a given frequency
+    if(.not.finalizeStats(iFreq)) cycle
+
+    ! loop through model variables
+    do iVar = 1,size(meta)
+
+      ! handle time first
+      if (meta(iVar)%varName=='time')then
+        ! Write the time step values
+        select type(dat)      ! forcStruc
+          class is (var_d)    ! x%var(:)
+            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
+
+      ! define the statistics index
+      iStat = meta(iVar)%statIndex(iFreq)
+
+      ! check that the variable is desired
+      if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle
+
+      ! stats output: only scalar variable type
+      if(meta(iVar)%varType==iLookVarType%scalarv) then
+        select type(stat)
+          class is (var_dlength)
+            select case(trim(structName))
+            case('forc')
+              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%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
+            case('diag')
+              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%gru(indxGRU)%hru(indxHRU)%var(map(iVar))%tim(iStep)%dat(iFreq) = stat%var(map(iVar))%dat(iFreq)
+            case('indx')
+              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
+          class default; err=20; message=trim(message)//'stats must be scalarv and of type gru_hru_doubleVec'; return
+        end select  ! stat
+
+        ! non-scalar variables: regular data structures
+      else
+
+        ! get the model layers
+        nSoil   = indx%var(iLookIndex%nSoil)%dat(1)
+        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%gru(indxGRU)%hru(indxHRU)%var(iLookIndex%nSnow)%tim(iStep)%dat(1) = nSnow
+        nLayers = indx%var(iLookIndex%nLayers)%dat(1)
+        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)
+          case(iLookVarType%wLength); datLength = maxSpectral
+          case(iLookVarType%midToto); datLength = nLayers
+          case(iLookVarType%midSnow); datLength = nSnow
+          case(iLookVarType%midSoil); datLength = nSoil
+          case(iLookVarType%ifcToto); datLength = nLayers+1
+          case(iLookVarType%ifcSnow); datLength = nSnow+1
+          case(iLookVarType%ifcSoil); datLength = nSoil+1
+          case default; cycle
+        end select ! vartype
+      
+        ! get the data vectors
+        select type (dat)
+          class is (var_dlength)
+            select case(trim(structName))
+              case('prog')
+                outputStructure(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = dat%var(iVar)%dat(:)
+              case('diag')
+                outputStructure(1)%diagStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = dat%var(iVar)%dat(:)
+              case('flux')
+                outputStructure(1)%fluxStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = 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%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = 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
+
+        ! get the maximum length of each data vector
+        select case (meta(iVar)%varType)
+          case(iLookVarType%wLength); maxLength = maxSpectral
+          case(iLookVarType%midToto); maxLength = maxLayers
+          case(iLookVarType%midSnow); maxLength = maxLayers-nSoil
+          case(iLookVarType%midSoil); maxLength = nSoil
+          case(iLookVarType%ifcToto); maxLength = maxLayers+1
+          case(iLookVarType%ifcSnow); maxLength = (maxLayers-nSoil)+1
+          case(iLookVarType%ifcSoil); maxLength = nSoil+1
+          case default; cycle
+        end select ! vartype
+
+      end if ! not scalarv
+
+        ! process error code
+      if (err/=0) message=trim(message)//trim(meta(iVar)%varName)//'_'//trim(get_statName(iStat))
+        call netcdf_err(err,message); if (err/=0) return
+
+    end do ! iVar
+  end do ! iFreq
+end subroutine writeData
+
+! **************************************************************************************
+! public subroutine writeBasin: write basin-average variables
+! **************************************************************************************
+subroutine writeBasin(indxGRU,indxHRU,iStep,finalizeStats,&
+                      outputTimestep,meta,stat,dat,map,err,message)
+ USE data_types,only:var_info                       ! metadata type
+ USE var_lookup,only:maxVarStat                     ! index into stats structure
+ USE var_lookup,only:iLookVarType                   ! index into type structure
+ USE globalData,only:outFreq                   ! output file information
+ USE get_ixName_module,only:get_varTypeName         ! to access type strings for error messages
+ USE get_ixName_module,only:get_statName            ! to access type strings for error messages
+ implicit none
+
+ ! declare dummy variables
+ integer(i4b)  ,intent(in)     :: indxGRU
+ integer(i4b)  ,intent(in)     :: indxHRU
+ integer(i4b)  ,intent(in)     :: iStep
+ logical(lgt)  ,intent(in)     :: finalizeStats(:)  ! flags to finalize statistics
+ integer(i4b)  ,intent(in)     :: outputTimestep(:) ! output time step
+ type(var_info),intent(in)     :: meta(:)           ! meta data
+ type(dlength) ,intent(in)     :: stat(:)           ! stats data
+ type(dlength) ,intent(in)     :: dat(:)            ! timestep data
+ integer(i4b)  ,intent(in)     :: map(:)            ! map into stats child struct
+ integer(i4b)  ,intent(out)    :: err               ! error code
+ character(*)  ,intent(out)    :: message           ! error message
+ ! local variables
+ integer(i4b)                  :: iVar              ! variable index
+ integer(i4b)                  :: iStat             ! statistics index
+ integer(i4b)                  :: iFreq             ! frequency index
+ ! initialize error control
+ err=0;message="hru_modelwrite.f90-writeBasin/"
+
+ ! loop through output frequencies
+ do iFreq=1,maxvarFreq
+
+  ! skip frequencies that are not needed
+  if(.not.outFreq(iFreq)) cycle
+
+  ! check that we have finalized statistics for a given frequency
+  if(.not.finalizeStats(iFreq)) cycle
+
+  ! loop through model variables
+  do iVar = 1,size(meta)
+
+   ! define the statistics index
+   iStat = meta(iVar)%statIndex(iFreq)
+
+   ! check that the variable is desired
+   if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle
+
+   ! stats/data output - select data type
+   select case (meta(iVar)%varType)
+
+    case (iLookVarType%scalarv)
+      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%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(iFreq) = dat(iVar)%dat(iFreq)
+     end if
+
+    case default
+     err=40; message=trim(message)//"unknownVariableType[name='"//trim(meta(iVar)%varName)//"';type='"//trim(get_varTypeName(meta(iVar)%varType))//    "']"; return
+   end select ! variable type
+
+   ! process error code
+   if (err.ne.0) message=trim(message)//trim(meta(iVar)%varName)//'_'//trim(get_statName(iStat))
+   if (err/=0) return
+
+  end do ! iVar
+ end do ! iFreq
+
+end subroutine writeBasin
+
+ ! **************************************************************************************
+ ! public subroutine writeTime: write current time to all files
+ ! **************************************************************************************
+subroutine writeTime(indxGRU,indxHRU,iStep,finalizeStats,meta,dat,err,message)
+ USE data_types,only:var_info                       ! metadata type
+ USE var_lookup,only:iLookStat                      ! index into stat structure
+ implicit none
+
+ ! declare dummy variables
+ integer(i4b)  ,intent(in)     :: indxGRU
+ integer(i4b)  ,intent(in)     :: indxHRU
+ integer(i4b)  ,intent(in)     :: iStep
+ logical(lgt)  ,intent(in)     :: finalizeStats(:)  ! flags to finalize statistics
+ type(var_info),intent(in)     :: meta(:)           ! meta data
+ integer       ,intent(in)     :: dat(:)            ! timestep data
+ integer(i4b)  ,intent(out)    :: err               ! error code
+ character(*)  ,intent(out)    :: message           ! error message
+ ! local variables
+ integer(i4b)                  :: iVar              ! variable index
+ integer(i4b)                  :: iFreq             ! frequency index
+ ! initialize error control
+ err=0;message="hru_modelwrite.f90-writeTime/"
+
+ ! loop through output frequencies
+ do iFreq=1,maxvarFreq
+
+  ! check that we have finalized statistics for a given frequency
+  if(.not.finalizeStats(iFreq)) cycle
+
+  ! loop through model variables
+  do iVar = 1,size(meta)
+
+   ! check instantaneous
+   if (meta(iVar)%statIndex(iFreq)/=iLookStat%inst) cycle
+
+   ! add to outputStructure
+   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
+
+  end do ! iVar
+ end do ! iFreq
+
+end subroutine writeTime
+
+! *********************************************************************************************************
+! public subroutine printRestartFile: print a re-start file
+! *********************************************************************************************************
+subroutine writeRestart(filename,         & ! intent(in): name of restart file
+                         nGRU,             & ! intent(in): number of GRUs
+                         nHRU,             & ! intent(in): number of HRUs
+                         prog_meta,        & ! intent(in): prognostics metadata
+                         prog_data,        & ! intent(in): prognostics data
+                         bvar_meta,        & ! intent(in): basin (gru) variable metadata
+                         bvar_data,        & ! intent(in): basin (gru) variable data
+                         maxLayers,        & ! intent(in): maximum number of layers
+                         maxSnowLayers,    & ! intent(in): maximum number of snow layers
+                         indx_meta,        & ! intent(in): index metadata
+                         indx_data,        & ! intent(in): index data
+                         err,message)        ! intent(out): error control
+ ! --------------------------------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------------------------------
+ ! access the derived types to define the data structures
+ USE data_types,only:var_info               ! metadata
+ ! access named variables defining elements in the data structures
+ USE var_lookup,only:iLookINDEX             ! named variables for structure elements
+ USE var_lookup,only:iLookVarType           ! named variables for structure elements
+ USE var_lookup,only:iLookBVAR              ! named variables for structure elements
+ ! constants
+ USE globalData,only:gru_struc              ! gru-hru mapping structures
+ ! external routines
+ USE netcdf_util_module,only:nc_file_close  ! close netcdf file
+ USE netcdf_util_module,only:nc_file_open   ! open netcdf file
+ USE globalData,only:nTimeDelay             ! number of timesteps in the time delay histogram
+ 
+ implicit none
+ ! --------------------------------------------------------------------------------------------------------
+ ! input
+ character(len=256),intent(in)      :: filename      ! name of the restart file
+ integer(i4b),intent(in)            :: nGRU          ! number of GRUs
+ integer(i4b),intent(in)            :: nHRU          ! number of HRUs
+ type(var_info),intent(in)          :: prog_meta(:)  ! prognostic variable metadata
+ type(var_dlength),intent(in)       :: prog_data     ! prognostic vars
+ type(var_info),intent(in)          :: bvar_meta(:)  ! basin variable metadata
+ type(var_dlength),intent(in)       :: bvar_data     ! basin variables
+ type(var_info),intent(in)          :: indx_meta(:)  ! metadata
+ type(var_ilength),intent(in)       :: indx_data     ! indexing vars
+ ! output: error control
+ integer(i4b),intent(out)           :: err           ! error code
+ character(*),intent(out)           :: message       ! error message
+ ! --------------------------------------------------------------------------------------------------------
+ ! dummy variables
+ integer(i4b), intent(in)           :: maxLayers     ! maximum number of total layers
+ integer(i4b), intent(in)           :: maxSnowLayers ! maximum number of snow layers
+
+ ! local variables
+ integer(i4b)                       :: ncid          ! netcdf file id
+ integer(i4b),allocatable           :: ncVarID(:)    ! netcdf variable id
+ integer(i4b)                       :: ncSnowID      ! index variable id
+ integer(i4b)                       :: ncSoilID      ! index variable id
+
+ integer(i4b)                       :: nSoil         ! number of soil layers
+ integer(i4b)                       :: nSnow         ! number of snow layers
+ integer(i4b)                       :: maxSnow       ! maximum number of snow layers
+ integer(i4b)                       :: maxSoil       ! maximum number of soil layers
+ integer(i4b)                       :: nLayers       ! number of total layers
+ integer(i4b),parameter             :: nSpectral=2   ! number of spectal bands
+ integer(i4b),parameter             :: nScalar=1     ! size of a scalar
+ integer(i4b)                       :: nProgVars     ! number of prognostic variables written to state file
+
+ integer(i4b)                       :: hruDimID      ! variable dimension ID
+ integer(i4b)                       :: gruDimID      ! variable dimension ID
+ integer(i4b)                       :: tdhDimID      ! variable dimension ID
+ integer(i4b)                       :: scalDimID     ! variable dimension ID
+ integer(i4b)                       :: specDimID     ! variable dimension ID
+ integer(i4b)                       :: midSnowDimID  ! variable dimension ID
+ integer(i4b)                       :: midSoilDimID  ! variable dimension ID
+ integer(i4b)                       :: midTotoDimID  ! variable dimension ID
+ integer(i4b)                       :: ifcSnowDimID  ! variable dimension ID
+ integer(i4b)                       :: ifcSoilDimID  ! variable dimension ID
+ integer(i4b)                       :: ifcTotoDimID  ! variable dimension ID
+
+ character(len=32),parameter        :: hruDimName    ='hru'      ! dimension name for HRUs
+ character(len=32),parameter        :: gruDimName    ='gru'      ! dimension name for GRUs
+ character(len=32),parameter        :: tdhDimName    ='tdh'      ! dimension name for time-delay basin variables
+ character(len=32),parameter        :: scalDimName   ='scalarv'  ! dimension name for scalar data
+ character(len=32),parameter        :: specDimName   ='spectral' ! dimension name for spectral bands
+ character(len=32),parameter        :: midSnowDimName='midSnow'  ! dimension name for snow-only layers
+ 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        :: ifcSnowDimName='ifcSnow'  ! dimension name for snow-only layers
+ character(len=32),parameter        :: ifcSoilDimName='ifcSoil'  ! dimension name for soil-only layers
+ character(len=32),parameter        :: ifcTotoDimName='ifcToto'  ! dimension name for layered variables
+
+ integer(i4b)                       :: cHRU          ! count of HRUs
+ integer(i4b)                       :: iHRU          ! index of HRUs
+ integer(i4b)                       :: iGRU          ! index of GRUs
+ integer(i4b)                       :: iVar          ! variable index
+ logical(lgt)                       :: okLength      ! flag to check if the vector length is OK
+ character(len=256)                 :: cmessage      ! downstream error message
+ ! --------------------------------------------------------------------------------------------------------
+
+ ! initialize error control
+ err=0; message='writeRestart/'
+
+ ! size of prognostic variable vector
+ nProgVars = size(prog_meta)
+ allocate(ncVarID(nProgVars+1))     ! include 1 additional basin variable in ID array (possibly more later)
+
+ ! maximum number of soil layers
+ maxSoil = gru_struc(1)%hruInfo(1)%nSoil
+
+ ! maximum number of snow layers
+ maxSnow = maxSnowLayers
+ 
+ ! create file
+ !TODO: Verify if this is needed
+ err = nf90_create(trim(filename),nf90_classic_model,ncid)
+ message='iCreate[create]'; call netcdf_err(err,message); if(err/=0)return
+
+ ! define dimensions
+ !TODO: Verify if this is needed
+                err = nf90_def_dim(ncid,trim(hruDimName)    ,nHRU       ,    hruDimID); message='iCreate[hru]'     ; call netcdf_err(err,message); if(err/=0)return
+                err = nf90_def_dim(ncid,trim(gruDimName)    ,nGRU       ,    gruDimID); message='iCreate[gru]'     ; call netcdf_err(err,message); if(err/=0)return
+                err = nf90_def_dim(ncid,trim(tdhDimName)    ,nTimeDelay ,    tdhDimID); message='iCreate[tdh]'     ; call netcdf_err(err,message); if(err/=0)return
+                err = nf90_def_dim(ncid,trim(scalDimName)   ,nScalar    ,   scalDimID); message='iCreate[scalar]'  ; call netcdf_err(err,message); if(err/=0)return
+                err = nf90_def_dim(ncid,trim(specDimName)   ,nSpectral  ,   specDimID); message='iCreate[spectral]'; call netcdf_err(err,message); if(err/=0)return
+                err = nf90_def_dim(ncid,trim(midSoilDimName),maxSoil    ,midSoilDimID); message='iCreate[ifcSoil]' ; call netcdf_err(err,message); if(err/=0)return
+                err = nf90_def_dim(ncid,trim(midTotoDimName),maxLayers  ,midTotoDimID); message='iCreate[midToto]' ; call netcdf_err(err,message); if(err/=0)return
+                err = nf90_def_dim(ncid,trim(ifcSoilDimName),maxSoil+1  ,ifcSoilDimID); message='iCreate[ifcSoil]' ; call netcdf_err(err,message); if(err/=0)return
+                err = nf90_def_dim(ncid,trim(ifcTotoDimName),maxLayers+1,ifcTotoDimID); message='iCreate[ifcToto]' ; call netcdf_err(err,message); if(err/=0)return
+ if (maxSnow>0) err = nf90_def_dim(ncid,trim(midSnowDimName),maxSnow    ,midSnowDimID); message='iCreate[ifcSnow]' ; call netcdf_err(err,message); if(err/=0)return
+ if (maxSnow>0) err = nf90_def_dim(ncid,trim(ifcSnowDimName),maxSnow+1  ,ifcSnowDimID); message='iCreate[ifcSnow]' ; call netcdf_err(err,message); if(err/=0)return
+ ! re-initialize error control
+ err=0; message='writeRestart/'
+
+ ! define prognostic variables
+ do iVar = 1,nProgVars
+  if (prog_meta(iVar)%varType==iLookvarType%unknown) cycle
+
+  ! define variable
+  select case(prog_meta(iVar)%varType)
+    !TODO: Verify if this is needed
+   case(iLookvarType%scalarv);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,  scalDimID /),ncVarID(iVar))
+   case(iLookvarType%wLength);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,  specDimID /),ncVarID(iVar))
+   case(iLookvarType%midSoil);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,midSoilDimID/),ncVarID(iVar))
+   case(iLookvarType%midToto);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,midTotoDimID/),ncVarID(iVar))
+   case(iLookvarType%ifcSoil);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,ifcSoilDimID/),ncVarID(iVar))
+   case(iLookvarType%ifcToto);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,ifcTotoDimID/),ncVarID(iVar))
+   case(iLookvarType%midSnow); if (maxSnow>0) err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,midSnowDimID/),ncVarID(iVar))
+   case(iLookvarType%ifcSnow); if (maxSnow>0) err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,ifcSnowDimID/),ncVarID(iVar))
+  end select
+
+  ! check errors
+  if(err/=0)then
+   message=trim(message)//trim(cmessage)//' [variable '//trim(prog_meta(iVar)%varName)//']'
+   return
+  end if
+
+  ! add parameter description
+  !TODO: Verify if this is needed
+  err = nf90_put_att(ncid,ncVarID(iVar),'long_name',trim(prog_meta(iVar)%vardesc))
+  call netcdf_err(err,message)
+
+  ! add parameter units
+  !TODO: Verify if this is needed
+  err = nf90_put_att(ncid,ncVarID(iVar),'units',trim(prog_meta(iVar)%varunit))
+  call netcdf_err(err,message)
+
+ end do ! iVar
+ 
+ ! define selected basin variables (derived) -- e.g., hillslope routing
+ !TODO: Verify if this is needed
+ err = nf90_def_var(ncid, trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varName), nf90_double, (/gruDimID, tdhDimID /), ncVarID(nProgVars+1))
+ err = nf90_put_att(ncid,ncVarID(nProgVars+1),'long_name',trim(bvar_meta(iLookBVAR%routingRunoffFuture)%vardesc));   call netcdf_err(err,message)
+ err = nf90_put_att(ncid,ncVarID(nProgVars+1),'units'    ,trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varunit));   call netcdf_err(err,message)
+
+ ! define index variables - snow
+ !TODO: Verify if this is needed
+ err = nf90_def_var(ncid,trim(indx_meta(iLookIndex%nSnow)%varName),nf90_int,(/hruDimID/),ncSnowID); call netcdf_err(err,message)
+ err = nf90_put_att(ncid,ncSnowID,'long_name',trim(indx_meta(iLookIndex%nSnow)%vardesc));           call netcdf_err(err,message)
+ err = nf90_put_att(ncid,ncSnowID,'units'    ,trim(indx_meta(iLookIndex%nSnow)%varunit));           call netcdf_err(err,message)
+
+ ! define index variables - soil
+ !TODO: Verify if this is needed
+ err = nf90_def_var(ncid,trim(indx_meta(iLookIndex%nSoil)%varName),nf90_int,(/hruDimID/),ncSoilID); call netcdf_err(err,message)
+ err = nf90_put_att(ncid,ncSoilID,'long_name',trim(indx_meta(iLookIndex%nSoil)%vardesc));           call netcdf_err(err,message)
+ err = nf90_put_att(ncid,ncSoilID,'units'    ,trim(indx_meta(iLookIndex%nSoil)%varunit));           call netcdf_err(err,message)
+
+ ! end definition phase
+ !TODO: Verify if this is needed
+ err = nf90_enddef(ncid); call netcdf_err(err,message); if (err/=0) return
+
+ ! write variables
+ do iGRU = 1,nGRU
+  do iHRU = 1,gru_struc(iGRU)%hruCount
+   cHRU = gru_struc(iGRU)%hruInfo(iHRU)%hru_ix
+   do iVar = 1,size(prog_meta)
+
+    ! excape if this variable is not used
+    if (prog_meta(iVar)%varType==iLookvarType%unknown) cycle
+
+    ! actual number of layers
+    nSnow = gru_struc(iGRU)%hruInfo(iHRU)%nSnow
+    nSoil = gru_struc(iGRU)%hruInfo(iHRU)%nSoil
+    nLayers = nSoil + nSnow
+
+    ! check size
+    ! NOTE: this may take time that we do not wish to use
+    okLength=.true.
+    select case (prog_meta(iVar)%varType)
+    ! case(iLookVarType%scalarv);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nScalar  )
+    ! case(iLookVarType%wlength);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSpectral)
+    ! case(iLookVarType%midSoil);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSoil    )
+    ! case(iLookVarType%midToto);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nLayers  )
+    ! case(iLookVarType%ifcSoil);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSoil+1  )
+    ! case(iLookVarType%ifcToto);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nLayers+1)
+    ! case(iLookVarType%midSnow); if (nSnow>0) okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSnow    )
+    ! case(iLookVarType%ifcSnow); if (nSnow>0) okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSnow+1  )
+     case(iLookVarType%scalarv);              okLength = (size(prog_data%var(iVar)%dat) == nScalar  )
+     case(iLookVarType%wlength);              okLength = (size(prog_data%var(iVar)%dat) == nSpectral)
+     case(iLookVarType%midSoil);              okLength = (size(prog_data%var(iVar)%dat) == nSoil    )
+     case(iLookVarType%midToto);              okLength = (size(prog_data%var(iVar)%dat) == nLayers  )
+     case(iLookVarType%ifcSoil);              okLength = (size(prog_data%var(iVar)%dat) == nSoil+1  )
+     case(iLookVarType%ifcToto);              okLength = (size(prog_data%var(iVar)%dat) == nLayers+1)
+     case(iLookVarType%midSnow); if (nSnow>0) okLength = (size(prog_data%var(iVar)%dat) == nSnow    )
+     case(iLookVarType%ifcSnow); if (nSnow>0) okLength = (size(prog_data%var(iVar)%dat) == nSnow+1  )
+     case default; err=20; message=trim(message)//'unknown var type'; return
+    end select
+
+    ! error check
+    if(.not.okLength)then
+     message=trim(message)//'bad vector length for variable '//trim(prog_meta(iVar)%varname)
+     err=20; return
+    endif
+
+    ! write data
+    select case (prog_meta(iVar)%varType)
+     case(iLookVarType%scalarv);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nScalar  /))
+     case(iLookVarType%wlength);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSpectral/))
+     case(iLookVarType%midSoil);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSoil    /))
+     case(iLookVarType%midToto);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nLayers  /))
+     case(iLookVarType%ifcSoil);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSoil+1  /))
+     case(iLookVarType%ifcToto);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nLayers+1/))
+     case(iLookVarType%midSnow); if (nSnow>0) err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSnow    /))
+     case(iLookVarType%ifcSnow); if (nSnow>0) err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSnow+1  /))
+     case default; err=20; message=trim(message)//'unknown var type'; return
+    end select
+
+    ! error check
+    if (err.ne.0) message=trim(message)//'writing variable:'//trim(prog_meta(iVar)%varName)
+    call netcdf_err(err,message); if (err/=0) return
+    err=0; message='writeRestart/'
+
+   end do ! iVar loop
+
+   ! write index variables
+   !TODO: Verify if this is needed
+   err=nf90_put_var(ncid,ncSnowID,(/indx_data%var(iLookIndex%nSnow)%dat/),start=(/cHRU/),count=(/1/))
+   err=nf90_put_var(ncid,ncSoilID,(/indx_data%var(iLookIndex%nSoil)%dat/),start=(/cHRU/),count=(/1/))
+
+  end do ! iHRU loop
+  
+  ! write selected basin variables
+ !TODO: Verify if this is needed
+  err=nf90_put_var(ncid,ncVarID(nProgVars+1),(/bvar_data%var(iLookBVAR%routingRunoffFuture)%dat/),  start=(/iGRU/),count=(/1,nTimeDelay/))
+
+ end do  ! iGRU loop
+
+ ! close file
+ call nc_file_close(ncid,err,cmessage)
+ if(err/=0)then;message=trim(message)//trim(cmessage);return;end if
+
+ ! cleanup
+ deallocate(ncVarID)
+
+end subroutine writeRestart
 
 end module HRUwriteoOutput_module
\ No newline at end of file