diff --git a/.gitignore b/.gitignore
index 6a659174fd57b1389fa05c05a13e77c2abd7431e..9c3347b6ba3e7492740051e15a04a7514db21a1b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -27,3 +27,4 @@ build/source/testing/containers/output_container/out.txt
 build/summa_old
 bin/config.json
 bin/summa_sundials
+build/cmake_build/
diff --git a/build/CMakeLists.txt b/build/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..195911908d169ed8a4f266b440646defb3495b62
--- /dev/null
+++ b/build/CMakeLists.txt
@@ -0,0 +1,329 @@
+cmake_minimum_required(VERSION 3.10 FATAL_ERROR)
+project(summa_actors LANGUAGES CXX Fortran)
+enable_language(C)
+SET (CMAKE_Fortran_COMPILER  gfortran)
+include(FortranCInterface)
+FortranCInterface_VERIFY(CXX)
+
+get_filename_component(PARENT_DIR "${CMAKE_CURRENT_SOURCE_DIR}" DIRECTORY)
+set(EXEC_DIR ${PARENT_DIR}/bin) # set the output directory for executables
+SET(F_MASTER ${PARENT_DIR}/build/summa)
+# Add options for build type
+set(CMAKE_CONFIGURATION_TYPES BE BE_Cluster BE_Cluster_Debug)
+
+# Set Compiler Options
+if(CMAKE_BUILD_TYPE MATCHES Debug)
+    message("\nSetting Debug Options\n")
+    add_compile_definitions(DEBUG)
+    set(FLAGS_NOAH -g -O0 -fbacktrace -fbounds-check -ffree-form -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors ${FLAGS_OPT})
+    set(FLAGS_ALL  -g -O0 -fbacktrace -fbounds-check -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors -cpp ${FLAGS_OPT})
+    set(FLAGS_CXX  -g -O0 -fbounds-check -Wfatal-errors -std=c++17 ${FLAGS_OPT})
+else()
+    message("\nSetting Release Options")
+    set(FLAGS_NOAH -O3 -ffree-form -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors ${FLAGS_OPT})
+    set(FLAGS_ALL  -O3 -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors -cpp ${FLAGS_OPT})
+    set(FLAGS_CXX  -O3 -Wfatal-errors -std=c++17 ${FLAGS_OPT})
+endif()
+
+find_package(OpenBLAS REQUIRED)
+set(EXEC_NAME summa_be)
+
+if (CMAKE_BUILD_TYPE MATCHES Cluster)    
+    # Set include directories
+    set(INCLUDES $ENV{EBROOTNETCDFMINFORTRAN}/include ${netCDF_INCLUDES} ${OpenBLAS_INCLUDES})
+    set(LIBRARIES SUMMA_NOAHMP  ${OpenBLAS_LIBRARIES} -lnetcdff)
+
+    set(INC_ACTORS 
+        ${CAF_INCLUDES} 
+        ${PARENT_DIR}/build/includes/global 
+        ${PARENT_DIR}/build/includes/summa_actor 
+        ${PARENT_DIR}/build/includes/gru_actor 
+        ${PARENT_DIR}/build/includes/job_actor 
+        ${PARENT_DIR}/build/includes/file_access_actor 
+        ${PARENT_DIR}/build/includes/hru_actor)
+    set(LIB_ACTORS 
+        ${CAF_LIBRARIES} 
+        -lcaf_core 
+        -lcaf_io)
+
+else()
+    set(CMAKE_BUILD_RPATH "/usr/local/lib")
+    set(SUMMA_INCLUDES
+        "/usr/include"
+        ${netCDF_INCLUDES}
+        ${OpenBLAS_INCLUDES})
+    
+    set(SUMMA_LIBS
+        -lnetcdff
+        -lopenblas
+        SUMMA_NOAHMP)
+
+    set(SUMMA_ACTORS_INCLUDES
+        "/usr/local/actor-framework/debug/include"
+        ${CAF_INCLUDES}
+        ${OpenBLAS_INCLUDES}
+        "${PARENT_DIR}/build/includes/global"
+        "${PARENT_DIR}/build/includes/summa_actor"
+        "${PARENT_DIR}/build/includes/gru_actor"
+        "${PARENT_DIR}/build/includes/job_actor"
+        "${PARENT_DIR}/build/includes/file_access_actor"
+        "${PARENT_DIR}/build/includes/hru_actor")
+    
+    link_directories("/usr/local/actor-framework/debug/lib")
+    set(SUMMA_ACTORS_LIBS   
+        -lopenblas
+        -lcaf_core
+        -lcaf_io
+        summa
+        -lnetcdff)
+endif()
+
+# Define directories that contains source code
+set(DRIVER_DIR ${F_MASTER}/build/source/driver)
+set(DSHARE_DIR ${F_MASTER}/build/source/dshare)
+set(ENGINE_DIR ${F_MASTER}/build/source/engine)
+set(HOOKUP_DIR ${F_MASTER}/build/source/hookup)
+set(NETCDF_DIR ${F_MASTER}/build/source/netcdf)
+set(NOAHMP_DIR ${F_MASTER}/build/source/noah-mp)
+
+# Define directories for source files that might be replaced by actors (identical if not using actors)
+set(SUB_DRIVER_DIR ${PARENT_DIR}/build/source/driver)
+set(SUB_DSHARE_DIR ${PARENT_DIR}/build/source/dshare)
+set(SUB_ENGINE_DIR ${PARENT_DIR}/build/source/engine)
+set(SUB_HOOKUP_DIR ${PARENT_DIR}/build/source/hookup)
+set(SUB_NETCDF_DIR ${PARENT_DIR}/build/source/netcdf)
+
+# Define Actors specific directories
+set(ACTORS_DIR ${PARENT_DIR}/build/source/actors)
+set(FILE_ACCESS_DIR ${ACTORS_DIR}/file_access_actor)
+set(JOB_ACTOR_DIR   ${ACTORS_DIR}/job_actor)
+set(HRU_ACTOR_DIR   ${ACTORS_DIR}/hru_actor)
+set(GRU_ACTOR_DIR   ${ACTORS_DIR}/gru_actor)
+
+# NOAHMP modules
+set(NOAHMP
+    ${NOAHMP_DIR}/module_model_constants.F
+    ${NOAHMP_DIR}/module_sf_noahutl.F
+    ${NOAHMP_DIR}/module_sf_noahlsm.F
+    ${NOAHMP_DIR}/module_sf_noahmplsm.F)
+
+# Free versions of numerical recipes utilities for NOAH-MP modules
+set(NRUTIL
+    ${ENGINE_DIR}/f2008funcs.f90
+    ${ENGINE_DIR}/nr_utility.f90
+    ${ENGINE_DIR}/nrtype.f90)
+
+# Free versions of numerical recipes procedures for SUMMA modules
+set(NRPROC
+    ${ENGINE_DIR}/expIntegral.f90
+    ${ENGINE_DIR}/spline_int.f90)
+
+# Hook-up modules
+set(HOOKUP
+    ${HOOKUP_DIR}/ascii_util.f90
+    ${HOOKUP_DIR}/summaFileManager.f90)
+
+# Data modules
+set(DATAMS
+    ${SUB_DSHARE_DIR}/data_types.f90
+    ${DSHARE_DIR}/flxMapping.f90
+    ${DSHARE_DIR}/get_ixname.f90
+    ${DSHARE_DIR}/globalData.f90
+    ${DSHARE_DIR}/multiconst.f90
+    ${DSHARE_DIR}/outpt_stat.f90
+    ${DSHARE_DIR}/popMetadat.f90
+    ${DSHARE_DIR}/var_lookup.f90)
+
+# Utility modules
+set(UTILMS
+    ${ENGINE_DIR}/matrixOper.f90
+    ${ENGINE_DIR}/mDecisions.f90
+    ${ENGINE_DIR}/snow_utils.f90
+    ${ENGINE_DIR}/soil_utils.f90
+    ${ENGINE_DIR}/time_utils.f90
+    ${ENGINE_DIR}/updatState.f90)
+
+# NetCDF routines
+set(NETCDF
+    ${SUB_NETCDF_DIR}/def_output.f90
+    ${SUB_NETCDF_DIR}/modelwrite.f90
+    ${NETCDF_DIR}/netcdf_util.f90
+    ${NETCDF_DIR}/read_icond.f90)
+
+# Preliminary modules
+set(PRELIM
+    ${ENGINE_DIR}/allocspace.f90
+    ${ENGINE_DIR}/check_icond.f90
+    ${ENGINE_DIR}/checkStruc.f90
+    ${ENGINE_DIR}/childStruc.f90
+    ${ENGINE_DIR}/conv_funcs.f90
+    ${ENGINE_DIR}/convE2Temp.f90
+    ${SUB_ENGINE_DIR}/ffile_info.f90
+    ${ENGINE_DIR}/read_pinit.f90
+    ${ENGINE_DIR}/read_attrb.f90
+    ${ENGINE_DIR}/paramCheck.f90
+    ${ENGINE_DIR}/pOverwrite.f90
+    ${ENGINE_DIR}/sunGeomtry.f90
+    ${ENGINE_DIR}/read_param.f90)
+
+    # Model run support modules
+set(MODRUN
+    ${ENGINE_DIR}/canopySnow.f90
+    ${ENGINE_DIR}/derivforce.f90
+    ${ENGINE_DIR}/getVectorz.f90
+    ${ENGINE_DIR}/indexState.f90
+    ${ENGINE_DIR}/layerMerge.f90
+    ${ENGINE_DIR}/layerDivide.f90
+    ${ENGINE_DIR}/qTimeDelay.f90
+    ${ENGINE_DIR}/snowAlbedo.f90
+    ${ENGINE_DIR}/snwCompact.f90
+    ${ENGINE_DIR}/tempAdjust.f90
+    ${ENGINE_DIR}/updateVars.f90
+    ${ENGINE_DIR}/var_derive.f90
+    ${ENGINE_DIR}/volicePack.f90)
+
+# Solver main modules
+set(SOLVER
+    ${ENGINE_DIR}/bigAquifer.f90
+    ${ENGINE_DIR}/computFlux.f90
+    ${ENGINE_DIR}/computJacob.f90
+    ${ENGINE_DIR}/computResid.f90
+    ${ENGINE_DIR}/coupled_em.f90
+    ${ENGINE_DIR}/diagn_evar.f90
+    ${ENGINE_DIR}/eval8summa.f90
+    ${ENGINE_DIR}/groundwatr.f90
+    ${ENGINE_DIR}/opSplittin.f90
+    ${ENGINE_DIR}/snowLiqFlx.f90
+    ${ENGINE_DIR}/soilLiqFlx.f90
+    ${ENGINE_DIR}/ssdNrgFlux.f90
+    ${ENGINE_DIR}/stomResist.f90
+    ${ENGINE_DIR}/summaSolve.f90
+    ${ENGINE_DIR}/systemSolv.f90
+    ${ENGINE_DIR}/varSubstep.f90
+    ${ENGINE_DIR}/vegLiqFlux.f90
+    ${ENGINE_DIR}/vegNrgFlux.f90
+    ${ENGINE_DIR}/vegPhenlgy.f90
+    ${ENGINE_DIR}/vegSWavRad.f90)
+
+set(DRIVER
+    ${DRIVER_DIR}/summa_alarms.f90
+    ${DRIVER_DIR}/summa_globalData.f90)
+
+
+# Actors interface modules
+set(INTERFACE
+    ${ACTORS_DIR}/global/cppwrap_auxiliary.f90
+    ${ACTORS_DIR}/global/cppwrap_datatypes.f90
+    ${ACTORS_DIR}/global/cppwrap_metadata.f90)
+set(FILE_ACCESS_INTERFACE
+    ${FILE_ACCESS_DIR}/fortran_code/cppwrap_fileAccess.f90
+    ${FILE_ACCESS_DIR}/fortran_code/output_structure.f90
+    ${FILE_ACCESS_DIR}/fortran_code/read_force.f90
+    ${FILE_ACCESS_DIR}/fortran_code/write_to_netcdf.f90
+    ${FILE_ACCESS_DIR}/fortran_code/writeOutputFromOutputStructure.f90)
+set(JOB_INTERFACE
+    ${JOB_ACTOR_DIR}/job_actor.f90)
+set(HRU_INTERFACE
+    ${HRU_ACTOR_DIR}/fortran_code/hru_actor.f90
+    ${HRU_ACTOR_DIR}/fortran_code/hru_init.f90
+    ${HRU_ACTOR_DIR}/fortran_code/hru_modelRun.f90
+    ${HRU_ACTOR_DIR}/fortran_code/hru_modelwrite.f90
+    ${HRU_ACTOR_DIR}/fortran_code/hru_restart.f90
+    ${HRU_ACTOR_DIR}/fortran_code/hru_setup.f90
+    ${HRU_ACTOR_DIR}/fortran_code/hru_writeOutput.f90)
+
+# Actors actual actor modules
+set(ACTORS_GLOBAL
+    ${ACTORS_DIR}/global/auxiliary.cpp
+    ${ACTORS_DIR}/global/global.cpp
+    ${ACTORS_DIR}/global/message_atoms.cpp
+    ${ACTORS_DIR}/global/settings_functions.cpp
+    ${ACTORS_DIR}/global/timing_info.cpp)
+set(SUMMA_ACTOR
+    ${ACTORS_DIR}/summa_actor/batch/batch.cpp
+    ${ACTORS_DIR}/summa_actor/batch/batch_container.cpp
+    ${ACTORS_DIR}/summa_actor/client/client.cpp
+    ${ACTORS_DIR}/summa_actor/client/client_container.cpp
+    ${ACTORS_DIR}/summa_actor/summa_actor.cpp
+    ${ACTORS_DIR}/summa_actor/summa_backup_server.cpp
+    ${ACTORS_DIR}/summa_actor/summa_client.cpp
+    ${ACTORS_DIR}/summa_actor/summa_server.cpp)
+set(FILE_ACCESS_ACTOR
+    ${ACTORS_DIR}/file_access_actor/cpp_code/file_access_actor.cpp
+    ${ACTORS_DIR}/file_access_actor/cpp_code/forcing_file_info.cpp
+    ${ACTORS_DIR}/file_access_actor/cpp_code/output_container.cpp)
+set(JOB_ACTOR
+    ${ACTORS_DIR}/job_actor/GRU.cpp
+    ${ACTORS_DIR}/job_actor/job_actor.cpp)
+set(HRU_ACTOR
+    ${ACTORS_DIR}/hru_actor/cpp_code/hru_actor.cpp)
+
+#=========================================================================================
+# COMPILE PART 3: Collect the subroutines into build groups depending on build type
+#=========================================================================================
+set(COMM_ALL
+    ${NRPROC}
+    ${HOOKUP}
+    ${DATAMS}
+    ${UTILMS}
+    ${INTERFACE})
+
+set(SUMMA_ALL
+    ${NETCDF}
+    ${PRELIM}
+    ${MODRUN}
+    ${SOLVER}
+    ${DRIVER})
+
+set(SUMMA_ALL 
+    ${SUMMA_ALL}
+    ${FILE_ACCESS_INTERFACE}
+    ${JOB_INTERFACE}
+    ${HRU_INTERFACE})
+
+set(MAIN_ACTOR ${ACTORS_DIR}/main.cpp)
+
+# Define version number, not working correctly
+set(VERSIONFILE     ${DRIVER_DIR}/summaversion.inc)
+execute_process(COMMAND "    ${GIT_EXECUTABLE} tag | tail -n 1" OUTPUT_VARIABLE VERSION)
+execute_process(COMMAND "date" OUTPUT_VARIABLE BULTTIM)
+execute_process(COMMAND "    ${GIT_EXECUTABLE} describe --long --all --always | sed -e's/heads\///'" OUTPUT_VARIABLE GITBRCH)
+execute_process(COMMAND "    ${GIT_EXECUTABLE} rev-parse HEAD" OUTPUT_VARIABLE GITHASH)
+
+
+#=========================================================================================
+# COMPILE PART 4: Do the compilation
+#=========================================================================================
+# update version information, not working correctly
+file(WRITE  ${VERSIONFILE} "character(len=64), parameter     :: summaVersion = '${VERSION}'\n")
+file(APPEND ${VERSIONFILE} "character(len=64), parameter     :: buildTime = ''\n")
+file(APPEND ${VERSIONFILE} "character(len=64), parameter     :: gitBranch = '${GITBRCH}'\n")
+file(APPEND ${VERSIONFILE} "character(len=64), parameter     :: gitHash = '${GITHASH}'")
+
+# Build SUMMA_NOAHMP Object
+add_library(SUMMA_NOAHMP OBJECT ${NOAHMP} ${NRUTIL})
+target_compile_options(SUMMA_NOAHMP PRIVATE ${FLAGS_NOAH})
+
+# Build SUMMA_COMM Object
+add_library(SUMMA_COMM OBJECT ${COMM_ALL})
+target_compile_options(SUMMA_COMM PRIVATE ${FLAGS_ALL})
+target_include_directories(SUMMA_COMM PRIVATE ${INCLUDES})
+target_link_libraries(SUMMA_COMM PUBLIC SUMMA_NOAHMP ${FLAGS_ALL}) # added flags to the link step
+
+add_library(summaactors SHARED ${SUMMA_ALL})
+target_compile_options(summaactors PRIVATE ${FLAGS_ALL})
+target_include_directories(summaactors PUBLIC ${INCLUDES})
+target_link_libraries(summaactors PUBLIC ${LIBRARIES} SUMMA_NOAHMP SUMMA_COMM)
+add_executable(${EXEC_NAME} 
+               ${MAIN_ACTOR}
+               ${ACTORS_GLOBAL}
+               ${FILE_ACCESS_ACTOR}
+               ${JOB_ACTOR}
+               ${HRU_ACTOR}
+               ${SUMMA_ACTOR}
+               ${SUMMA_CLIENT}
+               ${SUMMA_SERVER})
+set_property(TARGET ${EXEC_NAME} PROPERTY LINKER_LANGUAGE Fortran)
+target_compile_options(${EXEC_NAME} PUBLIC ${FLAGS_CXX})
+target_include_directories(${EXEC_NAME} PUBLIC ${INC_ACTORS})
+target_link_libraries( ${EXEC_NAME} ${LIB_ACTORS} summaactors)
\ No newline at end of file
diff --git a/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90 b/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90
index 1c67351a140525d8cbd5f3a333fbec3d8bfcfb66..16351f8c7afd4ec74bad4f7acf8f9cd0bb2c2b14 100644
--- a/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90
+++ b/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90
@@ -74,14 +74,16 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
   USE globalData,only:checkHRU                                ! index of the HRU for a single HRU run
   
   ! look-up values for the choice of heat capacity computation
