diff --git a/build/includes/file_access_actor/file_access_actor.hpp b/build/includes/file_access_actor/file_access_actor.hpp
index 32dc0cfed22f0744946ed75ffaba0bf8a4776185..068d1e5f9a4d34c12c18fce47e343d042eb33cdc 100644
--- a/build/includes/file_access_actor/file_access_actor.hpp
+++ b/build/includes/file_access_actor/file_access_actor.hpp
@@ -29,6 +29,9 @@ extern "C" {
   void writeOutput_fortran(void* handle_ncid, int* num_steps, int* start_gru, 
       int* max_gru, bool* writeParamFlag, int* err);
 
+  void writeRestart_fortran(void* handle_ncid, int* start_gru, int* max_gru, int* timestep, 
+      int* year, int* month, int* day, int* hour, int* err); 
+
   void read_forcingFile(void* forcFileInfo, int* currentFile, int* stepsInFile,
       int* startGRU, int* numGRU, int* err);
 
@@ -59,6 +62,7 @@ struct file_access_state {
   int numFiles;
   int filesLoaded;
   int num_output_steps;
+  int err = 0; // this is to make compiler happy
 
   Output_Container* output_container;
 
@@ -72,6 +76,12 @@ struct file_access_state {
 
 
   bool write_params_flag = true;
+
+  // Checkpointing variables
+  int completed_checkpoints = 1;  // 
+  std::vector<int> hru_checkpoints;
+  std::vector<int> hru_timesteps;
+
 };
 
 // called to spawn a file_access_actor
@@ -87,5 +97,8 @@ behavior file_access_actor(stateful_actor<file_access_state>* self,
 void writeOutput(stateful_actor<file_access_state>* self, 
     Output_Partition* partition);
 
+void writeRestart(stateful_actor<file_access_state>* self, Output_Partition* partition, 
+    int start_gru, int num_gru, int timestep, int year, int month, int day, int hour);
+
  
 } // end namespace
\ No newline at end of file
diff --git a/build/includes/global/message_atoms.hpp b/build/includes/global/message_atoms.hpp
index f0854ffc23c6db24b2b54ae092abbb6850344ff6..d08ce36f8b4692342454f2fd33fc3e7d518537d4 100644
--- a/build/includes/global/message_atoms.hpp
+++ b/build/includes/global/message_atoms.hpp
@@ -247,6 +247,10 @@ CAF_BEGIN_TYPE_ID_BLOCK(summa, first_custom_type_id)
     // Reciever:
     // Summary:
     CAF_ADD_ATOM(summa, write_output)
+    // Sender: HRU Actor
+    // Reciever: File Access Actor 
+    // Summary: Updates FAA when hru reaches restart checkpoint
+    CAF_ADD_ATOM(summa, write_restart)
     // Sender:
     // Reciever:
     // Summary:
diff --git a/build/includes/hru_actor/hru_actor.hpp b/build/includes/hru_actor/hru_actor.hpp
index 4467993c76d83dca9ac3cac70ce35fb91089a0d9..7240d5d38b2d22970d4078c643a6e1055570c0ac 100644
--- a/build/includes/hru_actor/hru_actor.hpp
+++ b/build/includes/hru_actor/hru_actor.hpp
@@ -9,6 +9,15 @@
 #include "message_atoms.hpp"
 #include "global.hpp"
 
+
+struct Date {
+  int y;
+  int m;
+  int d;
+  int h;
+};
+
+
 /*********************************************
  * HRU Actor Fortran Functions
  *********************************************/
@@ -27,9 +36,10 @@ extern "C" {
   void RunPhysics(int* id, int* stepIndex, void* hru_data, double* dt, 
       int* dt_int_factor, double* walltime_timestep, int* err);
   
-  void hru_writeOutput(int* index_hru, int* index_gru, int* timestep, 
-      int* output_step, void* hru_data, int* err);
-  
+  void hru_writeOutput(int* index_hru, int* index_gru, int* timestep, int* output_step, void* hru_data, int* y, int* m, int* d, int* h, int* err);
+    
+  int  hru_writeRestart(int* index_hru, int* index_gru, int* timestep, int* output_step, void* hru_data, int* err);
+
   void setTimeZoneOffset(int* iFile, void* hru_data, int* err);
 
   void HRU_readForcing(int* index_gru, int* iStep, int* iRead, int* iFile, 
@@ -89,7 +99,15 @@ struct hru_state {
   double rtol = -9999; // -9999 uses default
   double atol = -9999; // -9999 uses default
 
-  double walltime_timestep = 0.0; // walltime for the current timestep		        
+  double walltime_timestep = 0.0; // walltime for the current timestep		
+
+  // Checkpointing variables
+  Date startDate =   {0,0,0,0}; // will be initalized when hru finishes first timestep
+  Date currentDate = {0,0,0,0}; // will be initalized when hru finishes first timestep
+  int  checkpoint= 0;
+  int daysInMonth[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
+  int restartFrequency; // 0=never(default) 1=hour 2=day 3=week? 4=month 5=year 
+        
 
   // Settings
   HRU_Actor_Settings hru_actor_settings;
@@ -113,5 +131,8 @@ void Initialize_HRU(stateful_actor<hru_state>* self);
 /** Function runs all of the hru time_steps */
 int Run_HRU(stateful_actor<hru_state>* self);
 
+//** Function checks if the HRU is at a restart checkpoint **/
+bool isCheckpoint(stateful_actor<hru_state>* self);
+
 
 }
\ No newline at end of file
diff --git a/build/source/file_access_actor/cppwrap_fileAccess.f90 b/build/source/file_access_actor/cppwrap_fileAccess.f90
index 50b4b689398e7f7ed670ad967135a370adbbc9fe..f99ca0e861c02e9a0eb845fb656513d78b6de9e0 100644
--- a/build/source/file_access_actor/cppwrap_fileAccess.f90
+++ b/build/source/file_access_actor/cppwrap_fileAccess.f90
@@ -67,7 +67,7 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   USE var_lookup,only:iLookPARAM
   USE var_lookup,only:iLookATTR                               ! look-up values for model attributes
   USE var_lookup,only:iLookBVAR                               ! look-up values for basin-average variables
-  USE output_structure_module,only:outputStructure            ! output structure
+  USE output_structure_module,only:summa_struct            ! output structure
 
   USE globalData,only:iRunModeFull,iRunModeGRU,iRunModeHRU
   USE globalData,only:iRunMode                                ! define the current running mode
@@ -217,8 +217,8 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   ! *****************************************************************************
 
   attrFile = trim(SETTINGS_PATH)//trim(LOCAL_ATTRIBUTES)
-  call read_attrb(trim(attrFile),num_gru,outputStructure(1)%attrStruct,&
-      outputStructure(1)%typeStruct,outputStructure(1)%idStruct,err,message)
+  call read_attrb(trim(attrFile),num_gru,summa_struct(1)%attrStruct,&
+      summa_struct(1)%typeStruct,summa_struct(1)%idStruct,err,message)
   if(err/=0)then; print*,trim(message); return; endif
 
 
@@ -226,12 +226,12 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   do iGRU=1, num_gru
     do iHRU=1, gru_struc(iGRU)%hruCount
       ! set parmameters to their default value
-      outputStructure(1)%dparStruct%gru(iGRU)%hru(iHRU)%var(:) = localParFallback(:)%default_val         ! x%hru(:)%var(:)
+      summa_struct(1)%dparStruct%gru(iGRU)%hru(iHRU)%var(:) = localParFallback(:)%default_val         ! x%hru(:)%var(:)
 
       ! overwrite default model parameters with information from the Noah-MP tables
-      call pOverwrite(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex),  &  ! vegetation category
-                      outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%soilTypeIndex), &  ! soil category
-                      outputStructure(1)%dparStruct%gru(iGRU)%hru(iHRU)%var,                          &  ! default model parameters
+      call pOverwrite(summa_struct(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex),  &  ! vegetation category
+                      summa_struct(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%soilTypeIndex), &  ! soil category
+                      summa_struct(1)%dparStruct%gru(iGRU)%hru(iHRU)%var,                          &  ! default model parameters
                       err,message)                                                   ! error control
       if(err/=0)then; print*, trim(message); return; endif
 
@@ -239,13 +239,13 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
    ! copy over to the parameter structure
    ! NOTE: constant for the dat(:) dimension (normally depth)
       do ivar=1,size(localParFallback)
-        outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(ivar)%dat(:) = outputStructure(1)%dparStruct%gru(iGRU)%hru(iHRU)%var(ivar)
+        summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(ivar)%dat(:) = summa_struct(1)%dparStruct%gru(iGRU)%hru(iHRU)%var(ivar)
       end do  ! looping through variables
     
     end do  ! looping through HRUs
     
     ! set default for basin-average parameters
-    outputStructure(1)%bparStruct%gru(iGRU)%var(:) = basinParFallback(:)%default_val
+    summa_struct(1)%bparStruct%gru(iGRU)%var(:) = basinParFallback(:)%default_val
     
   end do  ! looping through GRUs
 
@@ -254,8 +254,8 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   ! *** Read Parameters
   ! *****************************************************************************
   checkHRU = integerMissing
-  call read_param(iRunMode,checkHRU,start_gru,num_hru,num_gru,outputStructure(1)%idStruct,&
-      outputStructure(1)%mparStruct,outputStructure(1)%bparStruct,err,message)
+  call read_param(iRunMode,checkHRU,start_gru,num_hru,num_gru,summa_struct(1)%idStruct,&
+      summa_struct(1)%mparStruct,summa_struct(1)%bparStruct,err,message)
   if(err/=0)then; print*,trim(message); return; endif
 
   ! *****************************************************************************
@@ -264,8 +264,8 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   ! ! loop through GRUs
   do iGRU=1,num_gru
     ! calculate the fraction of runoff in future time steps
-    call fracFuture(outputStructure(1)%bparStruct%gru(iGRU)%var,     &  ! vector of basin-average model parameters
-                    outputStructure(1)%bvarStruct_init%gru(iGRU),    &  ! data structure of basin-average variables
+    call fracFuture(summa_struct(1)%bparStruct%gru(iGRU)%var,     &  ! vector of basin-average model parameters
+                    summa_struct(1)%bvarStruct_init%gru(iGRU),    &  ! data structure of basin-average variables
                     err,message)                   ! error control
     if(err/=0)then; print*, trim(message); return; endif
 
@@ -276,7 +276,7 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
       kHRU=0
       ! check the network topology (only expect there to be one downslope HRU)
       do jHRU=1,gru_struc(iGRU)%hruCount
-      if(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%downHRUindex) == outputStructure(1)%idStruct%gru(iGRU)%hru(jHRU)%var(iLookID%hruId))then
+      if(summa_struct(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%downHRUindex) == summa_struct(1)%idStruct%gru(iGRU)%hru(jHRU)%var(iLookID%hruId))then
         if(kHRU==0)then  ! check there is a unique match
         kHRU=jHRU
         else
@@ -287,13 +287,13 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
 
 
       ! check that the parameters are consistent
-      call paramCheck(outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU),err,message)
+      call paramCheck(summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU),err,message)
       if(err/=0)then; print*, message; return; endif
 
 #ifdef V4_ACTIVE      
       ! calculate a look-up table for the temperature-enthalpy conversion of snow for future snow layer merging
       ! NOTE2: H is the mixture enthalpy of snow liquid and ice
-      call T2H_lookup_snow(outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU),err,message)
+      call T2H_lookup_snow(summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU),err,message)
       if(err/=0)then; print*, message; return; endif
 
       ! calculate a lookup table for the temperature-enthalpy conversion of soil
