diff --git a/build/source/dshare/get_ixname.f90 b/build/source/dshare/get_ixname.f90
index a7866a4b1d7f1fa3bcc96be8b038e689fc6f5060..e3b1249d4866404c5161a348ef18c9d2617ea6b4 100755
--- a/build/source/dshare/get_ixname.f90
+++ b/build/source/dshare/get_ixname.f90
@@ -95,7 +95,9 @@ contains
case('subRouting' ); get_ixdecisions=iLookDECISIONS%subRouting ! choice of method for sub-grid routing
case('snowDenNew' ); get_ixdecisions=iLookDECISIONS%snowDenNew ! choice of method for new snow density
case('snowUnload' ); get_ixdecisions=iLookDECISIONS%snowUnload ! choice of parameterization for snow unloading from canopy
- ! get to here if cannot find the variable
+ case('howHeatCap' ); get_ixdecisions=iLookDECISIONS%howHeatCap ! how to compute heat capacity in energy equation
+ case('diffEqSolv' ); get_ixdecisions=iLookDECISIONS%diffEqSolv ! how to solve the system of differential equations
+ ! get to here if cannot find the variable
case default
get_ixdecisions = integerMissing
end select
diff --git a/build/source/dshare/var_lookup.f90 b/build/source/dshare/var_lookup.f90
index f1657613547a75f7a31b354caf2d9d743be0b063..7e9e48e376a61c74f09ec60529871e23523c0773 100755
--- a/build/source/dshare/var_lookup.f90
+++ b/build/source/dshare/var_lookup.f90
@@ -29,7 +29,6 @@ MODULE var_lookup
integer(8),parameter :: ix8Val=2 ! an example 8 byte integer
integer(i4b),parameter :: iLength =storage_size(ixVal) ! size of the example 4 byte integer
integer(i4b),parameter :: i8Length=storage_size(ix8Val) ! size of the example 8 byte integer
-
! ***************************************************************************************
! (0) define model decisions
! ***************************************************************************************
@@ -72,6 +71,8 @@ MODULE var_lookup
integer(i4b) :: spatial_gw = integerMissing ! choice of method for spatial representation of groundwater
integer(i4b) :: subRouting = integerMissing ! choice of method for sub-grid routing
integer(i4b) :: snowDenNew = integerMissing ! choice of method for new snow density
+ integer(i4b) :: howHeatCap = integerMissing ! how to compute heat capacity in energy equation
+ integer(i4b) :: diffEqSolv = integerMissing ! how to solve the system of differential equations
endtype iLook_decision
! ***********************************************************************************************************
@@ -777,7 +778,7 @@ MODULE var_lookup
type(iLook_decision),public,parameter :: iLookDECISIONS=iLook_decision( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,&
11, 12, 13, 14, 15, 16, 17, 18, 19, 20,&
21, 22, 23, 24, 25, 26, 27, 28, 29, 30,&
- 31, 32, 33, 34, 35, 36, 37, 38)
+ 31, 32, 33, 34, 35, 36, 37, 38, 39, 40)
! named variables: model time
type(iLook_time), public,parameter :: iLookTIME =iLook_time ( 1, 2, 3, 4, 5, 6, 7)
diff --git a/build/source/engine/mDecisions.f90 b/build/source/engine/mDecisions.f90
index f04bd566f863bc2fba068e6566a5266339a28678..4a863a7dab36fc0615446b574a4a2f73dee87277 100755
--- a/build/source/engine/mDecisions.f90
+++ b/build/source/engine/mDecisions.f90
@@ -58,6 +58,8 @@ integer(i4b),parameter,public :: minFunc = 62 ! do not enable c
integer(i4b),parameter,public :: constantScaling = 71 ! constant scaling factor
integer(i4b),parameter,public :: laiScaling = 72 ! exponential function of LAI (Leuning, Plant Cell Env 1995: "Scaling from..." [eq 9])
! 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
@@ -148,6 +150,12 @@ integer(i4b),parameter,public :: pahaut_76 = 314 ! Pahaut 1976, wi
! look-up values for the choice of snow unloading from the canopy
integer(i4b),parameter,public :: meltDripUnload = 321 ! Hedstrom and Pomeroy (1998), Storck et al 2002 (snowUnloadingCoeff & ratioDrip2Unloading)
integer(i4b),parameter,public :: windUnload = 322 ! Roesch et al 2001, formulate unloading based on wind and temperature
+! 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
@@ -399,8 +407,16 @@ subroutine mDecisions(numSteps,err,message)
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
@@ -408,6 +424,32 @@ subroutine mDecisions(numSteps,err,message)
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
+ 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
+ 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
@@ -634,31 +676,6 @@ subroutine mDecisions(numSteps,err,message)
end select
- ! -----------------------------------------------------------------------------------------------------------------------------------------------
- ! check for consistency among options
- ! -----------------------------------------------------------------------------------------------------------------------------------------------
-
- ! check there is prescribedHead for soil hydrology when zeroFlux or prescribedTemp for thermodynamics
- !select case(model_decisions(iLookDECISIONS%bcUpprTdyn)%iDecision)
- ! case(prescribedTemp,zeroFlux)
- ! if(model_decisions(iLookDECISIONS%bcUpprSoiH)%iDecision /= prescribedHead)then
- ! message=trim(message)//'upper boundary condition for soil hydology must be presHead with presTemp and zeroFlux options for thermodynamics'
- ! err=20; return
- ! end if
- !end select
-
- ! check there is prescribedTemp or zeroFlux for thermodynamics when using prescribedHead for soil hydrology
- !select case(model_decisions(iLookDECISIONS%bcUpprSoiH)%iDecision)
- ! case(prescribedHead)
- ! ! check that upper boundary condition for thermodynamics is presTemp or zeroFlux
- ! select case(model_decisions(iLookDECISIONS%bcUpprTdyn)%iDecision)
- ! case(prescribedTemp,zeroFlux) ! do nothing: this is OK
- ! case default
- ! message=trim(message)//'upper boundary condition for thermodynamics must be presTemp or zeroFlux with presHead option for soil hydology'
- ! err=20; return
- ! end select
- !end select
-
! check zero flux lower boundary for topmodel baseflow option
select case(model_decisions(iLookDECISIONS%groundwatr)%iDecision)
case(qbaseTopmodel)
@@ -685,13 +702,6 @@ subroutine mDecisions(numSteps,err,message)
end if
end if
- ! ensure that the LAI seaonality option is switched off (this was a silly idea, in retrospect)
- !if(model_decisions(iLookDECISIONS%LAI_method)%iDecision == specified)then
- ! message=trim(message)//'parameterization of LAI in terms of seasonal cycle of green veg fraction was a silly idea '&
- ! //' -- the LAI_method option ["specified"] is no longer supported'
- ! err=20; return
- !end if
-
end subroutine mDecisions
diff --git a/build/source/engine/opSplittin.f90 b/build/source/engine/opSplittin.f90
index e329f626c7a29940f158837d23952a5658f11f6b..f499e8c890028f36facb53191a82358fdbea6a15 100755
--- a/build/source/engine/opSplittin.f90
+++ b/build/source/engine/opSplittin.f90
@@ -112,7 +112,9 @@ USE mDecisions_module,only: &
USE mDecisions_module,only: &
qbaseTopmodel, & ! TOPMODEL-ish baseflow parameterization
bigBucket, & ! a big bucket (lumped aquifer model)
- noExplicit ! no explicit groundwater parameterization
+ noExplicit, & ! no explicit groundwater parameterization
+ sundialIDA, & ! IDA solver from Sundials package
+ backwEuler ! backward Euler method
! safety: set private unless specified otherwise
implicit none
@@ -122,7 +124,6 @@ public::opSplittin
! named variables for the coupling method
integer(i4b),parameter :: fullyCoupled=1 ! 1st try: fully coupled solution
integer(i4b),parameter :: stateTypeSplit=2 ! 2nd try: separate solutions for each state type
-integer(i4b),parameter :: nCoupling=2 ! number of possible solutions
! named variables for the state variable split
integer(i4b),parameter :: nrgSplit=1 ! order in sequence for the energy operation
@@ -156,16 +157,16 @@ real(dp),parameter :: dx = 1.e-8_dp ! finite difference increm
contains
- ! **********************************************************************************************************
- ! public subroutine opSplittin: run the coupled energy-mass model for one timestep
- !
- ! The logic of the solver is as follows:
- ! (1) Attempt different solutions in the following order: (a) fully coupled; (b) split by state type and by
- ! domain type for a given energy and mass split (vegetation, snow, and soil); and (c) scalar solution
- ! for a given state type and domain subset.
- ! (2) For a given split, compute a variable number of substeps (in varSubstep).
- ! **********************************************************************************************************
- subroutine opSplittin(&
+! **********************************************************************************************************
+! public subroutine opSplittin: run the coupled energy-mass model for one timestep
+!
+! The logic of the solver is as follows:
+! (1) Attempt different solutions in the following order: (a) fully coupled; (b) split by state type and by
+! domain type for a given energy and mass split (vegetation, snow, and soil); and (c) scalar solution
+! for a given state type and domain subset.
+! (2) For a given split, compute a variable number of substeps (in varSubstep).
+! **********************************************************************************************************
+subroutine opSplittin(&
! input: model control
nSnow, & ! intent(in): number of snow layers
nSoil, & ! intent(in): number of soil layers
@@ -192,173 +193,188 @@ contains
stepFailure, & ! intent(out): flag to denote step failure
ixCoupling, & ! intent(out): coupling method used in this iteration
err,message) ! intent(out): error code and error message
- ! ---------------------------------------------------------------------------------------
- ! structure allocations
- USE allocspace4chm_module,only:allocLocal ! allocate local data structures
- ! simulation of fluxes and residuals given a trial state vector
- USE soil_utils_module,only:matricHead ! compute the matric head based on volumetric water content
- USE soil_utils_module,only:liquidHead ! compute the liquid water matric potential
- ! population/extraction of state vectors
- USE indexState_module,only:indexSplit ! get state indices
- USE varSubstep_module,only:varSubstep ! complete substeps for a given split
- ! identify name of variable type (for error message)
- USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages
- implicit none
- ! ---------------------------------------------------------------------------------------
- ! * dummy variables
- ! ---------------------------------------------------------------------------------------
- ! input: model control
- integer(i4b),intent(in) :: nSnow ! number of snow layers
- integer(i4b),intent(in) :: nSoil ! number of soil layers
- integer(i4b),intent(in) :: nLayers ! total number of layers
- integer(i4b),intent(in) :: nState ! total number of state variables
- real(dp),intent(inout) :: dt ! time step (seconds)
- logical(lgt),intent(in) :: firstSubStep ! flag to indicate if we are processing the first sub-step
- logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
- ! input/output: data structures
- type(var_i),intent(in) :: type_data ! type of vegetation and soil
- type(var_d),intent(in) :: attr_data ! spatial attributes
- type(var_d),intent(in) :: forc_data ! model forcing data
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- type(var_ilength),intent(inout) :: indx_data ! indices for a local HRU
- type(var_dlength),intent(inout) :: prog_data ! prognostic variables for a local HRU
- type(var_dlength),intent(inout) :: diag_data ! diagnostic variables for a local HRU
- type(var_dlength),intent(inout) :: flux_data ! model fluxes for a local HRU
- type(var_dlength),intent(in) :: bvar_data ! model variables for the local basin
- type(zLookup),intent(in) :: lookup_data ! lookup tables
- type(model_options),intent(in) :: model_decisions(:) ! model decisions
- ! output: model control
- real(dp),intent(out) :: dtMultiplier ! substep multiplier (-)
- logical(lgt),intent(out) :: tooMuchMelt ! flag to denote that ice is insufficient to support melt
- logical(lgt),intent(out) :: stepFailure ! flag to denote step failure
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! *********************************************************************************************************************************************************
- ! *********************************************************************************************************************************************************
- ! ---------------------------------------------------------------------------------------
- ! * general local variables
- ! ---------------------------------------------------------------------------------------
- character(LEN=256) :: cmessage ! error message of downwind routine
- integer(i4b) :: minLayer ! the minimum layer used in assigning flags for flux aggregations
- integer(i4b) :: iOffset ! offset to account for different indices in the soil domain
- integer(i4b) :: iMin(1),iMax(1) ! bounds of a given vector
- integer(i4b) :: iLayer,jLayer ! index of model layer
- integer(i4b) :: iSoil ! index of soil layer
- integer(i4b) :: iVar ! index of variables in data structures
- logical(lgt) :: firstSuccess ! flag to define the first success
- logical(lgt) :: firstFluxCall ! flag to define the first flux call
- logical(lgt) :: reduceCoupledStep ! flag to define the need to reduce the length of the coupled step
- type(var_dlength) :: prog_temp ! temporary model prognostic variables
- type(var_dlength) :: diag_temp ! temporary model diagnostic variables
- type(var_dlength) :: flux_temp ! temporary model fluxes
- type(var_dlength) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
- real(dp),dimension(nLayers) :: mLayerVolFracIceInit ! initial vector for volumetric fraction of ice (-)
- ! ------------------------------------------------------------------------------------------------------
- ! * operator splitting
- ! ------------------------------------------------------------------------------------------------------
- ! minimum timestep
- real(dp),parameter :: dtmin_coupled=1800._dp ! minimum time step for the fully coupled solution (seconds)
- real(dp),parameter :: dtmin_split=60._dp ! minimum time step for the fully split solution (seconds)
- real(dp),parameter :: dtmin_scalar=10._dp ! minimum time step for the scalar solution (seconds)
- real(dp) :: dt_min ! minimum time step (seconds)
- real(dp) :: dtInit ! initial time step (seconds)
- ! explicit error tolerance (depends on state type split, so defined here)
- real(dp),parameter :: errorTolLiqFlux=0.01_dp ! error tolerance in the explicit solution (liquid flux)
- real(dp),parameter :: errorTolNrgFlux=10._dp ! error tolerance in the explicit solution (energy flux)
- ! number of substeps taken for a given split
- integer(i4b) :: nSubsteps ! number of substeps taken for a given split
- ! named variables defining the coupling and solution method
- integer(i4b) :: ixCoupling ! index of coupling method (1,2)
- integer(i4b) :: ixSolution ! index of solution method (1,2)
- integer(i4b) :: ixStateThenDomain ! switch between the state and domain (1,2)
- integer(i4b) :: tryDomainSplit ! (0,1) - flag to try the domain split
- ! actual number of splits
- integer(i4b) :: nStateTypeSplit ! number of splits for the state type
- integer(i4b) :: nDomainSplit ! number of splits for the domain
- integer(i4b) :: nStateSplit ! number of splits for the states within a given domain
- ! indices for the state type and the domain split
- integer(i4b) :: iStateTypeSplit ! index of the state type split
- integer(i4b) :: iDomainSplit ! index of the domain split
- integer(i4b) :: iStateSplit ! index of the state split
- ! flux masks
- logical(lgt) :: neededFlux(nFlux) ! .true. if flux is needed at all
- logical(lgt) :: desiredFlux ! .true. if flux is desired for a given split
- type(var_ilength) :: fluxCount ! number of times each flux is updated (should equal nSubsteps)
- type(var_flagVec) :: fluxMask ! mask defining model fluxes
- ! state masks
- integer(i4b),dimension(nState) :: stateCheck ! number of times each state variable is updated (should equal 1)
- logical(lgt),dimension(nState) :: stateMask ! mask defining desired state variables
- integer(i4b) :: nSubset ! number of selected state variables for a given split
- ! flags
- logical(lgt) :: failure ! flag to denote failure of substepping
- logical(lgt) :: doAdjustTemp ! flag to adjust temperature after the mass split
- logical(lgt) :: failedMinimumStep ! flag to denote failure of substepping for a given split
- integer(i4b) :: ixSaturation ! index of the lowest saturated layer (NOTE: only computed on the first iteration)
- ! ---------------------------------------------------------------------------------------
- ! point to variables in the data structures
- ! ---------------------------------------------------------------------------------------
- globalVars: associate(&
- ! model decisions
- ixGroundwater => model_decisions(iLookDECISIONS%groundwatr)%iDecision ,& ! intent(in): [i4b] groundwater parameterization
- ixSpatialGroundwater => model_decisions(iLookDECISIONS%spatial_gw)%iDecision ,& ! intent(in): [i4b] spatial representation of groundwater (local-column or single-basin)
- ! domain boundary conditions
- airtemp => forc_data%var(iLookFORCE%airtemp) ,& ! intent(in): [dp] temperature of the upper boundary of the snow and soil domains (K)
- ! vector of energy and hydrology indices for the snow and soil domains
- ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow+soil domain
- ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow+soil domain
- nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain
- nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology state variables in the snow+soil domain
- ! indices of model state variables
- ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset (missing for values not in the subset)
- ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (ixNrgState...)
- ixNrgCanair => indx_data%var(iLookINDEX%ixNrgCanair)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in canopy air space domain
- ixNrgCanopy => indx_data%var(iLookINDEX%ixNrgCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the canopy domain
- ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the canopy domain
- ixNrgLayer => indx_data%var(iLookINDEX%ixNrgLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the snow+soil domain
- ixHydLayer => indx_data%var(iLookINDEX%ixHydLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the snow+soil domain
- ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable
- ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable
- ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass)
- ! numerix tracking
- numberStateSplit => indx_data%var(iLookINDEX%numberStateSplit )%dat(1) ,& ! intent(inout): [i4b] number of state splitting solutions (-)
- numberDomainSplitNrg => indx_data%var(iLookINDEX%numberDomainSplitNrg )%dat(1) ,& ! intent(inout): [i4b] number of domain splitting solutions for energy (-)
- numberDomainSplitMass => indx_data%var(iLookINDEX%numberDomainSplitMass)%dat(1) ,& ! intent(inout): [i4b] number of domain splitting solutions for mass (-)
- numberScalarSolutions => indx_data%var(iLookINDEX%numberScalarSolutions)%dat(1) ,& ! intent(inout): [i4b] number of scalar solutions (-)
- ! domain configuration
- canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp] canopy depth (m)
- mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
- ! snow parameters
- snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) ,& ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1)
- ! depth-varying soil parameters
- vGn_m => diag_data%var(iLookDIAG%scalarVGn_m)%dat ,& ! intent(in): [dp(:)] van Genutchen "m" parameter (-)
- vGn_n => mpar_data%var(iLookPARAM%vGn_n)%dat ,& ! intent(in): [dp(:)] van Genutchen "n" parameter (-)
- vGn_alpha => mpar_data%var(iLookPARAM%vGn_alpha)%dat ,& ! intent(in): [dp(:)] van Genutchen "alpha" parameter (m-1)
- theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat ,& ! intent(in): [dp(:)] soil porosity (-)
- theta_res => mpar_data%var(iLookPARAM%theta_res)%dat ,& ! intent(in): [dp(:)] soil residual volumetric water content (-)
- ! soil parameters
- specificStorage => mpar_data%var(iLookPARAM%specificStorage)%dat(1) ,& ! intent(in): [dp] specific storage coefficient (m-1)
- ! model diagnostic variables (fraction of liquid water)
- scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg)%dat(1) ,& ! intent(out): [dp] fraction of liquid water on vegetation (-)
- mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow)%dat ,& ! intent(out): [dp(:)] fraction of liquid water in each snow layer (-)
- mLayerMeltFreeze => diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat ,& ! intent(out): [dp(:)] melt-freeze in each snow and soil layer (kg m-3)
- ! model state variables (vegetation canopy)
- scalarCanairTemp => prog_data%var(iLookPROG%scalarCanairTemp)%dat(1) ,& ! intent(out): [dp] temperature of the canopy air space (K)
- scalarCanopyTemp => prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1) ,& ! intent(out): [dp] temperature of the vegetation canopy (K)
- scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! intent(out): [dp] mass of ice on the vegetation canopy (kg m-2)
- scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! intent(out): [dp] mass of liquid water on the vegetation canopy (kg m-2)
- scalarCanopyWat => prog_data%var(iLookPROG%scalarCanopyWat)%dat(1) ,& ! intent(out): [dp] mass of total water on the vegetation canopy (kg m-2)
- ! model state variables (snow and soil domains)
- mLayerTemp => prog_data%var(iLookPROG%mLayerTemp)%dat ,& ! intent(out): [dp(:)] temperature of each snow/soil layer (K)
- mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(out): [dp(:)] volumetric fraction of ice (-)
- mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat ,& ! intent(out): [dp(:)] volumetric fraction of liquid water (-)
- mLayerVolFracWat => prog_data%var(iLookPROG%mLayerVolFracWat)%dat ,& ! intent(out): [dp(:)] volumetric fraction of total water (-)
- mLayerMatricHead => prog_data%var(iLookPROG%mLayerMatricHead)%dat ,& ! intent(out): [dp(:)] matric head (m)
- mLayerMatricHeadLiq => diag_data%var(iLookDIAG%mLayerMatricHeadLiq)%dat & ! intent(out): [dp(:)] matric potential of liquid water (m)
- )
- ! ---------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message="opSplittin/"
-
+ ! ---------------------------------------------------------------------------------------
+ ! structure allocations
+ USE allocspace4chm_module,only:allocLocal ! allocate local data structures
+ ! simulation of fluxes and residuals given a trial state vector
+ USE soil_utils_module,only:matricHead ! compute the matric head based on volumetric water content
+ USE soil_utils_module,only:liquidHead ! compute the liquid water matric potential
+ ! population/extraction of state vectors
+ USE indexState_module,only:indexSplit ! get state indices
+ USE varSubstep_module,only:varSubstep ! complete substeps for a given split
+ ! identify name of variable type (for error message)
+ USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages
+ implicit none
+ ! ---------------------------------------------------------------------------------------
+ ! * dummy variables
+ ! ---------------------------------------------------------------------------------------
+ ! input: model control
+ integer(i4b),intent(in) :: nSnow ! number of snow layers
+ integer(i4b),intent(in) :: nSoil ! number of soil layers
+ integer(i4b),intent(in) :: nLayers ! total number of layers
+ integer(i4b),intent(in) :: nState ! total number of state variables
+ real(dp),intent(inout) :: dt ! time step (seconds)
+ logical(lgt),intent(in) :: firstSubStep ! flag to indicate if we are processing the first sub-step
+ logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
+ ! input/output: data structures
+ type(var_i),intent(in) :: type_data ! type of vegetation and soil
+ type(var_d),intent(in) :: attr_data ! spatial attributes
+ type(var_d),intent(in) :: forc_data ! model forcing data
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ type(var_ilength),intent(inout) :: indx_data ! indices for a local HRU
+ type(var_dlength),intent(inout) :: prog_data ! prognostic variables for a local HRU
+ type(var_dlength),intent(inout) :: diag_data ! diagnostic variables for a local HRU
+ type(var_dlength),intent(inout) :: flux_data ! model fluxes for a local HRU
+ type(var_dlength),intent(in) :: bvar_data ! model variables for the local basin
+ type(zLookup),intent(in) :: lookup_data ! lookup tables
+ type(model_options),intent(in) :: model_decisions(:) ! model decisions
+ ! output: model control
+ real(dp),intent(out) :: dtMultiplier ! substep multiplier (-)
+ logical(lgt),intent(out) :: tooMuchMelt ! flag to denote that ice is insufficient to support melt
+ logical(lgt),intent(out) :: stepFailure ! flag to denote step failure
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! *********************************************************************************************************************************************************
+ ! *********************************************************************************************************************************************************
+ ! ---------------------------------------------------------------------------------------
+ ! * general local variables
+ ! ---------------------------------------------------------------------------------------
+ character(LEN=256) :: cmessage ! error message of downwind routine
+ integer(i4b) :: minLayer ! the minimum layer used in assigning flags for flux aggregations
+ integer(i4b) :: iOffset ! offset to account for different indices in the soil domain
+ integer(i4b) :: iMin(1),iMax(1) ! bounds of a given vector
+ integer(i4b) :: iLayer,jLayer ! index of model layer
+ integer(i4b) :: iSoil ! index of soil layer
+ integer(i4b) :: iVar ! index of variables in data structures
+ logical(lgt) :: firstSuccess ! flag to define the first success
+ logical(lgt) :: firstFluxCall ! flag to define the first flux call
+ logical(lgt) :: reduceCoupledStep ! flag to define the need to reduce the length of the coupled step
+ type(var_dlength) :: prog_temp ! temporary model prognostic variables
+ type(var_dlength) :: diag_temp ! temporary model diagnostic variables
+ type(var_dlength) :: flux_temp ! temporary model fluxes
+ type(var_dlength) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
+ real(dp),dimension(nLayers) :: mLayerVolFracIceInit ! initial vector for volumetric fraction of ice (-)
+ ! ------------------------------------------------------------------------------------------------------
+ ! * operator splitting
+ ! ------------------------------------------------------------------------------------------------------
+ ! minimum timestep
+ real(dp),parameter :: dtmin_coupled=1800._dp ! minimum time step for the fully coupled solution (seconds)
+ real(dp),parameter :: dtmin_split=60._dp ! minimum time step for the fully split solution (seconds)
+ real(dp),parameter :: dtmin_scalar=10._dp ! minimum time step for the scalar solution (seconds)
+ real(dp) :: dt_min ! minimum time step (seconds)
+ real(dp) :: dtInit ! initial time step (seconds)
+ ! explicit error tolerance (depends on state type split, so defined here)
+ real(dp),parameter :: errorTolLiqFlux=0.01_dp ! error tolerance in the explicit solution (liquid flux)
+ real(dp),parameter :: errorTolNrgFlux=10._dp ! error tolerance in the explicit solution (energy flux)
+ ! number of substeps taken for a given split
+ integer(i4b) :: nSubsteps ! number of substeps taken for a given split
+ ! named variables defining the coupling and solution method
+ integer(i4b) :: ixCoupling ! index of coupling method (1,2)
+ integer(i4b) :: ixSolution ! index of solution method (1,2)
+ integer(i4b) :: ixStateThenDomain ! switch between the state and domain (1,2)
+ integer(i4b) :: tryDomainSplit ! (0,1) - flag to try the domain split
+ ! actual number of splits
+ integer(i4b) :: nStateTypeSplit ! number of splits for the state type
+ integer(i4b) :: nDomainSplit ! number of splits for the domain
+ integer(i4b) :: nStateSplit ! number of splits for the states within a given domain
+ ! indices for the state type and the domain split
+ integer(i4b) :: iStateTypeSplit ! index of the state type split
+ integer(i4b) :: iDomainSplit ! index of the domain split
+ integer(i4b) :: iStateSplit ! index of the state split
+ ! flux masks
+ logical(lgt) :: neededFlux(nFlux) ! .true. if flux is needed at all
+ logical(lgt) :: desiredFlux ! .true. if flux is desired for a given split
+ type(var_ilength) :: fluxCount ! number of times each flux is updated (should equal nSubsteps)
+ type(var_flagVec) :: fluxMask ! mask defining model fluxes
+ ! state masks
+ integer(i4b),dimension(nState) :: stateCheck ! number of times each state variable is updated (should equal 1)
+ logical(lgt),dimension(nState) :: stateMask ! mask defining desired state variables
+ integer(i4b) :: nSubset ! number of selected state variables for a given split
+ ! flags
+ logical(lgt) :: failure ! flag to denote failure of substepping
+ logical(lgt) :: doAdjustTemp ! flag to adjust temperature after the mass split
+ logical(lgt) :: failedMinimumStep ! flag to denote failure of substepping for a given split
+ integer(i4b) :: ixSaturation ! index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ integer(i4b) :: nCoupling ! number of possible solutions
+
+ ! ---------------------------------------------------------------------------------------
+ ! point to variables in the data structures
+ ! ---------------------------------------------------------------------------------------
+ globalVars: associate(&
+ ! model decisions
+ ixGroundwater => model_decisions(iLookDECISIONS%groundwatr)%iDecision ,& ! intent(in): [i4b] groundwater parameterization
+ ixSpatialGroundwater => model_decisions(iLookDECISIONS%spatial_gw)%iDecision ,& ! intent(in): [i4b] spatial representation of groundwater (local-column or single-basin)
+ ! domain boundary conditions
+ airtemp => forc_data%var(iLookFORCE%airtemp) ,& ! intent(in): [dp] temperature of the upper boundary of the snow and soil domains (K)
+ ! vector of energy and hydrology indices for the snow and soil domains
+ ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow+soil domain
+ ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow+soil domain
+ nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain
+ nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology state variables in the snow+soil domain
+ ! indices of model state variables
+ ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset (missing for values not in the subset)
+ ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (ixNrgState...)
+ ixNrgCanair => indx_data%var(iLookINDEX%ixNrgCanair)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in canopy air space domain
+ ixNrgCanopy => indx_data%var(iLookINDEX%ixNrgCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the canopy domain
+ ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the canopy domain
+ ixNrgLayer => indx_data%var(iLookINDEX%ixNrgLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the snow+soil domain
+ ixHydLayer => indx_data%var(iLookINDEX%ixHydLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the snow+soil domain
+ ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable
+ ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable
+ ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass)
+ ! numerix tracking
+ numberStateSplit => indx_data%var(iLookINDEX%numberStateSplit )%dat(1) ,& ! intent(inout): [i4b] number of state splitting solutions (-)
+ numberDomainSplitNrg => indx_data%var(iLookINDEX%numberDomainSplitNrg )%dat(1) ,& ! intent(inout): [i4b] number of domain splitting solutions for energy (-)
+ numberDomainSplitMass => indx_data%var(iLookINDEX%numberDomainSplitMass)%dat(1) ,& ! intent(inout): [i4b] number of domain splitting solutions for mass (-)
+ numberScalarSolutions => indx_data%var(iLookINDEX%numberScalarSolutions)%dat(1) ,& ! intent(inout): [i4b] number of scalar solutions (-)
+ ! domain configuration
+ canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp] canopy depth (m)
+ mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
+ ! snow parameters
+ snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) ,& ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1)
+ ! depth-varying soil parameters
+ vGn_m => diag_data%var(iLookDIAG%scalarVGn_m)%dat ,& ! intent(in): [dp(:)] van Genutchen "m" parameter (-)
+ vGn_n => mpar_data%var(iLookPARAM%vGn_n)%dat ,& ! intent(in): [dp(:)] van Genutchen "n" parameter (-)
+ vGn_alpha => mpar_data%var(iLookPARAM%vGn_alpha)%dat ,& ! intent(in): [dp(:)] van Genutchen "alpha" parameter (m-1)
+ theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat ,& ! intent(in): [dp(:)] soil porosity (-)
+ theta_res => mpar_data%var(iLookPARAM%theta_res)%dat ,& ! intent(in): [dp(:)] soil residual volumetric water content (-)
+ ! soil parameters
+ specificStorage => mpar_data%var(iLookPARAM%specificStorage)%dat(1) ,& ! intent(in): [dp] specific storage coefficient (m-1)
+ ! model diagnostic variables (fraction of liquid water)
+ scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg)%dat(1) ,& ! intent(out): [dp] fraction of liquid water on vegetation (-)
+ mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow)%dat ,& ! intent(out): [dp(:)] fraction of liquid water in each snow layer (-)
+ mLayerMeltFreeze => diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat ,& ! intent(out): [dp(:)] melt-freeze in each snow and soil layer (kg m-3)
+ ! model state variables (vegetation canopy)
+ scalarCanairTemp => prog_data%var(iLookPROG%scalarCanairTemp)%dat(1) ,& ! intent(out): [dp] temperature of the canopy air space (K)
+ scalarCanopyTemp => prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1) ,& ! intent(out): [dp] temperature of the vegetation canopy (K)
+ scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! intent(out): [dp] mass of ice on the vegetation canopy (kg m-2)
+ scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! intent(out): [dp] mass of liquid water on the vegetation canopy (kg m-2)
+ scalarCanopyWat => prog_data%var(iLookPROG%scalarCanopyWat)%dat(1) ,& ! intent(out): [dp] mass of total water on the vegetation canopy (kg m-2)
+ ! model state variables (snow and soil domains)
+ mLayerTemp => prog_data%var(iLookPROG%mLayerTemp)%dat ,& ! intent(out): [dp(:)] temperature of each snow/soil layer (K)
+ mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(out): [dp(:)] volumetric fraction of ice (-)
+ mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat ,& ! intent(out): [dp(:)] volumetric fraction of liquid water (-)
+ mLayerVolFracWat => prog_data%var(iLookPROG%mLayerVolFracWat)%dat ,& ! intent(out): [dp(:)] volumetric fraction of total water (-)
+ mLayerMatricHead => prog_data%var(iLookPROG%mLayerMatricHead)%dat ,& ! intent(out): [dp(:)] matric head (m)
+ mLayerMatricHeadLiq => diag_data%var(iLookDIAG%mLayerMatricHeadLiq)%dat & ! intent(out): [dp(:)] matric potential of liquid water (m)
+ )
+ ! ---------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message="opSplittin/"
+
+ ! we just solve the fully coupled problem by ida
+ select case(model_decisions(iLookDECISIONS%diffEqSolv)%iDecision)
+ case(sundialIDA); nCoupling = 1
+ case(backwEuler); nCoupling = 2
+ case default
+ err=20
+ message=trim(message)//'expect case to be sundialIDA or backwEuler'
+ print*, message
+ return
+ end select
+
+
+ print*, "nCoupling", nCoupling
! *****
! (0) PRELIMINARIES...
! ********************
diff --git a/utils/laugh_tests/celia1990/settings/summa_zDecisions_celia1990.txt b/utils/laugh_tests/celia1990/settings/summa_zDecisions_celia1990.txt
index 92822de7e65f51e81711c5c982a52914eee16fe8..78562a25a63e9d19961c8b7fd7d633dccc6abb9b 100644
--- a/utils/laugh_tests/celia1990/settings/summa_zDecisions_celia1990.txt
+++ b/utils/laugh_tests/celia1990/settings/summa_zDecisions_celia1990.txt
@@ -36,6 +36,8 @@ thCondSnow jrdn1991 ! (26) choice of thermal conduct
thCondSoil mixConstit ! (27) choice of thermal conductivity representation for soil
spatial_gw localColumn ! (28) choice of method for the spatial representation of groundwater
subRouting timeDlay ! (29) choice of method for sub-grid routing
+howHeatCap enthalpyFD
+diffEqSolv backwEuler
! ***********************************************************************************************
! ***** description of the options available -- nothing below this point is read ****************
! ***********************************************************************************************
diff --git a/utils/laugh_tests/colbeck1976/settings/summa_zDecisions_colbeck1976.txt b/utils/laugh_tests/colbeck1976/settings/summa_zDecisions_colbeck1976.txt
index f5ac1fa2017c484c54809ff8a92f921a4af24451..60f645f2c9f1463c3a4ce38bdda62d05f3405007 100644
--- a/utils/laugh_tests/colbeck1976/settings/summa_zDecisions_colbeck1976.txt
+++ b/utils/laugh_tests/colbeck1976/settings/summa_zDecisions_colbeck1976.txt
@@ -36,6 +36,8 @@ thCondSnow smnv2000 ! (26) choice of thermal conduct
thCondSoil mixConstit ! (27) choice of thermal conductivity representation for soil
spatial_gw localColumn ! (28) choice of method for the spatial representation of groundwater
subRouting timeDlay ! (29) choice of method for sub-grid routing
+howHeatCap enthalpyFD
+diffEqSolv backwEuler
! ***********************************************************************************************
! ***** description of the options available -- nothing below this point is read ****************
! ***********************************************************************************************
diff --git a/utils/laugh_tests/miller1998/settings/summa_zDecisions_miller1998.txt b/utils/laugh_tests/miller1998/settings/summa_zDecisions_miller1998.txt
index d17071ceff8b467be5120db09ec0df8c69d7a2d8..eb5fe85530a04a9c7c458e580185db6799c49916 100644
--- a/utils/laugh_tests/miller1998/settings/summa_zDecisions_miller1998.txt
+++ b/utils/laugh_tests/miller1998/settings/summa_zDecisions_miller1998.txt
@@ -36,6 +36,8 @@ thCondSnow jrdn1991 ! (26) choice of thermal conducti
thCondSoil mixConstit ! (27) choice of thermal conductivity representation for soil
spatial_gw localColumn ! (28) choice of method for the spatial representation of groundwater
subRouting timeDlay ! (29) choice of method for sub-grid routing
+howHeatCap enthalpyFD
+diffEqSolv backwEuler
! ***********************************************************************************************
! ***** description of the options available -- nothing below this point is read ****************
! ***********************************************************************************************
diff --git a/utils/laugh_tests/mizoguchi1990/settings/summa_zDecisions_mizoguchi.txt b/utils/laugh_tests/mizoguchi1990/settings/summa_zDecisions_mizoguchi.txt
index 32ebccac5a0c549881def40f2c71c35b005c1652..f320c97d24a2eddfea29333ed872fef2b74009b6 100644
--- a/utils/laugh_tests/mizoguchi1990/settings/summa_zDecisions_mizoguchi.txt
+++ b/utils/laugh_tests/mizoguchi1990/settings/summa_zDecisions_mizoguchi.txt
@@ -36,6 +36,8 @@ thCondSnow jrdn1991 ! (26) choice of thermal conduct
thCondSoil hanssonVZJ ! (27) choice of thermal conductivity representation for soil
spatial_gw localColumn ! (28) choice of method for the spatial representation of groundwater
subRouting timeDlay ! (29) choice of method for sub-grid routing
+howHeatCap enthalpyFD
+diffEqSolv backwEuler
! ***********************************************************************************************
! ***** description of the options available -- nothing below this point is read ****************
! ***********************************************************************************************
diff --git a/utils/laugh_tests/vanderborght2005/settings/summa_zDecisions_vanderborght2005.txt b/utils/laugh_tests/vanderborght2005/settings/summa_zDecisions_vanderborght2005.txt
index 6d50dc7f44b4b52a93dc23ff7a0e971910cb467e..0ba3b21c0f8a9007d0258cb50624577e5edefeae 100644
--- a/utils/laugh_tests/vanderborght2005/settings/summa_zDecisions_vanderborght2005.txt
+++ b/utils/laugh_tests/vanderborght2005/settings/summa_zDecisions_vanderborght2005.txt
@@ -36,6 +36,8 @@ thCondSnow jrdn1991 ! (26) choice of thermal conduct
thCondSoil mixConstit ! (27) choice of thermal conductivity representation for soil
spatial_gw localColumn ! (28) choice of method for the spatial representation of groundwater
subRouting timeDlay ! (29) choice of method for sub-grid routing
+howHeatCap enthalpyFD
+diffEqSolv backwEuler
! ***********************************************************************************************
! ***** description of the options available -- nothing below this point is read ****************
! ***********************************************************************************************
diff --git a/utils/laugh_tests/wigmosta1999/settings/summa_zDecisions.txt b/utils/laugh_tests/wigmosta1999/settings/summa_zDecisions.txt
index 3548441f737f8fe5c4036ae6092ea814099e63b8..74b289f5297f0f5c586c5762fb2981ef492a2008 100644
--- a/utils/laugh_tests/wigmosta1999/settings/summa_zDecisions.txt
+++ b/utils/laugh_tests/wigmosta1999/settings/summa_zDecisions.txt
@@ -36,6 +36,8 @@ thCondSnow jrdn1991 ! (26) choice of thermal conduct
thCondSoil mixConstit ! (27) choice of thermal conductivity representation for soil
spatial_gw localColumn ! (28) choice of method for the spatial representation of groundwater
subRouting timeDlay ! (29) choice of method for sub-grid routing
+howHeatCap enthalpyFD
+diffEqSolv backwEuler
! ***********************************************************************************************
! ***** description of the options available -- nothing below this point is read ****************
! ***********************************************************************************************