-  USE mDecisions_module,only:  &
-  enthalpyFD                             ! heat capacity using enthalpy
+#ifdef SUNDIALS_ACTIVE
+  USE mDecisions_module,only:enthalpyFD                       ! heat capacity using enthalpy
+  USE t2enthalpy_module,only:T2E_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
+#endif
   USE mDecisions_module,only:&
- monthlyTable,& ! LAI/SAI taken directly from a monthly table for different vegetation classes
- specified      ! LAI/SAI computed from green vegetation fraction and winterSAI and summerLAI parameters
+                        monthlyTable,& ! LAI/SAI taken directly from a monthly table for different vegetation classes
+                        specified      ! LAI/SAI computed from green vegetation fraction and winterSAI and summerLAI parameters
 
   USE ConvE2Temp_module,only:E2T_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
-  USE t2enthalpy_module,only:T2E_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
+
 
   USE NOAHMP_VEG_PARAMETERS,only:SAIM,LAIM                    ! 2-d tables for stem area index and leaf area index (vegType,month)
   USE NOAHMP_VEG_PARAMETERS,only:HVT,HVB                      ! height at the top and bottom of vegetation (vegType)
@@ -294,6 +296,7 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
       if(err/=0)then; print*,message; return; endif
 
       ! calculate a lookup table to compute enthalpy from temperature, only for enthalpyFD
