diff --git a/build/source/engine/mDecisions.f90 b/build/source/engine/mDecisions.f90
index 4a863a7dab36fc0615446b574a4a2f73dee87277..72258c04f57a72805cd4a04a246e06aa84957352 100755
--- a/build/source/engine/mDecisions.f90
+++ b/build/source/engine/mDecisions.f90
@@ -19,7 +19,6 @@
 ! along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 module mDecisions_module
-USE data_types,only: var_i
 USE nrtype
 USE var_lookup, only: maxvarDecisions  ! maximum number of decisions
 implicit none
@@ -60,9 +59,6 @@ integer(i4b),parameter,public :: laiScaling           =  72    ! exponential fun
 ! look-up values for the choice of numerical method
 integer(i4b),parameter,public :: bEuler               =  81    ! home-grown backward Euler solution with long time steps
 integer(i4b),parameter,public :: sundials             =  82    ! SUNDIALS/IDA solution
-integer(i4b),parameter,public :: iterative            =  81    ! iterative
-integer(i4b),parameter,public :: nonIterative         =  82    ! non-iterative
-integer(i4b),parameter,public :: iterSurfEnergyBal    =  83    ! iterate only on the surface energy balance
 ! look-up values for method used to compute derivative
 integer(i4b),parameter,public :: numerical            =  91    ! numerical solution
 integer(i4b),parameter,public :: analytical           =  92    ! analytical solution
@@ -153,27 +149,24 @@ integer(i4b),parameter,public :: windUnload           = 322    ! Roesch et al 20
 ! look-up values for the choice of energy equation
 integer(i4b),parameter,public :: enthalpyFD           =  323    ! enthalpyFD
 integer(i4b),parameter,public :: closedForm           =  324    ! closedForm
-! look-up values for the choice of DAE solver
-integer(i4b),parameter,public :: sundialIDA           =  325    ! IDA solver form Sundials package
-integer(i4b),parameter,public :: backwEuler           =  326    ! backward Euler method implemented by Martyn
 ! -----------------------------------------------------------------------------------------------------------
 
 contains
 
- ! ************************************************************************************************
- ! public subroutine mDecisions: save model decisions as named integers
- ! ************************************************************************************************
-subroutine mDecisions(numSteps,err,message)
+! ************************************************************************************************
+! public subroutine mDecisions: save model decisions as named integers
+! ************************************************************************************************
+subroutine mDecisions(err,message)
   ! model time structures
   USE multiconst,only:secprday               ! number of seconds in a day
   USE var_lookup,only:iLookTIME              ! named variables that identify indices in the time structures
+  USE globalData,only:refTime,refJulday      ! reference time
+  USE globalData,only:oldTime                ! time from the previous time step
+  USE globalData,only:startTime,finshTime    ! start/end time of simulation
   USE globalData,only:dJulianStart           ! julian day of start time of simulation
   USE globalData,only:dJulianFinsh           ! julian day of end time of simulation
-  USE globalData,only:startTime,finshTime,refTime
-  USE globalData,only:refJulday              ! julian day of start time of simulation
-  USE globalData,only:refJulday              ! julian day of end time of simulation
   USE globalData,only:data_step              ! length of data step (s)
-  USE globalData,only:numtim                 ! the number of timesteps in the model
+  USE globalData,only:numtim                 ! number of time steps in the simulation
   ! model decision structures
   USE globaldata,only:model_decisions        ! model decision structure
   USE var_lookup,only:iLookDECISIONS         ! named variables for elements of the decision structure
@@ -188,31 +181,25 @@ subroutine mDecisions(numSteps,err,message)
   USE time_utils_module,only:extractTime     ! extract time info from units string
   USE time_utils_module,only:compjulday      ! compute the julian day
   USE time_utils_module,only:fracDay         ! compute fractional day
-  USE summaActors_FileManager,only: SIM_START_TM, SIM_END_TM   ! time info from control file module
-
+  USE summaFileManager,only: SIM_START_TM, SIM_END_TM   ! time info from control file module
 
   implicit none
   ! define output
-  integer(i4b),intent(inout)           :: numSteps             
   integer(i4b),intent(out)             :: err            ! error code
   character(*),intent(out)             :: message        ! error message
   ! define local variables
   character(len=256)                   :: cmessage       ! error message for downwind routine
-  real(dp)                             :: dsec,dsec_tz   ! second
+  real(rkind)                             :: dsec,dsec_tz   ! second
   ! initialize error control
-  err=0;message='mDecisions/'
-
-  ! -------------------------------------------------------------------------------------------------
-  ! -------------------------------------------------------------------------------------------------
+  err=0; message='mDecisions/'
 
   ! read information from model decisions file, and populate model decisions structure
   call readoption(err,cmessage)
   if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
 
-  ! -------------------------------------------------------------------------------------------------
 
   ! put reference time information into the time structures
