From b55a84c07c5ef8fc073cff1a72a44f1543054b7e Mon Sep 17 00:00:00 2001 From: Kyle <kyle.c.klenk@gmail.com> Date: Mon, 29 Aug 2022 20:37:58 +0000 Subject: [PATCH] added model decisions for choosing backweuler and ida --- build/source/dshare/get_ixname.f90 | 4 +- build/source/dshare/var_lookup.f90 | 5 +- build/source/engine/mDecisions.f90 | 74 ++-- build/source/engine/opSplittin.f90 | 374 +++++++++--------- .../settings/summa_zDecisions_celia1990.txt | 2 + .../settings/summa_zDecisions_colbeck1976.txt | 2 + .../settings/summa_zDecisions_miller1998.txt | 2 + .../settings/summa_zDecisions_mizoguchi.txt | 2 + .../summa_zDecisions_vanderborght2005.txt | 2 + .../settings/summa_zDecisions.txt | 2 + 10 files changed, 255 insertions(+), 214 deletions(-) diff --git a/build/source/dshare/get_ixname.f90 b/build/source/dshare/get_ixname.f90 index a7866a4..e3b1249 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 f165761..7e9e48e 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 f04bd56..4a863a7 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 e329f62..f499e8c 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 92822de..78562a2 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 f5ac1fa..60f645f 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 d17071c..eb5fe85 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 32ebcca..f320c97 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 6d50dc7..0ba3b21 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 3548441..74b289f 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 **************** ! *********************************************************************************************** -- GitLab