+#ifdef SUNDIALS_ACTIVE      
       if(model_decisions(iLookDECISIONS%howHeatCap)%iDecision == enthalpyFD)then
         call T2E_lookup(gru_struc(iGRU)%hruInfo(iHRU)%nSoil,   &   ! intent(in):    number of soil layers
                         outputStructure(1)%mparStruct%gru(iGRU)%hru(iHRU),        &   ! intent(in):    parameter data structure
@@ -301,7 +304,7 @@ subroutine fileAccessActor_init_fortran(& ! Variables for forcing
                         err,message)                              ! intent(out):   error control
         if(err/=0)then; print*, message; return; endif
       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)
@@ -400,7 +403,7 @@ subroutine FileAccessActor_DeallocateStructures(handle_forcFileInfo, handle_ncid
   USE globalData,only:failedHRUs
   USE globalData,only:forcingDataStruct
   USE globalData,only:vectime
-  USE globalData,only:outputTimeStep
+  USE output_structure_module,only:outputTimeStep
   implicit none
   type(c_ptr),intent(in), value        :: handle_forcFileInfo
   type(c_ptr),intent(in), value        :: handle_ncid
diff --git a/build/source/actors/file_access_actor/fortran_code/output_structure.f90 b/build/source/actors/file_access_actor/fortran_code/output_structure.f90
index 9fd2c5ce422785299a2ecc242926d7c970c9756e..719993917817093b45a23656d5b87e6cdcb4ea89 100644
--- a/build/source/actors/file_access_actor/fortran_code/output_structure.f90
+++ b/build/source/actors/file_access_actor/fortran_code/output_structure.f90
@@ -96,19 +96,18 @@ module output_structure_module
     ! finalize stats structure
     type(gru_hru_time_flagVec)                        :: finalizeStats                 ! x%gru(:)%hru(:)%tim(:)%dat -- flags on when to write to file
 
-
     type(gru_d)                                       :: upArea
 
     integer(i4b)                                      :: nTimeSteps
   end type summa_output_type  
   
-  
+
   type(summa_output_type),allocatable,save,public :: outputStructure(:) ! summa_OutputStructure(1)%struc%var(:)%dat(nTimeSteps) 
+  type(ilength),allocatable,save,public           :: outputTimeStep(:)                 ! timestep in output files
   
   contains
 
 subroutine initOutputTimeStep(num_gru, err)
-  USE globalData,only:outputTimeStep
   USE var_lookup,only:maxvarFreq                ! maximum number of output files
   implicit none
   integer(i4b), intent(in)  :: num_gru
@@ -132,7 +131,9 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
   USE globalData,only:prog_meta,diag_meta,flux_meta,id_meta   ! metadata structures
   USE globalData,only:mpar_meta,indx_meta                     ! metadata structures
   USE globalData,only:bpar_meta,bvar_meta                     ! metadata structures
+#ifdef SUNDIALS_ACTIVE
   USE globalData,only:lookup_meta
+#endif
   USE globalData,only:statForc_meta                           ! child metadata for stats
   USE globalData,only:statProg_meta                           ! child metadata for stats
   USE globalData,only:statDiag_meta                           ! child metadata for stats
@@ -146,6 +147,7 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
   USE var_lookup,only:maxvarFreq                              ! maximum number of output files
 
   USE allocspace_module,only:allocGlobal                      ! module to allocate space for global data structures
+  USE globalData,only:maxSnowLayers
   
   implicit none
   type(file_info_array),intent(in)      :: forcFileInfo
@@ -248,7 +250,9 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
     case('bpar'); call allocGlobal(bpar_meta,outputStructure(1)%bparStruct    ,err, message);  ! basin-average params 
     case('bvar'); call allocGlobal(bvar_meta,outputStructure(1)%bvarStruct_init,err,message);  ! basin-average variables
     case('deriv'); cycle;
+#ifdef SUNDIALS_ACTIVE      
     case('lookup'); call allocGlobal(lookup_meta,outputStructure(1)%lookupStruct,err, message);
+#endif
     end select
   end do
   
@@ -274,10 +278,10 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
             case('forc')
               ! Structure
               call alloc_outputStruc(forc_meta,outputStructure(1)%forcStruct%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model forcing data
+                                      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), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model forcing data
+                                      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),&
               !                         nSteps=1,err=err,message=message)