-  call extractTime(forc_meta(iLookFORCE%time)%varunit,                    & ! date-time string
+  call extractTime(forc_meta(iLookFORCE%time)%varunit,                     & ! date-time string
                     refTime%var(iLookTIME%iyyy),                           & ! year
                     refTime%var(iLookTIME%im),                             & ! month
                     refTime%var(iLookTIME%id),                             & ! day
@@ -227,14 +214,14 @@ subroutine mDecisions(numSteps,err,message)
 
   ! compute the julian date (fraction of day) for the reference time
   call compjulday(&
-                  refTime%var(iLookTIME%iyyy),                           & ! year
-                  refTime%var(iLookTIME%im),                             & ! month
-                  refTime%var(iLookTIME%id),                             & ! day
-                  refTime%var(iLookTIME%ih),                             & ! hour
-                  refTime%var(iLookTIME%imin),                           & ! minute
-                  0._dp,                                                 & ! second
-                  refJulday,                                             & ! julian date for the start of the simulation
-                  err, cmessage)                                           ! error control
+                  refTime%var(iLookTIME%iyyy),                            & ! year
+                  refTime%var(iLookTIME%im),                              & ! month
+                  refTime%var(iLookTIME%id),                              & ! day
+                  refTime%var(iLookTIME%ih),                              & ! hour
+                  refTime%var(iLookTIME%imin),                            & ! minute
+                  0._rkind,                                               & ! second
+                  refJulday,                                              & ! julian date for the start of the simulation
+                  err, cmessage)                                            ! error control
   if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
 
   ! put simulation start time information into the time structures
@@ -253,18 +240,18 @@ subroutine mDecisions(numSteps,err,message)
 
   ! compute the julian date (fraction of day) for the start of the simulation
   call compjulday(&
-                  startTime%var(iLookTIME%iyyy),                         & ! year
-                  startTime%var(iLookTIME%im),                           & ! month
-                  startTime%var(iLookTIME%id),                           & ! day
-                  startTime%var(iLookTIME%ih),                           & ! hour
-                  startTime%var(iLookTIME%imin),                         & ! minute
-                  0._dp,                                                 & ! second
-                  dJulianStart,                                          & ! julian date for the start of the simulation
-                  err, cmessage)                                           ! error control
+                  startTime%var(iLookTIME%iyyy),                           & ! year
+                  startTime%var(iLookTIME%im),                             & ! month
+                  startTime%var(iLookTIME%id),                             & ! day
+                  startTime%var(iLookTIME%ih),                             & ! hour
+                  startTime%var(iLookTIME%imin),                           & ! minute
+                  0._rkind,                                                & ! second
+                  dJulianStart,                                            & ! julian date for the start of the simulation
+                  err, cmessage)                                             ! error control
   if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
 
   ! put simulation end time information into the time structures