@@ -301,45 +301,45 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
       !       multiply by Cp_liq*iden_water to get temperature component of enthalpy
       if(needLookup)then
         call T2L_lookup_soil(gru_struc(iGRU)%hruInfo(iHRU)%nSoil,   &   ! intent(in):    number of soil layers
-                        outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU),        &   ! intent(in):    parameter data structure
-                        outputStructure(1)%lookupStruct%gru(iGRU)%hru(iHRU),      &   ! intent(inout): lookup table data structure
+                        summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU),        &   ! intent(in):    parameter data structure
+                        summa_struct(1)%lookupStruct%gru(iGRU)%hru(iHRU),      &   ! intent(inout): lookup table data structure
                         err,message)                              ! intent(out):   error control
         if(err/=0)then; print*, message; return; endif
       endif
 else
       ! calculate a look-up table for the temperature-enthalpy conversion
-      call E2T_lookup(outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU),err,message)
+      call E2T_lookup(summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU),err,message)
       if(err/=0)then; message=trim(message); print*, message; return; endif
 
 #endif
       ! overwrite the vegetation height
-      HVT(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex)) = outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%heightCanopyTop)%dat(1)
-      HVB(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex)) = outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%heightCanopyBottom)%dat(1)
+      HVT(summa_struct(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex)) = summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%heightCanopyTop)%dat(1)
+      HVB(summa_struct(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex)) = summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%heightCanopyBottom)%dat(1)
          
       ! overwrite the tables for LAI and SAI
       if(model_decisions(iLookDECISIONS%LAI_method)%iDecision == specified)then
-        SAIM(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex),:) = outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%winterSAI)%dat(1)
-        LAIM(outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex),:) = outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%summerLAI)%dat(1)*greenVegFrac_monthly
+        SAIM(summa_struct(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex),:) = summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%winterSAI)%dat(1)
+        LAIM(summa_struct(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookTYPE%vegTypeIndex),:) = summa_struct(1)%mparStruct%gru(iGRU)%hru(iHRU)%var(iLookPARAM%summerLAI)%dat(1)*greenVegFrac_monthly
       endif
 
     end do ! HRU
     
     ! compute total area of the upstream HRUS that flow into each HRU
     do iHRU=1,gru_struc(iGRU)%hruCount
-      outputStructure(1)%upArea%gru(iGRU)%hru(iHRU) = 0._rkind
+      summa_struct(1)%upArea%gru(iGRU)%hru(iHRU) = 0._rkind
       do jHRU=1,gru_struc(iGRU)%hruCount
        ! check if jHRU flows into iHRU; assume no exchange between GRUs
-       if(outputStructure(1)%typeStruct%gru(iGRU)%hru(jHRU)%var(iLookTYPE%downHRUindex)==outputStructure(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookID%hruId))then
-        outputStructure(1)%upArea%gru(iGRU)%hru(iHRU) = outputStructure(1)%upArea%gru(iGRU)%hru(iHRU) + outputStructure(1)%attrStruct%gru(iGRU)%hru(jHRU)%var(iLookATTR%HRUarea)
+       if(summa_struct(1)%typeStruct%gru(iGRU)%hru(jHRU)%var(iLookTYPE%downHRUindex)==summa_struct(1)%typeStruct%gru(iGRU)%hru(iHRU)%var(iLookID%hruId))then
+        summa_struct(1)%upArea%gru(iGRU)%hru(iHRU) = summa_struct(1)%upArea%gru(iGRU)%hru(iHRU) + summa_struct(1)%attrStruct%gru(iGRU)%hru(jHRU)%var(iLookATTR%HRUarea)
        endif   ! (if jHRU is an upstream HRU)
       end do  ! jHRU
     end do  ! iHRU
   
     ! identify the total basin area for a GRU (m2)  
-    outputStructure(1)%bvarStruct_init%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) = 0._rkind
+    summa_struct(1)%bvarStruct_init%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) = 0._rkind
     do iHRU=1,gru_struc(iGRU)%hruCount
-      outputStructure(1)%bvarStruct_init%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) = &
-      outputStructure(1)%bvarStruct_init%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) + outputStructure(1)%attrStruct%gru(iGRU)%hru(iHRU)%var(iLookATTR%HRUarea)
+      summa_struct(1)%bvarStruct_init%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) = &
+      summa_struct(1)%bvarStruct_init%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) + summa_struct(1)%attrStruct%gru(iGRU)%hru(iHRU)%var(iLookATTR%HRUarea)
     end do
   
   end do ! GRU
@@ -362,17 +362,17 @@ else
  ! read initial conditions
   call read_icond(restartFile,                        & ! intent(in):    name of initial conditions file
                   num_gru,                            & ! intent(in):    number of response units
-                  outputStructure(1)%mparStruct,      & ! intent(in):    model parameters
-                  outputStructure(1)%progStruct_init, & ! intent(inout): model prognostic variables
-                  outputStructure(1)%bvarStruct_init, & ! intent(inout): model basin (GRU) variables
-                  outputStructure(1)%indxStruct_init, & ! intent(inout): model indices
+                  summa_struct(1)%mparStruct,      & ! intent(in):    model parameters
+                  summa_struct(1)%progStruct_init, & ! intent(inout): model prognostic variables
+                  summa_struct(1)%bvarStruct_init, & ! intent(inout): model basin (GRU) variables
+                  summa_struct(1)%indxStruct_init, & ! intent(inout): model indices
                   err,message)                          ! intent(out):   error control
   if(err/=0)then; print*, message; return; endif
 
   call check_icond(num_gru,                            &
-                   outputStructure(1)%progStruct_init, &  ! intent(inout): model prognostic variables
-                   outputStructure(1)%mparStruct,      & ! intent(in):    model parameters
-                   outputStructure(1)%indxStruct_init, & ! intent(inout): model indices
+                   summa_struct(1)%progStruct_init, &  ! intent(inout): model prognostic variables
+                   summa_struct(1)%mparStruct,      & ! intent(in):    model parameters
+                   summa_struct(1)%indxStruct_init, & ! intent(inout): model indices
                    err,message)                          ! intent(out):   error control
   if(err/=0)then; print*, message; return; endif  
 end subroutine fileAccessActor_init_fortran
diff --git a/build/source/file_access_actor/fileAccess_writeOutput.f90 b/build/source/file_access_actor/fileAccess_writeOutput.f90
index 3f21df782a05d7b35003e509bba79acdc631700d..d9745cf2dd0df5758dd91766779b41d41326e1f8 100644
--- a/build/source/file_access_actor/fileAccess_writeOutput.f90
+++ b/build/source/file_access_actor/fileAccess_writeOutput.f90
@@ -34,7 +34,7 @@ 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
+USE output_structure_module,only:summa_struct
 USE output_structure_module,only:outputTimeStep
 
 
@@ -127,14 +127,14 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, write
       do iGRU=start_gru, max_gru
         select case(trim(structInfo(iStruct)%structName))
         case('attr'); call writeParm(ncid,gru_struc(iGRU)%hruInfo(indxHRU)%hru_ix, &
-          outputStructure(1)%attrStruct%gru(iGRU)%hru(indxHRU),attr_meta,err,cmessage)
+          summa_struct(1)%attrStruct%gru(iGRU)%hru(indxHRU),attr_meta,err,cmessage)
         case('type'); call writeParm(ncid,gru_struc(iGRU)%hruInfo(indxHRU)%hru_ix, &
-          outputStructure(1)%typeStruct%gru(iGRU)%hru(indxHRU),type_meta,err,cmessage)
+          summa_struct(1)%typeStruct%gru(iGRU)%hru(indxHRU),type_meta,err,cmessage)
         case('mpar'); call writeParm(ncid,gru_struc(iGRU)%hruInfo(indxHRU)%hru_ix, &
-          outputStructure(1)%mparStruct%gru(iGRU)%hru(indxHRU),mpar_meta,err,cmessage)
+          summa_struct(1)%mparStruct%gru(iGRU)%hru(indxHRU),mpar_meta,err,cmessage)
         end select
         if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