@@ -287,31 +291,31 @@ subroutine initOutputStructure(forcFileInfo, maxSteps, num_gru, err)
             case('indx')
               ! Structure
               call alloc_outputStruc(indx_meta,outputStructure(1)%indxStruct%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,str_name='indx',message=message);    ! model variables
+                                      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), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! index vars
+                                      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=nSnow,nSoil=nSoil,err=err,message=message);    ! model prognostic (state) variables
+                                      nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model prognostic (state) variables
               ! Statistics
               call alloc_outputStruc(statProg_meta(:)%var_info,outputStructure(1)%progStat%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model prognostic 
+                                      nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model prognostic 
             case('diag')
               ! Structure
               call alloc_outputStruc(diag_meta,outputStructure(1)%diagStruct%gru(iGRU)%hru(iHRU), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model diagnostic variables
+                                      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), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model diagnostic
+                                      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), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model fluxes
+                                      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), &
-                                      nSteps=maxSteps,nSnow=nSnow,nSoil=nSoil,err=err,message=message);    ! model fluxes
+                                      nSteps=maxSteps,nSnow=maxSnowLayers,nSoil=nSoil,err=err,message=message);    ! model fluxes
             case('bpar'); cycle;
             case('bvar')
               ! Structure
diff --git a/build/source/actors/file_access_actor/fortran_code/writeOutputFromOutputStructure.f90 b/build/source/actors/file_access_actor/fortran_code/writeOutputFromOutputStructure.f90
index 364d14537077dde693c61e89a9237a351175c6ee..f5a3ed7ffba02fbaa06b387b2033ad8862ea04e8 100644
--- a/build/source/actors/file_access_actor/fortran_code/writeOutputFromOutputStructure.f90
+++ b/build/source/actors/file_access_actor/fortran_code/writeOutputFromOutputStructure.f90
@@ -35,6 +35,7 @@ USE globalData,only: integerMissing, realMissing
 USE globalData,only:gru_struc                             ! gru->hru mapping structure
 
 USE output_structure_module,only:outputStructure