-  call extractTime(trim(SIM_END_TM),                                      & ! date-time string
+  call extractTime(trim(SIM_END_TM),                                       & ! date-time string
                     finshTime%var(iLookTIME%iyyy),                         & ! year
                     finshTime%var(iLookTIME%im),                           & ! month
                     finshTime%var(iLookTIME%id),                           & ! day
@@ -284,27 +271,25 @@ subroutine mDecisions(numSteps,err,message)
                   finshTime%var(iLookTIME%id),                           & ! day
                   finshTime%var(iLookTIME%ih),                           & ! hour
                   finshTime%var(iLookTIME%imin),                         & ! minute
-                  0._dp,                                                 & ! second
+                  0._rkind,                                              & ! second
                   dJulianFinsh,                                          & ! julian date for the end of the simulation
                   err, cmessage)                                           ! error control
   if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
 
   ! check start and finish time
-  !  write(*,'(a,i4,1x,4(i2,1x))') 'startTime: iyyy, im, id, ih, imin = ', startTime%var(1:5)
-  !  write(*,'(a,i4,1x,4(i2,1x))') 'finshTime: iyyy, im, id, ih, imin = ', finshTime%var(1:5)
+  write(*,'(a,i4,1x,4(i2,1x))') 'startTime: iyyy, im, id, ih, imin = ', startTime%var(1:5)
+  write(*,'(a,i4,1x,4(i2,1x))') 'finshTime: iyyy, im, id, ih, imin = ', finshTime%var(1:5)
 
   ! check that simulation end time is > start time
   if(dJulianFinsh < dJulianStart)then; err=20; message=trim(message)//'end time of simulation occurs before start time'; return; end if
 
   ! initialize the old time vector (time from the previous time step)
+  oldTime%var(:) = startTime%var(:)
 
   ! compute the number of time steps
-  numSteps = nint( (dJulianFinsh - dJulianStart)*secprday/data_step ) + 1
-  numTim = numSteps
-
-  !  write(*,'(a,1x,i10)') 'number of time steps = ', numSteps
+  numtim = nint( (dJulianFinsh - dJulianStart)*secprday/data_step ) + 1
+  write(*,'(a,1x,i10)') 'number of time steps = ', numtim
 
-  ! -------------------------------------------------------------------------------------------------
 
   ! set Noah-MP options
   DVEG=3      ! option for dynamic vegetation
@@ -322,7 +307,7 @@ subroutine mDecisions(numSteps,err,message)
     case('CLM_Type'); model_decisions(iLookDECISIONS%soilStress)%iDecision = CLM_Type             ! thresholded linear function of matric head
     case('SiB_Type'); model_decisions(iLookDECISIONS%soilStress)%iDecision = SiB_Type             ! exponential of the log of matric head
     case default
-    err=10; message=trim(message)//"unknown soil moisture function [option="//trim(model_decisions(iLookDECISIONS%soilStress)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown soil moisture function [option="//trim(model_decisions(iLookDECISIONS%soilStress)%cDecision)//"]"; return
   end select
 
   ! identify the choice of function for stomatal resistance
@@ -333,129 +318,103 @@ subroutine mDecisions(numSteps,err,message)
     case('BallBerryFlex'      ); model_decisions(iLookDECISIONS%stomResist)%iDecision = BallBerryFlex       ! flexible Ball-Berry scheme
     case('BallBerryTest'      ); model_decisions(iLookDECISIONS%stomResist)%iDecision = BallBerryTest       ! flexible Ball-Berry scheme (testing)
     case default
-    err=10; message=trim(message)//"unknown stomatal resistance function [option="//trim(model_decisions(iLookDECISIONS%stomResist)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown stomatal resistance function [option="//trim(model_decisions(iLookDECISIONS%stomResist)%cDecision)//"]"; return
   end select
 
   ! identify the leaf temperature controls on photosynthesis + stomatal resistance
   if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then
     select case(trim(model_decisions(iLookDECISIONS%bbTempFunc)%cDecision))
-    case('q10Func'            ); model_decisions(iLookDECISIONS%bbTempFunc)%iDecision = q10Func
-    case('Arrhenius'          ); model_decisions(iLookDECISIONS%bbTempFunc)%iDecision = Arrhenius
-    case default
-      err=10; message=trim(message)//"unknown leaf temperature function [option="//trim(model_decisions(iLookDECISIONS%bbTempFunc)%cDecision)//"]"; return
+      case('q10Func'            ); model_decisions(iLookDECISIONS%bbTempFunc)%iDecision = q10Func
+      case('Arrhenius'          ); model_decisions(iLookDECISIONS%bbTempFunc)%iDecision = Arrhenius
+      case default
+        err=10; message=trim(message)//"unknown leaf temperature function [option="//trim(model_decisions(iLookDECISIONS%bbTempFunc)%cDecision)//"]"; return
     end select
   end if
 
   ! identify the humidity controls on stomatal resistance
   if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then
     select case(trim(model_decisions(iLookDECISIONS%bbHumdFunc)%cDecision))
-    case('humidLeafSurface'   ); model_decisions(iLookDECISIONS%bbHumdFunc)%iDecision = humidLeafSurface
-    case('scaledHyperbolic'   ); model_decisions(iLookDECISIONS%bbHumdFunc)%iDecision = scaledHyperbolic
-    case default
-      err=10; message=trim(message)//"unknown humidity function [option="//trim(model_decisions(iLookDECISIONS%bbHumdFunc)%cDecision)//"]"; return
+      case('humidLeafSurface'   ); model_decisions(iLookDECISIONS%bbHumdFunc)%iDecision = humidLeafSurface
+      case('scaledHyperbolic'   ); model_decisions(iLookDECISIONS%bbHumdFunc)%iDecision = scaledHyperbolic
+      case default
+        err=10; message=trim(message)//"unknown humidity function [option="//trim(model_decisions(iLookDECISIONS%bbHumdFunc)%cDecision)//"]"; return
     end select
   end if
 
   ! identify functions for electron transport function (dependence of photosynthesis on PAR)
   if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then
     select case(trim(model_decisions(iLookDECISIONS%bbElecFunc)%cDecision))
-    case('linear'             ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = linear
-    case('linearJmax'         ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = linearJmax
-    case('quadraticJmax'      ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = quadraticJmax
-    case default
-      err=10; message=trim(message)//"unknown electron transport function [option="//trim(model_decisions(iLookDECISIONS%bbElecFunc)%cDecision)//"]"; return
+      case('linear'             ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = linear
+      case('linearJmax'         ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = linearJmax
+      case('quadraticJmax'      ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = quadraticJmax
+      case default
+        err=10; message=trim(message)//"unknown electron transport function [option="//trim(model_decisions(iLookDECISIONS%bbElecFunc)%cDecision)//"]"; return
     end select
   end if
 
   ! identify the use of the co2 compensation point in the stomatal conductance calaculations
   if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then
     select case(trim(model_decisions(iLookDECISIONS%bbCO2point)%cDecision))
-    case('origBWB'            ); model_decisions(iLookDECISIONS%bbCO2point)%iDecision = origBWB
-    case('Leuning'            ); model_decisions(iLookDECISIONS%bbCO2point)%iDecision = Leuning
-    case default
-      err=10; message=trim(message)//"unknown option for the co2 compensation point [option="//trim(model_decisions(iLookDECISIONS%bbCO2point)%cDecision)//"]"; return
+      case('origBWB'            ); model_decisions(iLookDECISIONS%bbCO2point)%iDecision = origBWB
+      case('Leuning'            ); model_decisions(iLookDECISIONS%bbCO2point)%iDecision = Leuning
+      case default
+        err=10; message=trim(message)//"unknown option for the co2 compensation point [option="//trim(model_decisions(iLookDECISIONS%bbCO2point)%cDecision)//"]"; return
     end select
   end if
 
   ! identify the iterative numerical solution method used in the Ball-Berry stomatal resistance parameterization
   if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then
     select case(trim(model_decisions(iLookDECISIONS%bbNumerics)%cDecision))
-    case('NoahMPsolution'     ); model_decisions(iLookDECISIONS%bbNumerics)%iDecision = NoahMPsolution  ! the NoahMP solution (and CLM4): fixed point iteration; max 3 iterations
-    case('newtonRaphson'      ); model_decisions(iLookDECISIONS%bbNumerics)%iDecision = newtonRaphson   ! full Newton-Raphson iterative solution to convergence
-    case default
-      err=10; message=trim(message)//"unknown option for the Ball-Berry numerical solution [option="//trim(model_decisions(iLookDECISIONS%bbNumerics)%cDecision)//"]"; return
+      case('NoahMPsolution'     ); model_decisions(iLookDECISIONS%bbNumerics)%iDecision = NoahMPsolution  ! the NoahMP solution (and CLM4): fixed point iteration; max 3 iterations
+      case('newtonRaphson'      ); model_decisions(iLookDECISIONS%bbNumerics)%iDecision = newtonRaphson   ! full Newton-Raphson iterative solution to convergence
+      case default
+        err=10; message=trim(message)//"unknown option for the Ball-Berry numerical solution [option="//trim(model_decisions(iLookDECISIONS%bbNumerics)%cDecision)//"]"; return
     end select
   end if
 
   ! identify the controls on carbon assimilation
   if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then
     select case(trim(model_decisions(iLookDECISIONS%bbAssimFnc)%cDecision))
-    case('colimitation'       ); model_decisions(iLookDECISIONS%bbAssimFnc)%iDecision = colimitation    ! enable colimitation, as described by Collatz et al. (1991) and Sellers et al. (1996)
-    case('minFunc'            ); model_decisions(iLookDECISIONS%bbAssimFnc)%iDecision = minFunc         ! do not enable colimitation: use minimum of the three controls on carbon assimilation
-    case default
-      err=10; message=trim(message)//"unknown option for the controls on carbon assimilation [option="//trim(model_decisions(iLookDECISIONS%bbAssimFnc)%cDecision)//"]"; return
+      case('colimitation'       ); model_decisions(iLookDECISIONS%bbAssimFnc)%iDecision = colimitation    ! enable colimitation, as described by Collatz et al. (1991) and Sellers et al. (1996)
+      case('minFunc'            ); model_decisions(iLookDECISIONS%bbAssimFnc)%iDecision = minFunc         ! do not enable colimitation: use minimum of the three controls on carbon assimilation
+      case default
+        err=10; message=trim(message)//"unknown option for the controls on carbon assimilation [option="//trim(model_decisions(iLookDECISIONS%bbAssimFnc)%cDecision)//"]"; return
     end select
   end if
 
   ! identify the scaling of photosynthesis from the leaf to the canopy
   if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then
     select case(trim(model_decisions(iLookDECISIONS%bbCanIntg8)%cDecision))
-    case('constantScaling'    ); model_decisions(iLookDECISIONS%bbCanIntg8)%iDecision = constantScaling ! constant scaling factor
-    case('laiScaling'         ); model_decisions(iLookDECISIONS%bbCanIntg8)%iDecision = laiScaling      ! exponential function of LAI (Leuning, Plant Cell Env 1995: "Scaling from..." [eq 9])
-    case default
-      err=10; message=trim(message)//"unknown option for scaling of photosynthesis from the leaf to the canopy [option="//trim(model_decisions(iLookDECISIONS%bbCanIntg8)%cDecision)//"]"; return
+      case('constantScaling'    ); model_decisions(iLookDECISIONS%bbCanIntg8)%iDecision = constantScaling ! constant scaling factor
+      case('laiScaling'         ); model_decisions(iLookDECISIONS%bbCanIntg8)%iDecision = laiScaling      ! exponential function of LAI (Leuning, Plant Cell Env 1995: "Scaling from..." [eq 9])
+      case default
+        err=10; message=trim(message)//"unknown option for scaling of photosynthesis from the leaf to the canopy [option="//trim(model_decisions(iLookDECISIONS%bbCanIntg8)%cDecision)//"]"; return
     end select
   end if
 
-  ! ************************************
-  ! ************************************
-  ! ************************************
-  ! ************************************
-  ! SUNDIALS ADDITIONS
-
   ! identify the numerical method
   select case(trim(model_decisions(iLookDECISIONS%num_method)%cDecision))
-    case('bEuler');   model_decisions(iLookDECISIONS%num_method)%iDecision = bEuler 
-    case('sundials'); model_decisions(iLookDECISIONS%num_method)%iDecision = sundials
-    case('itertive'); model_decisions(iLookDECISIONS%num_method)%iDecision = iterative           ! iterative
-    case('non_iter'); model_decisions(iLookDECISIONS%num_method)%iDecision = nonIterative        ! non-iterative
-    case('itersurf'); model_decisions(iLookDECISIONS%num_method)%iDecision = iterSurfEnergyBal   ! iterate only on the surface energy balance
+    case('bEuler'   ); model_decisions(iLookDECISIONS%num_method)%iDecision = bEuler             ! home-grown backward Euler solution with long time steps
+    case('itertive' ); model_decisions(iLookDECISIONS%num_method)%iDecision = bEuler             ! home-grown backward Euler solution (included for backwards compatibility)
+    case('sundials' ); model_decisions(iLookDECISIONS%num_method)%iDecision = sundials           ! SUNDIALS/IDA solution
     case default
-    err=10; message=trim(message)//"unknown numerical method [option="//trim(model_decisions(iLookDECISIONS%num_method)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown numerical method [option="//trim(model_decisions(iLookDECISIONS%num_method)%cDecision)//"]"; return
   end select
 
-    ! how to compute heat capacity in energy equation
+  ! how to compute heat capacity in energy equation
   select case(trim(model_decisions(iLookDECISIONS%howHeatCap)%cDecision))
     case('enthalpyFD'); model_decisions(iLookDECISIONS%howHeatCap)%iDecision = enthalpyFD        ! enthalpyFD
     case('closedForm'); model_decisions(iLookDECISIONS%howHeatCap)%iDecision = closedForm        ! closedForm
     case default
-      ! TODO: after adding howHeatCap decision in corresponding file we should delete the next line
-      model_decisions(iLookDECISIONS%howHeatCap)%iDecision = closedForm
-      ! err=10; message=trim(message)//"unknown Cp computation [option="//trim(model_decisions(iLookDECISIONS%howHeatCap)%cDecision)//"]"; return
-  end select
-
-    ! how to solve the system of differential equations
-  select case(trim(model_decisions(iLookDECISIONS%diffEqSolv)%cDecision))
-    case('sundialIDA'); model_decisions(iLookDECISIONS%diffEqSolv)%iDecision = sundialIDA        ! enthalpyFD
-    case('backwEuler'); model_decisions(iLookDECISIONS%diffEqSolv)%iDecision = backwEuler        ! closedForm
-    case default
-      ! TODO: after adding diffEqSolv decision in corresponding file we should delete the next line
-      model_decisions(iLookDECISIONS%diffEqSolv)%iDecision = sundialIDA
-      ! err=10; message=trim(message)//"unknown DAE solver [option="//trim(model_decisions(iLookDECISIONS%diffEqSolv)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown Cp computation [option="//trim(model_decisions(iLookDECISIONS%howHeatCap)%cDecision)//"]"; return
   end select
 
-  ! ************************************
-  ! ************************************
-  ! ************************************
-  ! ************************************
-  ! END SUNDIALS ADDITTION
-
   ! identify the method used to calculate flux derivatives
   select case(trim(model_decisions(iLookDECISIONS%fDerivMeth)%cDecision))
     case('numericl'); model_decisions(iLookDECISIONS%fDerivMeth)%iDecision = numerical           ! numerical
     case('analytic'); model_decisions(iLookDECISIONS%fDerivMeth)%iDecision = analytical          ! analytical
     case default
-    err=10; message=trim(message)//"unknown method used to calculate flux derivatives [option="//trim(model_decisions(iLookDECISIONS%fDerivMeth)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown method used to calculate flux derivatives [option="//trim(model_decisions(iLookDECISIONS%fDerivMeth)%cDecision)//"]"; return
   end select
 
   ! identify the method used to determine LAI and SAI
@@ -463,7 +422,7 @@ subroutine mDecisions(numSteps,err,message)
     case('monTable');  model_decisions(iLookDECISIONS%LAI_method)%iDecision = monthlyTable       ! LAI/SAI taken directly from a monthly table for different vegetation classes
     case('specified'); model_decisions(iLookDECISIONS%LAI_method)%iDecision = specified          ! LAI/SAI computed from green vegetation fraction and winterSAI and summerLAI parameters
     case default
-    err=10; message=trim(message)//"unknown method to determine LAI and SAI [option="//trim(model_decisions(iLookDECISIONS%LAI_method)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown method to determine LAI and SAI [option="//trim(model_decisions(iLookDECISIONS%LAI_method)%cDecision)//"]"; return
   end select
 
   ! identify the canopy interception parameterization
@@ -472,7 +431,7 @@ subroutine mDecisions(numSteps,err,message)
     case('sparseCanopy');    model_decisions(iLookDECISIONS%cIntercept)%iDecision = sparseCanopy
     case('storageFunc');     model_decisions(iLookDECISIONS%cIntercept)%iDecision = storageFunc
     case default
-    err=10; message=trim(message)//"unknown canopy interception parameterization [option="//trim(model_decisions(iLookDECISIONS%cIntercept)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown canopy interception parameterization [option="//trim(model_decisions(iLookDECISIONS%cIntercept)%cDecision)//"]"; return
   end select
 
   ! identify the form of Richards' equation
@@ -480,7 +439,7 @@ subroutine mDecisions(numSteps,err,message)
     case('moisture'); model_decisions(iLookDECISIONS%f_Richards)%iDecision = moisture            ! moisture-based form
     case('mixdform'); model_decisions(iLookDECISIONS%f_Richards)%iDecision = mixdform            ! mixed form
     case default
-    err=10; message=trim(message)//"unknown form of Richards' equation [option="//trim(model_decisions(iLookDECISIONS%f_Richards)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown form of Richards' equation [option="//trim(model_decisions(iLookDECISIONS%f_Richards)%cDecision)//"]"; return
   end select
 
   ! identify the groundwater parameterization
@@ -489,7 +448,7 @@ subroutine mDecisions(numSteps,err,message)
     case('bigBuckt'); model_decisions(iLookDECISIONS%groundwatr)%iDecision = bigBucket           ! a big bucket (lumped aquifer model)
     case('noXplict'); model_decisions(iLookDECISIONS%groundwatr)%iDecision = noExplicit          ! no explicit groundwater parameterization
     case default
-    err=10; message=trim(message)//"unknown groundwater parameterization [option="//trim(model_decisions(iLookDECISIONS%groundwatr)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown groundwater parameterization [option="//trim(model_decisions(iLookDECISIONS%groundwatr)%cDecision)//"]"; return
   end select
 
   ! identify the hydraulic conductivity profile
@@ -497,7 +456,7 @@ subroutine mDecisions(numSteps,err,message)
     case('constant'); model_decisions(iLookDECISIONS%hc_profile)%iDecision = constant            ! constant hydraulic conductivity with depth
     case('pow_prof'); model_decisions(iLookDECISIONS%hc_profile)%iDecision = powerLaw_profile    ! power-law profile
     case default
-    err=10; message=trim(message)//"unknown hydraulic conductivity profile [option="//trim(model_decisions(iLookDECISIONS%hc_profile)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown hydraulic conductivity profile [option="//trim(model_decisions(iLookDECISIONS%hc_profile)%cDecision)//"]"; return
   end select
 
   ! identify the upper boundary conditions for thermodynamics
@@ -506,7 +465,7 @@ subroutine mDecisions(numSteps,err,message)
     case('nrg_flux'); model_decisions(iLookDECISIONS%bcUpprTdyn)%iDecision = energyFlux          ! energy flux
     case('zeroFlux'); model_decisions(iLookDECISIONS%bcUpprTdyn)%iDecision = zeroFlux            ! zero flux
     case default
-    err=10; message=trim(message)//"unknown upper boundary conditions for thermodynamics [option="//trim(model_decisions(iLookDECISIONS%bcUpprTdyn)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown upper boundary conditions for thermodynamics [option="//trim(model_decisions(iLookDECISIONS%bcUpprTdyn)%cDecision)//"]"; return
   end select
 
   ! identify the lower boundary conditions for thermodynamics
@@ -514,7 +473,7 @@ subroutine mDecisions(numSteps,err,message)
     case('presTemp'); model_decisions(iLookDECISIONS%bcLowrTdyn)%iDecision = prescribedTemp      ! prescribed temperature
     case('zeroFlux'); model_decisions(iLookDECISIONS%bcLowrTdyn)%iDecision = zeroFlux            ! zero flux
     case default
-    err=10; message=trim(message)//"unknown lower boundary conditions for thermodynamics [option="//trim(model_decisions(iLookDECISIONS%bcLowrTdyn)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown lower boundary conditions for thermodynamics [option="//trim(model_decisions(iLookDECISIONS%bcLowrTdyn)%cDecision)//"]"; return
   end select
 
   ! identify the upper boundary conditions for soil hydrology
@@ -522,7 +481,7 @@ subroutine mDecisions(numSteps,err,message)
     case('presHead'); model_decisions(iLookDECISIONS%bcUpprSoiH)%iDecision = prescribedHead      ! prescribed head (volumetric liquid water content for mixed form of Richards' eqn)
     case('liq_flux'); model_decisions(iLookDECISIONS%bcUpprSoiH)%iDecision = liquidFlux          ! liquid water flux
     case default
-    err=10; message=trim(message)//"unknown upper boundary conditions for soil hydrology [option="//trim(model_decisions(iLookDECISIONS%bcUpprSoiH)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown upper boundary conditions for soil hydrology [option="//trim(model_decisions(iLookDECISIONS%bcUpprSoiH)%cDecision)//"]"; return
   end select
 
   ! identify the lower boundary conditions for soil hydrology
@@ -532,7 +491,7 @@ subroutine mDecisions(numSteps,err,message)
     case('drainage'); model_decisions(iLookDECISIONS%bcLowrSoiH)%iDecision = freeDrainage        ! free drainage
     case('zeroFlux'); model_decisions(iLookDECISIONS%bcLowrSoiH)%iDecision = zeroFlux            ! zero flux
     case default
-    err=10; message=trim(message)//"unknown lower boundary conditions for soil hydrology [option="//trim(model_decisions(iLookDECISIONS%bcLowrSoiH)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown lower boundary conditions for soil hydrology [option="//trim(model_decisions(iLookDECISIONS%bcLowrSoiH)%cDecision)//"]"; return
   end select
 
   ! identify the choice of parameterization for vegetation roughness length and displacement height
@@ -541,7 +500,7 @@ subroutine mDecisions(numSteps,err,message)
     case('CM_QJRMS1988'   ); model_decisions(iLookDECISIONS%veg_traits)%iDecision = CM_QJRMS1988     ! Choudhury and Monteith (QJRMS 1998) "A four layer model for the heat budget..."
     case('vegTypeTable'   ); model_decisions(iLookDECISIONS%veg_traits)%iDecision = vegTypeTable     ! constant parameters dependent on the vegetation type
     case default
-    err=10; message=trim(message)//"unknown parameterization for vegetation roughness length and displacement height [option="//trim(model_decisions(iLookDECISIONS%veg_traits)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown parameterization for vegetation roughness length and displacement height [option="//trim(model_decisions(iLookDECISIONS%veg_traits)%cDecision)//"]"; return
   end select
 
   ! identify the choice of parameterization for the rooting profile
@@ -550,7 +509,7 @@ subroutine mDecisions(numSteps,err,message)
     case('powerLaw','notPopulatedYet');  model_decisions(iLookDECISIONS%rootProfil)%iDecision = powerLaw      ! simple power-law rooting profile
     case('doubleExp');                   model_decisions(iLookDECISIONS%rootProfil)%iDecision = doubleExp     ! the double exponential function of Xeng et al. (JHM 2001)
     case default
-    err=10; message=trim(message)//"unknown parameterization for rooting profile [option="//trim(model_decisions(iLookDECISIONS%rootProfil)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown parameterization for rooting profile [option="//trim(model_decisions(iLookDECISIONS%rootProfil)%cDecision)//"]"; return
   end select
 
   ! identify the choice of parameterization for canopy emissivity
@@ -558,7 +517,7 @@ subroutine mDecisions(numSteps,err,message)
     case('simplExp'); model_decisions(iLookDECISIONS%canopyEmis)%iDecision = simplExp            ! simple exponential function
     case('difTrans'); model_decisions(iLookDECISIONS%canopyEmis)%iDecision = difTrans            ! parameterized as a function of diffuse transmissivity
     case default
-    err=10; message=trim(message)//"unknown parameterization for canopy emissivity [option="//trim(model_decisions(iLookDECISIONS%canopyEmis)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown parameterization for canopy emissivity [option="//trim(model_decisions(iLookDECISIONS%canopyEmis)%cDecision)//"]"; return
   end select
 
   ! choice of parameterization for snow interception
@@ -566,7 +525,7 @@ subroutine mDecisions(numSteps,err,message)
     case('stickySnow'); model_decisions(iLookDECISIONS%snowIncept)%iDecision = stickySnow        ! maximum interception capacity an increasing function of temerature
     case('lightSnow' ); model_decisions(iLookDECISIONS%snowIncept)%iDecision = lightSnow         ! maximum interception capacity an inverse function of new snow density
     case default
-    err=10; message=trim(message)//"unknown option for snow interception capacity[option="//trim(model_decisions(iLookDECISIONS%snowIncept)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for snow interception capacity[option="//trim(model_decisions(iLookDECISIONS%snowIncept)%cDecision)//"]"; return
   end select
 
   ! identify the choice of wind profile
@@ -574,7 +533,7 @@ subroutine mDecisions(numSteps,err,message)
     case('exponential'   ); model_decisions(iLookDECISIONS%windPrfile)%iDecision = exponential      ! exponential wind profile extends to the surface
     case('logBelowCanopy'); model_decisions(iLookDECISIONS%windPrfile)%iDecision = logBelowCanopy   ! logarithmic profile below the vegetation canopy
     case default
-    err=10; message=trim(message)//"unknown option for choice of wind profile[option="//trim(model_decisions(iLookDECISIONS%windPrfile)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for choice of wind profile[option="//trim(model_decisions(iLookDECISIONS%windPrfile)%cDecision)//"]"; return
   end select
 
   ! identify the choice of atmospheric stability function
@@ -583,7 +542,7 @@ subroutine mDecisions(numSteps,err,message)
     case('louisinv'); model_decisions(iLookDECISIONS%astability)%iDecision = louisInversePower   ! Louis (1979) inverse power function
     case('mahrtexp'); model_decisions(iLookDECISIONS%astability)%iDecision = mahrtExponential    ! Mahrt (1987) exponential
     case default
-    err=10; message=trim(message)//"unknown stability function [option="//trim(model_decisions(iLookDECISIONS%astability)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown stability function [option="//trim(model_decisions(iLookDECISIONS%astability)%cDecision)//"]"; return
   end select
 
   ! choice of canopy shortwave radiation method
@@ -594,7 +553,7 @@ subroutine mDecisions(numSteps,err,message)
     case('NL_scatter' ); model_decisions(iLookDECISIONS%canopySrad)%iDecision = NL_scatter       ! Simplified method Nijssen and Lettenmaier (JGR 1999)
     case('BeersLaw'   ); model_decisions(iLookDECISIONS%canopySrad)%iDecision = BeersLaw         ! Beer's Law (as implemented in VIC)
     case default
-    err=10; message=trim(message)//"unknown canopy radiation method [option="//trim(model_decisions(iLookDECISIONS%canopySrad)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown canopy radiation method [option="//trim(model_decisions(iLookDECISIONS%canopySrad)%cDecision)//"]"; return
   end select
 
   ! choice of albedo representation
@@ -602,7 +561,7 @@ subroutine mDecisions(numSteps,err,message)
     case('conDecay'); model_decisions(iLookDECISIONS%alb_method)%iDecision = constantDecay       ! constant decay (e.g., VIC, CLASS)
     case('varDecay'); model_decisions(iLookDECISIONS%alb_method)%iDecision = variableDecay       ! variable decay (e.g., BATS approach, with destructive metamorphism + soot content)
     case default
-    err=10; message=trim(message)//"unknown option for snow albedo [option="//trim(model_decisions(iLookDECISIONS%alb_method)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for snow albedo [option="//trim(model_decisions(iLookDECISIONS%alb_method)%cDecision)//"]"; return
   end select
 
   ! choice of snow compaction routine
@@ -610,7 +569,7 @@ subroutine mDecisions(numSteps,err,message)
     case('consettl'); model_decisions(iLookDECISIONS%compaction)%iDecision = constantSettlement  ! constant settlement rate
     case('anderson'); model_decisions(iLookDECISIONS%compaction)%iDecision = andersonEmpirical   ! semi-empirical method of Anderson (1976)
     case default
-    err=10; message=trim(message)//"unknown option for snow compaction [option="//trim(model_decisions(iLookDECISIONS%compaction)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for snow compaction [option="//trim(model_decisions(iLookDECISIONS%compaction)%cDecision)//"]"; return
   end select
 
   ! choice of method to combine and sub-divide snow layers
@@ -618,7 +577,7 @@ subroutine mDecisions(numSteps,err,message)
     case('jrdn1991'); model_decisions(iLookDECISIONS%snowLayers)%iDecision = sameRulesAllLayers    ! SNTHERM option: same combination/sub-dividion rules applied to all layers
     case('CLM_2010'); model_decisions(iLookDECISIONS%snowLayers)%iDecision = rulesDependLayerIndex ! CLM option: combination/sub-dividion rules depend on layer index
     case default
-    err=10; message=trim(message)//"unknown option for combination/sub-division of snow layers [option="//trim(model_decisions(iLookDECISIONS%snowLayers)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for combination/sub-division of snow layers [option="//trim(model_decisions(iLookDECISIONS%snowLayers)%cDecision)//"]"; return
   end select
 
   ! choice of thermal conductivity representation for snow
@@ -628,7 +587,7 @@ subroutine mDecisions(numSteps,err,message)
     case('jrdn1991'); model_decisions(iLookDECISIONS%thCondSnow)%iDecision = Jordan1991          ! Jordan (1991)
     case('smnv2000'); model_decisions(iLookDECISIONS%thCondSnow)%iDecision = Smirnova2000        ! Smirnova et al. (2000)
     case default
-    err=10; message=trim(message)//"unknown option for thermal conductivity of snow [option="//trim(model_decisions(iLookDECISIONS%thCondSnow)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for thermal conductivity of snow [option="//trim(model_decisions(iLookDECISIONS%thCondSnow)%cDecision)//"]"; return
   end select
 
   ! choice of thermal conductivity representation for soil
@@ -637,7 +596,7 @@ subroutine mDecisions(numSteps,err,message)
     case('mixConstit' ); model_decisions(iLookDECISIONS%thCondSoil)%iDecision = mixConstit       ! mixture of constituents
     case('hanssonVZJ' ); model_decisions(iLookDECISIONS%thCondSoil)%iDecision = hanssonVZJ       ! test case for the mizoguchi lab experiment, Hansson et al. VZJ 2004
     case default
-    err=10; message=trim(message)//"unknown option for thermal conductivity of soil [option="//trim(model_decisions(iLookDECISIONS%thCondSoil)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for thermal conductivity of soil [option="//trim(model_decisions(iLookDECISIONS%thCondSoil)%cDecision)//"]"; return
   end select
 
   ! choice of method for the spatial representation of groundwater
@@ -645,7 +604,7 @@ subroutine mDecisions(numSteps,err,message)
     case('localColumn'); model_decisions(iLookDECISIONS%spatial_gw)%iDecision = localColumn       ! separate groundwater in each local soil column
     case('singleBasin'); model_decisions(iLookDECISIONS%spatial_gw)%iDecision = singleBasin       ! single groundwater store over the entire basin
     case default
-    err=10; message=trim(message)//"unknown option for spatial representation of groundwater [option="//trim(model_decisions(iLookDECISIONS%spatial_gw)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for spatial representation of groundwater [option="//trim(model_decisions(iLookDECISIONS%spatial_gw)%cDecision)//"]"; return
   end select
 
   ! choice of routing method
@@ -653,7 +612,7 @@ subroutine mDecisions(numSteps,err,message)
     case('timeDlay'); model_decisions(iLookDECISIONS%subRouting)%iDecision = timeDelay           ! time-delay histogram
     case('qInstant'); model_decisions(iLookDECISIONS%subRouting)%iDecision = qInstant            ! instantaneous routing
     case default
-    err=10; message=trim(message)//"unknown option for sub-grid routing [option="//trim(model_decisions(iLookDECISIONS%subRouting)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for sub-grid routing [option="//trim(model_decisions(iLookDECISIONS%subRouting)%cDecision)//"]"; return
   end select
 
   ! choice of new snow density
@@ -664,7 +623,7 @@ subroutine mDecisions(numSteps,err,message)
     case('pahaut_76');                   model_decisions(iLookDECISIONS%snowDenNew)%iDecision = pahaut_76           ! Pahaut 1976, wind speed dependent (derived from Col de Porte, French Alps)
     case('constDens');                   model_decisions(iLookDECISIONS%snowDenNew)%iDecision = constDens           ! Constant new snow density
     case default
-    err=10; message=trim(message)//"unknown option for new snow density [option="//trim(model_decisions(iLookDECISIONS%snowDenNew)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for new snow density [option="//trim(model_decisions(iLookDECISIONS%snowDenNew)%cDecision)//"]"; return
   end select
 
   ! choice of snow unloading from canopy
@@ -672,49 +631,52 @@ subroutine mDecisions(numSteps,err,message)
     case('meltDripUnload','notPopulatedYet'); model_decisions(iLookDECISIONS%snowUnload)%iDecision = meltDripUnload  ! Hedstrom and Pomeroy (1998), Storck et al 2002 (snowUnloadingCoeff & ratioDrip2Unloading)
     case('windUnload');                       model_decisions(iLookDECISIONS%snowUnload)%iDecision = windUnload          ! Roesch et al 2001, formulate unloading based on wind and temperature
     case default
-    err=10; message=trim(message)//"unknown option for snow unloading [option="//trim(model_decisions(iLookDECISIONS%snowUnload)%cDecision)//"]"; return
+      err=10; message=trim(message)//"unknown option for snow unloading [option="//trim(model_decisions(iLookDECISIONS%snowUnload)%cDecision)//"]"; return
   end select
 
 
+  ! -----------------------------------------------------------------------------------------------------------------------------------------------
+  ! check for consistency among options
+  ! -----------------------------------------------------------------------------------------------------------------------------------------------
   ! check zero flux lower boundary for topmodel baseflow option
   select case(model_decisions(iLookDECISIONS%groundwatr)%iDecision)
     case(qbaseTopmodel)
-    if(model_decisions(iLookDECISIONS%bcLowrSoiH)%iDecision /= zeroFlux)then
-      message=trim(message)//'lower boundary condition for soil hydology must be zeroFlux with qbaseTopmodel option for groundwater'
-      err=20; return
-    end if
+      if(model_decisions(iLookDECISIONS%bcLowrSoiH)%iDecision /= zeroFlux)then
+        message=trim(message)//'lower boundary condition for soil hydology must be zeroFlux with qbaseTopmodel option for groundwater'
+        err=20; return
+      end if
   end select
 
   ! check power-law profile is selected when using topmodel baseflow option
   select case(model_decisions(iLookDECISIONS%groundwatr)%iDecision)
     case(qbaseTopmodel)
-    if(model_decisions(iLookDECISIONS%hc_profile)%iDecision /= powerLaw_profile)then
-      message=trim(message)//'power-law transmissivity profile must be selected when using topmodel baseflow option'
-      err=20; return
-    end if
+      if(model_decisions(iLookDECISIONS%hc_profile)%iDecision /= powerLaw_profile)then
+        message=trim(message)//'power-law transmissivity profile must be selected when using topmodel baseflow option'
+        err=20; return
+      end if
   end select
 
   ! check bigBucket groundwater option is used when for spatial groundwater is singleBasin
   if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == singleBasin)then
     if(model_decisions(iLookDECISIONS%groundwatr)%iDecision /= bigBucket)then
-    message=trim(message)//'groundwater parameterization must be bigBucket when using singleBasin for spatial_gw'
-    err=20; return
+      message=trim(message)//'groundwater parameterization must be bigBucket when using singleBasin for spatial_gw'
+      err=20; return
     end if
   end if
 
 end subroutine mDecisions
 
 
- ! ************************************************************************************************
- ! private subroutine readoption: read information from model decisions file
- ! ************************************************************************************************
+! ************************************************************************************************
+! private subroutine readoption: read information from model decisions file
+! ************************************************************************************************
 subroutine readoption(err,message)
   ! used to read information from model decisions file
   USE ascii_util_module,only:file_open       ! open file
   USE ascii_util_module,only:linewidth       ! max character number for one line
   USE ascii_util_module,only:get_vlines      ! get a vector of non-comment lines
-  USE summaActors_FileManager,only:SETTINGS_PATH    ! path for metadata files
-  USE summaActors_FileManager,only:M_DECISIONS      ! definition of modeling options
+  USE summaFileManager,only:SETTINGS_PATH    ! path for metadata files
+  USE summaFileManager,only:M_DECISIONS      ! definition of modeling options
   USE get_ixname_module,only:get_ixdecisions ! identify index of named variable
   USE globalData,only:model_decisions        ! model decision structure
   implicit none
@@ -735,7 +697,7 @@ subroutine readoption(err,message)
   err=0; message='readoption/'
   ! build filename
   infile = trim(SETTINGS_PATH)//trim(M_DECISIONS)
-  ! write(*,'(2(a,1x))') 'decisions file = ', trim(infile)
+  write(*,'(2(a,1x))') 'decisions file = ', trim(infile)
   ! open file
   call file_open(trim(infile),unt,err,cmessage)
   if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
@@ -753,13 +715,11 @@ subroutine readoption(err,message)
     if (err/=0) then; err=30; message=trim(message)//"errorReadLine"; return; end if
     ! get the index of the decision in the data structure
     iVar = get_ixdecisions(trim(option))
-    ! write(*,'(i4,1x,a)') iDecision, trim(option)//': '//trim(decision)
+    write(*,'(i4,1x,a)') iDecision, trim(option)//': '//trim(decision)
     if(iVar<=0)then; err=40; message=trim(message)//"cannotFindDecisionIndex[name='"//trim(option)//"']"; return; end if
     ! populate the model decisions structure
     model_decisions(iVar)%cOption   = trim(option)
     model_decisions(iVar)%cDecision = trim(decision)
   end do
 end subroutine readoption
-
-
 end module mDecisions_module