-        call writeParm(ncid,iGRU,outputStructure(1)%bparStruct%gru(iGRU),bpar_meta,err,cmessage)
+        call writeParm(ncid,iGRU,summa_struct(1)%bparStruct%gru(iGRU),bpar_meta,err,cmessage)
         if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
       end do ! GRU
     end do ! structInfo
@@ -150,7 +150,7 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, write
     stepCounter(:) = outputTimeStep(iGRU)%dat(:) ! We want to avoid updating outputTimeStep
     do iStep=1, num_steps
       call writeTime(ncid,outputTimeStep(iGRU)%dat(:),iStep,time_meta, &
-              outputStructure(1)%timeStruct%gru(iGRU)%hru(indxHRU)%var,err,cmessage)
+              summa_struct(1)%timeStruct%gru(iGRU)%hru(indxHRU)%var,err,cmessage)
     end do ! istep
   end do ! iGRU
 
@@ -161,7 +161,7 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, write
   ! ****************************************************************************
   call writeBasin(ncid,outputTimeStep(start_gru)%dat(:),outputTimeStepUpdate,num_steps,&
                   start_gru, max_gru, numGRU, &
-                  bvar_meta,outputStructure(1)%bvarStat,outputStructure(1)%bvarStruct, &
+                  bvar_meta,summa_struct(1)%bvarStat,summa_struct(1)%bvarStruct, &
                   bvarChild_map,err,cmessage)
 
   ! ****************************************************************************
@@ -172,28 +172,28 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, write
       case('forc')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, & 
-                        forc_meta,outputStructure(1)%forcStat,outputStructure(1)%forcStruct,'forc', &
-                        forcChild_map,outputStructure(1)%indxStruct,err,cmessage)
+                        forc_meta,summa_struct(1)%forcStat,summa_struct(1)%forcStruct,'forc', &
+                        forcChild_map,summa_struct(1)%indxStruct,err,cmessage)
       case('prog')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, &
-                        prog_meta,outputStructure(1)%progStat,outputStructure(1)%progStruct,'prog', &
-                        progChild_map,outputStructure(1)%indxStruct,err,cmessage)
+                        prog_meta,summa_struct(1)%progStat,summa_struct(1)%progStruct,'prog', &
+                        progChild_map,summa_struct(1)%indxStruct,err,cmessage)
       case('diag')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, &
-                        diag_meta,outputStructure(1)%diagStat,outputStructure(1)%diagStruct,'diag', &
-                        diagChild_map,outputStructure(1)%indxStruct,err,cmessage)
+                        diag_meta,summa_struct(1)%diagStat,summa_struct(1)%diagStruct,'diag', &
+                        diagChild_map,summa_struct(1)%indxStruct,err,cmessage)
       case('flux')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, &
-                        flux_meta,outputStructure(1)%fluxStat,outputStructure(1)%fluxStruct,'flux', &
-                        fluxChild_map,outputStructure(1)%indxStruct,err,cmessage)
+                        flux_meta,summa_struct(1)%fluxStat,summa_struct(1)%fluxStruct,'flux', &
+                        fluxChild_map,summa_struct(1)%indxStruct,err,cmessage)
       case('indx')
         call writeData(ncid,outputTimeStep(start_gru)%dat(:),outputTimestepUpdate,maxLayers,num_steps,&
                         start_gru, max_gru, numGRU, &
-                        indx_meta,outputStructure(1)%indxStat,outputStructure(1)%indxStruct,'indx', &
-                        indxChild_map,outputStructure(1)%indxStruct,err,cmessage)
+                        indx_meta,summa_struct(1)%indxStat,summa_struct(1)%indxStruct,'indx', &
+                        indxChild_map,summa_struct(1)%indxStruct,err,cmessage)
     end select
     if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
   end do  ! (looping through structures)
@@ -206,6 +206,89 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, write
 end subroutine writeOutput_fortran
 
 
+subroutine writeRestart_fortran(handle_ncid,  start_gru, num_gru, checkpoint, year, month, day, hour, err) bind(C, name="writeRestart_fortran")
+  USE var_lookup,only:maxVarFreq                               ! # of available output frequencies
+  USE globalData,only:structInfo
+  USE globalData,only:bvarChild_map,forcChild_map,progChild_map,diagChild_map,fluxChild_map,indxChild_map             ! index of the child data structure: stats bvar
+  USE globalData,only:attr_meta,bvar_meta,type_meta,time_meta,forc_meta,prog_meta,diag_meta,flux_meta,indx_meta,bpar_meta,mpar_meta
+  USE globalData,only:maxLayers                               ! maximum number of layers
+  USE globalData,only:maxSnowLayers                           ! maximum number of snow layers
+
+  USE summaFileManager,only:OUTPUT_PATH,OUTPUT_PREFIX         ! define output file
+  USE summaFileManager,only:STATE_PATH                        ! optional path to state output files (defaults to OUTPUT_PATH)
+
+  ! USE fileAccess_writeRestart,only:writeRestart_fortran
+  implicit none
+  ! dummy variables
+  type(c_ptr),intent(in), value        :: handle_ncid       ! ncid of the output file
+  ! integer(c_int),intent(in)            :: num_steps         ! number of steps to write
+  integer(c_int),intent(in)            :: start_gru         ! index of GRU we are currently writing for
+  integer(c_int),intent(in)            :: num_gru           ! index of HRU we are currently writing for
+  
+  integer(c_int),intent(in)            :: checkpoint           ! slowest timestep of all grus in job
+  integer(c_int),intent(in)            :: year 
+  integer(c_int),intent(in)            :: month
+  integer(c_int),intent(in)            :: day
+  integer(c_int),intent(in)            :: hour
+  ! logical(c_bool),intent(in)           :: write_parm_flag   ! flag to write parameters
+  integer(c_int),intent(out)           :: err               ! Error code
+  ! local variables
+  type(var_i),pointer                  :: ncid
+  integer(i4b)                         :: iGRU              ! loop through GRUs
+  integer(i4b)                         :: iStep             ! loop through time steps
+  integer(i4b)                         :: iFreq             ! loop through output frequencies
+  integer(i4b)                         :: indxHRU=1         ! index of HRU to write
+  integer(i4b), dimension(maxVarFreq)  :: outputTimestepUpdate
+  integer(i4b), dimension(maxVarFreq)  :: stepCounter
+  character(LEN=256)                   :: message
+  character(LEN=256)                   :: cmessage
+  
+  character (len = 11) :: output_fileSuffix
+  character(len=256)                    :: restartFile       ! restart file name
+  character(len=256)                    :: timeString        ! portion of restart file name that contains the write-out time
+  integer(i4b)                          :: restart_flag
+
+
+  integer(i4b)                         :: iStruct
+  integer(i4b)                         :: numGRU
+  
+  ! Change the C pointer to a fortran pointer
+  call c_f_pointer(handle_ncid, ncid)
+
+! *****************************************************************************
+  ! *** write restart file
+  ! *****************************************************************************
+  restart_flag = 1 ! temp
+  ! print a restart file if requested
+  if( restart_flag == 1 )then ! temp bare bones check
+    write(timeString,        '(I4.4,I2.2,I2.2,I2.2)') year,month,day,hour
+    write(output_fileSuffix, '(I5.5,"-",I5.5)') start_gru, start_gru + num_gru - 1
+
+    if(STATE_PATH == '') then
+      restartFile=trim(OUTPUT_PATH)//trim(OUTPUT_PREFIX)//'_restart_'//trim(timeString)//"_G"//output_fileSuffix//'.nc'
+    else
+      restartFile= trim(STATE_PATH)//trim(OUTPUT_PREFIX)//'_restart_'//trim(timeString)//"_G"//trim(output_fileSuffix)//'.nc'
+    endif
+
+    call writeRestart(restartFile,                   &  ! filename
+                      num_gru,                       &  ! nHRU
+                      checkpoint,                    &  ! checkpoint
+                      prog_meta,                     &  ! prog_meta
+                      summa_struct(1)%progStruct, &  ! prog_data
+                      bvar_meta,                     &  ! bvar_meta
+                      summa_struct(1)%bvarStruct, &  ! bvar_data
+                      maxLayers,                     &  ! maxLayers
+                      maxSnowLayers,                 &  ! maxSnowLayers
+                      indx_meta,                     &  ! indx_meta
+                      summa_struct(1)%indxStruct, &  ! indx_data
+                      err,                           &  ! err
+                      cmessage)                         ! message 
+    if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+  end if
+
+end subroutine writeRestart_fortran
+
+
 ! **********************************************************************************************************
 ! public subroutine writeParm: write model parameters
 ! **********************************************************************************************************
@@ -322,7 +405,6 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,nSteps,
   do iFreq=1,maxvarFreq
     ! skip frequencies that are not needed
     if(.not.outFreq(iFreq)) cycle
-
     ! loop through model variables
     do iVar = 1,size(meta)
       stepCounter = 0
@@ -333,15 +415,12 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,nSteps,
         call netcdf_err(err,message); if (err/=0) return
         
         do iStep = 1, nSteps