+USE output_structure_module,only:outputTimeStep
 
 USE data_types,only:var_i
 
@@ -93,7 +94,6 @@ subroutine writeOutput_fortran(handle_ncid, num_steps, start_gru, max_gru, err)
   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:outputTimeStep
   USE globalData,only:bvar_meta,time_meta,forc_meta,prog_meta,diag_meta,flux_meta,indx_meta
   USE globalData,only:maxLayers
   implicit none
@@ -492,11 +492,11 @@ subroutine writeVector(ncid, outputTimestep, maxLayers, nSteps, minGRU, maxGRU,
       select type (dat)
           class is (gru_hru_time_doubleVec)
               if(.not.outputStructure(1)%finalizeStats%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
-              realArray(gruCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(:)
+              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
-              intArray(gruCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(:)
+              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
 
@@ -520,11 +520,11 @@ subroutine writeVector(ncid, outputTimestep, maxLayers, nSteps, minGRU, maxGRU,
       case(ixReal)
         err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realArray(1:numGRU,1:maxLength),start=(/minGRU,1,stepCounter/),count=(/numGRU,maxLength,1/))
         if(err/=0)then; print*, "ERROR: with nf90_put_var in data vector (ixReal)"; return; endif
-
+        realArray(:,:) = realMissing ! reset the realArray
       case(ixInteger)
         err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),intArray(1:numGRU,1:maxLength),start=(/minGRU,1,stepCounter/),count=(/numGRU,maxLength,1/))
         if(err/=0)then; print*, "ERROR: with nf90_put_var in data vector (ixInteger)"; return; endif
-
+        intArray(:,:) = integerMissing ! reset the intArray
       case default; err=20; message=trim(message)//'data must be of type integer or real'; return
     end select ! data type
     stepCounter = stepCounter + 1
diff --git a/build/source/actors/hru_actor/fortran_code/hru_actor.f90 b/build/source/actors/hru_actor/fortran_code/hru_actor.f90
index b04ee25fb0c5e54eb46521b3d17bbecac26a1a0a..4e08e85645f3051addfb7f469b208c10be8d7bb3 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_actor.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_actor.f90
@@ -330,7 +330,8 @@ subroutine setIDATolerances(handle_hru_data,    &
   type(hru_type),pointer                  :: hru_data          !  model time data
 
   call c_f_pointer(handle_hru_data, hru_data)
-
+  
+#ifdef SUNDIALS_ACTIVE
   hru_data%mparStruct%var(iLookPARAM%relTolTempCas)%dat(1)       = relTolTempCas 
   hru_data%mparStruct%var(iLookPARAM%absTolTempCas)%dat(1)       = absTolTempCas
   hru_data%mparStruct%var(iLookPARAM%relTolTempVeg)%dat(1)       = relTolTempVeg
@@ -345,6 +346,7 @@ subroutine setIDATolerances(handle_hru_data,    &
   hru_data%mparStruct%var(iLookPARAM%absTolMatric)%dat(1)        = absTolMatric
   hru_data%mparStruct%var(iLookPARAM%relTolAquifr)%dat(1)        = relTolAquifr
   hru_data%mparStruct%var(iLookPARAM%absTolAquifr)%dat(1)        = absTolAquifr
-end subroutine setIDATolerances
+#endif
+  end subroutine setIDATolerances
 
 end module hru_actor
\ No newline at end of file
diff --git a/build/source/actors/hru_actor/fortran_code/hru_init.f90 b/build/source/actors/hru_actor/fortran_code/hru_init.f90
index 480eda55fdfcb838a8fba354442fd787b2209bea..7b7ed505b9a5805beb4aa44ea6169488221a2352 100755
--- a/build/source/actors/hru_actor/fortran_code/hru_init.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_init.f90
@@ -23,7 +23,7 @@ USE globalData,only:prog_meta,diag_meta,flux_meta,id_meta   ! metadata structure
 USE globalData,only:mpar_meta,indx_meta                     ! metadata structures
 USE globalData,only:bpar_meta,bvar_meta                     ! metadata structures
 USE globalData,only:averageFlux_meta                        ! metadata for time-step average fluxes
-USE globalData,only:lookup_meta 
+! USE globalData,only:lookup_meta 
 
 ! statistics metadata structures
 USE globalData,only:statForc_meta                           ! child metadata for stats
diff --git a/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90 b/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90
index d035d43129a72175af189249ed6f5e852361de08..76b02b27492d4a5088f8a720636eb520aae91d10 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_modelRun.f90
@@ -114,7 +114,7 @@ subroutine runPhysics(&
   ! ---------------------------------------------------------------------------------------
   ! Dummy Variables
   ! ---------------------------------------------------------------------------------------
-  integer(c_int),intent(in)                :: indxHRU                ! id of HRU                   
+  integer(c_long),intent(in)               :: indxHRU                ! id of HRU                   
   integer(c_int),intent(in)                :: modelTimeStep          ! time step index
   type(c_ptr), intent(in), value           :: handle_hru_data         ! c_ptr to -- hru data
   real(c_double),intent(inout)             :: fracJulDay                    ! fractional julian days since the start of year
@@ -160,19 +160,19 @@ subroutine runPhysics(&
       ! (compute the exposed LAI and SAI and whether veg is buried by snow)
       call vegPhenlgy(&
                       ! model control
-                      fracJulDay,                     & ! intent(in):    fractional julian days since the start of year
-                      yearLength,                     & ! intent(in):    number of days in the current year
+                      fracJulDay,             & ! intent(in):    fractional julian days since the start of year
+                      yearLength,             & ! intent(in):    number of days in the current year
                       ! input/output: data structures
-                      model_decisions,                & ! intent(in):    model decisions
-                      hru_data%typeStruct,                     & ! intent(in):    type of vegetation and soil
-                      hru_data%attrStruct,                     & ! intent(in):    spatial attributes
-                      hru_data%mparStruct,                     & ! intent(in):    model parameters
-                      hru_data%progStruct,                     & ! intent(in):    model prognostic variables for a local HRU
-                      hru_data%diagStruct,                     & ! intent(inout): model diagnostic variables for a local HRU
+                      model_decisions,        & ! intent(in):    model decisions
+                      hru_data%typeStruct,    & ! intent(in):    type of vegetation and soil
+                      hru_data%attrStruct,    & ! intent(in):    spatial attributes
+                      hru_data%mparStruct,    & ! intent(in):    model parameters
+                      hru_data%progStruct,    & ! intent(in):    model prognostic variables for a local HRU
+                      hru_data%diagStruct,    & ! intent(inout): model diagnostic variables for a local HRU
                       ! output
-                      computeVegFluxFlag,             & ! intent(out): flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
-                      notUsed_canopyDepth,            & ! intent(out): NOT USED: canopy depth (m)
-                      notUsed_exposedVAI,             & ! intent(out): NOT USED: exposed vegetation area index (m2 m-2)
+                      computeVegFluxFlag,     & ! intent(out): flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
+                      notUsed_canopyDepth,    & ! intent(out): NOT USED: canopy depth (m)
+                      notUsed_exposedVAI,     & ! intent(out): NOT USED: exposed vegetation area index (m2 m-2)
                       err,cmessage)                     ! intent(out): error control
       if(err/=0)then;message=trim(message)//trim(cmessage); print*, message; return; endif
 
@@ -267,8 +267,8 @@ subroutine runPhysics(&
         hru_data%progStruct,         & ! data structure of model prognostic variables
         hru_data%diagStruct,         & ! data structure of model diagnostic variables
         hru_data%fluxStruct,         & ! data structure of model fluxes
-        tmZoneOffsetFracDay,& ! time zone offset in fractional days
-        err,cmessage)       ! error control
+        tmZoneOffsetFracDay,         & ! time zone offset in fractional days
+        err,cmessage)                  ! error control
   if(err/=0)then;err=20; message=trim(message)//cmessage; print*, message; return; endif
  
   ! initialize the number of flux calls
@@ -277,19 +277,21 @@ subroutine runPhysics(&
   ! run the model for a single HRU
   call coupled_em(&
                   ! model control
-                  indxHRU,            & ! intent(in):    hruID
-                  dt_init,            & ! intent(inout): initial time step
-                  dt_init_factor,     & ! Used to adjust the length of the timestep in the event of a failure
-                  computeVegFluxFlag, & ! intent(inout): flag to indicate if we are computing fluxes over vegetation
-                  fracJulDay,        & ! intent(in):    fractional julian days since the start of year
-                  yearLength,        & ! intent(in):    number of days in the current year
+                  indxHRU,                     & ! intent(in):    hruID
+                  dt_init,                     & ! intent(inout): initial time step
+                  dt_init_factor,              & ! Used to adjust the length of the timestep in the event of a failure
+                  computeVegFluxFlag,          & ! intent(inout): flag to indicate if we are computing fluxes over vegetation
+                  fracJulDay,                  & ! intent(in):    fractional julian days since the start of year
+                  yearLength,                  & ! intent(in):    number of days in the current year
                   ! data structures (input)
                   hru_data%typeStruct,         & ! intent(in):    local classification of soil veg etc. for each HRU
                   hru_data%attrStruct,         & ! intent(in):    local attributes for each HRU
                   hru_data%forcStruct,         & ! intent(in):    model forcing data
                   hru_data%mparStruct,         & ! intent(in):    model parameters
                   hru_data%bvarStruct,         & ! intent(in):    basin-average model variables
+#ifdef SUNDIALS_ACTIVE                  
                   hru_data%lookupStruct,       &
+#endif
                   ! data structures (input-output)
                   hru_data%indxStruct,         & ! intent(inout): model indices
                   hru_data%progStruct,         & ! intent(inout): model prognostic variables for a local HRU
diff --git a/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90 b/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90
index 96c95688dc159e3527bf898de4f190044ca7278d..f562fd11e99abe10f19aa2cfd9ec76f6061c584e 100755
--- a/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_modelwrite.f90
@@ -273,16 +273,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(:) = dat%var(iVar)%dat(:)
+                outputStructure(1)%progStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = dat%var(iVar)%dat(:)
               case('diag')
-                outputStructure(1)%diagStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
+                outputStructure(1)%diagStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = dat%var(iVar)%dat(:)
               case('flux')
-                outputStructure(1)%fluxStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
+                outputStructure(1)%fluxStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = dat%var(iVar)%dat(:)
               case default
                 err=21; message=trim(message)//'data structure not found for output'
             end select
           class is (var_ilength) 
-            outputStructure(1)%indxStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(:) = dat%var(iVar)%dat(:)
+            outputStructure(1)%indxStruct%gru(indxGRU)%hru(indxHRU)%var(iVar)%tim(iStep)%dat(1:datLength) = dat%var(iVar)%dat(:)
           class default; err=20; message=trim(message)//'data must not be scalarv and either of type var_dlength or var_ilength'; return
         end select
 
diff --git a/build/source/actors/hru_actor/fortran_code/hru_setup.f90 b/build/source/actors/hru_actor/fortran_code/hru_setup.f90
index 2a43598b673e202596a73cac30c40295b4a9085b..e014839a39e8b9457d3df3d9f7a889aad2a40826 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_setup.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_setup.f90
@@ -83,7 +83,9 @@ subroutine setupHRUParam(&
   USE paramCheck_module,only:paramCheck                       ! module to check consistency of model parameters
   USE pOverwrite_module,only:pOverwrite                       ! module to overwrite default parameter values with info from the Noah tables
   USE ConvE2Temp_module,only:E2T_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
+#ifdef SUNDIALS_ACTIVE  
   USE t2enthalpy_module,only:T2E_lookup                       ! module to calculate a look-up table for the temperature-enthalpy conversion
+#endif
   USE var_derive_module,only:fracFuture                       ! module to calculate the fraction of runoff in future time steps (time delay histogram)
   USE module_sf_noahmplsm,only:read_mp_veg_parameters         ! module to read NOAH vegetation tables
   ! global data structures
@@ -150,6 +152,7 @@ subroutine setupHRUParam(&
     hru_data%bvarStruct%var(ivar)%dat(:) = outputStructure(1)%bvarStruct_init%gru(indxGRU)%var(ivar)%dat(:)
   enddo
   ! Copy the lookup Struct if its allocated
+#ifdef SUNDIALS_ACTIVE
   if (allocated(outputStructure(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z)) then
     do i_z=1, size(outputStructure(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z(:))
       do iVar=1, size(outputStructure(1)%lookupStruct%gru(indxGRU)%hru(indxHRU)%z(i_z)%var(:))
@@ -157,6 +160,7 @@ subroutine setupHRUParam(&
       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(:)
diff --git a/build/source/netcdf/def_output.f90 b/build/source/netcdf/def_output.f90
index 73df9a12f8ce3d9a852e4ea79c3d5a71b9546876..1b8a9d3cfd5332cbdf4ccab0ea3638d1979ceb36 100755
--- a/build/source/netcdf/def_output.f90
+++ b/build/source/netcdf/def_output.f90
@@ -79,7 +79,6 @@ subroutine def_output(ncid,startGRU,nGRU,nHRU,actor_info,err,message)
   USE globalData,only:bpar_meta,bvar_meta,time_meta            ! metaData structures
   USE globalData,only:model_decisions                          ! model decisions
   USE globalData,only:outFreq                                  ! output frequencies
-  USE globalData,only:fname
   ! Some global variabels required in the writing process
   USE globalData,only:nHRUrun
   USE globalData,only:nGRUrun
@@ -113,6 +112,7 @@ subroutine def_output(ncid,startGRU,nGRU,nHRU,actor_info,err,message)
   integer(i4b)                         :: iGRU
   character(LEN=256)                   :: startGRUString    ! String Variable to convert startGRU
   character(LEN=256)                   :: numGRUString      ! String Varaible to convert numGRU
+  character(len=1024)                  :: fname                         ! temporary filename
 
 
   ! initialize errors
diff --git a/utils/netcdf/OutputVerification/checkOutput.py b/utils/netcdf/OutputVerification/checkOutput.py
index 6bc6809f5d9cbe299507a36a914f1887e6543116..268524640924a2ac91a96bbf634678aeadf6f5e9 100755
--- a/utils/netcdf/OutputVerification/checkOutput.py
+++ b/utils/netcdf/OutputVerification/checkOutput.py
@@ -6,15 +6,7 @@ import numpy as np
 import sys
 
 
-def check_variable_length(hru_from_dataset_1, hru_from_dataset_2, variable):
-    if len(hru_from_dataset_1[variable].values) != len(hru_from_dataset_2[variable].values):
-        print("ERROR: output variable", variable, "does not contain the same amount of data")
-        print("     hru_from_dataset_1 = ", len(hru_from_dataset_1[variable].values))
-        print("     hru_from_dataset_2 = ", len(hru_from_dataset_2[variable].values))
-        return False
-    else:
-        return True
-
+# Get data in usable format
 def extract_variable_data(hru_dataset, var):
     hru_variable_data_from_dataset = []
 
@@ -29,29 +21,49 @@ def extract_variable_data(hru_dataset, var):
             hru_variable_data_from_dataset.append(data)
     return hru_variable_data_from_dataset
 
+
+# Check Functions
+def check_variable_length(hru_from_dataset_1, hru_from_dataset_2, variable):
+    if len(hru_from_dataset_1[variable].values) != len(hru_from_dataset_2[variable].values):
+        print("ERROR: output variable", variable, "does not contain the same amount of data")
+        print("     hru_from_dataset_1 = ", len(hru_from_dataset_1[variable].values))
+        print("     hru_from_dataset_2 = ", len(hru_from_dataset_2[variable].values))
+        return False
+    else:
+        return True
 def check_data_for_errors(dataset_1, dataset_2, tolerance):
     error_counter = 0
     for i in range(0, len(dataset_1)):
-        if abs(dataset_1[i] - dataset_2[i]) > tolerance:
+        if abs(dataset_1[i] - dataset_2[i]) > tolerance: 
             error_counter += 1
+            # Open an error file and append the error to it
+            error_file = open("error_file.txt", "a")
+            error_file.write("     dataset_1 = " + str(dataset_1[i]) + "\n")
+            error_file.write("     dataset_2 = " + str(dataset_2[i]) + "\n")
+            error_file.write("     " + str(i))
+            error_file.write("\n")
+            error_file.close()
+
+
     return error_counter
 
 
 
 def verify_data(dataset_1, dataset_2, num_hru, output_variables):
-   
     dataset_1 = xr.open_dataset(dataset_1)
     dataset_2 = xr.open_dataset(dataset_2)
-
     total_errors = 0
     for i_hru in range(0, num_hru):
         hru_from_dataset_1 = dataset_1.isel(hru=i_hru).copy()
         hru_from_dataset_2 = dataset_2.isel(hru=i_hru).copy()
-        
-        # print("\nHRU - hru_dataset_1", hru_from_dataset_1["hruId"].values)
-        # print("HRU - hru_dataset_2", hru_from_dataset_2["hruId"].values, "\n")
+        print(hru_from_dataset_1)
+        # print("CHECKING HRU", i_hru, "OF", num_hru)
 
         for var in output_variables:
+            # Open an error file and append the error to it
+            error_file = open("error_file.txt", "a")
+            error_file.write("Checking: " + str(var) + "\n")
+            error_file.close()
             if not check_variable_length(hru_from_dataset_1, hru_from_dataset_2, var):
                 print("ERROR: output variable", var, "does not contain the same amount of data")
 
@@ -81,14 +93,14 @@ def get_output_vars(model_output_file):
 
 
 
-num_hru = 50
+num_hru = 1
 print("Checking output for", num_hru, "HRUs")
 
 
 dataset_1 = sys.argv[1]
 dataset_2 = sys.argv[2]
 
-model_output_file = "/project/gwf/gwf_cmt/kck540/domain_NorthAmerica/Summa-Projects/input_data/summa_actors_input/outputControl.txt"
+model_output_file = "/project/gwf/gwf_cmt/kck540/domain_NorthAmerica/Summa-Projects/input_data/summa_actors_input/outputControl_state_vars.txt"
 
 output_vars = get_output_vars(model_output_file)
 verify_data(dataset_1, dataset_2, num_hru, output_vars)