diff --git a/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp b/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp
index 7bf60e01d296b4a6f0f09fc6bc2db43cc11f7ee8..c04458dc2bdf99c574429cb3c1783a42ecedf0d4 100644
--- a/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp
+++ b/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp
@@ -48,4 +48,10 @@ extern "C" {
     void* handle_diag_struct, void* handle_flux_stat, void* handle_flux_struct,
     void* handle_indx_stat, void* handle_indx_struct, void* handle_output_timestep,
     int* err);
+
+  void writeBasinToNetCDF(void* handle_ncid, int* index_gru, void* handle_finalize_stats,
+    void* handle_output_timestep, void* handle_bvar_stat, void* handle_bvar_struct, int* err);
+
+  void writeTimeToNetCDF(void* handle_ncid, void* handle_finalize_stats, void* handle_output_timestep,
+    void* handle_time_struct, int* err);
 }
diff --git a/build/source/actors/file_access_actor/file_access_actor.cpp b/build/source/actors/file_access_actor/file_access_actor.cpp
index 3e583b433213a68dae8e48fb8d2cafec351e04d4..6456b56413ad3435b12f5f8cd0566ea2825ed891 100644
--- a/build/source/actors/file_access_actor/file_access_actor.cpp
+++ b/build/source/actors/file_access_actor/file_access_actor.cpp
@@ -136,46 +136,82 @@ behavior file_access_actor(stateful_actor<file_access_state>* self, int startGRU
             }
         },
 