-          ! Find HRU that is not missing or NaN
-          ! check if we want this timestep
-          do iGRU = minGRU, maxGRU
-            if(.not.outputStructure(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
-            val = outputStructure(1)%forcStruct%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)
-            exit
-          end do
+
+          if(.not.summa_struct(1)%finalizeStats%gru(minGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+
           stepCounter = stepCounter+1
-          timeVec(stepCounter) = val
+          timeVec(stepCounter) = summa_struct(1)%forcStruct%gru(minGRU)%hru(1)%var(iVar)%tim(iStep)
+
         end do ! iStep
         err = nf90_put_var(ncid%var(iFreq),ncVarID,timeVec(1:stepCounter),start=(/outputTimestep(iFreq)/),count=(/stepCounter/))
         call netcdf_err(err,message); if (err/=0)then; print*, "err"; return; endif
@@ -411,13 +490,14 @@ subroutine writeScalar(ncid, outputTimestep, outputTimestepUpdate, nSteps, minGR
         stepCounter = 0
         gruCounter = gruCounter + 1
         do iStep = 1, nSteps
-          if(.not.outputStructure(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) then
-            val = realMissing
-          else
-            val = stat%gru(iGRU)%hru(1)%var(map(iVar))%tim(iStep)%dat(iFreq)
-          end if
+          
+          if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+
+
           stepCounter = stepCounter + 1
-          realVec(gruCounter, stepCounter) = val
+
+          realVec(gruCounter, stepCounter) = stat%gru(iGRU)%hru(1)%var(map(iVar))%tim(iStep)%dat(iFreq)
+
           outputTimeStepUpdate(iFreq) = stepCounter
         end do ! iStep
       end do ! iGRU 
@@ -515,11 +595,11 @@ subroutine writeVector(ncid, outputTimestep, maxLayers, nSteps, minGRU, maxGRU,
       ! get the data vectors
       select type (dat)
           class is (gru_hru_time_doubleVec)
-              if(.not.outputStructure(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+              if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
               realArray(gruCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(1:datLength)
 
           class is (gru_hru_time_intVec)
-              if(.not.outputStructure(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+              if(.not.summa_struct(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
               intArray(gruCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(1:datLength)
           class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return
       end select
@@ -656,7 +736,7 @@ subroutine writeTime(ncid,outputTimestep,iStep,meta,dat,err,message)
   do iFreq=1,maxvarFreq
 
     ! check that we have finalized statistics for a given frequency
-    if(.not.outputStructure(1)%finalizeStats%gru(1)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+    if(.not.summa_struct(1)%finalizeStats%gru(1)%hru(1)%tim(iStep)%dat(iFreq)) cycle
 
     ! loop through model variables
     do iVar = 1,size(meta)
@@ -678,4 +758,280 @@ subroutine writeTime(ncid,outputTimestep,iStep,meta,dat,err,message)
 
 
 end subroutine writeTime   
+
+subroutine writeRestart(filename,          & ! intent(in): name of restart file
+  ! minGRU, &
+  ! maxGRU, &
+  ! nGRU,             & ! intent(in): number of GRUs
+  nGRU,             & ! intent(in): number of HRUs
+  checkpoint, &
+  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
+!  integer(i4b)  ,intent(in)        :: minGRU            ! minGRU index to write
+!  integer(i4b)  ,intent(in)        :: maxGRU            ! maxGRU index to write - probably not needed
+integer(i4b),intent(in)            :: checkpoint      ! checkpoint the restart file is based on
+type(var_info),intent(in)          :: prog_meta(:)  ! prognostic variable metadata
+type(gru_hru_time_doubleVec),intent(in) :: prog_data     ! prognostic vars
+type(var_info),intent(in)          :: bvar_meta(:)  ! basin variable metadata
+type(gru_hru_time_doubleVec),intent(in)     :: bvar_data     ! basin variables
+type(var_info),intent(in)          :: indx_meta(:)  ! metadata
+type(gru_hru_time_intVec),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)                       :: scalar_val
+integer(i4b)                       :: varID
+
+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 = gru_struc(1)%hruInfo(1)%nSnow
+
+! create file
+err = nf90_create(trim(filename),nf90_classic_model,ncid)
+message='iCreate[create]'; call netcdf_err(err,message); if(err/=0)return
+
+
+! define dimensions
+err = nf90_def_dim(ncid,trim(hruDimName)    ,nGRU       ,    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)
+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
+err = nf90_put_att(ncid,ncVarID(iVar),'long_name',trim(prog_meta(iVar)%vardesc))
+call netcdf_err(err,message)
+
+! add parameter units
+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
+err = nf90_def_var(ncid, trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varName), nf90_double, (/hruDimID, 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
+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
+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
+err = nf90_enddef(ncid); call netcdf_err(err,message); if (err/=0) return
+
+! write variables
+do iGRU = 1,nGRU
+do iHRU = 1,1
+cHRU = iGRU!gru_struc(iGRU)%hruInfo(iHRU)%hru_ix
+do iVar = 1,nProgVars
+! excape if this variable is not used
+if (prog_meta(iVar)%varType==iLookvarType%unknown) cycle
+
+! actual number of layers
+nSnow = 5! TEMP gru_struc(iGRU)%hruInfo(iHRU)%nSnow
+nSoil = gru_struc(iHRU)%hruInfo(1)%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) ! should this be prog_meta?
+case(iLookVarType%scalarv);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat) == nScalar  )
+case(iLookVarType%wlength);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat) == nSpectral)
+case(iLookVarType%midSoil);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat) == nSoil    )
+case(iLookVarType%midToto);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat) == nLayers  )
+case(iLookVarType%ifcSoil);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat) == nSoil+1  )
+case(iLookVarType%ifcToto);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat) == nLayers+1)
+case(iLookVarType%midSnow); if (nSnow>0) okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat) == nSnow    )
+case(iLookVarType%ifcSnow); if (nSnow>0) okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%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%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat/),start=(/cHRU,1/),count=(/1,nScalar  /))
+case(iLookVarType%wlength);              
+err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat/),start=(/cHRU,1/),count=(/1,nSpectral/))
+case(iLookVarType%midSoil);              
+err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat/),start=(/cHRU,1/),count=(/1,nSoil    /))
+case(iLookVarType%midToto);              
+err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat/),start=(/cHRU,1/),count=(/1,nLayers - nSnow  /))
+case(iLookVarType%ifcSoil);              
+err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat/),start=(/cHRU,1/),count=(/1,nSoil+1  /))
+case(iLookVarType%ifcToto);              
+err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat/),start=(/cHRU,1/),count=(/1,nLayers - nSnow +1/))
+case(iLookVarType%midSnow); 
+if (nSnow>0) err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%dat/),start=(/cHRU,1/),count=(/1,nSnow    /))
+case(iLookVarType%ifcSnow); 
+if (nSnow>0) err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%tim(checkpoint)%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) then 
+return
+end if
+err=0; message='writeRestart/'
+
+end do ! iVar loop
+
+! write index variables
+err=nf90_put_var(ncid,ncSnowID,(/gru_struc(iGRU)%hruInfo(iHRU)%nSnow/),start=(/cHRU/),count=(/1/))
+err=nf90_put_var(ncid,ncSoilID,(/gru_struc(iGRU)%hruInfo(iHRU)%nSoil/),start=(/cHRU/),count=(/1/))
+
+! write selected basin variables
+err=nf90_put_var(ncid,ncVarID(nProgVars+1),(/bvar_data%gru(iGRU)%hru(iHRU)%var(iLookBVAR%routingRunoffFuture)%tim(checkpoint)%dat/),  start=(/iGRU/),count=(/1,nTimeDelay/))
+end do ! iHRU loop
+
+
+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 fileAccess_writeOutput
\ No newline at end of file
diff --git a/build/source/file_access_actor/file_access_actor.cpp b/build/source/file_access_actor/file_access_actor.cpp
index e715bcaab4cf1abc74a82c54ae5e5be2c31cbfe3..0b2ededc077e37daedb6a79a44ab59cecf87305e 100644
--- a/build/source/file_access_actor/file_access_actor.cpp
+++ b/build/source/file_access_actor/file_access_actor.cpp
@@ -71,6 +71,13 @@ behavior file_access_actor(stateful_actor<file_access_state>* self,
         fa_settings.num_timesteps_in_output_buffer, self->state.num_steps);
   }
 
