diff --git a/build/source/engine/sundials/updateVarsSundials.f90 b/build/source/engine/sundials/updateVarsSundials.f90
index 0f54d15b1ef86edb8253a05526a000e70610e498..3a33e4e4ce3688cc276edc2a8dd158ff6851784e 100644
--- a/build/source/engine/sundials/updateVarsSundials.f90
+++ b/build/source/engine/sundials/updateVarsSundials.f90
@@ -18,831 +18,1087 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.
-module updateVarsSundials_module
+module varSubstepSundials_module
! data types
USE nrtype
- ! missing values
+ ! access missing values
USE globalData,only:integerMissing ! missing integer
- USE globalData,only:realMissing ! missing real number
+ USE globalData,only:realMissing ! missing double precision number
+ USE globalData,only:quadMissing ! missing quadruple precision number
! access the global print flag
USE globalData,only:globalPrintFlag
! domain types
- USE globalData,only:iname_cas ! named variables for canopy air space
- USE globalData,only:iname_veg ! named variables for vegetation canopy
+ USE globalData,only:iname_cas ! named variables for the canopy air space
+ USE globalData,only:iname_veg ! named variables for vegetation
USE globalData,only:iname_snow ! named variables for snow
USE globalData,only:iname_soil ! named variables for soil
- USE globalData,only:iname_aquifer ! named variables for the aquifer
-
- ! named variables to describe the state variable type
- USE globalData,only:iname_nrgCanair ! named variable defining the energy of the canopy air space
- USE globalData,only:iname_nrgCanopy ! named variable defining the energy of the vegetation canopy
- USE globalData,only:iname_watCanopy ! named variable defining the mass of total water on the vegetation canopy
- USE globalData,only:iname_liqCanopy ! named variable defining the mass of liquid water on the vegetation canopy
- USE globalData,only:iname_nrgLayer ! named variable defining the energy state variable for snow+soil layers
- USE globalData,only:iname_watLayer ! named variable defining the total water state variable for snow+soil layers
- USE globalData,only:iname_liqLayer ! named variable defining the liquid water state variable for snow+soil layers
- USE globalData,only:iname_matLayer ! named variable defining the matric head state variable for soil layers
- USE globalData,only:iname_lmpLayer ! named variable defining the liquid matric potential state variable for soil layers
-
- ! metadata for information in the data structures
- USE globalData,only:indx_meta ! metadata for the variables in the index structure
- ! constants
- USE multiconst,only:&
- gravity, & ! acceleration of gravity (m s-2)
- Tfreeze, & ! temperature at freezing (K)
- Cp_air, & ! specific heat of air (J kg-1 K-1)
- Cp_ice, & ! specific heat of ice (J kg-1 K-1)
- Cp_water, & ! specific heat of liquid water (J kg-1 K-1)
- LH_fus, & ! latent heat of fusion (J kg-1)
- iden_air, & ! intrinsic density of air (kg m-3)
- iden_ice, & ! intrinsic density of ice (kg m-3)
- iden_water ! intrinsic density of liquid water (kg m-3)
+ ! global metadata
+ USE globalData,only:flux_meta ! metadata on the model fluxes
- ! provide access to the derived types to define the data structures
+ ! derived types to define the data structures
USE data_types,only:&
var_i, & ! data vector (i4b)
var_d, & ! data vector (rkind)
+ var_flagVec, & ! data vector with variable length dimension (i4b)
var_ilength, & ! data vector with variable length dimension (i4b)
- var_dlength ! data vector with variable length dimension (rkind)
+ var_dlength, & ! data vector with variable length dimension (rkind)
+ zLookup, & ! data vector with variable length dimension (rkind)
+ model_options ! defines the model decisions
! provide access to indices that define elements of the data structures
- USE var_lookup,only:iLookDIAG ! named variables for structure elements
- USE var_lookup,only:iLookPROG ! named variables for structure elements
- USE var_lookup,only:iLookDERIV ! named variables for structure elements
- USE var_lookup,only:iLookPARAM ! named variables for structure elements
- USE var_lookup,only:iLookINDEX ! named variables for structure elements
-
- ! provide access to routines to update states
- USE updatStateSundials_module,only:updateSnowSundials ! update snow states
- USE updatStateSundials_module,only:updateSoilSundials ! update soil states
-
- ! provide access to functions for the constitutive functions and derivatives
- USE snow_utils_module,only:fracliquid ! compute the fraction of liquid water (snow)
- USE snow_utils_module,only:dFracLiq_dTk ! differentiate the freezing curve w.r.t. temperature (snow)
- USE soil_utils_module,only:dTheta_dTk ! differentiate the freezing curve w.r.t. temperature (soil)
- USE soil_utils_module,only:dTheta_dPsi ! derivative in the soil water characteristic (soil)
- USE soil_utils_module,only:dPsi_dTheta ! derivative in the soil water characteristic (soil)
- USE soil_utils_module,only:matricHead ! compute the matric head based on volumetric water content
- USE soil_utils_module,only:volFracLiq ! compute volumetric fraction of liquid water
- USE soil_utils_module,only:crit_soilT ! compute critical temperature below which ice exists
- USE soil_utilsAddSundials_module,only:liquidHeadSundials ! compute the liquid water matric potential
- USE soil_utilsAddSundials_module,only:d2Theta_dPsi2
- USE soil_utilsAddSundials_module,only:d2Theta_dTk2
-
- ! IEEE checks
- USE, intrinsic :: ieee_arithmetic ! check values (NaN, etc.)
+ USE var_lookup,only:iLookFLUX ! named variables for structure elements
+ USE var_lookup,only:iLookPROG ! named variables for structure elements
+ USE var_lookup,only:iLookDIAG ! named variables for structure elements
+ USE var_lookup,only:iLookPARAM ! named variables for structure elements
+ USE var_lookup,only:iLookINDEX ! named variables for structure elements
+ ! look up structure for variable types
+ USE var_lookup,only:iLookVarType
+
+ ! constants
+ USE multiconst,only:&
+ Tfreeze, & ! freezing temperature (K)
+ LH_fus, & ! latent heat of fusion (J kg-1)
+ LH_vap, & ! latent heat of vaporization (J kg-1)
+ iden_ice, & ! intrinsic density of ice (kg m-3)
+ iden_water, & ! intrinsic density of liquid water (kg m-3)
+ ! specific heat
+ Cp_air, & ! specific heat of air (J kg-1 K-1)
+ Cp_ice, & ! specific heat of ice (J kg-1 K-1)
+ Cp_soil, & ! specific heat of soil (J kg-1 K-1)
+ Cp_water ! specific heat of liquid water (J kg-1 K-1)
+
+ ! safety: set private unless specified otherwise
implicit none
private
- public::updateVarsSundials
+ public::varSubstepSundials
+
+ ! algorithmic parameters
+ real(rkind),parameter :: verySmall=1.e-6_rkind ! used as an additive constant to check if substantial difference among real numbers
contains
+
! **********************************************************************************************************
- ! public subroutine updateVarsSundials: compute diagnostic variables and derivatives for Sundials Jacobian
+ ! public subroutine varSubstepSundials: run the model for a collection of substeps for a given state subset
! **********************************************************************************************************
- subroutine updateVarsSundials(&
- ! input
- dt, & ! intent(in): time step
- computJac, & ! intent(in): logical flag if computing Jacobian for Sundials solver
- do_adjustTemp, & ! intent(in): logical flag to adjust temperature to account for the energy used in melt+freeze
- mpar_data, & ! intent(in): model parameters for a local HRU
- indx_data, & ! intent(in): indices defining model states and layers
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- mLayerVolFracWatPrev, & ! intent(in): previous vector of total water matric potential (m)
- mLayerMatricHeadPrev, & ! intent(in): previous vector of volumetric total water content (-)
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables
- ! output: variables for the vegetation canopy
- scalarCanopyTempTrial, & ! intent(inout): trial value of canopy temperature (K)
- scalarCanopyWatTrial, & ! intent(inout): trial value of canopy total water (kg m-2)
- scalarCanopyLiqTrial, & ! intent(inout): trial value of canopy liquid water (kg m-2)
- scalarCanopyIceTrial, & ! intent(inout): trial value of canopy ice content (kg m-2)
- scalarCanopyTempPrime, & ! intent(inout): trial value of time derivative canopy temperature (K)
- scalarCanopyWatPrime, & ! intent(inout): trial value of time derivative canopy total water (kg m-2)
- scalarCanopyLiqPrime, & ! intent(inout): trial value of time derivative canopy liquid water (kg m-2)
- scalarCanopyIcePrime, & ! intent(inout): trial value of time derivative canopy ice content (kg m-2)
- ! output: variables for the snow-soil domain
- mLayerTempTrial, & ! intent(inout): trial vector of layer temperature (K)
- mLayerVolFracWatTrial, & ! intent(inout): trial vector of volumetric total water content (-)
- mLayerVolFracLiqTrial, & ! intent(inout): trial vector of volumetric liquid water content (-)
- mLayerVolFracIceTrial, & ! intent(inout): trial vector of volumetric ice water content (-)
- mLayerMatricHeadTrial, & ! intent(inout): trial vector of total water matric potential (m)
- mLayerMatricHeadLiqTrial, & ! intent(inout): trial vector of liquid water matric potential (m)
- mLayerTempPrime, & ! intent(inout): trial vector of time derivative layer temperature (K)
- mLayerVolFracWatPrime, & ! intent(inout): trial vector of time derivative volumetric total water content (-)
- mLayerVolFracLiqPrime, & ! intent(inout): trial vector of time derivative volumetric liquid water content (-)
- mLayerVolFracIcePrime, & ! intent(inout): trial vector of time derivative volumetric ice water content (-)
- mLayerMatricHeadPrime, & ! intent(inout): trial vector of time derivative total water matric potential (m)
- mLayerMatricHeadLiqPrime, & ! intent(inout): trial vector of time derivative liquid water matric potential (m)
- ! output: error control
- err,message) ! intent(out): error control
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------------------------------------------------------
+ subroutine varSubstepSundials(&
+ ! input: model control
+ dt, & ! intent(in) : time step (s)
+ dtInit, & ! intent(in) : initial time step (seconds)
+ dt_min, & ! intent(in) : minimum time step (seconds)
+ nState, & ! intent(in) : total number of state variables
+ doAdjustTemp, & ! intent(in) : flag to indicate if we adjust the temperature
+ firstSubStep, & ! intent(in) : flag to denote first sub-step
+ firstFluxCall, & ! intent(inout) : flag to indicate if we are processing the first flux call
+ computeVegFlux, & ! intent(in) : flag to denote if computing energy flux over vegetation
+ scalarSolution, & ! intent(in) : flag to denote implementing the scalar solution
+ iStateSplit, & ! intent(in) : index of the state in the splitting operation
+ fluxMask, & ! intent(in) : mask for the fluxes used in this given state subset
+ fluxCount, & ! intent(inout) : number of times that fluxes are updated (should equal nSubsteps)
+ ! input/output: data structures
+ model_decisions, & ! intent(in) : model decisions
+ lookup_data, & ! intent(in) : lookup tables
+ type_data, & ! intent(in) : type of vegetation and soil
+ attr_data, & ! intent(in) : spatial attributes
+ forc_data, & ! intent(in) : model forcing data
+ mpar_data, & ! intent(in) : model parameters
+ indx_data, & ! intent(inout) : index data
+ prog_data, & ! intent(inout) : model prognostic variables for a local HRU
+ diag_data, & ! intent(inout) : model diagnostic variables for a local HRU
+ flux_data, & ! intent(inout) : model fluxes for a local HRU
+ deriv_data, & ! intent(inout) : derivatives in model fluxes w.r.t. relevant state variables
+ bvar_data, & ! intent(in) : model variables for the local basin
+ ! output: model control
+ ixSaturation, & ! intent(inout) : index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ dtMultiplier, & ! intent(out) : substep multiplier (-)
+ nSubsteps, & ! intent(out) : number of substeps taken for a given split
+ failedMinimumStep, & ! intent(out) : flag to denote success of substepping for a given split
+ reduceCoupledStep, & ! intent(out) : flag to denote need to reduce the length of the coupled step
+ tooMuchMelt, & ! intent(out) : flag to denote that ice is insufficient to support melt
+ dt_out, & ! intent(out)
+ err,message) ! intent(out) : error code and error message
+ ! ---------------------------------------------------------------------------------------
+ ! structure allocations
+ USE allocspace_module,only:allocLocal ! allocate local data structures
+ ! simulation of fluxes and residuals given a trial state vector
+ USE systemSolv_module,only:systemSolv ! solve the system of equations for one time step
+ USE getVectorz_module,only:popStateVec ! populate the state vector
+ USE updateVarsSundials_module,only:updateVarsSundials ! update prognostic variables
+ USE getVectorzAddSundials_module,only:varExtractSundials
+ ! identify name of variable type (for error message)
+ USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages
+ USE systemSolvSundials_module,only:systemSolvSundials
implicit none
- ! input
- real(rkind) ,intent(in) :: dt ! time step
- logical(lgt) ,intent(in) :: computJac ! flag if computing Jacobian for Sundials solver
- logical(lgt) ,intent(in) :: do_adjustTemp ! flag to adjust temperature to account for the energy used in melt+freeze
- type(var_dlength),intent(in) :: mpar_data ! model parameters for a local HRU
- type(var_ilength),intent(in) :: indx_data ! indices defining model states and layers
- type(var_dlength),intent(in) :: prog_data ! prognostic variables for a local HRU
- real(rkind),intent(in) :: mLayerVolFracWatPrev(:) ! previous vector of total water matric potential (m)
- real(rkind),intent(in) :: mLayerMatricHeadPrev(:) ! previous vector of volumetric total water content (-)
- type(var_dlength),intent(inout) :: diag_data ! diagnostic variables for a local HRU
- type(var_dlength),intent(inout) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
- ! output: variables for the vegetation canopy
- real(rkind),intent(inout) :: scalarCanopyTempTrial ! trial value of canopy temperature (K)
- real(rkind),intent(inout) :: scalarCanopyWatTrial ! trial value of canopy total water (kg m-2)
- real(rkind),intent(inout) :: scalarCanopyLiqTrial ! trial value of canopy liquid water (kg m-2)
- real(rkind),intent(inout) :: scalarCanopyIceTrial ! trial value of canopy ice content (kg m-2)
- real(rkind),intent(inout) :: scalarCanopyTempPrime ! trial value of time derivative canopy temperature (K)
- real(rkind),intent(inout) :: scalarCanopyWatPrime ! trial value of time derivative canopy total water (kg m-2)
- real(rkind),intent(inout) :: scalarCanopyLiqPrime ! trial value of time derivative canopy liquid water (kg m-2)
- real(rkind),intent(inout) :: scalarCanopyIcePrime ! trial value of time derivative canopy ice content (kg m-2)
- ! output: variables for the snow-soil domain
- real(rkind),intent(inout) :: mLayerTempTrial(:) ! trial vector of layer temperature (K)
- real(rkind),intent(inout) :: mLayerVolFracWatTrial(:) ! trial vector of volumetric total water content (-)
- real(rkind),intent(inout) :: mLayerVolFracLiqTrial(:) ! trial vector of volumetric liquid water content (-)
- real(rkind),intent(inout) :: mLayerVolFracIceTrial(:) ! trial vector of volumetric ice water content (-)
- real(rkind),intent(inout) :: mLayerMatricHeadTrial(:) ! trial vector of total water matric potential (m)
- real(rkind),intent(inout) :: mLayerMatricHeadLiqTrial(:) ! trial vector of liquid water matric potential (m)
- real(rkind),intent(inout) :: mLayerTempPrime(:) ! trial vector of time derivative layer temperature (K)
- real(rkind),intent(inout) :: mLayerVolFracWatPrime(:) ! trial vector of time derivative volumetric total water content (-)
- real(rkind),intent(inout) :: mLayerVolFracLiqPrime(:) ! trial vector of time derivative volumetric liquid water content (-)
- real(rkind),intent(inout) :: mLayerVolFracIcePrime(:) ! trial vector of time derivative volumetric ice water content (-)
- real(rkind),intent(inout) :: mLayerMatricHeadPrime(:) ! trial vector of time derivative total water matric potential (m)
- real(rkind),intent(inout) :: mLayerMatricHeadLiqPrime(:) ! trial vector of time derivative liquid water matric potential (m)
- ! output: error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! --------------------------------------------------------------------------------------------------------------------------------
+ ! ---------------------------------------------------------------------------------------
+ ! * dummy variables
+ ! ---------------------------------------------------------------------------------------
+ ! input: model control
+ real(rkind),intent(in) :: dt ! time step (seconds)
+ real(rkind),intent(in) :: dtInit ! initial time step (seconds)
+ real(rkind),intent(in) :: dt_min ! minimum time step (seconds)
+ integer(i4b),intent(in) :: nState ! total number of state variables
+ logical(lgt),intent(in) :: doAdjustTemp ! flag to indicate if we adjust the temperature
+ logical(lgt),intent(in) :: firstSubStep ! flag to indicate if we are processing the first sub-step
+ logical(lgt),intent(inout) :: firstFluxCall ! flag to define the first flux call
+ logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
+ logical(lgt),intent(in) :: scalarSolution ! flag to denote implementing the scalar solution
+ integer(i4b),intent(in) :: iStateSplit ! index of the state in the splitting operation
+ type(var_flagVec),intent(in) :: fluxMask ! flags to denote if the flux is calculated in the given state subset
+ type(var_ilength),intent(inout) :: fluxCount ! number of times that the flux is updated (should equal nSubsteps)
+ ! input/output: data structures
+ type(model_options),intent(in) :: model_decisions(:) ! model decisions
+ type(zLookup),intent(in) :: lookup_data ! lookup tables
+ 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(inout) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
+ type(var_dlength),intent(in) :: bvar_data ! model variables for the local basin
+ ! output: model control
+ integer(i4b),intent(inout) :: ixSaturation ! index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ real(rkind),intent(out) :: dtMultiplier ! substep multiplier (-)
+ integer(i4b),intent(out) :: nSubsteps ! number of substeps taken for a given split
+ logical(lgt),intent(out) :: failedMinimumStep ! flag to denote success of substepping for a given split
+ logical(lgt),intent(out) :: reduceCoupledStep ! flag to denote need to reduce the length of the coupled step
+ logical(lgt),intent(out) :: tooMuchMelt ! flag to denote that ice is insufficient to support melt
+ real(qp),intent(out) :: dt_out
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! ---------------------------------------------------------------------------------------
+ ! * general local variables
+ ! ---------------------------------------------------------------------------------------
+ ! error control
+ character(LEN=256) :: cmessage ! error message of downwind routine
! general local variables
- integer(i4b) :: iState ! index of model state variable
- integer(i4b) :: iLayer ! index of layer within the snow+soil domain
- integer(i4b) :: ixFullVector ! index within full state vector
- integer(i4b) :: ixDomainType ! name of a given model domain
- integer(i4b) :: ixControlIndex ! index within a given model domain
- integer(i4b) :: ixOther,ixOtherLocal ! index of the coupled state variable within the (full, local) vector
- logical(lgt) :: isCoupled ! .true. if a given variable shared another state variable in the same control volume
- logical(lgt) :: isNrgState ! .true. if a given variable is an energy state
- logical(lgt),allocatable :: computedCoupling(:) ! .true. if computed the coupling for a given state variable
- real(rkind) :: scalarVolFracLiq ! volumetric fraction of liquid water (-)
- real(rkind) :: scalarVolFracIce ! volumetric fraction of ice (-)
- real(rkind) :: scalarVolFracLiqPrime ! time derivative volumetric fraction of liquid water (-)
- real(rkind) :: scalarVolFracIcePrime ! time derivative volumetric fraction of ice (-)
- real(rkind) :: Tcrit ! critical soil temperature below which ice exists (K)
- real(rkind) :: xTemp ! temporary temperature (K)
- real(rkind) :: fLiq ! fraction of liquid water (-)
- real(rkind) :: effSat ! effective saturation (-)
- real(rkind) :: avPore ! available pore space (-)
- character(len=256) :: cMessage ! error message of downwind routine
- logical(lgt),parameter :: printFlag=.false. ! flag to turn on printing
- ! iterative solution for temperature
- real(rkind) :: meltNrg ! energy for melt+freeze (J m-3)
- real(rkind) :: residual ! residual in the energy equation (J m-3)
- real(rkind) :: derivative ! derivative in the energy equation (J m-3 K-1)
- real(rkind) :: tempInc ! iteration increment (K)
- integer(i4b) :: iter ! iteration index
- integer(i4b) :: niter ! number of iterations
- integer(i4b),parameter :: maxiter=100 ! maximum number of iterations
- real(rkind),parameter :: nrgConvTol=1.e-4_rkind ! convergence tolerance for energy (J m-3)
- real(rkind),parameter :: tempConvTol=1.e-6_rkind ! convergence tolerance for temperature (K)
- real(rkind) :: critDiff ! temperature difference from critical (K)
- real(rkind) :: tempMin ! minimum bracket for temperature (K)
- real(rkind) :: tempMax ! maximum bracket for temperature (K)
- logical(lgt) :: bFlag ! flag to denote that iteration increment was constrained using bi-section
- real(rkind),parameter :: epsT=1.e-7_rkind ! small interval above/below critical temperature (K)
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! make association with variables in the data structures
- associate(&
- ! number of model layers, and layer type
- nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1) ,& ! intent(in): [i4b] total number of snow layers
- nSoil => indx_data%var(iLookINDEX%nSoil)%dat(1) ,& ! intent(in): [i4b] total number of soil layers
- nLayers => indx_data%var(iLookINDEX%nLayers)%dat(1) ,& ! intent(in): [i4b] total number of snow and soil layers
- mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
- ! indices defining model states and layers
- 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)
- ! indices in the full vector for specific domains
- 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
- ! mapping between the full state vector and the state subset
- ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset for each state in the full state vector
- ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] [state subset] list of indices of the full state vector in the state subset
- ! type of domain, type of state variable, and index of control volume within domain
- ixDomainType_subset => indx_data%var(iLookINDEX%ixDomainType_subset)%dat ,& ! intent(in): [i4b(:)] [state subset] id of domain for desired model state variables
- ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of the control volume for different domains (veg, snow, soil)
- ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (iname_nrgLayer...)
- ! 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 model 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 (-)
- ! model diagnostic variables (heat capacity)
- canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp ] canopy depth (m)
- scalarBulkVolHeatCapVeg => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg)%dat(1),& ! intent(in): [dp ] volumetric heat capacity of the vegetation (J m-3 K-1)
- mLayerVolHtCapBulk => diag_data%var(iLookDIAG%mLayerVolHtCapBulk)%dat ,& ! intent(in): [dp(:)] volumetric heat capacity in each layer (J m-3 K-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 (-)
- ! model states for the vegetation canopy
- scalarCanairTemp => prog_data%var(iLookPROG%scalarCanairTemp)%dat(1) ,& ! intent(in): [dp] temperature of the canopy air space (K)
- scalarCanopyTemp => prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1) ,& ! intent(in): [dp] temperature of the vegetation canopy (K)
- scalarCanopyWat => prog_data%var(iLookPROG%scalarCanopyWat)%dat(1) ,& ! intent(in): [dp] mass of total water on the vegetation canopy (kg m-2)
- ! model state variable vectors for the snow-soil layers
- mLayerTemp => prog_data%var(iLookPROG%mLayerTemp)%dat ,& ! intent(in): [dp(:)] temperature of each snow/soil layer (K)
- mLayerVolFracWat => prog_data%var(iLookPROG%mLayerVolFracWat)%dat ,& ! intent(in): [dp(:)] volumetric fraction of total water (-)
- mLayerMatricHead => prog_data%var(iLookPROG%mLayerMatricHead)%dat ,& ! intent(in): [dp(:)] total water matric potential (m)
- mLayerMatricHeadLiq => diag_data%var(iLookDIAG%mLayerMatricHeadLiq)%dat ,& ! intent(in): [dp(:)] liquid water matric potential (m)
- ! model diagnostic variables from a previous solution
- scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! intent(in): [dp(:)] mass of liquid water on the vegetation canopy (kg m-2)
- scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! intent(in): [dp(:)] mass of ice on the vegetation canopy (kg m-2)
- mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat ,& ! intent(in): [dp(:)] volumetric fraction of liquid water (-)
- mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(in): [dp(:)] volumetric fraction of ice (-)
- ! derivatives
- dVolTot_dPsi0 => deriv_data%var(iLookDERIV%dVolTot_dPsi0 )%dat ,& ! intent(out): [dp(:)] derivative in total water content w.r.t. total water matric potential
- dPsiLiq_dPsi0 => deriv_data%var(iLookDERIV%dPsiLiq_dPsi0 )%dat ,& ! intent(out): [dp(:)] derivative in liquid water matric pot w.r.t. the total water matric pot (-)
- dPsiLiq_dTemp => deriv_data%var(iLookDERIV%dPsiLiq_dTemp )%dat ,& ! intent(out): [dp(:)] derivative in the liquid water matric potential w.r.t. temperature
- mLayerdTheta_dTk => deriv_data%var(iLookDERIV%mLayerdTheta_dTk)%dat ,& ! intent(out): [dp(:)] derivative of volumetric liquid water content w.r.t. temperature
- dTheta_dTkCanopy => deriv_data%var(iLookDERIV%dTheta_dTkCanopy)%dat(1) ,& ! intent(out): [dp] derivative of volumetric liquid water content w.r.t. temperature
- dFracLiqSnow_dTk => deriv_data%var(iLookDERIV%dFracLiqSnow_dTk )%dat ,& ! intent(out): [dp(:)] derivative in fraction of liquid snow w.r.t. temperature
- dFracLiqVeg_dTkCanopy => deriv_data%var(iLookDERIV%dFracLiqVeg_dTkCanopy )%dat(1) ,& ! intent(out): [dp ] derivative in fraction of (throughfall + drainage) w.r.t. temperature
- ! derivatives inside solver for Jacobian only
- dVolHtCapBulk_dPsi0 => deriv_data%var(iLookDERIV%dVolHtCapBulk_dPsi0 )%dat ,& ! intent(out): [dp(:)] derivative in bulk heat capacity w.r.t. matric potential
- dVolHtCapBulk_dTheta => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTheta )%dat ,& ! intent(out): [dp(:)] derivative in bulk heat capacity w.r.t. volumetric water content
- dVolHtCapBulk_dCanWat => deriv_data%var(iLookDERIV%dVolHtCapBulk_dCanWat)%dat(1) ,& ! intent(out): [dp ] derivative in bulk heat capacity w.r.t. volumetric water content
- dVolHtCapBulk_dTk => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTk )%dat ,& ! intent(out): [dp(:)] derivative in bulk heat capacity w.r.t. temperature
- dVolHtCapBulk_dTkCanopy => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTkCanopy)%dat(1) ,& ! intent(out): [dp ] derivative in bulk heat capacity w.r.t. temperature
- d2VolTot_d2Psi0 => deriv_data%var(iLookDERIV%d2VolTot_d2Psi0 )%dat ,& ! intent(out): [dp(:)] second derivative in total water content w.r.t. total water matric potential
- mLayerd2Theta_dTk2 => deriv_data%var(iLookDERIV%mLayerd2Theta_dTk2 )%dat ,& ! intent(out): [dp(:)] second derivative of volumetric liquid water content w.r.t. temperature
- d2Theta_dTkCanopy2 => deriv_data%var(iLookDERIV%d2Theta_dTkCanopy2 )%dat(1) & ! intent(out): [dp ] second derivative of volumetric liquid water content w.r.t. temperature
- ) ! association with variables in the data structures
-
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------------------------------------------------------
+ integer(i4b) :: iVar ! index of variables in data structures
+ integer(i4b) :: iSoil ! index of soil layers
+ integer(i4b) :: ixLayer ! index in a given domain
+ integer(i4b), dimension(1) :: ixMin,ixMax ! bounds of a given flux vector
+ ! time stepping
+ real(rkind) :: dtSum ! sum of time from successful steps (seconds)
+ real(rkind) :: dt_wght ! weight given to a given flux calculation
+ real(rkind) :: dtSubstep ! length of a substep (s)
+ ! adaptive sub-stepping for the explicit solution
+ logical(lgt) :: failedSubstep ! flag to denote success of substepping for a given split
+ real(rkind),parameter :: safety=0.85_rkind ! safety factor in adaptive sub-stepping
+ real(rkind),parameter :: reduceMin=0.1_rkind ! mimimum factor that time step is reduced
+ real(rkind),parameter :: increaseMax=4.0_rkind ! maximum factor that time step is increased
+ ! adaptive sub-stepping for the implicit solution
+ integer(i4b),parameter :: n_inc=5 ! minimum number of iterations to increase time step
+ integer(i4b),parameter :: n_dec=15 ! maximum number of iterations to decrease time step
+ real(rkind),parameter :: F_inc = 1.25_rkind ! factor used to increase time step
+ real(rkind),parameter :: F_dec = 0.90_rkind ! factor used to decrease time step
+ ! state and flux vectors
+ real(rkind) :: untappedMelt(nState) ! un-tapped melt energy (J m-3 s-1)
+ real(rkind) :: stateVecInit(nState) ! initial state vector (mixed units)
+ real(rkind) :: stateVecTrial(nState) ! trial state vector (mixed units)
+ real(rkind) :: stateVecPrime(nState) ! trial state vector (mixed units)
+ type(var_dlength) :: flux_temp ! temporary model fluxes
+ ! flags
+ logical(lgt) :: firstSplitOper ! flag to indicate if we are processing the first flux call in a splitting operation
+ logical(lgt) :: checkMassBalance ! flag to check the mass balance
+ logical(lgt) :: checkNrgBalance
+ logical(lgt) :: waterBalanceError ! flag to denote that there is a water balance error
+ logical(lgt) :: nrgFluxModified ! flag to denote that the energy fluxes were modified
+ ! energy fluxes
+ real(rkind) :: sumCanopyEvaporation ! sum of canopy evaporation/condensation (kg m-2 s-1)
+ real(rkind) :: sumLatHeatCanopyEvap ! sum of latent heat flux for evaporation from the canopy to the canopy air space (W m-2)
+ real(rkind) :: sumSenHeatCanopy ! sum of sensible heat flux from the canopy to the canopy air space (W m-2)
+ real(rkind) :: sumSoilCompress
+ real(rkind),allocatable :: sumLayerCompress(:)
+ ! ---------------------------------------------------------------------------------------
+ ! point to variables in the data structures
+ ! ---------------------------------------------------------------------------------------
+ globalVars: associate(&
+ ! number of layers
+ nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1) ,& ! intent(in): [i4b] number of snow layers
+ nSoil => indx_data%var(iLookINDEX%nSoil)%dat(1) ,& ! intent(in): [i4b] number of soil layers
+ nLayers => indx_data%var(iLookINDEX%nLayers)%dat(1) ,& ! intent(in): [i4b] total number of layers
+ nSoilOnlyHyd => indx_data%var(iLookINDEX%nSoilOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the soil domain
+ mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
+ ! mapping between state vectors and control volumes
+ ixLayerActive => indx_data%var(iLookINDEX%ixLayerActive)%dat ,& ! intent(in): [i4b(:)] list of indices for all active layers (inactive=integerMissing)
+ ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] mapping of full state vector to the state subset
+ ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of control volume for different domains (veg, snow, soil)
+ ! model state variables (vegetation canopy)
+ scalarCanairTemp => prog_data%var(iLookPROG%scalarCanairTemp)%dat(1) ,& ! intent(inout): [dp] temperature of the canopy air space (K)
+ scalarCanopyTemp => prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1) ,& ! intent(inout): [dp] temperature of the vegetation canopy (K)
+ scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! intent(inout): [dp] mass of ice on the vegetation canopy (kg m-2)
+ scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! intent(inout): [dp] mass of liquid water on the vegetation canopy (kg m-2)
+ scalarCanopyWat => prog_data%var(iLookPROG%scalarCanopyWat)%dat(1) ,& ! intent(inout): [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(inout): [dp(:)] temperature of each snow/soil layer (K)
+ mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(inout): [dp(:)] volumetric fraction of ice (-)
+ mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat ,& ! intent(inout): [dp(:)] volumetric fraction of liquid water (-)
+ mLayerVolFracWat => prog_data%var(iLookPROG%mLayerVolFracWat)%dat ,& ! intent(inout): [dp(:)] volumetric fraction of total water (-)
+ mLayerMatricHead => prog_data%var(iLookPROG%mLayerMatricHead)%dat ,& ! intent(inout): [dp(:)] matric head (m)
+ mLayerMatricHeadLiq => diag_data%var(iLookDIAG%mLayerMatricHeadLiq)%dat & ! intent(inout): [dp(:)] matric potential of liquid water (m)
+ ) ! end association with variables in the data structures
+ ! *********************************************************************************************************************************************************
+ ! *********************************************************************************************************************************************************
+ ! Procedure starts here
! initialize error control
- err=0; message='updateVarsSundials/'
+ err=0; message='varSubstepSundials/'
+
+ ! initialize flag for the success of the substepping
+ failedMinimumStep=.false.
+
+ ! initialize the length of the substep
+ dtSubstep = dtInit
- ! allocate space and assign values to the flag vector
- allocate(computedCoupling(size(ixMapSubset2Full)),stat=err) ! .true. if computed the coupling for a given state variable
- if(err/=0)then; message=trim(message)//'problem allocating computedCoupling'; return; endif
- computedCoupling(:)=.false.
+ ! allocate space for the temporary model flux structure
+ call allocLocal(flux_meta(:),flux_temp,nSnow,nSoil,err,cmessage)
+ if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
- ! loop through model state variables
- do iState=1,size(ixMapSubset2Full)
+ ! initialize the model fluxes (some model fluxes are not computed in the iterations)
+ do iVar=1,size(flux_data%var)
+ flux_temp%var(iVar)%dat(:) = flux_data%var(iVar)%dat(:)
+ end do
- ! check the need for the computations
- if(computedCoupling(iState)) cycle
+ ! initialize the total energy fluxes (modified in updateProgSundials)
+ sumCanopyEvaporation = 0._rkind ! canopy evaporation/condensation (kg m-2 s-1)
+ sumLatHeatCanopyEvap = 0._rkind ! latent heat flux for evaporation from the canopy to the canopy air space (W m-2)
+ sumSenHeatCanopy = 0._rkind ! sensible heat flux from the canopy to the canopy air space (W m-2)
+ sumSoilCompress = 0._rkind ! total soil compression
+ allocate(sumLayerCompress(nSoil)); sumLayerCompress = 0._rkind ! soil compression by layer
+
+ ! define the first flux call in a splitting operation
+ firstSplitOper = (.not.scalarSolution .or. iStateSplit==1)
+
+ ! initialize subStep
+ dtSum = 0._rkind ! keep track of the portion of the time step that is completed
+ nSubsteps = 0
+
+ ! loop through substeps
+ ! NOTE: continuous do statement with exit clause
+ substeps: do
+
+ ! initialize error control
+ err=0; message='varSubstepSundials/'
+
+ !write(*,'(a,1x,3(f13.2,1x))') '***** new subStep: dtSubstep, dtSum, dt = ', dtSubstep, dtSum, dt
+ !print*, 'scalarCanopyIce = ', prog_data%var(iLookPROG%scalarCanopyIce)%dat(1)
+ !print*, 'scalarCanopyTemp = ', prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1)
+
+ ! -----
+ ! * populate state vectors...
+ ! ---------------------------
+
+ ! initialize state vectors
+ call popStateVec(&
+ ! input
+ nState, & ! intent(in): number of desired state variables
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ diag_data, & ! intent(in): model diagnostic variables for a local HRU
+ indx_data, & ! intent(in): indices defining model states and layers
+ ! output
+ stateVecInit, & ! intent(out): initial model state vector (mixed units)
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; endif ! (check for errors)
! -----
- ! - compute indices...
- ! --------------------
-
- ! get domain type, and index of the control volume within the domain
- ixFullVector = ixMapSubset2Full(iState) ! index within full state vector
- ixDomainType = ixDomainType_subset(iState) ! named variables defining the domain (iname_cas, iname_veg, etc.)
- ixControlIndex = ixControlVolume(ixFullVector) ! index within a given domain
-
- ! get the layer index
- select case(ixDomainType)
- case(iname_cas); cycle ! canopy air space: do nothing
- case(iname_veg); iLayer = 0
- case(iname_snow); iLayer = ixControlIndex
- case(iname_soil); iLayer = ixControlIndex + nSnow
- case(iname_aquifer); cycle ! aquifer: do nothing
- case default; err=20; message=trim(message)//'expect case to be iname_cas, iname_veg, iname_snow, iname_soil, iname_aquifer'; return
- end select
-
- ! get the index of the other (energy or mass) state variable within the full state vector
- select case(ixDomainType)
- case(iname_veg) ; ixOther = merge(ixHydCanopy(1), ixNrgCanopy(1), ixStateType(ixFullVector)==iname_nrgCanopy)
- case(iname_snow, iname_soil); ixOther = merge(ixHydLayer(iLayer),ixNrgLayer(iLayer),ixStateType(ixFullVector)==iname_nrgLayer)
- case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
- end select
-
- ! get the index in the local state vector
- ixOtherLocal = ixMapFull2Subset(ixOther) ! ixOtherLocal could equal integerMissing
- if(ixOtherLocal/=integerMissing) computedCoupling(ixOtherLocal)=.true.
-
- ! check if we have a coupled solution
- isCoupled = (ixOtherLocal/=integerMissing)
-
- ! check if we are an energy state
- isNrgState = (ixStateType(ixFullVector)==iname_nrgCanopy .or. ixStateType(ixFullVector)==iname_nrgLayer)
-
- if(printFlag)then
- print*, 'iState = ', iState, size(ixMapSubset2Full)
- print*, 'ixFullVector = ', ixFullVector
- print*, 'ixDomainType = ', ixDomainType
- print*, 'ixControlIndex = ', ixControlIndex
- print*, 'ixOther = ', ixOther
- print*, 'ixOtherLocal = ', ixOtherLocal
- print*, 'do_adjustTemp = ', do_adjustTemp
- print*, 'isCoupled = ', isCoupled
- print*, 'isNrgState = ', isNrgState
+ ! * iterative solution...
+ ! -----------------------
+ ! solve the system of equations for a given state subset
+ call systemSolvSundials(&
+ ! input: model control
+ dtSubstep, & ! intent(in): time step (s)
+ nState, & ! intent(in): total number of state variables
+ firstSubStep, & ! intent(in): flag to denote first sub-step
+ firstFluxCall, & ! intent(inout): flag to indicate if we are processing the first flux call
+ firstSplitOper, & ! intent(in): flag to indicate if we are processing the first flux call in a splitting operation
+ computeVegFlux, & ! intent(in): flag to denote if computing energy flux over vegetation
+ scalarSolution, & ! intent(in): flag to denote if implementing the scalar solution
+ ! input/output: data structures
+ lookup_data, & ! intent(in): lookup tables
+ type_data, & ! intent(in): type of vegetation and soil
+ attr_data, & ! intent(in): spatial attributes
+ forc_data, & ! intent(in): model forcing data
+ mpar_data, & ! intent(in): model parameters
+ indx_data, & ! intent(inout): index data
+ prog_data, & ! intent(inout): model prognostic variables for a local HRU
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ flux_temp, & ! intent(inout): model fluxes for a local HRU
+ bvar_data, & ! intent(in): model variables for the local basin
+ model_decisions, & ! intent(in): model decisions
+ stateVecInit, & ! intent(in): initial state vector
+ ! output: model control
+ deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables
+ ixSaturation, & ! intent(inout): index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ stateVecTrial, & ! intent(out): updated state vector
+ stateVecPrime, & ! intent(out): updated state vector
+ reduceCoupledStep, & ! intent(out): flag to reduce the length of the coupled step
+ tooMuchMelt, & ! intent(out): flag to denote that ice is insufficient to support melt
+ dt_out, & ! intent(out): time step (s)
+ err,cmessage) ! intent(out): error code and error message
+
+ if(err/=0)then
+ message=trim(message)//trim(cmessage)
+ if(err>0) return
endif
- ! =======================================================================================================================================
- ! =======================================================================================================================================
- ! =======================================================================================================================================
- ! =======================================================================================================================================
- ! =======================================================================================================================================
- ! =======================================================================================================================================
-
- ! update hydrology state variables for the uncoupled solution
- if(.not.isNrgState .and. .not.isCoupled)then
-
- if(.not.computJac) stop 1 ! this does not work yet? FIX
-
- ! update the total water from volumetric liquid water
- if(ixStateType(ixFullVector)==iname_liqCanopy .or. ixStateType(ixFullVector)==iname_liqLayer)then
- select case(ixDomainType)
- case(iname_veg)
- scalarCanopyWatTrial = scalarCanopyLiqTrial + scalarCanopyIceTrial
- scalarCanopyWatPrime = scalarCanopyLiqPrime + scalarCanopyIcePrime
- case(iname_snow)
- mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer)*iden_ice/iden_water
- mLayerVolFracWatPrime(iLayer) = mLayerVolFracLiqPrime(iLayer) + mLayerVolFracIcePrime(iLayer)*iden_ice/iden_water
- case(iname_soil)
- mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer) ! no volume expansion
- mLayerVolFracWatPrime(iLayer) = mLayerVolFracLiqPrime(iLayer) + mLayerVolFracIcePrime(iLayer)
- case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, or iname_soil'; return
- end select
- endif
+ ! set untapped melt energy to zero
+ untappedMelt(:) = 0._rkind
+
+ ! if too much melt or need to reduce length of the coupled step then return
+ ! NOTE: need to go all the way back to coupled_em and merge snow layers, as all splitting operations need to occur with the same layer geometry
+ if(tooMuchMelt .or. reduceCoupledStep) return
+
+ ! identify failure
+ failedSubstep = (err<0)
+
+ ! reduce step based on failure
+ if(failedSubstep)then
+ err=0; message='varSubstepSundials/' ! recover from failed convergence
+ dtMultiplier = 0.5_rkind ! system failure: step halving
+ else
+
+ endif ! switch between failure and success
+
+ ! check if we failed the substep
+ if(failedSubstep)then
+
+ ! check that the substep is greater than the minimum step
+ if(dtSubstep*dtMultiplier<dt_min)then
+ ! --> exit, and either (1) try another solution method; or (2) reduce coupled step
+ failedMinimumStep=.true.
+ exit subSteps
- ! update the total water and the total water matric potential
- if(ixDomainType==iname_soil)then
- select case( ixStateType(ixFullVector) )
- ! --> update the total water from the liquid water matric potential
- case(iname_lmpLayer)
-
- effSat = volFracLiq(mLayerMatricHeadLiqTrial(ixControlIndex),vGn_alpha(ixControlIndex),0._rkind,1._rkind,vGn_n(ixControlIndex),vGn_m(ixControlIndex)) ! effective saturation
- avPore = theta_sat(ixControlIndex) - mLayerVolFracIceTrial(iLayer) - theta_res(ixControlIndex) ! available pore space
- mLayerVolFracLiqTrial(iLayer) = effSat*avPore + theta_res(ixControlIndex)
- mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer) ! no volume expansion
- mLayerVolFracWatPrime(iLayer) = mLayerVolFracLiqPrime(iLayer) + mLayerVolFracIcePrime(iLayer)
- mLayerMatricHeadTrial(ixControlIndex) = matricHead(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
- mLayerMatricHeadPrime(ixControlIndex) = dPsi_dTheta(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) * mLayerVolFracWatPrime(iLayer)
- !write(*,'(a,1x,i4,1x,3(f20.10,1x))') 'mLayerVolFracLiqTrial(iLayer) 1 = ', iLayer, mLayerVolFracLiqTrial(iLayer), mLayerVolFracIceTrial(iLayer), mLayerVolFracWatTrial(iLayer)
- ! --> update the total water from the total water matric potential
- case(iname_matLayer)
-
- mLayerVolFracWatTrial(iLayer) = volFracLiq(mLayerMatricHeadTrial(ixControlIndex),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
- mLayerVolFracWatPrime(iLayer) = dTheta_dPsi(mLayerMatricHeadTrial(ixControlIndex),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) *mLayerMatricHeadPrime(ixControlIndex)
- ! --> update the total water matric potential (assume already have mLayerVolFracWatTrial given block above)
- case(iname_liqLayer, iname_watLayer)
-
- mLayerMatricHeadTrial(ixControlIndex) = matricHead(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
- mLayerMatricHeadPrime(ixControlIndex) = dPsi_dTheta(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) * mLayerVolFracWatPrime(iLayer)
- case default; err=20; message=trim(message)//'expect iname_lmpLayer, iname_matLayer, iname_liqLayer, or iname_watLayer'; return
- end select
- endif ! if in the soil domain
-
- endif ! if hydrology state variable or uncoupled solution
-
- ! compute the critical soil temperature below which ice exists
- select case(ixDomainType)
- case(iname_veg, iname_snow); Tcrit = Tfreeze
- case(iname_soil); Tcrit = crit_soilT( mLayerMatricHeadTrial(ixControlIndex) )
- case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
- end select
-
- ! initialize temperature
- select case(ixDomainType)
- case(iname_veg); xTemp = scalarCanopyTempTrial
- case(iname_snow, iname_soil); xTemp = mLayerTempTrial(iLayer)
- case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
- end select
-
- ! define brackets for the root
- ! NOTE: start with an enormous range; updated quickly in the iterations
- tempMin = xTemp - 10._rkind
- tempMax = xTemp + 10._rkind
-
- ! get iterations (set to maximum iterations if adjusting the temperature)
- niter = merge(maxiter, 1, do_adjustTemp)
-
- ! iterate
- iterations: do iter=1,niter
-
- ! restrict temperature
- if(xTemp <= tempMin .or. xTemp >= tempMax)then
- xTemp = 0.5_rkind*(tempMin + tempMax) ! new value
- bFlag = .true.
+ else ! step is still OK
+ dtSubstep = dtSubstep*dtMultiplier
+ cycle subSteps
+ endif ! if step is less than the minimum
+
+ endif ! if failed the substep
+
+ ! -----
+ ! * update model fluxes...
+ ! ------------------------
+
+ ! NOTE: if we get to here then we are accepting the step
+
+ ! NOTE: we get to here if iterations are successful
+ if(err/=0)then
+ message=trim(message)//'expect err=0 if updating fluxes'
+ return
+ endif
+
+ ! identify the need to check the mass balance
+ checkMassBalance = .true. ! (.not.scalarSolution)
+ checkNrgBalance = .true.
+
+ ! update prognostic variables
+ call updateProgSundials(dt_out,nSnow,nSoil,nLayers,doAdjustTemp,computeVegFlux,untappedMelt,stateVecTrial,stateVecPrime,checkMassBalance, checkNrgBalance, & ! input: model control
+ lookup_data,mpar_data,indx_data,flux_temp,prog_data,diag_data,deriv_data, & ! input-output: data structures
+ waterBalanceError,nrgFluxModified,err,cmessage) ! output: flags and error control
+ if(err/=0)then
+ message=trim(message)//trim(cmessage)
+ if(err>0) return
+ endif
+
+ ! if water balance error then reduce the length of the coupled step
+ if(waterBalanceError)then
+ message=trim(message)//'water balance error'
+ reduceCoupledStep=.true.
+ err=-20; return
+ endif
+
+ if(globalPrintFlag)&
+ print*, trim(cmessage)//': dt = ', dtSubstep
+
+ ! recover from errors in prognostic update
+ if(err<0)then
+
+ ! modify step
+ err=0 ! error recovery
+ dtSubstep = dtSubstep/2._rkind
+
+ ! check minimum: fail minimum step if there is an error in the update
+ if(dtSubstep<dt_min)then
+ failedMinimumStep=.true.
+ exit subSteps
+ ! minimum OK -- try again
else
- bFlag = .false.
+ cycle substeps
endif
- ! -----
- ! - compute derivatives...
- ! ------------------------
-
- ! compute the derivative in bulk heat capacity w.r.t. total water content or water matric potential (m-1)
- ! compute the derivative in total water content w.r.t. total water matric potential (m-1)
- ! NOTE 1: valid for frozen and unfrozen conditions
- ! NOTE 2: for case "iname_lmpLayer", dVolTot_dPsi0 = dVolLiq_dPsi, dVolHtCapBulk_dPsi0 may be wrong
- select case(ixDomainType)
- case(iname_veg)
- if(computJac)then
- fLiq = fracLiquid(xTemp,snowfrz_scale)
- dVolHtCapBulk_dCanWat = ( -Cp_ice*( fLiq-1._rkind ) + Cp_water*fLiq )/canopyDepth !this is iden_water/(iden_water*canopyDepth)
- endif
- case(iname_snow)
- if(computJac)then
- fLiq = fracLiquid(xTemp,snowfrz_scale)
- dVolHtCapBulk_dTheta(iLayer) = iden_water * ( -Cp_ice*( fLiq-1._rkind ) + Cp_water*fLiq ) + iden_air * ( ( fLiq-1._rkind )*iden_water/iden_ice - fLiq ) * Cp_air
- endif
- case(iname_soil)
- select case( ixStateType(ixFullVector) )
- case(iname_lmpLayer)
- dVolTot_dPsi0(ixControlIndex) = dTheta_dPsi(mLayerMatricHeadLiqTrial(ixControlIndex),vGn_alpha(ixControlIndex),0._rkind,1._rkind,vGn_n(ixControlIndex),vGn_m(ixControlIndex))*avPore
- if(computJac) d2VolTot_d2Psi0(ixControlIndex) = d2Theta_dPsi2(mLayerMatricHeadLiqTrial(ixControlIndex),vGn_alpha(ixControlIndex),0._rkind,1._rkind,vGn_n(ixControlIndex),vGn_m(ixControlIndex))*avPore
- case default
- dVolTot_dPsi0(ixControlIndex) = dTheta_dPsi(mLayerMatricHeadTrial(ixControlIndex),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
- if(computJac) d2VolTot_d2Psi0(ixControlIndex) = d2Theta_dPsi2(mLayerMatricHeadTrial(ixControlIndex),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),&
- vGn_n(ixControlIndex),vGn_m(ixControlIndex))
- end select
- ! dVolHtCapBulk_dPsi0 uses the derivative in water retention curve above critical temp w.r.t.state variable dVolTot_dPsi0
- if(computJac)then
- dVolHtCapBulk_dTheta(iLayer) = realMissing ! do not use
- if(xTemp< Tcrit) dVolHtCapBulk_dPsi0(ixControlIndex) = (iden_ice * Cp_ice - iden_air * Cp_air) * dVolTot_dPsi0(ixControlIndex)
- if(xTemp>=Tcrit) dVolHtCapBulk_dPsi0(ixControlIndex) = (iden_water * Cp_water - iden_air * Cp_air) * dVolTot_dPsi0(ixControlIndex)
- endif
- end select
+ endif ! if errors in prognostic update
+
+ ! get the total energy fluxes (modified in updateProgSundials)
+ if(nrgFluxModified .or. indx_data%var(iLookINDEX%ixVegNrg)%dat(1)/=integerMissing)then
+ sumCanopyEvaporation = sumCanopyEvaporation + dt_out*flux_temp%var(iLookFLUX%scalarCanopyEvaporation)%dat(1) ! canopy evaporation/condensation (kg m-2 s-1)
+ sumLatHeatCanopyEvap = sumLatHeatCanopyEvap + dt_out*flux_temp%var(iLookFLUX%scalarLatHeatCanopyEvap)%dat(1) ! latent heat flux for evaporation from the canopy to the canopy air space (W m-2)
+ sumSenHeatCanopy = sumSenHeatCanopy + dt_out*flux_temp%var(iLookFLUX%scalarSenHeatCanopy)%dat(1) ! sensible heat flux from the canopy to the canopy air space (W m-2)
+ else
+ sumCanopyEvaporation = sumCanopyEvaporation + dt_out*flux_data%var(iLookFLUX%scalarCanopyEvaporation)%dat(1) ! canopy evaporation/condensation (kg m-2 s-1)
+ sumLatHeatCanopyEvap = sumLatHeatCanopyEvap + dt_out*flux_data%var(iLookFLUX%scalarLatHeatCanopyEvap)%dat(1) ! latent heat flux for evaporation from the canopy to the canopy air space (W m-2)
+ sumSenHeatCanopy = sumSenHeatCanopy + dt_out*flux_data%var(iLookFLUX%scalarSenHeatCanopy)%dat(1) ! sensible heat flux from the canopy to the canopy air space (W m-2)
+ endif ! if energy fluxes were modified
+
+ ! get the total soil compression
+ if (count(indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat/=integerMissing)>0) then
+ ! scalar compression
+ if(.not.scalarSolution .or. iStateSplit==nSoil)&
+ sumSoilCompress = sumSoilCompress + diag_data%var(iLookDIAG%scalarSoilCompress)%dat(1) ! total soil compression
+ ! vector compression
+ do iSoil=1,nSoil
+ if(indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat(iSoil)/=integerMissing)&
+ sumLayerCompress(iSoil) = sumLayerCompress(iSoil) + diag_data%var(iLookDIAG%mLayerCompress)%dat(iSoil) ! soil compression in layers
+ end do
+ endif
- ! compute the derivative in liquid water content w.r.t. temperature
- ! --> partially frozen: dependence of liquid water on temperature
- ! compute the derivative in bulk heat capacity w.r.t. temperature
- if(xTemp<Tcrit)then
- select case(ixDomainType)
- case(iname_veg)
- dFracLiqVeg_dTkCanopy = dFracLiq_dTk(xTemp,snowfrz_scale)
- dTheta_dTkCanopy = dFracLiqVeg_dTkCanopy * scalarCanopyWatTrial/(iden_water*canopyDepth)
- if(computJac)then
- fLiq = fracLiquid(xTemp,snowfrz_scale)
- d2Theta_dTkCanopy2 = 2._rkind * snowfrz_scale**2._rkind * ( (Tfreeze - xTemp) * 2._rkind * fLiq * dFracLiqVeg_dTkCanopy - fLiq**2._rkind ) * scalarCanopyWatTrial/(iden_water*canopyDepth)
- dVolHtCapBulk_dTkCanopy = iden_water * (-Cp_ice + Cp_water) * dTheta_dTkCanopy !same as snow but there is no derivative in air
- endif
- case(iname_snow)
- dFracLiqSnow_dTk(iLayer) = dFracLiq_dTk(xTemp,snowfrz_scale)
- mLayerdTheta_dTk(iLayer) = dFracLiqSnow_dTk(iLayer) * mLayerVolFracWatTrial(iLayer)
- if(computJac)then
- fLiq = fracLiquid(xTemp,snowfrz_scale)
- mLayerd2Theta_dTk2(iLayer) = 2._rkind * snowfrz_scale**2._rkind * ( (Tfreeze - xTemp) * 2._rkind * fLiq * dFracLiqSnow_dTk(iLayer) - fLiq**2._rkind ) * mLayerVolFracWatTrial(iLayer)
- dVolHtCapBulk_dTk(iLayer) = ( iden_water * (-Cp_ice + Cp_water) + iden_air * (iden_water/iden_ice - 1._rkind) * Cp_air ) * mLayerdTheta_dTk(iLayer)
- endif
- case(iname_soil)
- dFracLiqSnow_dTk(iLayer) = 0._rkind !dTheta_dTk(xTemp,theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))/ mLayerVolFracWatTrial(iLayer)
- mLayerdTheta_dTk(iLayer) = dTheta_dTk(xTemp,theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
- if(computJac)then
- mLayerd2Theta_dTk2(iLayer) = d2Theta_dTk2(xTemp,theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
- dVolHtCapBulk_dTk(iLayer) = (-iden_ice * Cp_ice + iden_water * Cp_water) * mLayerdTheta_dTk(iLayer)
+ ! print progress
+ if(globalPrintFlag)&
+ write(*,'(a,1x,3(f13.2,1x))') 'updating: dtSubstep, dtSum, dt = ', dtSubstep, dtSum, dt
+
+ ! increment fluxes
+ dt_wght = 1._qp !dt_out/dt ! (define weight applied to each splitting operation)
+ do iVar=1,size(flux_meta)
+ if(count(fluxMask%var(iVar)%dat)>0) then
+
+ !print*, flux_meta(iVar)%varname, fluxMask%var(iVar)%dat
+
+ ! ** no domain splitting
+ if(count(ixLayerActive/=integerMissing)==nLayers)then
+ flux_data%var(iVar)%dat(:) = flux_data%var(iVar)%dat(:) + flux_temp%var(iVar)%dat(:)*dt_wght
+ fluxCount%var(iVar)%dat(:) = fluxCount%var(iVar)%dat(:) + 1
+
+ ! ** domain splitting
+ else
+ ixMin=lbound(flux_data%var(iVar)%dat)
+ ixMax=ubound(flux_data%var(iVar)%dat)
+ do ixLayer=ixMin(1),ixMax(1)
+ if(fluxMask%var(iVar)%dat(ixLayer)) then
+ flux_data%var(iVar)%dat(ixLayer) = flux_data%var(iVar)%dat(ixLayer) + flux_temp%var(iVar)%dat(ixLayer)*dt_wght
+ fluxCount%var(iVar)%dat(ixLayer) = fluxCount%var(iVar)%dat(ixLayer) + 1
endif
- case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
- end select ! domain type
+ end do
+ endif ! (domain splitting)
- ! --> unfrozen: no dependence of liquid water on temperature
- else
- select case(ixDomainType)
- case(iname_veg); dTheta_dTkCanopy = 0._rkind; d2Theta_dTkCanopy2 = 0._rkind; dFracLiqVeg_dTkCanopy = 0._rkind; dVolHtCapBulk_dTkCanopy = 0._rkind
- case(iname_snow, iname_soil); mLayerdTheta_dTk(iLayer) = 0._rkind; mLayerd2Theta_dTk2(iLayer) = 0._rkind; dFracLiqSnow_dTk(iLayer) = 0._rkind; dVolHtCapBulk_dTk(iLayer) = 0._rkind
- end select ! domain type
- endif
+ endif ! (if the flux is desired)
+ end do ! (loop through fluxes)
- ! -----
- ! - update volumetric fraction of liquid water and ice...
- ! => case of hydrology state uncoupled with energy (and when not adjusting the temperature)...
- ! -----------------------------------------------------------------------------------------------
-
- ! case of hydrology state uncoupled with energy (and when not adjusting the temperature)
- if(.not.do_adjustTemp .and. .not.isNrgState .and. .not.isCoupled)then
-
- ! compute the fraction of snow
- select case(ixDomainType)
- case(iname_veg); scalarFracLiqVeg = fracliquid(xTemp,snowfrz_scale)
- case(iname_snow); mLayerFracLiqSnow(iLayer) = fracliquid(xTemp,snowfrz_scale)
- case(iname_soil) ! do nothing
- case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
- end select ! domain type
-
- ! -----
- ! - update volumetric fraction of liquid water and ice...
- ! => case of energy state or coupled solution (or adjusting the temperature)...
- ! --------------------------------------------------------------------------------
-
- ! case of energy state OR coupled solution (or adjusting the temperature)
- elseif(do_adjustTemp .or. ( (isNrgState .or. isCoupled) ) )then
-
- ! identify domain type
- select case(ixDomainType)
-
- ! *** vegetation canopy
- case(iname_veg)
-
- ! compute volumetric fraction of liquid water and ice
- call updateSnowSundials(&
- xTemp, & ! intent(in) : temperature (K)
- scalarCanopyWatTrial/(iden_water*canopyDepth),& ! intent(in) : volumetric fraction of total water (-)
- snowfrz_scale, & ! intent(in) : scaling parameter for the snow freezing curve (K-1)
- scalarCanopyTempPrime, & ! intent(in) : canopy temperature time derivative (K/s)
- scalarCanopyWatPrime/(iden_water*canopyDepth),& ! intent(in) : volumetric fraction of total water time derivative (-)
- scalarVolFracLiq, & ! intent(out) : trial canopy liquid water (-)
- scalarVolFracIce, & ! intent(out) : trial volumetric canopy ice (-)
- scalarVolFracLiqPrime, & ! intent(out) : trial volumetric canopy liquid water (-)
- scalarVolFracIcePrime, & ! intent(out) : trial volumetric canopy ice (-)
- scalarFracLiqVeg, & ! intent(out) : fraction of liquid water (-)
- err,cmessage) ! intent(out) : error control
- if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+ ! ------------------------------------------------------
+ ! ------------------------------------------------------
- ! compute mass of water on the canopy
- ! NOTE: possibilities for speed-up here
- scalarCanopyLiqTrial = scalarFracLiqVeg *scalarCanopyWatTrial !(kg m-2), scalarVolFracLiq*iden_water*canopyDepth
- scalarCanopyLiqPrime = scalarVolFracLiqPrime*iden_water*canopyDepth
- scalarCanopyIceTrial = (1._rkind - scalarFracLiqVeg)*scalarCanopyWatTrial !(kg m-2), scalarVolFracIce* iden_ice *canopyDepth
- scalarCanopyIcePrime = scalarVolFracIcePrime* iden_ice *canopyDepth
-
- ! *** snow layers
- case(iname_snow)
-
- ! compute volumetric fraction of liquid water and ice
- call updateSnowSundials(&
- xTemp, & ! intent(in) : temperature (K)
- mLayerVolFracWatTrial(iLayer), & ! intent(in) : mass state variable = trial volumetric fraction of water (-)
- snowfrz_scale, & ! intent(in) : scaling parameter for the snow freezing curve (K-1)
- mLayerTempPrime(iLayer), & ! intent(in) : temperature time derivative (K/s)
- mLayerVolFracWatPrime(iLayer), & ! intent(in) : volumetric fraction of total water time derivative (-)
- mLayerVolFracLiqTrial(iLayer), & ! intent(out) : trial volumetric fraction of liquid water (-)
- mLayerVolFracIceTrial(iLayer), & ! intent(out) : trial volumetric fraction if ice (-)
- mLayerVolFracLiqPrime(iLayer), & ! intent(out) : volumetric fraction of liquid water time derivative (-)
- mLayerVolFracIcePrime(iLayer), & ! intent(out) : volumetric fraction of ice time derivative (-)
- mLayerFracLiqSnow(iLayer), & ! intent(out) : fraction of liquid water (-)
- err,cmessage) ! intent(out) : error control
- if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+ ! increment the number of substeps
+ nSubsteps = nSubsteps+1
+
+ ! increment the sub-step legth
+ dtSum = dtSum + dtSubstep
+ !print*, 'dtSum, dtSubstep, dt, nSubsteps = ', dtSum, dtSubstep, dt, nSubsteps
- ! *** soil layers
- case(iname_soil)
-
- ! compute volumetric fraction of liquid water and ice
- call updateSoilSundials(&
- dt, & ! intent(in) : time step
- computJac, & ! intent(in) : logical flag if inside Sundials solver
- xTemp, & ! intent(in) : temperature (K)
- mLayerMatricHeadTrial(ixControlIndex), & ! intent(in) : total water matric potential (m)
- mLayerMatricHeadPrev(ixControlIndex), & ! intent(in) : previous values, will be same as current if computJac
- mLayerVolFracWatPrev(iLayer), & ! intent(in) : previous values, will be same as current if computJac
- mLayerTempPrime(iLayer), & ! intent(in) : temperature time derivative (K/s)
- mLayerMatricHeadPrime(ixControlIndex), & ! intent(in) : total water matric potential time derivative (m/s)
- vGn_alpha(ixControlIndex), & ! intent(in) : van Genutchen "alpha" parameter
- vGn_n(ixControlIndex), & ! intent(in) : van Genutchen "n" parameter
- theta_sat(ixControlIndex), & ! intent(in) : soil porosity (-)
- theta_res(ixControlIndex), & ! intent(in) : soil residual volumetric water content (-)
- vGn_m(ixControlIndex), & ! intent(in) : van Genutchen "m" parameter (-)
- mLayerVolFracWatTrial(iLayer), & ! intent(in) : mass state variable = trial volumetric fraction of water (-)
- mLayerVolFracLiqTrial(iLayer), & ! intent(out) : trial volumetric fraction of liquid water (-)
- mLayerVolFracIceTrial(iLayer), & ! intent(out) : trial volumetric fraction if ice (-)
- mLayerVolFracWatPrime(iLayer), & ! intent(out) : volumetric fraction of total water time derivative (-)
- mLayerVolFracLiqPrime(iLayer), & ! intent(out) : volumetric fraction of liquid water time derivative (-)
- mLayerVolFracIcePrime(iLayer), & ! intent(out) : volumetric fraction of ice time derivative (-)
- err,cmessage) ! intent(out) : error control
+ ! check that we have completed the sub-step
+ if(dtSum >= dt-verySmall)then
+ failedMinimumStep=.false.
+ exit subSteps
+ endif
+
+ ! adjust length of the sub-step (make sure that we don't exceed the step)
+ dtSubstep = min(dt - dtSum, max(dtSubstep*dtMultiplier, dt_min) )
+
+ end do substeps ! time steps for variable-dependent sub-stepping
+
+ ! save the energy fluxes
+ flux_data%var(iLookFLUX%scalarCanopyEvaporation)%dat(1) = sumCanopyEvaporation /dt_out ! canopy evaporation/condensation (kg m-2 s-1)
+ flux_data%var(iLookFLUX%scalarLatHeatCanopyEvap)%dat(1) = sumLatHeatCanopyEvap /dt_out ! latent heat flux for evaporation from the canopy to the canopy air space (W m-2)
+ flux_data%var(iLookFLUX%scalarSenHeatCanopy)%dat(1) = sumSenHeatCanopy /dt_out ! sensible heat flux from the canopy to the canopy air space (W m-2)
+
+ ! save the soil compression diagnostics
+ diag_data%var(iLookDIAG%scalarSoilCompress)%dat(1) = sumSoilCompress
+ do iSoil=1,nSoil
+ if(indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat(iSoil)/=integerMissing)&
+ diag_data%var(iLookDIAG%mLayerCompress)%dat(iSoil) = sumLayerCompress(iSoil)
+ end do
+ deallocate(sumLayerCompress)
+
+ ! end associate statements
+ end associate globalVars
+
+ ! update error codes
+ if(failedMinimumStep)then
+ err=-20 ! negative = recoverable error
+ message=trim(message)//'failed minimum step'
+ endif
+
+
+ end subroutine varSubstepSundials
+
+
+ ! **********************************************************************************************************
+ ! private subroutine updateProgSundials: update prognostic variables
+ ! **********************************************************************************************************
+ subroutine updateProgSundials(dt,nSnow,nSoil,nLayers,doAdjustTemp,computeVegFlux,untappedMelt,stateVecTrial,stateVecPrime,checkMassBalance, checkNrgBalance, & ! input: model control
+ lookup_data,mpar_data,indx_data,flux_data,prog_data,diag_data,deriv_data, & ! input-output: data structures
+ waterBalanceError,nrgFluxModified,err,message) ! output: flags and error control
+ USE getVectorz_module,only:varExtract ! extract variables from the state vector
+ USE updateVarsSundials_module,only:updateVarsSundials ! update prognostic variables
+ USE getVectorzAddSundials_module, only:varExtractSundials
+ USE computEnthalpy_module,only:computEnthalpy
+ USE t2enthalpy_module, only:t2enthalpy ! compute enthalpy
+ implicit none
+ ! model control
+ real(rkind) ,intent(in) :: dt ! time step (s)
+ 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
+ logical(lgt) ,intent(in) :: doAdjustTemp ! flag to indicate if we adjust the temperature
+ logical(lgt) ,intent(in) :: computeVegFlux ! flag to compute the vegetation flux
+ real(rkind) ,intent(in) :: untappedMelt(:) ! un-tapped melt energy (J m-3 s-1)
+ real(rkind) ,intent(in) :: stateVecTrial(:) ! trial state vector (mixed units)
+ real(rkind) ,intent(in) :: stateVecPrime(:) ! trial state vector (mixed units)
+ logical(lgt) ,intent(in) :: checkMassBalance ! flag to check the mass balance
+ logical(lgt) ,intent(in) :: checkNrgBalance ! flag to check the energy balance
+ ! data structures
+ type(zLookup),intent(in) :: lookup_data ! lookup tables
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ type(var_ilength),intent(in) :: indx_data ! indices for a local HRU
+ type(var_dlength),intent(inout) :: flux_data ! model fluxes 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) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
+ ! flags and error control
+ logical(lgt) ,intent(out) :: waterBalanceError ! flag to denote that there is a water balance error
+ logical(lgt) ,intent(out) :: nrgFluxModified ! flag to denote that the energy fluxes were modified
+ integer(i4b) ,intent(out) :: err ! error code
+ character(*) ,intent(out) :: message ! error message
+ ! ==================================================================================================================
+ ! general
+ integer(i4b) :: iState ! index of model state variable
+ integer(i4b) :: ixSubset ! index within the state subset
+ integer(i4b) :: ixFullVector ! index within full state vector
+ integer(i4b) :: ixControlIndex ! index within a given domain
+ real(rkind) :: volMelt ! volumetric melt (kg m-3)
+ real(rkind),parameter :: verySmall=epsilon(1._rkind)*2._rkind ! a very small number (deal with precision issues)
+ ! mass balance
+ real(rkind) :: canopyBalance0,canopyBalance1 ! canopy storage at start/end of time step
+ real(rkind) :: soilBalance0,soilBalance1 ! soil storage at start/end of time step
+ real(rkind) :: vertFlux ! change in storage due to vertical fluxes
+ real(rkind) :: tranSink,baseSink,compSink ! change in storage due to sink terms
+ real(rkind) :: liqError ! water balance error
+ real(rkind) :: fluxNet ! net water fluxes (kg m-2 s-1)
+ real(rkind) :: superflousWat ! superflous water used for evaporation (kg m-2 s-1)
+ real(rkind) :: superflousNrg ! superflous energy that cannot be used for evaporation (W m-2 [J m-2 s-1])
+ character(LEN=256) :: cmessage ! error message of downwind routine
+ ! trial state variables
+ real(rkind) :: scalarCanairTempTrial ! trial value for temperature of the canopy air space (K)
+ real(rkind) :: scalarCanopyTempTrial ! trial value for temperature of the vegetation canopy (K)
+ real(rkind) :: scalarCanopyWatTrial ! trial value for liquid water storage in the canopy (kg m-2)
+ real(rkind),dimension(nLayers) :: mLayerTempTrial ! trial vector for temperature of layers in the snow and soil domains (K)
+ real(rkind),dimension(nLayers) :: mLayerVolFracWatTrial ! trial vector for volumetric fraction of total water (-)
+ real(rkind),dimension(nSoil) :: mLayerMatricHeadTrial ! trial vector for total water matric potential (m)
+ real(rkind),dimension(nSoil) :: mLayerMatricHeadLiqTrial ! trial vector for liquid water matric potential (m)
+ real(rkind) :: scalarAquiferStorageTrial ! trial value for storage of water in the aquifer (m)
+ ! diagnostic variables
+ real(rkind) :: scalarCanopyLiqTrial ! trial value for mass of liquid water on the vegetation canopy (kg m-2)
+ real(rkind) :: scalarCanopyIceTrial ! trial value for mass of ice on the vegetation canopy (kg m-2)
+ real(rkind),dimension(nLayers) :: mLayerVolFracLiqTrial ! trial vector for volumetric fraction of liquid water (-)
+ real(rkind),dimension(nLayers) :: mLayerVolFracIceTrial ! trial vector for volumetric fraction of ice (-)
+ real(rkind) :: scalarCanairEnthalpyTrial ! enthalpy of the canopy air space (J m-3)
+ real(rkind) :: scalarCanopyEnthalpyTrial ! enthalpy of the vegetation canopy (J m-3)
+ real(rkind),dimension(nLayers) :: mLayerEnthalpyTrial ! enthalpy of snow + soil (J m-3)
+ ! derivative of state variables
+ real(rkind) :: scalarCanairTempPrime ! trial value for temperature of the canopy air space (K)
+ real(rkind) :: scalarCanopyTempPrime ! trial value for temperature of the vegetation canopy (K)
+ real(rkind) :: scalarCanopyWatPrime ! trial value for liquid water storage in the canopy (kg m-2)
+ real(rkind),dimension(nLayers) :: mLayerTempPrime ! trial vector for temperature of layers in the snow and soil domains (K)
+ real(rkind),dimension(nLayers) :: mLayerVolFracWatPrime ! trial vector for volumetric fraction of total water (-)
+ real(rkind),dimension(nSoil) :: mLayerMatricHeadPrime ! trial vector for total water matric potential (m)
+ real(rkind),dimension(nSoil) :: mLayerMatricHeadLiqPrime ! trial vector for liquid water matric potential (m)
+ real(rkind) :: scalarAquiferStoragePrime ! trial value for storage of water in the aquifer (m)
+ ! derivative of diagnostic variables
+ real(rkind) :: scalarCanopyLiqPrime ! trial value for mass of liquid water on the vegetation canopy (kg m-2)
+ real(rkind) :: scalarCanopyIcePrime ! trial value for mass of ice on the vegetation canopy (kg m-2)
+ real(rkind),dimension(nLayers) :: mLayerVolFracLiqPrime ! trial vector for volumetric fraction of liquid water (-)
+ real(rkind),dimension(nLayers) :: mLayerVolFracIcePrime ! trial vector for volumetric fraction of ice (-)
+ ! -------------------------------------------------------------------------------------------------------------------
+
+ ! -------------------------------------------------------------------------------------------------------------------
+ ! point to flux variables in the data structure
+ associate(&
+ ! get indices for mass balance
+ ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in) : [i4b] index of canopy hydrology state variable (mass)
+ ixSoilOnlyHyd => indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat ,& ! intent(in) : [i4b(:)] index in the state subset for hydrology state variables in the soil domain
+ ! get indices for the un-tapped melt
+ ixNrgOnly => indx_data%var(iLookINDEX%ixNrgOnly)%dat ,& ! intent(in) : [i4b(:)] list of indices for all energy states
+ ixDomainType => indx_data%var(iLookINDEX%ixDomainType)%dat ,& ! intent(in) : [i4b(:)] indices defining the domain of the state (iname_veg, iname_snow, iname_soil)
+ ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in) : [i4b(:)] index of the control volume for different domains (veg, snow, soil)
+ ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in) : [i4b(:)] [state subset] list of indices of the full state vector in the state subset
+ ! water fluxes
+ scalarRainfall => flux_data%var(iLookFLUX%scalarRainfall)%dat(1) ,& ! intent(in) : [dp] rainfall rate (kg m-2 s-1)
+ scalarThroughfallRain => flux_data%var(iLookFLUX%scalarThroughfallRain)%dat(1) ,& ! intent(in) : [dp] rain reaches ground without touching the canopy (kg m-2 s-1)
+ scalarCanopyEvaporation => flux_data%var(iLookFLUX%scalarCanopyEvaporation)%dat(1) ,& ! intent(in) : [dp] canopy evaporation/condensation (kg m-2 s-1)
+ scalarCanopyTranspiration => flux_data%var(iLookFLUX%scalarCanopyTranspiration)%dat(1) ,& ! intent(in) : [dp] canopy transpiration (kg m-2 s-1)
+ scalarCanopyLiqDrainage => flux_data%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) ,& ! intent(in) : [dp] drainage liquid water from vegetation canopy (kg m-2 s-1)
+ iLayerLiqFluxSoil => flux_data%var(iLookFLUX%iLayerLiqFluxSoil)%dat ,& ! intent(in) : [dp(0:)] vertical liquid water flux at soil layer interfaces (-)
+ iLayerNrgFlux => flux_data%var(iLookFLUX%iLayerNrgFlux)%dat ,& ! intent(in) :
+ mLayerNrgFlux => flux_data%var(iLookFLUX%mLayerNrgFlux)%dat ,& ! intent(out): [dp] net energy flux for each layer within the snow+soil domain (J m-3 s-1)
+ mLayerTranspire => flux_data%var(iLookFLUX%mLayerTranspire)%dat ,& ! intent(in) : [dp(:)] transpiration loss from each soil layer (m s-1)
+ mLayerBaseflow => flux_data%var(iLookFLUX%mLayerBaseflow)%dat ,& ! intent(in) : [dp(:)] baseflow from each soil layer (m s-1)
+ mLayerCompress => diag_data%var(iLookDIAG%mLayerCompress)%dat ,& ! intent(in) : [dp(:)] change in storage associated with compression of the soil matrix (-)
+ scalarCanopySublimation => flux_data%var(iLookFLUX%scalarCanopySublimation)%dat(1) ,& ! intent(in) : [dp] sublimation of ice from the vegetation canopy (kg m-2 s-1)
+ scalarSnowSublimation => flux_data%var(iLookFLUX%scalarSnowSublimation)%dat(1) ,& ! intent(in) : [dp] sublimation of ice from the snow surface (kg m-2 s-1)
+ ! energy fluxes
+ scalarLatHeatCanopyEvap => flux_data%var(iLookFLUX%scalarLatHeatCanopyEvap)%dat(1) ,& ! intent(in) : [dp] latent heat flux for evaporation from the canopy to the canopy air space (W m-2)
+ scalarSenHeatCanopy => flux_data%var(iLookFLUX%scalarSenHeatCanopy)%dat(1) ,& ! intent(in) : [dp] sensible heat flux from the canopy to the canopy air space (W m-2)
+ ! domain depth
+ 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)
+ ! model state variables (vegetation canopy)
+ scalarCanairTemp => prog_data%var(iLookPROG%scalarCanairTemp)%dat(1) ,& ! intent(inout) : [dp] temperature of the canopy air space (K)
+ scalarCanopyTemp => prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1) ,& ! intent(inout) : [dp] temperature of the vegetation canopy (K)
+ scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! intent(inout) : [dp] mass of ice on the vegetation canopy (kg m-2)
+ scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! intent(inout) : [dp] mass of liquid water on the vegetation canopy (kg m-2)
+ scalarCanopyWat => prog_data%var(iLookPROG%scalarCanopyWat)%dat(1) ,& ! intent(inout) : [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(inout) : [dp(:)] temperature of each snow/soil layer (K)
+ mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(inout) : [dp(:)] volumetric fraction of ice (-)
+ mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat ,& ! intent(inout) : [dp(:)] volumetric fraction of liquid water (-)
+ mLayerVolFracWat => prog_data%var(iLookPROG%mLayerVolFracWat)%dat ,& ! intent(inout) : [dp(:)] volumetric fraction of total water (-)
+ mLayerMatricHead => prog_data%var(iLookPROG%mLayerMatricHead)%dat ,& ! intent(inout) : [dp(:)] matric head (m)
+ mLayerMatricHeadLiq => diag_data%var(iLookDIAG%mLayerMatricHeadLiq)%dat ,& ! intent(inout) : [dp(:)] matric potential of liquid water (m)
+ ! enthalpy
+ scalarCanairEnthalpy => diag_data%var(iLookDIAG%scalarCanairEnthalpy)%dat(1) ,& ! intent(inout): [dp] enthalpy of the canopy air space (J m-3)
+ scalarCanopyEnthalpy => diag_data%var(iLookDIAG%scalarCanopyEnthalpy)%dat(1) ,& ! intent(inout): [dp] enthalpy of the vegetation canopy (J m-3)
+ mLayerEnthalpy => diag_data%var(iLookDIAG%mLayerEnthalpy)%dat ,& ! intent(inout): [dp(:)] enthalpy of the snow+soil layers (J m-3)
+ ! model state variables (aquifer)
+ scalarAquiferStorage => prog_data%var(iLookPROG%scalarAquiferStorage)%dat(1) ,& ! intent(inout) : [dp(:)] storage of water in the aquifer (m)
+ ! error tolerance
+ absConvTol_liquid => mpar_data%var(iLookPARAM%absConvTol_liquid)%dat(1) & ! intent(in) : [dp] absolute convergence tolerance for vol frac liq water (-)
+ ) ! associating flux variables in the data structure
+ ! -------------------------------------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message='updateProgSundials/'
+
+ ! initialize water balancmLayerVolFracWatTriale error
+ waterBalanceError=.false.
+
+ ! get storage at the start of the step
+ canopyBalance0 = merge(scalarCanopyWat, realMissing, computeVegFlux)
+ soilBalance0 = sum( (mLayerVolFracLiq(nSnow+1:nLayers) + mLayerVolFracIce(nSnow+1:nLayers) )*mLayerDepth(nSnow+1:nLayers) )
+
+ ! -----
+ ! * update states...
+ ! ------------------
+ ! these will need to be initialized as they do not have updated prognostic structures in Sundials
+ ! should all be set to previous values if splits, but for now operator splitting is not hooked up
+ scalarCanairTempPrime = realMissing
+ scalarCanopyTempPrime = realMissing
+ scalarCanopyWatPrime = realMissing
+ scalarCanopyLiqPrime = realMissing
+ scalarCanopyIcePrime = realMissing
+ mLayerTempPrime = realMissing
+ mLayerVolFracWatPrime = realMissing
+ mLayerVolFracLiqPrime = realMissing
+ mLayerVolFracIcePrime = realMissing
+ mLayerMatricHeadPrime = realMissing
+ mLayerMatricHeadLiqPrime = realMissing
+ scalarAquiferStoragePrime= realMissing
+ ! set to previous value from prognostic structure, correct because outside Sundials
+ scalarCanairTempTrial = scalarCanairTemp
+ scalarCanopyTempTrial = scalarCanopyTemp
+ scalarCanopyWatTrial = scalarCanopyWat
+ scalarCanopyLiqTrial = scalarCanopyLiq
+ scalarCanopyIceTrial = scalarCanopyIce
+ mLayerTempTrial = mLayerTemp
+ mLayerVolFracWatTrial = mLayerVolFracWat
+ mLayerVolFracLiqTrial = mLayerVolFracLiq
+ mLayerVolFracIceTrial = mLayerVolFracIce
+ mLayerMatricHeadTrial = mLayerMatricHead
+ mLayerMatricHeadLiqTrial = mLayerMatricHeadLiq
+ scalarAquiferStorageTrial= scalarAquiferStorage
+
+ ! extract variables from the model state vector
+ call varExtractSundials(&
+ ! input
+ stateVecTrial, & ! intent(in): model state vector (mixed units)
+ stateVecPrime, & ! intent(in): model state vector (mixed units)
+ diag_data, & ! intent(in): model diagnostic variables for a local HRU
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ indx_data, & ! intent(in): indices defining model states and layers
+ ! output: variables for the vegetation canopy
+ scalarCanairTempTrial, & ! intent(inout): trial value of canopy air temperature (K)
+ scalarCanopyTempTrial, & ! intent(inout): trial value of canopy temperature (K)
+ scalarCanopyWatTrial, & ! intent(inout): trial value of canopy total water (kg m-2)
+ scalarCanopyLiqTrial, & ! intent(inout): trial value of canopy liquid water (kg m-2)
+ scalarCanairTempPrime, & ! intent(inout): derivative of canopy air temperature (K)
+ scalarCanopyTempPrime, & ! intent(inout): derivative of canopy temperature (K)
+ scalarCanopyWatPrime, & ! intent(inout): derivative of canopy total water (kg m-2)
+ scalarCanopyLiqPrime, & ! intent(inout): derivative of canopy liquid water (kg m-2)
+ ! output: variables for the snow-soil domain
+ mLayerTempTrial, & ! intent(inout): trial vector of layer temperature (K)
+ mLayerVolFracWatTrial, & ! intent(inout): trial vector of volumetric total water content (-)
+ mLayerVolFracLiqTrial, & ! intent(inout): trial vector of volumetric liquid water content (-)
+ mLayerMatricHeadTrial, & ! intent(inout): trial vector of total water matric potential (m)
+ mLayerMatricHeadLiqTrial, & ! intent(inout): trial vector of liquid water matric potential (m)
+ mLayerTempPrime, & ! intent(inout): derivative of layer temperature (K)
+ mLayerVolFracWatPrime, & ! intent(inout): derivative of volumetric total water content (-)
+ mLayerVolFracLiqPrime, & ! intent(inout): derivative of volumetric liquid water content (-)
+ mLayerMatricHeadPrime, & ! intent(inout): derivative of total water matric potential (m)
+ mLayerMatricHeadLiqPrime, & ! intent(inout): derivative of liquid water matric potential (m)
+ ! output: variables for the aquifer
+ scalarAquiferStorageTrial,& ! intent(inout): trial value of storage of water in the aquifer (m)
+ scalarAquiferStoragePrime,& ! intent(inout): derivative of storage of water in the aquifer (m)
+ ! output: error control
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors)
+
+
+ ! update diagnostic variables
+ call updateVarsSundials(&
+ ! input
+ dt, &
+ .false., & ! intent(in): logical flag if computing Jacobian for Sundials solver
+ doAdjustTemp, & ! intent(in): logical flag to adjust temperature to account for the energy used in melt+freeze
+ mpar_data, & ! intent(in): model parameters for a local HRU
+ indx_data, & ! intent(in): indices defining model states and layers
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ mLayerVolFracWatTrial, & ! intent(in): use current vector for prev vector of volumetric total water content (-)
+ mLayerMatricHeadTrial, & ! intent(in): use current vector for prev vector of total water matric potential (m)
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables
+ ! output: variables for the vegetation canopy
+ scalarCanopyTempTrial, & ! intent(inout): trial value of canopy temperature (K)
+ scalarCanopyWatTrial, & ! intent(inout): trial value of canopy total water (kg m-2)
+ scalarCanopyLiqTrial, & ! intent(inout): trial value of canopy liquid water (kg m-2)
+ scalarCanopyIceTrial, & ! intent(inout): trial value of canopy ice content (kg m-2)
+ scalarCanopyTempPrime, & ! intent(inout): trial value of canopy temperature (K)
+ scalarCanopyWatPrime, & ! intent(inout): trial value of canopy total water (kg m-2)
+ scalarCanopyLiqPrime, & ! intent(inout): trial value of canopy liquid water (kg m-2)
+ scalarCanopyIcePrime, & ! intent(inout): trial value of canopy ice content (kg m-2)
+ ! output: variables for the snow-soil domain
+ mLayerTempTrial, & ! intent(inout): trial vector of layer temperature (K)
+ mLayerVolFracWatTrial, & ! intent(inout): trial vector of volumetric total water content (-)
+ mLayerVolFracLiqTrial, & ! intent(inout): trial vector of volumetric liquid water content (-)
+ mLayerVolFracIceTrial, & ! intent(inout): trial vector of volumetric ice water content (-)
+ mLayerMatricHeadTrial, & ! intent(inout): trial vector of total water matric potential (m)
+ mLayerMatricHeadLiqTrial, & ! intent(inout): trial vector of liquid water matric potential (m)
+ mLayerTempPrime, & !
+ mLayerVolFracWatPrime, & ! intent(inout): Prime vector of volumetric total water content (-)
+ mLayerVolFracLiqPrime, & ! intent(inout): Prime vector of volumetric liquid water content (-)
+ mLayerVolFracIcePrime, & !
+ mLayerMatricHeadPrime, & ! intent(inout): Prime vector of total water matric potential (m)
+ mLayerMatricHeadLiqPrime, & ! intent(inout): Prime vector of liquid water matric potential (m)
+ ! output: error control
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors)
+
+ ! ----
+ ! * check energy balance
+ !------------------------
+ ! NOTE: for now, we just compute enthalpy
+ if(checkNrgBalance)then
+ ! compute enthalpy at t_{n+1}
+ call t2enthalpy(&
+ ! input: data structures
+ diag_data, & ! intent(in): model diagnostic variables for a local HRU
+ mpar_data, & ! intent(in): parameter data structure
+ indx_data, & ! intent(in): model indices
+ lookup_data, & ! intent(in): lookup table data structure
+ ! input: state variables for the vegetation canopy
+ scalarCanairTempTrial, & ! intent(in): trial value of canopy air temperature (K)
+ scalarCanopyTempTrial, & ! intent(in): trial value of canopy temperature (K)
+ scalarCanopyWatTrial, & ! intent(in): trial value of canopy total water (kg m-2)
+ scalarCanopyIceTrial, & ! intent(in): trial value of canopy ice content (kg m-2)
+ ! input: variables for the snow-soil domain
+ mLayerTempTrial, & ! intent(in): trial vector of layer temperature (K)
+ mLayerVolFracWatTrial, & ! intent(in): trial vector of volumetric total water content (-)
+ mLayerMatricHeadTrial, & ! intent(in): trial vector of total water matric potential (m)
+ mLayerVolFracIceTrial, & ! intent(in): trial vector of volumetric fraction of ice (-)
+ ! output: enthalpy
+ scalarCanairEnthalpyTrial, & ! intent(out): enthalpy of the canopy air space (J m-3)
+ scalarCanopyEnthalpyTrial, & ! intent(out): enthalpy of the vegetation canopy (J m-3)
+ mLayerEnthalpyTrial, & ! intent(out): enthalpy of each snow+soil layer (J m-3)
+ ! output: error control
+ err,cmessage) ! intent(out): error control
if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
- ! check
- case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
+ endif
- end select ! domain type
+ ! -----
+ ! * check mass balance...
+ ! -----------------------
- ! final check
- else
+ ! NOTE: should not need to do this, since mass balance is checked in the solver
+ if(checkMassBalance)then
- ! do nothing (input = output) -- and check that we got here correctly
- if( (isNrgState .or. isCoupled) )then
- scalarVolFracLiq = realMissing
- scalarVolFracIce = realMissing
- else
- message=trim(message)//'unexpected else branch'
- err=20; return
+ ! check mass balance for the canopy
+ if(ixVegHyd/=integerMissing)then
+
+ ! handle cases where fluxes empty the canopy
+ fluxNet = scalarRainfall + scalarCanopyEvaporation - scalarThroughfallRain - scalarCanopyLiqDrainage
+ if(-fluxNet*dt > canopyBalance0)then
+
+ ! --> first add water
+ canopyBalance1 = canopyBalance0 + (scalarRainfall - scalarThroughfallRain)*dt
+
+ ! --> next, remove canopy evaporation -- put the unsatisfied evap into sensible heat
+ canopyBalance1 = canopyBalance1 + scalarCanopyEvaporation*dt
+ if(canopyBalance1 < 0._rkind)then
+ ! * get superfluous water and energy
+ superflousWat = -canopyBalance1/dt ! kg m-2 s-1
+ superflousNrg = superflousWat*LH_vap ! W m-2 (J m-2 s-1)
+ ! * update fluxes and states
+ canopyBalance1 = 0._rkind
+ scalarCanopyEvaporation = scalarCanopyEvaporation + superflousWat
+ scalarLatHeatCanopyEvap = scalarLatHeatCanopyEvap + superflousNrg
+ scalarSenHeatCanopy = scalarSenHeatCanopy - superflousNrg
endif
- endif ! if energy state or solution is coupled
-
- ! -----
- ! ------------------------
-
- ! check the need to adjust temperature (will always be false if inside solver)
- ! can be true if inside varSubstepSundials, outside solver, but currently will not work so turn off
- if(do_adjustTemp .and. computJac)then
-
- ! get the melt energy
- meltNrg = merge(LH_fus*iden_ice, LH_fus*iden_water, ixDomainType==iname_snow)
-
- ! compute the residual and the derivative
- select case(ixDomainType)
-
- ! * vegetation
- case(iname_veg)
- call xTempSolve(&
- ! constant over iterations
- meltNrg = meltNrg ,& ! intent(in) : energy for melt+freeze (J m-3)
- heatCap = scalarBulkVolHeatCapVeg ,& ! intent(in) : volumetric heat capacity (J m-3 K-1)
- tempInit = scalarCanopyTemp ,& ! intent(in) : initial temperature (K)
- volFracIceInit = scalarCanopyIce/(iden_water*canopyDepth),& ! intent(in) : initial volumetric fraction of ice (-)
- ! trial values
- xTemp = xTemp ,& ! intent(inout) : trial value of temperature
- dLiq_dT = dTheta_dTkCanopy ,& ! intent(in) : derivative in liquid water content w.r.t. temperature (K-1)
- volFracIceTrial = scalarVolFracIce ,& ! intent(in) : trial value for volumetric fraction of ice
- ! residual and derivative
- residual = residual ,& ! intent(out) : residual (J m-3)
- derivative = derivative ) ! intent(out) : derivative (J m-3 K-1)
-
- ! * snow and soil
- case(iname_snow, iname_soil)
- call xTempSolve(&
- ! constant over iterations
- meltNrg = meltNrg ,& ! intent(in) : energy for melt+freeze (J m-3)
- heatCap = mLayerVolHtCapBulk(iLayer) ,& ! intent(in) : volumetric heat capacity (J m-3 K-1)
- tempInit = mLayerTemp(iLayer) ,& ! intent(in) : initial temperature (K)
- volFracIceInit = mLayerVolFracIce(iLayer) ,& ! intent(in) : initial volumetric fraction of ice (-)
- ! trial values
- xTemp = xTemp ,& ! intent(inout) : trial value of temperature
- dLiq_dT = mLayerdTheta_dTk(iLayer) ,& ! intent(in) : derivative in liquid water content w.r.t. temperature (K-1)
- volFracIceTrial = mLayerVolFracIceTrial(iLayer) ,& ! intent(in) : trial value for volumetric fraction of ice
- ! residual and derivative
- residual = residual ,& ! intent(out) : residual (J m-3)
- derivative = derivative ) ! intent(out) : derivative (J m-3 K-1)
-
- ! * check
- case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
-
- end select ! domain type
-
- ! check validity of residual
- if( ieee_is_nan(residual) )then
- message=trim(message)//'residual is not valid'
- err=20; return
+ ! --> next, remove canopy drainage
+ canopyBalance1 = canopyBalance1 - scalarCanopyLiqDrainage*dt
+ if(canopyBalance1 < 0._rkind)then
+ superflousWat = -canopyBalance1/dt ! kg m-2 s-1
+ canopyBalance1 = 0._rkind
+ scalarCanopyLiqDrainage = scalarCanopyLiqDrainage + superflousWat
endif
- ! update bracket
- if(residual < 0._rkind)then
- tempMax = min(xTemp,tempMax)
+ ! update the trial state
+ scalarCanopyWatTrial = canopyBalance1
+
+ ! set the modification flag
+ nrgFluxModified = .true.
+
+ else
+ canopyBalance1 = canopyBalance0 + fluxNet*dt
+ nrgFluxModified = .false.
+ endif ! cases where fluxes empty the canopy
+
+ ! check the mass balance
+ fluxNet = scalarRainfall + scalarCanopyEvaporation - scalarThroughfallRain - scalarCanopyLiqDrainage
+ liqError = (canopyBalance0 + fluxNet*dt) - scalarCanopyWatTrial
+ if(abs(liqError) > absConvTol_liquid*10._rkind*iden_water)then ! *10 because of precision issues
+ !write(*,'(a,1x,f20.10)') 'dt = ', dt
+ !write(*,'(a,1x,f20.10)') 'scalarCanopyWatTrial = ', scalarCanopyWatTrial
+ !write(*,'(a,1x,f20.10)') 'canopyBalance0 = ', canopyBalance0
+ !write(*,'(a,1x,f20.10)') 'canopyBalance1 = ', canopyBalance1
+ !write(*,'(a,1x,f20.10)') 'scalarRainfall*dt = ', scalarRainfall*dt
+ !write(*,'(a,1x,f20.10)') 'scalarCanopyLiqDrainage*dt = ', scalarCanopyLiqDrainage*dt
+ !write(*,'(a,1x,f20.10)') 'scalarCanopyEvaporation*dt = ', scalarCanopyEvaporation*dt
+ !write(*,'(a,1x,f20.10)') 'scalarThroughfallRain*dt = ', scalarThroughfallRain*dt
+ !write(*,'(a,1x,f20.10)') 'liqError = ', liqError
+ waterBalanceError = .true.
+ return
+ endif ! if there is a water balance error
+ endif ! if veg canopy
+
+ ! check mass balance for soil
+ ! NOTE: fatal errors, though possible to recover using negative error codes
+ if(count(ixSoilOnlyHyd/=integerMissing)==nSoil)then
+ soilBalance1 = sum( (mLayerVolFracLiqTrial(nSnow+1:nLayers) + mLayerVolFracIceTrial(nSnow+1:nLayers) )*mLayerDepth(nSnow+1:nLayers) )
+ vertFlux = -(iLayerLiqFluxSoil(nSoil) - iLayerLiqFluxSoil(0))*dt ! m s-1 --> m
+ tranSink = sum(mLayerTranspire)*dt ! m s-1 --> m
+ baseSink = sum(mLayerBaseflow)*dt ! m s-1 --> m
+ compSink = sum(mLayerCompress(1:nSoil) * mLayerDepth(nSnow+1:nLayers) ) ! dimensionless --> m
+ liqError = soilBalance1 - (soilBalance0 + vertFlux + tranSink - baseSink - compSink)
+ if(abs(liqError) > absConvTol_liquid*10._rkind)then ! *10 because of precision issues
+ !write(*,'(a,1x,f20.10)') 'dt = ', dt
+ !write(*,'(a,1x,f20.10)') 'soilBalance0 = ', soilBalance0
+ !write(*,'(a,1x,f20.10)') 'soilBalance1 = ', soilBalance1
+ !write(*,'(a,1x,f20.10)') 'vertFlux = ', vertFlux
+ !write(*,'(a,1x,f20.10)') 'tranSink = ', tranSink
+ !write(*,'(a,1x,f20.10)') 'baseSink = ', baseSink
+ !write(*,'(a,1x,f20.10)') 'compSink = ', compSink
+ !write(*,'(a,1x,f20.10)') 'liqError = ', liqError
+ !write(*,'(a,1x,f20.10)') 'absConvTol_liquid = ', absConvTol_liquid
+ waterBalanceError = .true.
+ return
+ endif ! if there is a water balance error
+ endif ! if hydrology states exist in the soil domain
+
+ endif ! if checking the mass balance
+
+ ! -----
+ ! * remove untapped melt energy...
+ ! --------------------------------
+
+ ! only work with energy state variables
+ if(size(ixNrgOnly)>0)then ! energy state variables exist
+
+ ! loop through energy state variables
+ do iState=1,size(ixNrgOnly)
+
+ ! get index of the control volume within the domain
+ ixSubset = ixNrgOnly(iState) ! index within the state subset
+ ixFullVector = ixMapSubset2Full(ixSubset) ! index within full state vector
+ ixControlIndex = ixControlVolume(ixFullVector) ! index within a given domain
+
+ ! compute volumetric melt (kg m-3)
+ volMelt = dt*untappedMelt(ixSubset)/LH_fus ! (kg m-3)
+
+ ! update ice content
+ select case( ixDomainType(ixFullVector) )
+ case(iname_cas); cycle ! do nothing, since there is no snow stored in the canopy air space
+ case(iname_veg); scalarCanopyIceTrial = scalarCanopyIceTrial - volMelt*canopyDepth ! (kg m-2)
+ case(iname_snow); mLayerVolFracIceTrial(ixControlIndex) = mLayerVolFracIceTrial(ixControlIndex) - volMelt/iden_ice ! (-)
+ case(iname_soil); mLayerVolFracIceTrial(ixControlIndex+nSnow) = mLayerVolFracIceTrial(ixControlIndex+nSnow) - volMelt/iden_water ! (-)
+ case default; err=20; message=trim(message)//'unable to identify domain type [remove untapped melt energy]'; return
+ end select
+
+ ! update liquid water content
+ select case( ixDomainType(ixFullVector) )
+ case(iname_cas); cycle ! do nothing, since there is no snow stored in the canopy air space
+ case(iname_veg); scalarCanopyLiqTrial = scalarCanopyLiqTrial + volMelt*canopyDepth ! (kg m-2)
+ case(iname_snow); mLayerVolFracLiqTrial(ixControlIndex) = mLayerVolFracLiqTrial(ixControlIndex) + volMelt/iden_water ! (-)
+ case(iname_soil); mLayerVolFracLiqTrial(ixControlIndex+nSnow) = mLayerVolFracLiqTrial(ixControlIndex+nSnow) + volMelt/iden_water ! (-)
+ case default; err=20; message=trim(message)//'unable to identify domain type [remove untapped melt energy]'; return
+ end select
+
+ end do ! looping through energy variables
+
+ ! ========================================================================================================
+
+ ! *** ice
+
+ ! --> check if we removed too much water
+ if(scalarCanopyIceTrial < 0._rkind .or. any(mLayerVolFracIceTrial < 0._rkind) )then
+
+ ! **
+ ! canopy within numerical precision
+ if(scalarCanopyIceTrial < 0._rkind)then
+
+ if(scalarCanopyIceTrial > -verySmall)then
+ scalarCanopyLiqTrial = scalarCanopyLiqTrial - scalarCanopyIceTrial
+ scalarCanopyIceTrial = 0._rkind
+
+ ! encountered an inconsistency: spit the dummy
else
- tempMin = max(tempMin,xTemp)
- end if
-
- ! compute iteration increment
- tempInc = residual/derivative ! K
-
- ! check
- if(globalPrintFlag)&
- write(*,'(i4,1x,e20.10,1x,5(f20.10,1x),L1)') iter, residual, xTemp-Tcrit, tempInc, Tcrit, tempMin, tempMax, bFlag
-
- ! check convergence
- if(abs(residual) < nrgConvTol .or. abs(tempInc) < tempConvTol) exit iterations
-
- ! add constraints for snow temperature
- if(ixDomainType==iname_veg .or. ixDomainType==iname_snow)then
- if(tempInc > Tcrit - xTemp) tempInc=(Tcrit - xTemp)*0.5_rkind ! simple bi-section method
- endif ! if the domain is vegetation or snow
-
- ! deal with the discontinuity between partially frozen and unfrozen soil
- if(ixDomainType==iname_soil)then
- ! difference from the temperature below which ice exists
- critDiff = Tcrit - xTemp
- ! --> initially frozen (T < Tcrit)
- if(critDiff > 0._rkind)then
- if(tempInc > critDiff) tempInc = critDiff + epsT ! set iteration increment to slightly above critical temperature
- ! --> initially unfrozen (T > Tcrit)
+ print*, 'dt = ', dt
+ print*, 'untappedMelt = ', untappedMelt
+ print*, 'untappedMelt*dt = ', untappedMelt*dt
+ print*, 'scalarCanopyiceTrial = ', scalarCanopyIceTrial
+ message=trim(message)//'melted more than the available water'
+ err=20; return
+ endif ! (inconsistency)
+
+ endif ! if checking the canopy
+ ! **
+ ! snow+soil within numerical precision
+ do iState=1,size(mLayerVolFracIceTrial)
+
+ ! snow layer within numerical precision
+ if(mLayerVolFracIceTrial(iState) < 0._rkind)then
+
+ if(mLayerVolFracIceTrial(iState) > -verySmall)then
+ mLayerVolFracLiqTrial(iState) = mLayerVolFracLiqTrial(iState) - mLayerVolFracIceTrial(iState)
+ mLayerVolFracIceTrial(iState) = 0._rkind
+
+ ! encountered an inconsistency: spit the dummy
else
- if(tempInc < critDiff) tempInc = critDiff - epsT ! set iteration increment to slightly below critical temperature
- endif
- endif ! if the domain is soil
+ print*, 'dt = ', dt
+ print*, 'untappedMelt = ', untappedMelt
+ print*, 'untappedMelt*dt = ', untappedMelt*dt
+ print*, 'mLayerVolFracIceTrial = ', mLayerVolFracIceTrial
+ message=trim(message)//'melted more than the available water'
+ err=20; return
+ endif ! (inconsistency)
- ! update the temperature trial
- xTemp = xTemp + tempInc
+ endif ! if checking a snow layer
- ! check failed convergence
- if(iter==maxiter)then
- message=trim(message)//'failed to converge'
- err=-20; return ! negative error code = try to recover
- endif
+ end do ! (looping through state variables)
- endif ! if adjusting the temperature
+ endif ! (if we removed too much water)
- end do iterations ! iterating
+ ! ========================================================================================================
- ! save temperature
- select case(ixDomainType)
- case(iname_veg); scalarCanopyTempTrial = xTemp
- case(iname_snow, iname_soil); mLayerTempTrial(iLayer) = xTemp
- end select
+ ! *** liquid water
- ! =======================================================================================================================================
- ! =======================================================================================================================================
+ ! --> check if we removed too much water
+ if(scalarCanopyLiqTrial < 0._rkind .or. any(mLayerVolFracLiqTrial < 0._rkind) )then
- ! -----
- ! - compute the liquid water matric potential (and necessay derivatives)...
- ! -------------------------------------------------------------------------
+ ! **
+ ! canopy within numerical precision
+ if(scalarCanopyLiqTrial < 0._rkind)then
- ! only for soil
- if(ixDomainType==iname_soil)then
+ if(scalarCanopyLiqTrial > -verySmall)then
+ scalarCanopyIceTrial = scalarCanopyIceTrial - scalarCanopyLiqTrial
+ scalarCanopyLiqTrial = 0._rkind
- ! check liquid water
- if(mLayerVolFracLiqTrial(iLayer) > theta_sat(ixControlIndex) )then
- message=trim(message)//'liquid water greater than porosity'
- err=20; return
- endif
- ! case of hydrology state uncoupled with energy
- if(.not.isNrgState .and. .not.isCoupled)then
+ ! encountered an inconsistency: spit the dummy
+ else
+ print*, 'dt = ', dt
+ print*, 'untappedMelt = ', untappedMelt
+ print*, 'untappedMelt*dt = ', untappedMelt*dt
+ print*, 'scalarCanopyLiqTrial = ', scalarCanopyLiqTrial
+ message=trim(message)//'frozen more than the available water'
+ err=20; return
+ endif ! (inconsistency)
- ! derivatives relating liquid water matric potential to total water matric potential and temperature
- dPsiLiq_dPsi0(ixControlIndex) = 1._rkind ! exact correspondence (psiLiq=psi0)
- dPsiLiq_dTemp(ixControlIndex) = 0._rkind ! no relationship between liquid water matric potential and temperature
+ endif ! checking the canopy
- ! case of energy state or coupled solution
- else
- ! compute the liquid matric potential (and the derivatives w.r.t. total matric potential and temperature)
- call liquidHeadSundials(&
- ! input
- mLayerMatricHeadTrial(ixControlIndex) ,& ! intent(in) : total water matric potential (m)
- mLayerMatricHeadPrime(ixControlIndex) ,& !
- mLayerVolFracLiqTrial(iLayer) ,& ! intent(in) : volumetric fraction of liquid water (-)
- mLayerVolFracIceTrial(iLayer) ,& ! intent(in) : volumetric fraction of ice (-)
- vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),theta_sat(ixControlIndex),theta_res(ixControlIndex),vGn_m(ixControlIndex), & ! intent(in) : soil parameters
- dVolTot_dPsi0(ixControlIndex) ,& ! intent(in) : derivative in the soil water characteristic (m-1)
- mLayerdTheta_dTk(iLayer) ,& ! intent(in) : derivative in volumetric total water w.r.t. temperature (K-1)
- mLayerTempPrime(ixControlIndex) ,&
- mLayerVolFracLiqPrime(iLayer) ,&
- mLayerVolFracIcePrime(iLayer) ,&
- ! output
- mLayerMatricHeadLiqTrial(ixControlIndex) ,& ! intent(out): liquid water matric potential (m)
- mLayerMatricHeadLiqPrime(ixControlIndex) ,& !
- dPsiLiq_dPsi0(ixControlIndex) ,& ! intent(out): derivative in the liquid water matric potential w.r.t. the total water matric potential (-)
- dPsiLiq_dTemp(ixControlIndex) ,& ! intent(out): derivative in the liquid water matric potential w.r.t. temperature (m K-1)
- err,cmessage) ! intent(out): error control
- if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
-
- endif ! switch between hydrology and energy state
-
- endif ! if domain is soil
-
- end do ! looping through state variables
-
- ! deallocate space
- deallocate(computedCoupling,stat=err) ! .true. if computed the coupling for a given state variable
- if(err/=0)then; message=trim(message)//'problem deallocating computedCoupling'; return; endif
-
- ! end association to the variables in the data structures
- end associate
+ ! **
+ ! snow+soil within numerical precision
+ do iState=1,size(mLayerVolFracLiqTrial)
- end subroutine updateVarsSundials
+ ! snow layer within numerical precision
+ if(mLayerVolFracLiqTrial(iState) < 0._rkind)then
+ if(mLayerVolFracLiqTrial(iState) > -verySmall)then
+ mLayerVolFracIceTrial(iState) = mLayerVolFracIceTrial(iState) - mLayerVolFracLiqTrial(iState)
+ mLayerVolFracLiqTrial(iState) = 0._rkind
- ! **********************************************************************************************************
- ! private subroutine xTempSolve: compute residual and derivative for temperature
- ! **********************************************************************************************************
- subroutine xTempSolve(&
- ! input: constant over iterations
- meltNrg ,& ! intent(in) : energy for melt+freeze (J m-3)
- heatCap ,& ! intent(in) : volumetric heat capacity (J m-3 K-1)
- tempInit ,& ! intent(in) : initial temperature (K)
- volFracIceInit ,& ! intent(in) : initial volumetric fraction of ice (-)
- ! input-output: trial values
- xTemp ,& ! intent(inout) : trial value of temperature
- dLiq_dT ,& ! intent(in) : derivative in liquid water content w.r.t. temperature (K-1)
- volFracIceTrial ,& ! intent(in) : trial value for volumetric fraction of ice
- ! output: residual and derivative
- residual ,& ! intent(out) : residual (J m-3)
- derivative ) ! intent(out) : derivative (J m-3 K-1)
- implicit none
- ! input: constant over iterations
- real(rkind),intent(in) :: meltNrg ! energy for melt+freeze (J m-3)
- real(rkind),intent(in) :: heatCap ! volumetric heat capacity (J m-3 K-1)
- real(rkind),intent(in) :: tempInit ! initial temperature (K)
- real(rkind),intent(in) :: volFracIceInit ! initial volumetric fraction of ice (-)
- ! input-output: trial values
- real(rkind),intent(inout) :: xTemp ! trial value for temperature
- real(rkind),intent(in) :: dLiq_dT ! derivative in liquid water content w.r.t. temperature (K-1)
- real(rkind),intent(in) :: volFracIceTrial ! trial value for the volumetric fraction of ice (-)
- ! output: residual and derivative
- real(rkind),intent(out) :: residual ! residual (J m-3)
- real(rkind),intent(out) :: derivative ! derivative (J m-3 K-1)
- ! subroutine starts here
- residual = -heatCap*(xTemp - tempInit) + meltNrg*(volFracIceTrial - volFracIceInit) ! J m-3
- derivative = heatCap + LH_fus*iden_water*dLiq_dT ! J m-3 K-1
- end subroutine xTempSolve
-
- end module updateVarsSundials_module
+ ! encountered an inconsistency: spit the dummy
+ else
+ print*, 'dt = ', dt
+ print*, 'untappedMelt = ', untappedMelt
+ print*, 'untappedMelt*dt = ', untappedMelt*dt
+ print*, 'mLayerVolFracLiqTrial = ', mLayerVolFracLiqTrial
+ message=trim(message)//'frozen more than the available water'
+ err=20; return
+ endif ! (inconsistency)
+
+ endif ! checking a snow layer
+
+ end do ! (looping through state variables)
+
+ endif ! (if we removed too much water)
+
+ endif ! (if energy state variables exist)
+
+ ! -----
+ ! * update enthalpy as a diagnostic variable...
+ ! --------------------------------
+ scalarCanairEnthalpy = scalarCanairEnthalpyTrial
+ scalarCanopyEnthalpy = scalarCanopyEnthalpyTrial
+ mLayerEnthalpy = mLayerEnthalpyTrial
+
+ ! -----
+ ! * update prognostic variables...
+ ! --------------------------------
+ ! update state variables for the vegetation canopy
+ scalarCanairTemp = scalarCanairTempTrial ! trial value of canopy air temperature (K)
+ scalarCanopyTemp = scalarCanopyTempTrial ! trial value of canopy temperature (K)
+ scalarCanopyWat = scalarCanopyWatTrial ! trial value of canopy total water (kg m-2)
+ scalarCanopyLiq = scalarCanopyLiqTrial ! trial value of canopy liquid water (kg m-2)
+ scalarCanopyIce = scalarCanopyIceTrial ! trial value of canopy ice content (kg m-2)
+
+ ! update state variables for the snow+soil domain
+ mLayerTemp = mLayerTempTrial ! trial vector of layer temperature (K)
+ mLayerVolFracWat = mLayerVolFracWatTrial ! trial vector of volumetric total water content (-)
+ mLayerVolFracLiq = mLayerVolFracLiqTrial ! trial vector of volumetric liquid water content (-)
+ mLayerVolFracIce = mLayerVolFracIceTrial ! trial vector of volumetric ice water content (-)
+ mLayerMatricHead = mLayerMatricHeadTrial ! trial vector of matric head (m)
+ mLayerMatricHeadLiq = mLayerMatricHeadLiqTrial ! trial vector of matric head (m)
+
+ ! update state variables for the aquifer
+ scalarAquiferStorage = scalarAquiferStorageTrial
+
+ ! end associations to info in the data structures
+ end associate
+
+ end subroutine updateProgSundials
+
+ end module varSubstepSundials_module
\ No newline at end of file