-        [=](write_output, int index_gru, int index_hru, std::vector<int> finalize_stats, 
-            std::vector<std::vector<double>> forc_stat, std::vector<double> forc_struct,
-            std::vector<std::vector<double>> prog_stat, std::vector<std::vector<double>> prog_struct,
-            std::vector<std::vector<double>> diag_stat, std::vector<std::vector<double>> diag_struct,
-            std::vector<std::vector<double>> flux_stat, std::vector<std::vector<double>> flux_struct,
-            std::vector<std::vector<double>> indx_stat, std::vector<std::vector<int>> indx_struct,
-            std::vector<int> output_timestep) {
+        [=](write_output, int index_gru, int index_hru, 
+            // statistic structures
+            std::vector<std::vector<double>> forc_stat, std::vector<std::vector<double>> prog_stat, std::vector<std::vector<double>> diag_stat,
+            std::vector<std::vector<double>> flux_stat, std::vector<std::vector<double>> indx_stat, std::vector<std::vector<double>> bvar_stat,
+            // primary data structures (scalars)
+            std::vector<int> time_struct, std::vector<double> forc_struct, std::vector<double> attr_struct,
+            std::vector<int> type_struct, std::vector<long int> id_struct,  
+            // primary data structures (variable length vectors)
+            std::vector<std::vector<int>> indx_struct, std::vector<std::vector<double>> mpar_struct, std::vector<std::vector<double>> prog_struct,
+            std::vector<std::vector<double>> diag_struct, std::vector<std::vector<double>> flux_struct,
+             // basin-average structures
+            std::vector<double> bpar_struct, std::vector<std::vector<double>> bvar_struct,
+            // ancillary data structures
+            std::vector<double> dpar_struct, std::vector<int> finalize_stats, std::vector<int> output_timestep ) {
             
             int err = 0;
-            void* handle_finalize_stats   = new_handle_var_i();
+            // statistic structures
             void* handle_forc_stat        = new_handle_var_dlength();
-            void* handle_forc_struct      = new_handle_var_d();
             void* handle_prog_stat        = new_handle_var_dlength();
-            void* handle_prog_struct       = new_handle_var_dlength();
             void* handle_diag_stat        = new_handle_var_dlength();
-            void* handle_diag_struct      = new_handle_var_dlength();
             void* handle_flux_stat        = new_handle_var_dlength();
-            void* handle_flux_struct      = new_handle_var_dlength();
             void* handle_indx_stat        = new_handle_var_dlength();
-            void* handle_indx_struct      = new_handle_var_ilength();
-            void* handle_output_time_step = new_handle_var_i();
-            
-            set_var_i(finalize_stats, handle_finalize_stats);
+            void* handle_bvar_stat        = new_handle_var_dlength();
             set_var_dlength(forc_stat, handle_forc_stat);
-            set_var_d(forc_struct, handle_forc_struct);
             set_var_dlength(prog_stat, handle_prog_stat);
-            set_var_dlength(prog_struct, handle_prog_struct);
             set_var_dlength(diag_stat, handle_diag_stat);
-            set_var_dlength(diag_struct, handle_diag_struct);
             set_var_dlength(flux_stat, handle_flux_stat);
-            set_var_dlength(flux_struct, handle_flux_struct);
             set_var_dlength(indx_stat, handle_indx_stat);
+            set_var_dlength(bvar_stat, handle_bvar_stat);
+            // primary data structures (scalars)
+            void* handle_time_struct      = new_handle_var_i();
+            void* handle_forc_struct      = new_handle_var_d();
+            void* handle_attr_struct      = new_handle_var_d();
+            void* handle_type_struct      = new_handle_var_i();
+            void* handle_id_struct        = new_handle_var_i8();
+            set_var_i(time_struct, handle_time_struct);
+            set_var_d(forc_struct, handle_forc_struct);
+            set_var_d(attr_struct, handle_attr_struct);
+            set_var_i(type_struct, handle_type_struct);
+            set_var_i8(id_struct, handle_id_struct);
+            // primary data structures (variable length vectors)
+            void* handle_indx_struct      = new_handle_var_ilength();
+            void* handle_mpar_struct      = new_handle_var_dlength();
+            void* handle_prog_struct      = new_handle_var_dlength();
+            void* handle_diag_struct      = new_handle_var_dlength();
+            void* handle_flux_struct      = new_handle_var_dlength();
             set_var_ilength(indx_struct, handle_indx_struct);
-            set_var_i(output_timestep, handle_output_time_step);
+            set_var_dlength(mpar_struct, handle_mpar_struct);
+            set_var_dlength(prog_struct, handle_prog_struct);
+            set_var_dlength(diag_struct, handle_diag_struct);
+            set_var_dlength(flux_struct, handle_flux_struct);
+            // basin-average structures
+            void* handle_bpar_struct      = new_handle_var_d();
+            void* handle_bvar_struct      = new_handle_var_dlength();
+            set_var_d(bpar_struct,handle_bpar_struct);
+            set_var_dlength(bvar_struct,handle_bvar_struct);
+            // ancillary data structures
+            void* handle_dpar_struct      = new_handle_var_d();
+            void* handle_finalize_stats   = new_handle_var_i();
+            void* handle_output_timestep = new_handle_var_i();
+            set_var_d(dpar_struct, handle_dpar_struct);
+            set_var_i(finalize_stats, handle_finalize_stats);
+            set_var_i(output_timestep, handle_output_timestep);
+
+            writeBasinToNetCDF(self->state.handle_ncid, &index_gru,
+                handle_finalize_stats, handle_output_timestep, handle_bvar_stat,
+                handle_bvar_struct, &err);
+        
+            writeTimeToNetCDF(self->state.handle_ncid,
+                handle_finalize_stats, handle_output_timestep, handle_time_struct, &err);
 
             writeDataToNetCDF(self->state.handle_ncid, &index_gru, &index_hru,
                 handle_finalize_stats, handle_forc_stat, handle_forc_struct,
                 handle_prog_stat, handle_prog_struct, handle_diag_stat, 
                 handle_diag_struct, handle_flux_stat, handle_flux_struct,
-                handle_indx_stat, handle_indx_struct, handle_output_time_step,
+                handle_indx_stat, handle_indx_struct, handle_output_timestep,
                 &err);
 
         },
diff --git a/build/source/actors/file_access_actor/file_access_actor.f90 b/build/source/actors/file_access_actor/file_access_actor.f90
index 0582bc0382687fb1051a212f1e4f78ae74074a06..dffef44269cf5eb930646097b164868f0fa2227d 100644
--- a/build/source/actors/file_access_actor/file_access_actor.f90
+++ b/build/source/actors/file_access_actor/file_access_actor.f90
@@ -195,7 +195,7 @@ subroutine writeDataToNetCDF(handle_ncid,          &
   end do  ! (looping through structures)
 end subroutine writeDataToNetCDF
 
-subroutine writeBasinToNetCDF(handle_ncid, index_gru, index_hru, handle_finalize_stats, &
+subroutine writeBasinToNetCDF(handle_ncid, index_gru, handle_finalize_stats, &
   handle_output_timestep, handle_bvar_stat, handle_bvar_struct, err) bind(C, name="writeBasinToNetCDF")
   USE writeOutput_module,only:writeBasin
   USE globalData,only:bvar_meta                 ! metadata on basin-average variables
@@ -205,7 +205,6 @@ subroutine writeBasinToNetCDF(handle_ncid, index_gru, index_hru, handle_finalize
   ! dummy variables
   type(c_ptr),    intent(in), value  :: handle_ncid
   integer(c_int), intent(in)         :: index_gru
-  integer(c_int), intent(in)         :: index_hru
   type(c_ptr),    intent(in), value  :: handle_finalize_stats
   type(c_ptr),    intent(in), value  :: handle_output_timestep
   type(c_ptr),    intent(in), value  :: handle_bvar_stat
diff --git a/build/source/actors/hru_actor/hru_actor.cpp b/build/source/actors/hru_actor/hru_actor.cpp
index 8ddea1b8966cce8bf995246997e614e648b14288..0b56832d2fb261e164ed9c0e004c9991ab5f64cd 100644
--- a/build/source/actors/hru_actor/hru_actor.cpp
+++ b/build/source/actors/hru_actor/hru_actor.cpp
@@ -90,6 +90,7 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
                 err = Run_HRU(self); // Simulate a Timestep
 
                 // Get Data from fortran
+                // statistic structures
                 std::vector<std::vector<double>> forc_stat_array    = get_var_dlength(self->state.handle_forcStat);
                 std::vector<std::vector<double>> prog_stat_array    = get_var_dlength(self->state.handle_progStat);
                 std::vector<std::vector<double>> diag_stat_array    = get_var_dlength(self->state.handle_diagStat);
@@ -97,39 +98,53 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
                 std::vector<std::vector<double>> indx_stat_array    = get_var_dlength(self->state.handle_indxStat);
                 std::vector<std::vector<double>> bvar_stat_array    = get_var_dlength(self->state.handle_bvarStat);
                 // primary data structures (scalars)
-                std::vector<int> time_struct_array                  = get_var_i(self->state.handle_timeStruct);
-                std::vector<double> forc_struct_array               = get_var_d(self->state.handle_forcStruct);
-                std::vector<double> attr_struct_array               = get_var_d(self->state.handle_attrStruct); 
-                std::vector<int> type_struct_array                  = get_var_i(self->state.handle_typeStruct);
+                std::vector<int>      time_struct_array             = get_var_i(self->state.handle_timeStruct);
+                std::vector<double>   forc_struct_array             = get_var_d(self->state.handle_forcStruct);
+                std::vector<double>   attr_struct_array             = get_var_d(self->state.handle_attrStruct); 
+                std::vector<int>      type_struct_array             = get_var_i(self->state.handle_typeStruct);
                 std::vector<long int> id_struct_array               = get_var_i8(self->state.handle_idStruct);
                 // primary data structures (variable length vectors)
-                std::vector<std::vector<int>> indx_struct_array     = get_var_ilength(self->state.handle_indxStruct);
+                std::vector<std::vector<int>>    indx_struct_array  = get_var_ilength(self->state.handle_indxStruct);
                 std::vector<std::vector<double>> mpar_struct_array  = get_var_dlength(self->state.handle_mparStruct);
-                std::vector<std::vector<double>> prog_struc_array   = get_var_dlength(self->state.handle_progStruct);
+                std::vector<std::vector<double>> prog_struct_array   = get_var_dlength(self->state.handle_progStruct);
                 std::vector<std::vector<double>> diag_struct_array  = get_var_dlength(self->state.handle_diagStruct);
                 std::vector<std::vector<double>> flux_struct_array  = get_var_dlength(self->state.handle_fluxStruct);
                 // basin-average structures
-                std::vector<double> bpar_struct_array               = get_var_d(self->state.handle_bparStruct);
+                std::vector<double>              bpar_struct_array  = get_var_d(self->state.handle_bparStruct);
                 std::vector<std::vector<double>> bvar_struct_array  = get_var_dlength(self->state.handle_bvarStruct);
                 // ancillary data structures
-                std::vector<double> dpar_struct_array               = get_var_d(self->state.handle_dparStruct);
-                std::vector<int> finalize_stats_array               = get_flagVec(self->state.handle_finalizeStats);
-                std::vector<int> output_time_step_array             = get_var_i(self->state.handle_outputTimeStep);
+                std::vector<double>   dpar_struct_array             = get_var_d(self->state.handle_dparStruct);
+                std::vector<int>      finalize_stats_array          = get_flagVec(self->state.handle_finalizeStats);
+                std::vector<int>      output_time_step_array        = get_var_i(self->state.handle_outputTimeStep);
 
                 self->send(self->state.file_access_actor, write_output_v,
                     self->state.indxGRU,
                     self->state.indxHRU,
-                    finalize_stats_array,
+                    // statistic structures
                     forc_stat_array,
-                    forc_struct_array,
                     prog_stat_array,
-                    prog_struc_array,
                     diag_stat_array,
-                    diag_struct_array,
                     flux_stat_array,
-                    flux_struct_array,
                     indx_stat_array,
+                    bvar_stat_array,
+                    // primary data structures (scalars)
+                    time_struct_array,
+                    forc_struct_array,
+                    attr_struct_array,
+                    type_struct_array,
+                    id_struct_array,
+                    // primary data structures (variable length vectors)
                     indx_struct_array,
+                    mpar_struct_array,
+                    prog_struct_array,
+                    diag_struct_array,
+                    flux_struct_array,
+                    // basin-average structures
+                    bpar_struct_array,
+                    bvar_struct_array,
+                    // ancillary data structures
+                    dpar_struct_array,
+                    finalize_stats_array,
                     output_time_step_array);
 
                 updateCounters(self->state.handle_timeStruct, self->state.handle_statCounter, self->state.handle_outputTimeStep,
diff --git a/build/source/actors/hru_actor/hru_actor.f90 b/build/source/actors/hru_actor/hru_actor.f90
index 60ad9a336a3b418d693791044d8580de56c12208..c4c1d9b27d77fabac4e67c8cee40dfeddb0f72c2 100644
--- a/build/source/actors/hru_actor/hru_actor.f90
+++ b/build/source/actors/hru_actor/hru_actor.f90
@@ -197,6 +197,10 @@ subroutine prepareOutput(&
                             err, cmessage)                                  ! error control
   if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
 
+  print*, "iLookPROG%mLayerHeight = ", iLookPROG%mLayerHeight
+  print*, "mLayerHeight = ", progStruct%var(iLookPROG%mLayerHeight)%dat(100)
+  print*, "mLayerHeight = ", progStruct%var(iLookPROG%mLayerHeight)%dat(103)
+
  ! ****************************************************************************
  ! *** calculate output statistics
  ! ****************************************************************************
diff --git a/build/source/netcdf/writeOutput.f90 b/build/source/netcdf/writeOutput.f90
index 354522efdeaa40227b398bffe0878f0c02bce62c..c94e9d20eae8138a8c1b4e8ad3c4e47ba9e79625 100644
--- a/build/source/netcdf/writeOutput.f90
+++ b/build/source/netcdf/writeOutput.f90
@@ -186,6 +186,12 @@ subroutine writeDataNew(ncid, finalize_stats, output_timestep, max_layers, index
   ! output array
   integer(i4b)                       :: datLength         ! length of each data vector
   integer(i4b)                       :: maxLength         ! maximum length of each data vector
+  real(rkind)                        :: realVec(num_gru)  ! real vector for all HRUs in the run domain
+  real(rkind)                        :: realArray(num_gru,max_layers+1)  ! real array for all HRUs in the run domain
+  integer(i4b)                       :: intArray(num_gru,max_layers+1)  ! real array for all HRUs in the run domain
+  integer(i4b)                       :: dataType          ! type of data
+  integer(i4b),parameter             :: ixInteger=1001    ! named variable for integer
+  integer(i4b),parameter             :: ixReal=1002       ! named variable for real
   
   err=0;message="writeOutput.f90 - writeDataNew/"
   ! loop through output frequencies
@@ -239,8 +245,8 @@ subroutine writeDataNew(ncid, finalize_stats, output_timestep, max_layers, index
       if(meta(iVar)%varType==iLookVarType%scalarv) then
         select type(stat)
           class is (var_dlength)
-            ! realVec(gru_struc(index_gru)%hruInfo(index_hru)%hru_ix) = stat%var(map(iVar))%dat(iFreq)
-            err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),(stat%var(map(iVar))%dat(iFreq)),start=(/index_gru,output_timestep(iFreq)/))
+            realVec(index_gru) = stat%var(map(iVar))%dat(iFreq)
+            err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realVec,start=(/index_gru,output_timestep(iFreq)/),count=(/num_gru,1/))
             if (err/=0) then
               print*, message
               return
@@ -253,6 +259,19 @@ subroutine writeDataNew(ncid, finalize_stats, output_timestep, max_layers, index
         end select ! stat
 
       else
+
+         ! Write the data
+        select type(dat)
+          class is (var_dlength)
+            realArray(:,:) = realMissing;    dataType=ixReal
+          class is (var_ilength)
+            intArray(:,:) = integerMissing; dataType=ixInteger
+          class default
+          err=20
+          message=trim(message)//'data must be of type integer or real'
+          print*, message
+          return
+      end select
         ! vector
         nSoil   = indx%var(iLookIndex%nSoil)%dat(1)
         nSnow   = indx%var(iLookIndex%nSnow)%dat(1)
@@ -270,6 +289,13 @@ subroutine writeDataNew(ncid, finalize_stats, output_timestep, max_layers, index
           case default; cycle
         end select ! vartype
 
+                ! get the data vectors
+        select type (dat)
+          class is (var_dlength); realArray(index_gru,1:datLength) = dat%var(iVar)%dat(:)
+          class is (var_ilength);  intArray(index_gru,1:datLength) = dat%var(iVar)%dat(:)
+          class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return
+        end select
+
         ! get the maximum length of each data vector
         select case (meta(iVar)%varType)
           case(iLookVarType%wLength); maxLength = maxSpectral
@@ -283,12 +309,12 @@ subroutine writeDataNew(ncid, finalize_stats, output_timestep, max_layers, index
         end select ! vartype
 
         ! Write the data
-        select type(dat)
-          class is (var_dlength)
-            err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),(dat%var(iVar)%dat(:)),start=(/index_gru,1,output_timestep(iFreq)/),count=(/num_gru,maxLength,1/))
-          class is (var_ilength)
-            err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),(dat%var(iVar)%dat(:)),start=(/index_gru,1,output_timestep(iFreq)/),count=(/num_gru,maxLength,1/))
-          class default
+        select case(dataType)
+          case(ixReal)
+            err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realArray(1:num_gru,1:maxLength),start=(/index_gru,1,output_timestep(iFreq)/),count=(/num_gru,maxLength,1/))
+          case(ixInteger)
+            err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realArray(1:num_gru,1:maxLength),start=(/index_gru,1,output_timestep(iFreq)/),count=(/num_gru,maxLength,1/))
+          case default
             err=20
             message=trim(message)//'data must be of type integer or real'
             print*, message