+  // initialize vecs for keeping track of each hru's timestep and checkpoint
+  for (int i = 1; i < self->state.num_gru +1; ++i) {
+      self->state.hru_timesteps.push_back(0);
+      self->state.hru_checkpoints.push_back(0);
+  } 
+
+
   self->state.file_access_timing.updateEndPoint("init_duration");
 
   return {
@@ -171,6 +178,27 @@ behavior file_access_actor(stateful_actor<file_access_state>* self,
       self->state.file_access_timing.updateEndPoint("write_duration");
     },
 
+    [=] (write_restart, int gru, int gru_timestep, int gru_checkpoint, int output_stucture_index, int year, int month, int day, int hour){
+    // update hru progress vecs 
+    int gru_index = abs(self->state.start_gru - gru); 
+    self->state.hru_timesteps[gru_index] = gru_timestep;
+    self->state.hru_checkpoints[gru_index] = gru_checkpoint;
+
+    // find slowest time step of all hrus in job, stored in self->state.hru_timesteps
+    int slowest_timestep = *std::min_element(self->state.hru_timesteps.begin(), self->state.hru_timesteps.end());  
+    int slowest_checkpoint = *std::min_element(self->state.hru_checkpoints.begin(), self->state.hru_checkpoints.end());  
+
+    // if the slowest hru is past the ith checkpoint (current threshold)            
+    if ( slowest_checkpoint >= (self->state.completed_checkpoints)){// temp for dubuging
+        Output_Partition *output_partition = self->state.output_container->getOutputPartition(gru-1);
+        writeRestart(self, output_partition, self->state.start_gru, self->state.hru_timesteps.size(), output_stucture_index,
+                     year, month, day, hour);
+        // update checkpint counter
+        self->state.completed_checkpoints++;
+    }
+},
+
+
     // Write message from the job actor TODO: This could be async
     [=](write_output, int steps_to_write, int start_gru, int max_gru) {
       self->state.file_access_timing.updateStartPoint("write_duration");
@@ -257,4 +285,13 @@ void writeOutput(stateful_actor<file_access_state>* self,
   partition->resetReadyToWriteList();
 }
 
+void writeRestart(stateful_actor<file_access_state>* self , Output_Partition* partition, 
+  int start_gru, int num_gru, int timestep,
+  int year, int month, int day, int hour){  
+  
+  writeRestart_fortran(self->state.handle_ncid, &start_gru, &num_gru, &timestep, 
+    &year, &month, &day, &hour, &self->state.err);
+  }
+
+
 } // end namespace
\ No newline at end of file
diff --git a/build/source/file_access_actor/output_structure.f90 b/build/source/file_access_actor/output_structure.f90
index 7bef7cc4df5acf1d33d3cdd69385b2ed8e40f8fa..c7da30bd56e529048b4ff4cfa509de4f4b3e977d 100644
--- a/build/source/file_access_actor/output_structure.f90
+++ b/build/source/file_access_actor/output_structure.f90
@@ -105,7 +105,7 @@ module output_structure_module
   end type summa_output_type  
   
 
-  type(summa_output_type),allocatable,save,public     :: outputStructure(:)            ! summa_OutputStructure(1)%struc%var(:)%dat(nTimeSteps) 
+  type(summa_output_type),allocatable,save,public     :: summa_struct(:)            ! summa_OutputStructure(1)%struc%var(:)%dat(nTimeSteps) 
   type(ilength),allocatable,save,public               :: outputTimeStep(:)             ! timestep in output files
   
   contains
@@ -172,8 +172,8 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
 
 
   ! Allocate structure to hold output files
-  if (.not.allocated(outputStructure))then
-    allocate(outputStructure(1))
+  if (.not.allocated(summa_struct))then
+    allocate(summa_struct(1))
   else
     print*, "Already Allocated"; return;
   end if
@@ -181,58 +181,58 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
   ! LookupStructure
 
   ! Statistics Structures
-  allocate(outputStructure(1)%forcStat%gru(num_gru))
-  allocate(outputStructure(1)%progStat%gru(num_gru))
-  allocate(outputStructure(1)%diagStat%gru(num_gru))
-  allocate(outputStructure(1)%fluxStat%gru(num_gru))
-  allocate(outputStructure(1)%indxStat%gru(num_gru))
-  allocate(outputStructure(1)%bvarStat%gru(num_gru))
+  allocate(summa_struct(1)%forcStat%gru(num_gru))
+  allocate(summa_struct(1)%progStat%gru(num_gru))
+  allocate(summa_struct(1)%diagStat%gru(num_gru))
+  allocate(summa_struct(1)%fluxStat%gru(num_gru))
+  allocate(summa_struct(1)%indxStat%gru(num_gru))
+  allocate(summa_struct(1)%bvarStat%gru(num_gru))
   ! Primary Data Structures (scalars)
-  allocate(outputStructure(1)%timeStruct%gru(num_gru))
-  allocate(outputStructure(1)%forcStruct%gru(num_gru))
+  allocate(summa_struct(1)%timeStruct%gru(num_gru))
+  allocate(summa_struct(1)%forcStruct%gru(num_gru))
   ! Primary Data Structures (variable length vectors)
-  allocate(outputStructure(1)%indxStruct%gru(num_gru))
-  allocate(outputStructure(1)%progStruct%gru(num_gru))
-  allocate(outputStructure(1)%diagStruct%gru(num_gru))
-  allocate(outputStructure(1)%fluxStruct%gru(num_gru))
+  allocate(summa_struct(1)%indxStruct%gru(num_gru))
+  allocate(summa_struct(1)%progStruct%gru(num_gru))
+  allocate(summa_struct(1)%diagStruct%gru(num_gru))
+  allocate(summa_struct(1)%fluxStruct%gru(num_gru))
   ! Basin-Average structures
-  allocate(outputStructure(1)%bvarStruct%gru(num_gru))
+  allocate(summa_struct(1)%bvarStruct%gru(num_gru))
   ! Finalize Stats for writing
-  allocate(outputStructure(1)%finalizeStats%gru(num_gru))
+  allocate(summa_struct(1)%finalizeStats%gru(num_gru))
   ! Extras
-  allocate(outputStructure(1)%upArea%gru(num_gru))
+  allocate(summa_struct(1)%upArea%gru(num_gru))
   
   
   do iGRU = 1, num_gru
     num_hru = gru_struc(iGRU)%hruCount
     ! Statistics Structures
-    allocate(outputStructure(1)%forcStat%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%progStat%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%diagStat%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%fluxStat%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%indxStat%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%bvarStat%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%forcStat%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%progStat%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%diagStat%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%fluxStat%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%indxStat%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%bvarStat%gru(iGRU)%hru(num_hru))
 
     ! Primary Data Structures (scalars)
-    allocate(outputStructure(1)%timeStruct%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%forcStruct%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%timeStruct%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%forcStruct%gru(iGRU)%hru(num_hru))
 
     ! Primary Data Structures (variable length vectors)
-    allocate(outputStructure(1)%indxStruct%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%progStruct%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%diagStruct%gru(iGRU)%hru(num_hru))
-    allocate(outputStructure(1)%fluxStruct%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%indxStruct%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%progStruct%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%diagStruct%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%fluxStruct%gru(iGRU)%hru(num_hru))
   
     ! Basin-Average structures
-    allocate(outputStructure(1)%bvarStruct%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%bvarStruct%gru(iGRU)%hru(num_hru))
 
 
    ! define the ancillary data structures
 
     ! Finalize Stats for writing
-    allocate(outputStructure(1)%finalizeStats%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%finalizeStats%gru(iGRU)%hru(num_hru))
 
-    allocate(outputStructure(1)%upArea%gru(iGRU)%hru(num_hru))
+    allocate(summa_struct(1)%upArea%gru(iGRU)%hru(num_hru))
 
   end do
 
@@ -242,25 +242,25 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
     select case(trim(structInfo(iStruct)%structName))
     case('time'); cycle;
     case('forc'); cycle;
-    case('attr'); call allocGlobal(attr_meta,  outputStructure(1)%attrStruct,  err, message)
-    case('type'); call allocGlobal(type_meta,  outputStructure(1)%typeStruct,  err, message)
-    case('id'  ); call allocGlobal(id_meta,    outputStructure(1)%idStruct,    err, message)
-    case('mpar'); call allocGlobal(mpar_meta,  outputStructure(1)%mparStruct,  err, message); 
-    case('indx'); call allocGlobal(indx_meta,  outputStructure(1)%indxStruct_init, err, message);
-    case('prog'); call allocGlobal(prog_meta,  outputStructure(1)%progStruct_init, err, message);
+    case('attr'); call allocGlobal(attr_meta,  summa_struct(1)%attrStruct,  err, message)
+    case('type'); call allocGlobal(type_meta,  summa_struct(1)%typeStruct,  err, message)
+    case('id'  ); call allocGlobal(id_meta,    summa_struct(1)%idStruct,    err, message)
+    case('mpar'); call allocGlobal(mpar_meta,  summa_struct(1)%mparStruct,  err, message); 
+    case('indx'); call allocGlobal(indx_meta,  summa_struct(1)%indxStruct_init, err, message);
+    case('prog'); call allocGlobal(prog_meta,  summa_struct(1)%progStruct_init, err, message);
     case('diag'); cycle;
     case('flux'); cycle;
-    case('bpar'); call allocGlobal(bpar_meta,outputStructure(1)%bparStruct    ,err, message);  ! basin-average params 
-    case('bvar'); call allocGlobal(bvar_meta,outputStructure(1)%bvarStruct_init,err,message);  ! basin-average variables
+    case('bpar'); call allocGlobal(bpar_meta,summa_struct(1)%bparStruct    ,err, message);  ! basin-average params 
+    case('bvar'); call allocGlobal(bvar_meta,summa_struct(1)%bvarStruct_init,err,message);  ! basin-average variables
     case('deriv'); cycle;
 #ifdef V4_ACTIVE     
-    case('lookup'); call allocGlobal(lookup_meta,outputStructure(1)%lookupStruct,err, message);
+    case('lookup'); call allocGlobal(lookup_meta,summa_struct(1)%lookupStruct,err, message);
 #endif
     end select
   end do
   
 
-  call allocGlobal(mpar_meta,outputStructure(1)%dparStruct,err,message);
+  call allocGlobal(mpar_meta,summa_struct(1)%dparStruct,err,message);
 
 
 
@@ -276,57 +276,57 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
         ! allocate space structures
           select case(trim(structInfo(iStruct)%structName))    
             case('time')
-              call alloc_outputStruc(time_meta,outputStructure(1)%timeStruct%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(time_meta,summa_struct(1)%timeStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,err=err,message=message)     ! model forcing data
             case('forc')
               ! Structure
-              call alloc_outputStruc(forc_meta,outputStructure(1)%forcStruct%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(forc_meta,summa_struct(1)%forcStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model forcing data
               ! Statistics
-              call alloc_outputStruc(statForc_meta(:)%var_info,outputStructure(1)%forcStat%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(statForc_meta(:)%var_info,summa_struct(1)%forcStat%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model forcing data
             case('attr'); cycle;
-              ! call alloc_outputStruc(attr_meta, outputStructure(1)%attrStruct(1)%gru(iGRU)%hru(iHRU),&
+              ! call alloc_outputStruc(attr_meta, summa_struct(1)%attrStruct(1)%gru(iGRU)%hru(iHRU),&
               !                         nSteps=1,err=err,message=message)
             case('type'); cycle;
             case('id'  ); cycle;
             case('mpar'); cycle;
             case('indx')
               ! Structure
-              call alloc_outputStruc(indx_meta,outputStructure(1)%indxStruct%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(indx_meta,summa_struct(1)%indxStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,str_name='indx',message=message);    ! model variables
               ! Statistics
-              call alloc_outputStruc(statIndx_meta(:)%var_info,outputStructure(1)%indxStat%gru(iGRU)%hru(1), &
+              call alloc_outputStruc(statIndx_meta(:)%var_info,summa_struct(1)%indxStat%gru(iGRU)%hru(1), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! index vars
             case('prog')
               ! Structure
-              call alloc_outputStruc(prog_meta,outputStructure(1)%progStruct%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model prognostic (state) variables
+              call alloc_outputStruc(prog_meta,summa_struct(1)%progStruct%gru(iGRU)%hru(iHRU), &
+                                      nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,str_name='prog',message=message);    ! model prognostic (state) variables
               ! Statistics
-              call alloc_outputStruc(statProg_meta(:)%var_info,outputStructure(1)%progStat%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model prognostic 
+              call alloc_outputStruc(statProg_meta(:)%var_info,summa_struct(1)%progStat%gru(iGRU)%hru(iHRU), &
+                                      nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,str_name='prog',message=message);    ! model prognostic 
             case('diag')
               ! Structure
-              call alloc_outputStruc(diag_meta,outputStructure(1)%diagStruct%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(diag_meta,summa_struct(1)%diagStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model diagnostic variables
               ! Statistics
-              call alloc_outputStruc(statDiag_meta(:)%var_info,outputStructure(1)%diagStat%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(statDiag_meta(:)%var_info,summa_struct(1)%diagStat%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model diagnostic
             case('flux')
               ! Structure
-              call alloc_outputStruc(flux_meta,outputStructure(1)%fluxStruct%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(flux_meta,summa_struct(1)%fluxStruct%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model fluxes
               ! Statistics
-              call alloc_outputStruc(statFlux_meta(:)%var_info,outputStructure(1)%fluxStat%gru(iGRU)%hru(iHRU), &
+              call alloc_outputStruc(statFlux_meta(:)%var_info,summa_struct(1)%fluxStat%gru(iGRU)%hru(iHRU), &
                                       nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model fluxes
             case('bpar'); cycle;
             case('bvar')
               ! Structure
-              call alloc_outputStruc(bvar_meta,outputStructure(1)%bvarStruct%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=0,nSoil=0,err=err,message=message);  ! basin-average variables
+              call alloc_outputStruc(bvar_meta,summa_struct(1)%bvarStruct%gru(iGRU)%hru(iHRU), &
+                                      nSteps=maxSteps,nSnow=0,nSoil=0,err=err,str_name='bvar',message=message);  ! basin-average variables
               ! Statistics
-              call alloc_outputStruc(statBvar_meta(:)%var_info,outputStructure(1)%bvarStat%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=0,nSoil=0,err=err,message=message);  ! basin-average variables
+              call alloc_outputStruc(statBvar_meta(:)%var_info,summa_struct(1)%bvarStat%gru(iGRU)%hru(iHRU), &
+                                      nSteps=maxSteps,nSnow=0,nSoil=0,err=err,str_name='bvar',message=message);  ! basin-average variables
             case('deriv');  cycle
             case('lookup'); cycle  ! lookup tables
             case default; err=20; message='unable to find structure name: '//trim(structInfo(iStruct)%structName)
@@ -341,9 +341,9 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
       end do  ! looping through data structures
     
       ! Finalize stats structure for writing to output file
-      allocate(outputStructure(1)%finalizeStats%gru(iGRU)%hru(iHRU)%tim(maxSteps))
+      allocate(summa_struct(1)%finalizeStats%gru(iGRU)%hru(iHRU)%tim(maxSteps))
       do iStep = 1, maxSteps
-        allocate(outputStructure(1)%finalizeStats%gru(iGRU)%hru(iHRU)%tim(iStep)%dat(1:maxVarFreq))
+        allocate(summa_struct(1)%finalizeStats%gru(iGRU)%hru(iHRU)%tim(iStep)%dat(1:maxVarFreq))
       end do ! timeSteps
     end do ! Looping through GRUs
   end do
@@ -357,7 +357,7 @@ subroutine deallocateOutputStructure(err) bind(C, name="deallocateOutputStructur
   integer(i4b), intent(inout)   :: err
 
   err = 0
-  deallocate(outputStructure)
+  deallocate(summa_struct)
 
 end subroutine deallocateOutputStructure
 
@@ -519,6 +519,7 @@ subroutine alloc_outputStruc(metaStruct,dataStruct,nSteps,nSnow,nSoil,str_name,e
   character(*),intent(out)             :: message        ! error message
   ! local
   logical(lgt)                         :: check          ! .true. if the variables are allocated
+  logical(lgt)                         :: allocAllFlag   ! .true. if struct is to have all timesteps allocated
   integer(i4b)                         :: nVars          ! number of variables in the metadata structure
   integer(i4b)                         :: nLayers        ! total number of layers
   integer(i4b)                         :: iVar
@@ -527,6 +528,11 @@ subroutine alloc_outputStruc(metaStruct,dataStruct,nSteps,nSnow,nSoil,str_name,e
   ! initalize error control
   message='alloc_outputStruc'
 
+  allocAllFlag = .false.
+  if (present(str_name)) then
+    allocAllFlag = .true.
+  end if
+
   nVars = size(metaStruct)
   if(present(nSnow) .or. present(nSoil))then
     ! check both are present
@@ -554,7 +560,7 @@ subroutine alloc_outputStruc(metaStruct,dataStruct,nSteps,nSnow,nSoil,str_name,e
       end if
       do iVar=1, nVars
         ! Check if this variable is desired within any timeframe
-        if(is_var_desired(metaStruct,iVar))then
+        if(is_var_desired(metaStruct,iVar) .or. allocAllFlag)then
           allocate(dataStruct%var(iVar)%tim(nSteps))
         end if
       end do
@@ -568,7 +574,7 @@ subroutine alloc_outputStruc(metaStruct,dataStruct,nSteps,nSnow,nSoil,str_name,e
       end if 
       do iVar=1, nVars
         ! Check if this variable is desired within any timeframe
-        if(is_var_desired(metaStruct,iVar))then
+        if(is_var_desired(metaStruct,iVar) .or. allocAllFlag)then
           allocate(dataStruct%var(iVar)%tim(nSteps))
         end if
       end do
@@ -582,7 +588,7 @@ subroutine alloc_outputStruc(metaStruct,dataStruct,nSteps,nSnow,nSoil,str_name,e
       end if
       do iVar=1, nVars
         ! Check if this variable is desired within any timeframe
-        if(is_var_desired(metaStruct,iVar))then
+        if(is_var_desired(metaStruct,iVar) .or. allocAllFlag)then
           allocate(dataStruct%var(iVar)%tim(nSteps))
         end if
       end do
@@ -628,7 +634,7 @@ subroutine alloc_outputStruc(metaStruct,dataStruct,nSteps,nSnow,nSoil,str_name,e
       end if
       do iVar=1, nVars
         ! Check if this variable is desired within any timeframe
-        if(is_var_desired(metaStruct,iVar) .or. (present(str_name) .and. &
+        if(is_var_desired(metaStruct,iVar) .or. allocAllFlag .or. (present(str_name) .and. &
          ((iVar == iLookINDEX%nLayers) .or. (iVar == iLookINDEX%nSnow) .or. (iVar == iLookINDEX%nSoil)) ))then
         allocate(dataStruct%var(iVar)%tim(nSteps))
           call allocateDat_int(metaStruct,dataStruct,nSnow,nSoil,nSteps,iVar,err,cmessage)
@@ -643,7 +649,7 @@ subroutine alloc_outputStruc(metaStruct,dataStruct,nSteps,nSnow,nSoil,str_name,e
       end if
       do iVar=1, nVars
         ! Check if this variable is desired within any timeframe
-        if(is_var_desired(metaStruct,iVar))then
+        if(is_var_desired(metaStruct,iVar) .or. allocAllFlag)then
           allocate(dataStruct%var(iVar)%tim(nSteps))
           call allocateDat_rkind_nSteps(metaStruct,dataStruct,nSnow,nSoil,nSteps,iVar,err,cmessage)
         end if
diff --git a/build/source/hru_actor/hru_actor.cpp b/build/source/hru_actor/hru_actor.cpp
index 37bd152084cf409db21d6c12c656e664f73f6ae8..061e313c773775396f4bb45cc9f592d1f974564a 100644
--- a/build/source/hru_actor/hru_actor.cpp
+++ b/build/source/hru_actor/hru_actor.cpp
@@ -28,6 +28,10 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
   self->state.hru_actor_settings = hru_actor_settings;
   self->state.dt_init_factor = hru_actor_settings.dt_init_factor;
 
+  // Set the restart frequency
+  self->state.restartFrequency = 2; // TODO: obtain this value from command line arg
+
+
   return {
     [=](init_hru) {
       Initialize_HRU(self);
@@ -407,20 +411,93 @@ int Run_HRU(stateful_actor<hru_state>* self) {
     return 20;
   }
 
-  hru_writeOutput(&self->state.indxHRU, &self->state.indxGRU,
-      &self->state.timestep, &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;
+  if (self->state.timestep != 0){
+    // write variables for output to summa_struct, extract y,m,d,h from 
+    // fortran side ald save it to the hru's state
+    int y,m,d,h;
+
+    hru_writeOutput(&self->state.indxHRU, &self->state.indxGRU,
+        &self->state.timestep, &self->state.output_structure_step_index,
+        self->state.hru_data, &y, &m, &d, &h, &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;
+    }
+
+    self->state.currentDate.y = y;
+    self->state.currentDate.m = m;
+    self->state.currentDate.d = d;
+    self->state.currentDate.h = h;
+
+    // collect date infomation on first timestep to find starting date
+    if (self->state.timestep == 1 ){
+        self->state.startDate = self->state.currentDate;
+    }
+
+    // check if hru reached a checkpoint, if so it will write restart data to outputsructure and
+    // send an update to the FAA
+    if (isCheckpoint(self)){
+        self->state.checkpoint++;
+
+        hru_writeRestart(&self->state.indxHRU, 
+                        &self->state.indxGRU,
+                        &self->state.output_structure_step_index,
+                        &self->state.output_structure_step_index, //unused
+                        self->state.hru_data,
+                        &self->state.err);
+                  
+        self->send(self->state.file_access_actor, 
+            write_restart_v,
+            self->state.refGRU,
+            self->state.timestep,
+            self->state.checkpoint,
+            self->state.output_structure_step_index,
+            self->state.currentDate.y,
+            self->state.currentDate.m,
+            self->state.currentDate.d,
+            self->state.currentDate.h);
+    }
+
+
   }
 
   return 0;      
 }
 
+// given a hru_actor with a state, compared current date with starting date to deterimine if hru is on a checkpoint
+bool isCheckpoint(stateful_actor<hru_state>* self){
+    switch(self->state.restartFrequency){
+        case 0: // restart not enabled
+            break;
+        case 1: // every timestep
+            return true;
+        case 2: // daily
+            if (self->state.startDate.h == self->state.currentDate.h){
+                return true;
+            }
+            break;
+        case 3: // weekly not supported
+            break;
+        case 4: // monthly
+            if (self->state.startDate.d == self->state.currentDate.d &&
+                self->state.startDate.h == self->state.currentDate.h){
+                return true;
+            }    
+            break;   
+        case 5: // yearly
+            if (self->state.startDate.m == self->state.currentDate.m &&
+                self->state.startDate.d == self->state.currentDate.d &&
+                self->state.startDate.h == self->state.currentDate.h){
+                return true;
+            }
+            break;
+    }
+    return false;
+}
+
 
 }
\ No newline at end of file
diff --git a/build/source/hru_actor/hru_init.f90 b/build/source/hru_actor/hru_init.f90
index 2adc3280eb04a967414c55b55a0224ae48769205..e1e99fee294d70178954a414b547a4a2b6f061bb 100755
--- a/build/source/hru_actor/hru_init.f90
+++ b/build/source/hru_actor/hru_init.f90
@@ -249,7 +249,7 @@ subroutine setupHRUParam(indxGRU,                 & ! ID of hru
   ! * desired modules
   ! ---------------------------------------------------------------------------------------
   USE nrtype                                                  ! variable types, etc.
-  USE output_structure_module,only:outputStructure
+  USE output_structure_module,only:summa_struct
   ! subroutines and functions
   use time_utils_module,only:elapsedSec                       ! calculate the elapsed time
   USE mDecisions_module,only:mDecisions                       ! module to read model decisions
@@ -304,39 +304,39 @@ subroutine setupHRUParam(indxGRU,                 & ! ID of hru
   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(:)
+  hru_data%attrStruct%var(:) = summa_struct(1)%attrStruct%gru(indxGRU)%hru(indxHRU)%var(:)
   ! Copy the typeStruct
-  hru_data%typeStruct%var(:) = outputStructure(1)%typeStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  hru_data%typeStruct%var(:) = summa_struct(1)%typeStruct%gru(indxGRU)%hru(indxHRU)%var(:)
   ! Copy the idStruct
-  hru_data%idStruct%var(:) = outputStructure(1)%idStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  hru_data%idStruct%var(:) = summa_struct(1)%idStruct%gru(indxGRU)%hru(indxHRU)%var(:)
 
   ! Copy the mparStruct
-  hru_data%mparStruct%var(:) = outputStructure(1)%mparStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  hru_data%mparStruct%var(:) = summa_struct(1)%mparStruct%gru(indxGRU)%hru(indxHRU)%var(:)
   ! Copy the bparStruct
-  hru_data%bparStruct%var(:) = outputStructure(1)%bparStruct%gru(indxGRU)%var(:)
+  hru_data%bparStruct%var(:) = summa_struct(1)%bparStruct%gru(indxGRU)%var(:)
   ! Copy the dparStruct
-  hru_data%dparStruct%var(:) = outputStructure(1)%dparStruct%gru(indxGRU)%hru(indxHRU)%var(:)
+  hru_data%dparStruct%var(:) = summa_struct(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(:)
+  do ivar=1, size(summa_struct(1)%bvarStruct_init%gru(indxGRU)%var(:))
+    hru_data%bvarStruct%var(ivar)%dat(:) = summa_struct(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(:)
+  if (allocated(summa_struct(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z)) then
+    do i_z=1, size(summa_struct(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z(:))
+      do iVar=1, size(summa_struct(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z(i_z)%var(:))
+        hru_data%lookupStruct%z(i_z)%var(ivar)%lookup(:) = summa_struct(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(:)
+  do ivar=1, size(summa_struct(1)%progStruct_init%gru(indxGRU)%hru(indxHRU)%var(:))
+    hru_data%progStruct%var(ivar)%dat(:) = summa_struct(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(:)
+  do ivar=1, size(summa_struct(1)%indxStruct_init%gru(indxGRU)%hru(indxHRU)%var(:))
+    hru_data%indxStruct%var(ivar)%dat(:) = summa_struct(1)%indxStruct_init%gru(indxGRU)%hru(indxHRU)%var(ivar)%dat(:)
   enddo
 end subroutine setupHRUParam
 
diff --git a/build/source/hru_actor/hru_writeOutput.f90 b/build/source/hru_actor/hru_writeOutput.f90
index 78a5c6915d95418fb48aef9b1e675bd713fda1c2..4e479d8127a4c3f4efaac10ac583e1d2cef190fd 100644
--- a/build/source/hru_actor/hru_writeOutput.f90
+++ b/build/source/hru_actor/hru_writeOutput.f90
@@ -57,7 +57,7 @@ USE var_lookup,only:iLookBVAR                 ! named variables for basin parame
 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
+USE output_structure_module,only:summa_struct
 
 
 implicit none
@@ -76,7 +76,7 @@ subroutine hru_writeOutput(&
                             indxGRU,                   &
                             timestep,                  & ! model timestep
                             outputStep,                & ! index into the output Struc
-                            handle_hru_data,           & ! local HRU data  
+                            handle_hru_data, y,m,d,h,  & ! local HRU data  
                             err) bind(C, name="hru_writeOutput") 
   USE nrtype
   USE globalData,only:structInfo
@@ -95,7 +95,7 @@ subroutine hru_writeOutput(&
   USE output_stats,only:calcStats                             ! module for compiling output statistics
   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
+  USE output_structure_module,only:summa_struct
   USE netcdf_util_module,only:nc_file_close                   ! close netcdf file
   USE netcdf_util_module,only:nc_file_open                    ! open netcdf file
   USE var_lookup,only:maxvarFreq                              ! maximum number of output files
@@ -107,6 +107,11 @@ subroutine hru_writeOutput(&
   integer(c_int),intent(in)             :: outputStep            ! index into the output Struc
   type(c_ptr),intent(in),value          :: handle_hru_data       ! local HRU data
   integer(c_int),intent(out)            :: err
+  integer(c_int),intent(out)            :: y
+  integer(c_int),intent(out)            :: m
+  integer(c_int),intent(out)            :: d
+  integer(c_int),intent(out)            :: h
+
 
   ! local pointers
   type(hru_type), pointer               :: hru_data              ! local HRU data
@@ -125,6 +130,13 @@ subroutine hru_writeOutput(&
   err=0; message='summa_manageOutputFiles/'
   ! identify the start of the writing
 
+  ! updating date variables to be passed back to the actors
+  y = hru_data%timeStruct%var(iLookTIME%iyyy)
+  m = hru_data%timeStruct%var(iLookTIME%im)
+  d = hru_data%timeStruct%var(iLookTIME%id)
+  h = hru_data%timeStruct%var(iLookTIME%ih)
+
+
   ! Many variables get there values from summa4chm_util.f90:getCommandArguments()
   call summa_setWriteAlarms(hru_data%oldTime_hru%var, hru_data%timeStruct%var, hru_data%finishTime_hru%var,  &   ! time vectors
                             newOutputFile,  defNewOutputFile,            &
@@ -152,7 +164,7 @@ subroutine hru_writeOutput(&
 
 
  ! If we do not do this looping we segfault - I am not sure why
-  outputStructure(1)%finalizeStats%gru(indxGRU)%hru(indxHRU)%tim(outputStep)%dat(:) = hru_data%finalizeStats%dat(:)
+  summa_struct(1)%finalizeStats%gru(indxGRU)%hru(indxHRU)%tim(outputStep)%dat(:) = hru_data%finalizeStats%dat(:)
 
  ! ****************************************************************************
  ! *** calculate output statistics
@@ -224,6 +236,109 @@ subroutine hru_writeOutput(&
 
 end subroutine hru_writeOutput
 
+
+subroutine hru_writeRestart(&
+  indxHRU,                   &
+  indxGRU,                   &
+  checkPoint,                & ! model checkPoint, index into the output Struc
+  outputStep,                & ! unused :(
+  handle_hru_data,           & ! local HRU data  
+  err) bind(C, name="hru_writeRestart")
+  USE nrtype
+  USE globalData,only:structInfo
+  USE globalData,only:startWrite,endWrite
+  USE globalData,only:maxLayers                               ! maximum number of layers
+  USE globalData,only:maxSnowLayers                           ! maximum number of snow layers
+
+  USE globalData,only:ixProgress                              ! define frequency to write progress
+  USE globalData,only:ixRestart                               ! define frequency to write restart files
+  USE globalData,only:gru_struc
+
+  USE globalData,only:newOutputFile                           ! define option for new output files
+  USE summa_alarms,only:summa_setWriteAlarms
+
+  USE summaFileManager,only:OUTPUT_PATH,OUTPUT_PREFIX         ! define output file
+  USE summaFileManager,only:STATE_PATH                        ! optional path to state output files (defaults to OUTPUT_PATH)
+
+
+  USE globalData,only:forc_meta,attr_meta,type_meta           ! metaData structures
+  USE output_stats,only:calcStats                             ! module for compiling output statistics
+  USE time_utils_module,only:elapsedSec                       ! calculate the elapsed time
+  USE globalData,only:elapsedWrite                            ! elapsed time to write data
+  USE output_structure_module,only:summa_struct
+  USE netcdf_util_module,only:nc_file_close                   ! close netcdf file
+  USE netcdf_util_module,only:nc_file_open                    ! open netcdf file
+  USE var_lookup,only:maxvarFreq                              ! maximum number of output files
+  USE var_lookup,only:iLookVarType           ! named variables for structure elements
+
+
+  implicit none
+  integer(c_int),intent(in)             :: indxHRU               ! index of hru in GRU
+  integer(c_int),intent(in)             :: indxGRU               ! index of the GRU
+  integer(c_int),intent(in)             :: checkPoint              ! model checkPoint
+  integer(c_int),intent(in)             :: outputStep            ! index into the output Struc
+  type(c_ptr),intent(in),value          :: handle_hru_data       ! local HRU data
+  integer(c_int),intent(out)            :: err
+
+  ! local pointers
+  type(hru_type), pointer               :: hru_data              ! local HRU data
+  ! local variables
+  character(len=256)                    :: cmessage
+  character(len=256)                    :: message 
+  logical(lgt)                          :: defNewOutputFile=.false.
+  logical(lgt)                          :: printRestart=.false.
+  logical(lgt)                          :: printProgress=.false.
+  character(len=256)                    :: restartFile       ! restart file name
+  character(len=256)                    :: timeString        ! portion of restart file name that contains the write-out time
+  character (len = 5) :: output_fileSuffix
+
+  integer(i4b)                          :: iStruct           ! index of model structure
+  integer(i4b)                          :: iFreq             ! index of the output frequency
+  integer(i4b)                          :: iVar 
+  integer(i4b)                          :: iDat
+  
+  ! convert the C pointers to Fortran pointers
+  call c_f_pointer(handle_hru_data, hru_data)
+  err=0; message='summa_manageOutputFiles/'
+  
+  ! ****************************************************************************
+  ! *** write restart data
+  ! ****************************************************************************
+  
+  ! write prog vars
+  do iVar = 1, size(hru_data%progstruct%var(:))
+    select case (prog_meta(iVar)%varType)
+    case(iLookVarType%scalarv);
+      summa_struct(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(checkPoint)%dat(:) = hru_data%progstruct%var(iVar)%dat(:) 
+    case(iLookVarType%wlength);              
+      summa_struct(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(checkPoint)%dat(:) = hru_data%progstruct%var(iVar)%dat(:) 
+    case(iLookVarType%midSoil);              
+      summa_struct(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(checkPoint)%dat(:) = hru_data%progstruct%var(iVar)%dat(:) 
+    case(iLookVarType%midToto);              
+      do iDat = 1, size(hru_data%progstruct%var(iVar)%dat(:))
+        summa_struct(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(checkPoint)%dat(iDat) = hru_data%progstruct%var(iVar)%dat(iDat)
+      end do ! iDat 
+    case(iLookVarType%ifcSoil);              
+      summa_struct(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(checkPoint)%dat(:) = hru_data%progstruct%var(iVar)%dat(:) 
+    case(iLookVarType%ifcToto);              
+      do iDat = 0, size(hru_data%progstruct%var(iVar)%dat(:)) -1 !vartype 8 in hru_data begins its index at 0 instead of the default of 1. this case accomodates that
+        summa_struct(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(checkPoint)%dat(iDat+1) = hru_data%progstruct%var(iVar)%dat(iDat)
+      end do ! iDat 
+    case(iLookVarType%midSnow);
+      summa_struct(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(checkPoint)%dat(:) = hru_data%progstruct%var(iVar)%dat(:) 
+    case(iLookVarType%ifcSnow); 
+      summa_struct(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(checkPoint)%dat(:) = hru_data%progstruct%var(iVar)%dat(:) 
+    case default; err=20; message=trim(message)//'unknown var type'; return
+   end select
+  end do ! iVar
+
+  ! write Basin var
+  summa_struct(1)%bvarStruct%gru(indxGRU)%hru(indxHRU)%var(iLookBVAR%routingRunoffFuture)%tim(checkPoint)%dat(:) = hru_data%bvarstruct%var(iLookBVAR%routingRunoffFuture)%dat(:)
+
+  
+end subroutine hru_writeRestart
+
+
 ! **********************************************************************************************************
 ! public subroutine writeParm: write model parameters
 ! **********************************************************************************************************
@@ -263,17 +378,17 @@ subroutine writeParm(indxGRU,indxHRU,ispatial,struct,meta,structName,err,message
       select type (struct)
         class is (var_i)
         if (structName == "type")then
-          outputStructure(1)%typeStruct%gru(indxGRU)%hru(indxHRU)%var(iVar) = struct%var(iVar)
+          summa_struct(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)
+          summa_struct(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)
+          summa_struct(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
@@ -285,7 +400,7 @@ subroutine writeParm(indxGRU,indxHRU,ispatial,struct,meta,structName,err,message
       select type (struct)
         class is (var_d)
         if (structName == "bpar")then
-          outputStructure(1)%bparStruct%gru(indxGRU)%var(iVar) = struct%var(iVar) ! this will overwrite data
+          summa_struct(1)%bparStruct%gru(indxGRU)%var(iVar) = struct%var(iVar) ! this will overwrite data
           print*, "bpar"
         end if
         class is (var_i8)
@@ -360,7 +475,7 @@ subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
         ! Write the time step values
         select type(dat)      ! forcStruc
           class is (var_d)    ! x%var(:)
-            outputStructure(1)%forcStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep) = dat%var(iVar)
+            summa_struct(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
@@ -377,15 +492,15 @@ subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
           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)
+              summa_struct(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)
+              summa_struct(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)
+              summa_struct(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)
+              summa_struct(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)
+              summa_struct(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
@@ -397,11 +512,11 @@ subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
 
         ! 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
+        summa_struct(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
+        summa_struct(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
+        summa_struct(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)
@@ -420,16 +535,16 @@ subroutine writeData(indxGRU,indxHRU,iStep,structName,finalizeStats, &
           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(:)
+                summa_struct(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(:)
+                summa_struct(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(:)
+                summa_struct(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(:)
+            summa_struct(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
 
@@ -509,10 +624,10 @@ subroutine writeBasin(indxGRU,indxHRU,iStep,finalizeStats,&
    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)
+      summa_struct(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)
+      summa_struct(1)%bvarStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(iFreq) = dat(iVar)%dat(iFreq)
      end if
 
     case default
@@ -563,8 +678,8 @@ subroutine writeTime(indxGRU,indxHRU,iStep,finalizeStats,meta,dat,err,message)
    ! 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)
+   ! add to summa_struct
+   summa_struct(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
 
@@ -851,15 +966,15 @@ end subroutine writeRestart
 ! ****************************************************************************** 
 subroutine setFinalizeStatsFalse(indx_gru) & 
   bind(C, name='setFinalizeStatsFalse')
-  USE output_structure_module,only:outputStructure
+  USE output_structure_module,only:summa_struct
   implicit none
   integer(c_int), intent(in)        :: indx_gru
 
   integer(i4b)                      :: iStep
 
   ! set finalizeStats to false
-  do iStep=1, size(outputStructure(1)%finalizeStats%gru(indx_gru)%hru(1)%tim)
-    outputStructure(1)%finalizeStats%gru(indx_gru)%hru(1)%tim(iStep)%dat = .false.
+  do iStep=1, size(summa_struct(1)%finalizeStats%gru(indx_gru)%hru(1)%tim)
+    summa_struct(1)%finalizeStats%gru(indx_gru)%hru(1)%tim(iStep)%dat = .false.
   end do
 end subroutine setFinalizeStatsFalse