From 0de7f3197d96a76fc183addf3dfd5462cb002c0c Mon Sep 17 00:00:00 2001
From: Kyle <kyle.c.klenk@gmail.com>
Date: Fri, 9 Sep 2022 19:21:16 +0000
Subject: [PATCH] applied changes to updateVarsSundials

---
 .../engine/sundials/updateVarsSundials.f90    | 1449 ++++++++++-------
 1 file changed, 825 insertions(+), 624 deletions(-)

diff --git a/build/source/engine/sundials/updateVarsSundials.f90 b/build/source/engine/sundials/updateVarsSundials.f90
index 8e5f9fb..0f54d15 100644
--- a/build/source/engine/sundials/updateVarsSundials.f90
+++ b/build/source/engine/sundials/updateVarsSundials.f90
@@ -20,628 +20,829 @@
 
 module updateVarsSundials_module
 
-! data types
-USE nrtype
-
-! missing values
-USE globalData,only:integerMissing  ! missing integer
-USE globalData,only:realMissing     ! missing real 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_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)
-                    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)
-
-! provide access to the derived types to define the data structures
-USE data_types,only:&
-                    var_i,        & ! data vector (i4b)
-                    var_d,        & ! data vector (rkind)
-                    var_ilength,  & ! data vector with variable length dimension (i4b)
-                    var_dlength     ! data vector with variable length dimension (rkind)
-
-! 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:updateVegSundials     ! update snow 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_utilsSundials_module,only:liquidHeadSundials     ! compute the liquid water matric potential
-
-! IEEE checks
-USE, intrinsic :: ieee_arithmetic            ! check values (NaN, etc.)
-
-implicit none
-private
-public::updateVarsSundials
-
-contains
-
- ! **********************************************************************************************************
- ! public subroutine updateVarsSundials: compute diagnostic variables
- ! **********************************************************************************************************
- subroutine updateVarsSundials(&
+   ! data types
+   USE nrtype
+   
+   ! missing values
+   USE globalData,only:integerMissing  ! missing integer
+   USE globalData,only:realMissing     ! missing real 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_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)
+   
+   ! provide access to the derived types to define the data structures
+   USE data_types,only:&
+                       var_i,        & ! data vector (i4b)
+                       var_d,        & ! data vector (rkind)
+                       var_ilength,  & ! data vector with variable length dimension (i4b)
+                       var_dlength     ! data vector with variable length dimension (rkind)
+   
+   ! 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.)
+   
+   implicit none
+   private
+   public::updateVarsSundials
+   
+   contains
+   
+    ! **********************************************************************************************************
+    ! public subroutine updateVarsSundials: compute diagnostic variables and derivatives for Sundials Jacobian
+    ! **********************************************************************************************************
+    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
+    ! --------------------------------------------------------------------------------------------------------------------------------
+    ! --------------------------------------------------------------------------------------------------------------------------------
+    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
+    ! --------------------------------------------------------------------------------------------------------------------------------
+    ! 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
+   
+    ! --------------------------------------------------------------------------------------------------------------------------------
+    ! --------------------------------------------------------------------------------------------------------------------------------
+   
+    ! initialize error control
+    err=0; message='updateVarsSundials/'
+   
+    ! 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.
+   
+    ! loop through model state variables
+    do iState=1,size(ixMapSubset2Full)
+   
+     ! check the need for the computations
+     if(computedCoupling(iState)) cycle
+   
+     ! -----
+     ! - 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
+     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
+   
+      ! 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
+       bFlag = .false.
+      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
+   
+      ! 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)
+         endif
+        case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
+       end select  ! domain type
+   
+      ! --> 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
+   
+      ! -----
+      ! - 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
+   
+        ! *** 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
+         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
+   
+       end select  ! domain type
+   
+      ! final check
+      else
+   
+       ! 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
+       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
+       endif
+   
+       ! update bracket
+       if(residual < 0._rkind)then
+        tempMax = min(xTemp,tempMax)
+       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)
+        else
+         if(tempInc < critDiff) tempInc = critDiff - epsT  ! set iteration increment to slightly below critical temperature
+        endif
+       endif  ! if the domain is soil
+   
+       ! update the temperature trial
+       xTemp = xTemp + tempInc
+   
+       ! check failed convergence
+       if(iter==maxiter)then
+        message=trim(message)//'failed to converge'
+        err=-20; return ! negative error code = try to recover
+       endif
+   
+      endif   ! if adjusting the temperature
+   
+     end do iterations ! iterating
+   
+     ! save temperature
+     select case(ixDomainType)
+      case(iname_veg);              scalarCanopyTempTrial   = xTemp
+      case(iname_snow, iname_soil); mLayerTempTrial(iLayer) = xTemp
+     end select
+   
+     ! =======================================================================================================================================
+     ! =======================================================================================================================================
+   
+     ! -----
+     ! - compute the liquid water matric potential (and necessay derivatives)...
+     ! -------------------------------------------------------------------------
+   
+     ! only for soil
+     if(ixDomainType==iname_soil)then
+   
+      ! 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
+   
+       ! 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
+   
+      ! 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
-                       dt_cur,                                    &
-                       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)
-                       mLayerMatricHeadPrev,                      & ! intent(in)
-                       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,                           & ! reza
-                       mLayerVolFracWatPrime,                     & ! reza
-                       mLayerVolFracLiqPrime,                     & ! reza
-                       mLayerVolFracIcePrime,                     & ! reza
-                       mLayerMatricHeadPrime,                     & ! reza
-                       mLayerMatricHeadLiqPrime,                  & ! reza
-                       ! output: error control
-                       err,message)                                 ! intent(out):   error control
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------------------------------------------------------
- implicit none
- ! input
- real(rkind),intent(in)             :: dt_cur
- 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(:)
- real(rkind),intent(in)             :: mLayerMatricHeadPrev(:)
- 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 canopy temperature (K)
- real(rkind),intent(inout)          :: scalarCanopyWatPrime            ! trial value of canopy total water (kg m-2)
- real(rkind),intent(inout)          :: scalarCanopyLiqPrime            ! trial value of canopy liquid water (kg m-2)
- real(rkind),intent(inout)          :: scalarCanopyIcePrime            ! trial value of 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(:)
- real(rkind),intent(inout)          :: mLayerVolFracWatPrime(:)        ! reza
- real(rkind),intent(inout)          :: mLayerVolFracLiqPrime(:)        ! reza
- real(rkind),intent(inout)          :: mLayerVolFracIcePrime(:)        ! reza
- real(rkind),intent(inout)          :: mLayerMatricHeadPrime(:)        ! reza
- real(rkind),intent(inout)          :: mLayerMatricHeadLiqPrime(:)     ! reza
-
- ! output: error control
- integer(i4b),intent(out)        :: err                             ! error code
- character(*),intent(out)        :: message                         ! error message
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! 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)                        :: Tcrit                           ! critical soil temperature below which ice exists (K)
- real(rkind)                        :: xTemp                           ! temporary temperature (K)
- real(rkind)                        :: effSat                          ! effective saturation (-)
- real(rkind)                        :: avPore=-9999                    ! 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
- 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)                        :: 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
- ) ! association with variables in the data structures
-
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------------------------------------------------------
-
- ! initialize error control
- err=0; message='updateVarsSundials/'
-
- ! 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.
-
- ! loop through model state variables
- do iState=1,size(ixMapSubset2Full)
-
-  ! check the need for the computations
-  if(computedCoupling(iState)) cycle
-
-  ! -----
-  ! - 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
-  endif
-
-
-
-  ! =======================================================================================================================================
-  ! =======================================================================================================================================
-  ! =======================================================================================================================================
-  ! =======================================================================================================================================
-  ! =======================================================================================================================================
-  ! =======================================================================================================================================
-
-  ! update hydrology state variables for the uncoupled solution
-  if(.not.isNrgState .and. .not.isCoupled)then
-
-  stop 1
-
-   ! 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
-
-   ! 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
-    bFlag = .false.
-   endif
-
-   ! -----
-   ! - compute derivatives...
-   ! ------------------------
-
-   ! 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
-   if(ixDomainType==iname_soil)then
-    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
-     case default;          dVolTot_dPsi0(ixControlIndex) = dTheta_dPsi(mLayerMatricHeadTrial(ixControlIndex),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
-    end select
-   endif
-
-   ! compute the derivative in liquid water content w.r.t. temperature
-   ! --> partially frozen: dependence of liquid water on temperature
-   if(xTemp<Tcrit)then
-    select case(ixDomainType)
-     case(iname_veg);   dTheta_dTkCanopy         = dFracLiq_dTk(xTemp,snowfrz_scale)*scalarCanopyWat/(iden_water*canopyDepth)
-     case(iname_snow);  mLayerdTheta_dTk(iLayer) = dFracLiq_dTk(xTemp,snowfrz_scale)*mLayerVolFracWatTrial(iLayer)
-     case(iname_soil);  mLayerdTheta_dTk(iLayer) = dTheta_dTk(xTemp,theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
-     case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
-    end select  ! domain type
-
-   ! --> unfrozen: no dependence of liquid water on temperature
-   else
-    select case(ixDomainType)
-     case(iname_veg);             dTheta_dTkCanopy         = 0._rkind
-     case(iname_snow, iname_soil);   mLayerdTheta_dTk(iLayer) = 0._rkind
-     case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return
-    end select  ! domain type
-   endif
-
-
-
-
-   ! -----
-   ! - 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 mass of liquid water and ice
-      call updateVegSundials(&
-                      xTemp,                                        & ! intent(in)   : temperature (K)
-                      scalarCanopyWatTrial,                         & ! intent(in)   : mass of total water (-)
-                      snowfrz_scale,                                & ! intent(in)   : scaling parameter for the snow freezing curve (K-1)
-                      scalarCanopyTempPrime,                        & ! intent(in)
-                      scalarCanopyWatPrime,                         & ! intent(in)   : mass of total water (-)
-                      scalarCanopyLiqTrial,                         & ! intent(out)  : trial mass of liquid water (-)
-                      scalarCanopyIceTrial,                         & ! intent(out)  : trial mass of ice (-)
-                      scalarCanopyLiqPrime,                         & ! intent(out)  : trial mass of liquid water (-)
-                      scalarCanopyIcePrime,                         & ! intent(out)  : trial mass of 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
-
-     ! *** snow layers
-     case(iname_snow)
-
-      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),                      & !
-                      mLayerVolFracWatPrime(iLayer),                & ! intent(in)
-                      mLayerVolFracLiqTrial(iLayer),                & ! intent(out)  : trial volumetric fraction of liquid water (-)
-                      mLayerVolFracIceTrial(iLayer),                & ! intent(out)  : trial volumetric fraction if ice (-)
-                      mLayerVolFracLiqPrime(iLayer),                & ! intent(out)
-                      mLayerVolFracIcePrime(iLayer),                & ! intent(out)
-                      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
-
-     ! *** soil layers
-     case(iname_soil)
-
-      ! compute volumetric fraction of liquid water and ice, step size dt_cur changes here
-      call updateSoilSundials(&
-                      dt_cur,                                            &
-                      xTemp,                                             & ! intent(in)   : temperature (K)
-                      mLayerMatricHeadTrial(ixControlIndex),             & ! intent(in)   : total water matric potential (m)
-                      mLayerMatricHeadPrev(ixControlIndex),              & ! intent(in)
-                      mLayerVolFracWatPrev(iLayer),                      & ! intent(in)
-                      mLayerTempPrime(iLayer),                           &
-                      mLayerMatricHeadPrime(ixControlIndex),             &
-                     ! intent(in)   : soil parameters
-                      vGn_alpha(ixControlIndex),                         &
-                      vGn_n(ixControlIndex),                             &
-                      theta_sat(ixControlIndex),                         &
-                      theta_res(ixControlIndex),                         &
-                      vGn_m(ixControlIndex),                             &
-                      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),                     &
-                      mLayerVolFracLiqPrime(iLayer),                     &
-                      mLayerVolFracIcePrime(iLayer),                     &
-                      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
-
-    end select  ! domain type
-
-   ! final check
-   else
-
-    ! 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
-    endif
-
-   endif  ! if energy state or solution is coupled
-
-   ! -----
-   ! - update temperatures...
-   ! ------------------------
-
-
-  end do iterations ! iterating
-
-  ! save temperature
-  select case(ixDomainType)
-   case(iname_veg);              scalarCanopyTempTrial   = xTemp
-   case(iname_snow, iname_soil); mLayerTempTrial(iLayer) = xTemp
-  end select
-
-  ! =======================================================================================================================================
-  ! =======================================================================================================================================
-
-  ! -----
-  ! - compute the liquid water matric potential (and necessay derivatives)...
-  ! -------------------------------------------------------------------------
-
-  ! only for soil
-  if(ixDomainType==iname_soil)then
-
-   ! 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
-
-    ! 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
-
-   ! 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)
-                    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
-
- end subroutine updateVarsSundials
-
-
-end module updateVarsSundials_module
+                       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
+   
+    end subroutine updateVarsSundials
+   
+   
+    ! **********************************************************************************************************
+    ! 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
+   
\ No newline at end of file
-- 
GitLab