diff --git a/build/source/engine/computFlux.f90 b/build/source/engine/computFlux.f90
index 9d443f370006e5729d674ef5d31a5e2230ba3df1..73ed736c9b979dd1fe5fdbdb6d8088c878d2ea01 100755
--- a/build/source/engine/computFlux.f90
+++ b/build/source/engine/computFlux.f90
@@ -113,7 +113,7 @@ subroutine computFlux(&
firstSplitOper, & ! intent(in): flag to indicate if we are processing the first flux call in a splitting operation
computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation
scalarSolution, & ! intent(in): flag to indicate the scalar solution
- requireLWBal, & ! intent(in): flag to indicate if we need longwave to be balanced
+ insideIDA, & ! intent(in): flag if inside Sundials solver
drainageMeltPond, & ! intent(in): drainage from the surface melt pond (kg m-2 s-1)
! input: state variables
scalarCanairTempTrial, & ! intent(in): trial value for the temperature of the canopy air space (K)
@@ -167,7 +167,7 @@ subroutine computFlux(&
logical(lgt),intent(in) :: firstSplitOper ! flag to indicate if we are processing the first flux call in a splitting operation
logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if computing fluxes over vegetation
logical(lgt),intent(in) :: scalarSolution ! flag to denote if implementing the scalar solution
- logical(lgt),intent(in) :: requireLWBal ! flag to indicate if we need longwave to be balanced
+ logical(lgt),intent(in) :: insideIDA ! flag if inside Sundials solver
real(dp),intent(in) :: drainageMeltPond ! drainage from the surface melt pond (kg m-2 s-1)
! input: state variables
real(dp),intent(in) :: scalarCanairTempTrial ! trial value for temperature of the canopy air space (K)
@@ -441,7 +441,7 @@ subroutine computFlux(&
firstSubStep, & ! intent(in): flag to indicate if we are processing the first sub-step
firstFluxCall, & ! intent(in): flag to indicate if we are processing the first flux call
computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation
- requireLWBal, & ! intent(in): flag to indicate if we need longwave to be balanced
+ insideIDA, & ! intent(in): flag to indicate if we need longwave to be balanced
! input: model state variables
upperBoundTemp, & ! intent(in): temperature of the upper boundary (K) --> NOTE: use air temperature
scalarCanairTempTrial, & ! intent(in): trial value of the canopy air space temperature (K)
diff --git a/build/source/engine/computJacob.f90 b/build/source/engine/computJacob.f90
index 4073f566e8f316de439c54be25e956f06c0b31b0..6161682edbe8885907048e701531d14117a96f09 100755
--- a/build/source/engine/computJacob.f90
+++ b/build/source/engine/computJacob.f90
@@ -20,772 +20,994 @@
module computJacob_module
-! data types
-USE nrtype
-
-! derived types to define the data structures
-USE data_types,only:&
- var_ilength, & ! data vector with variable length dimension (i4b)
- var_dlength ! data vector with variable length dimension (dp)
-
-! named variables for structure elements
-USE var_lookup,only:iLookPROG ! named variables for structure elements
-USE var_lookup,only:iLookDIAG ! named variables for structure elements
-USE var_lookup,only:iLookINDEX ! named variables for structure elements
-USE var_lookup,only:iLookDERIV ! named variables for structure elements
-
-! access the global print flag
-USE globalData,only:globalPrintFlag
-
-! access missing values
-USE globalData,only:integerMissing ! missing integer
-USE globalData,only:realMissing ! missing real number
-
-! domain types
-USE globalData,only:iname_veg ! named variables for vegetation
-USE globalData,only:iname_snow ! named variables for snow
-USE globalData,only:iname_soil ! named variables for soil
-
-! 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 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
-
-! access named variables to describe the form and structure of the matrices used in the numerical solver
-USE globalData,only: ku ! number of super-diagonal bands
-USE globalData,only: kl ! number of sub-diagonal bands
-USE globalData,only: ixDiag ! index for the diagonal band
-USE globalData,only: nBands ! length of the leading dimension of the band diagonal matrix
-USE globalData,only: ixFullMatrix ! named variable for the full Jacobian matrix
-USE globalData,only: ixBandMatrix ! named variable for the band diagonal matrix
-USE globalData,only: iJac1 ! first layer of the Jacobian to print
-USE globalData,only: iJac2 ! last layer of the Jacobian to print
-
-! constants
-USE multiconst,only:&
- LH_fus, & ! latent heat of fusion (J kg-1)
- iden_ice, & ! intrinsic density of ice (kg m-3)
- iden_water ! intrinsic density of liquid water (kg m-3)
-
-implicit none
-! define constants
-real(dp),parameter :: verySmall=tiny(1.0_dp) ! a very small number
-integer(i4b),parameter :: ixBandOffset=kl+ku+1 ! offset in the band Jacobian matrix
-
-private
-public::computJacob
-contains
-
- ! **********************************************************************************************************
- ! public subroutine computJacob: compute the Jacobian matrix
- ! **********************************************************************************************************
- subroutine computJacob(&
- ! input: model control
- dt, & ! intent(in): length of the time step (seconds)
- nSnow, & ! intent(in): number of snow layers
- nSoil, & ! intent(in): number of soil layers
- nLayers, & ! intent(in): total number of layers
- computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation
- computeBaseflow, & ! intent(in): flag to indicate if we need to compute baseflow
- ixMatrix, & ! intent(in): form of the Jacobian matrix
- ! input: data structures
- indx_data, & ! intent(in): index data
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- diag_data, & ! intent(in): model diagnostic variables for a local HRU
- deriv_data, & ! intent(in): derivatives in model fluxes w.r.t. relevant state variables
- dBaseflow_dMatric, & ! intent(in): derivative in baseflow w.r.t. matric head (s-1)
- ! input-output: Jacobian and its diagonal
- dMat, & ! intent(inout): diagonal of the Jacobian matrix
- aJac, & ! intent(out): Jacobian matrix
- ! output: error control
- err,message) ! intent(out): error code and error message
- ! -----------------------------------------------------------------------------------------------------------------
- implicit none
- ! input: model control
- real(dp),intent(in) :: dt ! length of the time step (seconds)
- integer(i4b),intent(in) :: nSnow ! number of snow layers
- integer(i4b),intent(in) :: nSoil ! number of soil layers
- integer(i4b),intent(in) :: nLayers ! total number of layers in the snow+soil domain
- logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if computing fluxes over vegetation
- logical(lgt),intent(in) :: computeBaseflow ! flag to indicate if computing baseflow
- integer(i4b),intent(in) :: ixMatrix ! form of the Jacobian matrix
- ! input: data structures
- 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
- type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU
- type(var_dlength),intent(in) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
- real(dp),intent(in) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1)
- ! input-output: Jacobian and its diagonal
- real(dp),intent(inout) :: dMat(:) ! diagonal of the Jacobian matrix
- real(dp),intent(out) :: aJac(:,:) ! Jacobian matrix
- ! output variables
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! --------------------------------------------------------------
- ! * local variables
- ! --------------------------------------------------------------
- ! indices of model state variables
- integer(i4b) :: jState ! index of state within the state subset
- integer(i4b) :: qState ! index of cross-derivative state variable for baseflow
- integer(i4b) :: nrgState ! energy state variable
- integer(i4b) :: watState ! hydrology state variable
- integer(i4b) :: nState ! number of state variables
- ! indices of model layers
- integer(i4b) :: iLayer ! index of model layer
- integer(i4b) :: jLayer ! index of model layer within the full state vector (hydrology)
- integer(i4b) :: pLayer ! indices of soil layers (used for the baseflow derivatives)
- ! conversion factors
- real(dp) :: convLiq2tot ! factor to convert liquid water derivative to total water derivative
- ! --------------------------------------------------------------
- ! associate variables from data structures
- associate(&
- ! indices of model state variables
- ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable
- ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable
- ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass)
- ixTopNrg => indx_data%var(iLookINDEX%ixTopNrg)%dat(1) ,& ! intent(in): [i4b] index of upper-most energy state in the snow+soil subdomain
- ixTopHyd => indx_data%var(iLookINDEX%ixTopHyd)%dat(1) ,& ! intent(in): [i4b] index of upper-most hydrology state in the snow+soil subdomain
- ixAqWat => indx_data%var(iLookINDEX%ixAqWat)%dat(1) ,& ! intent(in): [i4b] index of water storage in the aquifer
- ! vectors of indices for specfic state types within specific sub-domains IN THE FULL STATE VECTOR
- 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
- ! vector of energy indices for the snow and soil domains
- ! NOTE: states not in the subset are equal to integerMissing
- ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow+soil domain
- ixSnowOnlyNrg => indx_data%var(iLookINDEX%ixSnowOnlyNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow domain
- ixSoilOnlyNrg => indx_data%var(iLookINDEX%ixSoilOnlyNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the soil domain
- ! vector of hydrology indices for the snow and soil domains
- ! NOTE: states not in the subset are equal to integerMissing
- ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow+soil domain
- ixSnowOnlyHyd => indx_data%var(iLookINDEX%ixSnowOnlyHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow domain
- ixSoilOnlyHyd => indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the soil domain
- ! number of state variables of a specific type
- nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain
- nSnowOnlyNrg => indx_data%var(iLookINDEX%nSnowOnlyNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow domain
- nSoilOnlyNrg => indx_data%var(iLookINDEX%nSoilOnlyNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the soil domain
- nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow+soil domain
- nSnowOnlyHyd => indx_data%var(iLookINDEX%nSnowOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow domain
- nSoilOnlyHyd => indx_data%var(iLookINDEX%nSoilOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the soil domain
- ! type and index of model control volume
- ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] index of the type of hydrology states in snow+soil domain
- ixDomainType => indx_data%var(iLookINDEX%ixDomainType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the domain (iname_veg, iname_snow, iname_soil)
- ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of the control volume for specific model domains
- ! mapping between states and model layers
- ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] list of indices in the full state vector that are in the state subset
- ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset in each element of the full state vector
- ! derivatives in net vegetation energy fluxes w.r.t. relevant state variables
- dCanairNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. canopy air temperature
- dCanairNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. canopy temperature
- dCanairNetFlux_dGroundTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dGroundTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. ground temperature
- dCanopyNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. canopy air temperature
- dCanopyNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. canopy temperature
- dCanopyNetFlux_dGroundTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dGroundTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. ground temperature
- dCanopyNetFlux_dCanWat => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in net canopy fluxes w.r.t. canopy liquid water content
- dGroundNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net ground flux w.r.t. canopy air temperature
- dGroundNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net ground flux w.r.t. canopy temperature
- dGroundNetFlux_dCanWat => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in net ground fluxes w.r.t. canopy liquid water content
- ! derivatives in evaporative fluxes w.r.t. relevant state variables
- dCanopyEvaporation_dTCanair => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTCanair )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy air temperature
- dCanopyEvaporation_dTCanopy => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTCanopy )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy temperature
- dCanopyEvaporation_dTGround => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTGround )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. ground temperature
- dCanopyEvaporation_dCanWat => deriv_data%var(iLookDERIV%dCanopyEvaporation_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy liquid water content
- dGroundEvaporation_dTCanair => deriv_data%var(iLookDERIV%dGroundEvaporation_dTCanair )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy air temperature
- dGroundEvaporation_dTCanopy => deriv_data%var(iLookDERIV%dGroundEvaporation_dTCanopy )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy temperature
- dGroundEvaporation_dTGround => deriv_data%var(iLookDERIV%dGroundEvaporation_dTGround )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. ground temperature
- dGroundEvaporation_dCanWat => deriv_data%var(iLookDERIV%dGroundEvaporation_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy liquid water content
- ! derivatives in canopy water w.r.t canopy temperature
- dCanLiq_dTcanopy => deriv_data%var(iLookDERIV%dCanLiq_dTcanopy )%dat(1) ,& ! intent(in): [dp] derivative of canopy liquid storage w.r.t. temperature
- dTheta_dTkCanopy => deriv_data%var(iLookDERIV%dTheta_dTkCanopy )%dat(1) ,& ! intent(in): [dp] derivative of volumetric liquid water content w.r.t. temperature
- ! derivatives in canopy liquid fluxes w.r.t. canopy water
- scalarCanopyLiqDeriv => deriv_data%var(iLookDERIV%scalarCanopyLiqDeriv )%dat(1) ,& ! intent(in): [dp] derivative in (throughfall + drainage) w.r.t. canopy liquid water
- ! derivatives in energy fluxes at the interface of snow+soil layers w.r.t. temperature in layers above and below
- dNrgFlux_dTempAbove => deriv_data%var(iLookDERIV%dNrgFlux_dTempAbove )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. temperature in the layer above
- dNrgFlux_dTempBelow => deriv_data%var(iLookDERIV%dNrgFlux_dTempBelow )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. temperature in the layer below
- ! derivative in liquid water fluxes at the interface of snow layers w.r.t. volumetric liquid water content in the layer above
- iLayerLiqFluxSnowDeriv => deriv_data%var(iLookDERIV%iLayerLiqFluxSnowDeriv )%dat ,& ! intent(in): [dp(:)] derivative in vertical liquid water flux at layer interfaces
- ! derivative in liquid water fluxes for the soil domain w.r.t hydrology state variables
- dVolTot_dPsi0 => deriv_data%var(iLookDERIV%dVolTot_dPsi0 )%dat ,& ! intent(in): [dp(:)] derivative in total water content w.r.t. total water matric potential
- dq_dHydStateAbove => deriv_data%var(iLookDERIV%dq_dHydStateAbove )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer above
- dq_dHydStateBelow => deriv_data%var(iLookDERIV%dq_dHydStateBelow )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer below
- dCompress_dPsi => deriv_data%var(iLookDERIV%dCompress_dPsi )%dat ,& ! intent(in): [dp(:)] derivative in compressibility w.r.t matric head
- ! derivative in baseflow flux w.r.t. aquifer storage
- dBaseflow_dAquifer => deriv_data%var(iLookDERIV%dBaseflow_dAquifer )%dat(1) ,& ! intent(out): [dp(:)] erivative in baseflow flux w.r.t. aquifer storage (s-1)
- ! derivative in liquid water fluxes for the soil domain w.r.t energy state variables
- dq_dNrgStateAbove => deriv_data%var(iLookDERIV%dq_dNrgStateAbove )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer above
- dq_dNrgStateBelow => deriv_data%var(iLookDERIV%dq_dNrgStateBelow )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer below
- mLayerdTheta_dTk => deriv_data%var(iLookDERIV%mLayerdTheta_dTk )%dat ,& ! intent(in): [dp(:)] derivative of volumetric liquid water content w.r.t. temperature
- ! diagnostic variables
- scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg)%dat(1) ,& ! intent(in): [dp] fraction of liquid water on vegetation (-)
- scalarBulkVolHeatCapVeg => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg)%dat(1) ,& ! intent(in): [dp] bulk volumetric heat capacity of vegetation (J m-3 K-1)
- mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow)%dat ,& ! intent(in): [dp(:)] fraction of liquid water in each snow layer (-)
- mLayerVolHtCapBulk => diag_data%var(iLookDIAG%mLayerVolHtCapBulk)%dat ,& ! intent(in): [dp(:)] bulk volumetric heat capacity in each snow and soil layer (J m-3 K-1)
- scalarSoilControl => diag_data%var(iLookDIAG%scalarSoilControl)%dat(1) ,& ! intent(in): [dp] soil control on infiltration, zero or one
- ! canopy and layer depth
- canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp ] canopy depth (m)
- mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat & ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
- ) ! making association with data in structures
- ! --------------------------------------------------------------
- ! initialize error control
- err=0; message='computJacob/'
-
- ! *********************************************************************************************************************************************************
- ! *********************************************************************************************************************************************************
- ! * PART 0: PRELIMINARIES (INITIALIZE JACOBIAN AND COMPUTE TIME-VARIABLE DIAGONAL TERMS)
- ! *********************************************************************************************************************************************************
- ! *********************************************************************************************************************************************************
-
- ! get the number of state variables
- nState = size(dMat)
-
- ! initialize the Jacobian
- ! NOTE: this needs to be done every time, since Jacobian matrix is modified in the solver
- aJac(:,:) = 0._dp ! analytical Jacobian matrix
-
- ! compute terms in the Jacobian for vegetation (excluding fluxes)
- ! NOTE: energy for vegetation is computed *within* the iteration loop as it includes phase change
- if(ixVegNrg/=integerMissing) dMat(ixVegNrg) = scalarBulkVolHeatCapVeg + LH_fus*iden_water*dTheta_dTkCanopy ! volumetric heat capacity of the vegetation (J m-3 K-1)
-
- ! compute additional terms for the Jacobian for the snow-soil domain (excluding fluxes)
- ! NOTE: energy for snow+soil is computed *within* the iteration loop as it includes phase change
- do iLayer=1,nLayers
- if(ixSnowSoilNrg(iLayer)/=integerMissing) dMat(ixSnowSoilNrg(iLayer)) = mLayerVolHtCapBulk(iLayer) + LH_fus*iden_water*mLayerdTheta_dTk(iLayer)
- end do
-
- ! compute additional terms for the Jacobian for the soil domain (excluding fluxes)
- do iLayer=1,nSoil
- if(ixSoilOnlyHyd(iLayer)/=integerMissing) dMat(ixSoilOnlyHyd(iLayer)) = dVolTot_dPsi0(iLayer) + dCompress_dPsi(iLayer)
- end do
-
- ! define the form of the matrix
- select case(ixMatrix)
-
- ! *********************************************************************************************************************************************************
- ! *********************************************************************************************************************************************************
- ! * PART 1: BAND MATRIX
- ! *********************************************************************************************************************************************************
- ! *********************************************************************************************************************************************************
- case(ixBandMatrix)
-
- ! check
- if(size(aJac,1)/=nBands .or. size(aJac,2)/=size(dMat))then
- message=trim(message)//'unexpected shape of the Jacobian matrix: expect aJac(nBands,nState)'
- err=20; return
- end if
-
- ! -----
- ! * energy and liquid fluxes over vegetation...
- ! ---------------------------------------------
- if(computeVegFlux)then ! (derivatives only defined when vegetation protrudes over the surface)
-
- ! * diagonal elements for the vegetation canopy (-)
- if(ixCasNrg/=integerMissing) aJac(ixDiag,ixCasNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanairTemp) + dMat(ixCasNrg)
- if(ixVegNrg/=integerMissing) aJac(ixDiag,ixVegNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanopyTemp) + dMat(ixVegNrg)
- if(ixVegHyd/=integerMissing) aJac(ixDiag,ixVegHyd) = -scalarFracLiqVeg*(dCanopyEvaporation_dCanWat - scalarCanopyLiqDeriv)*dt + 1._dp ! ixVegHyd: CORRECT
-
- ! * cross-derivative terms w.r.t. canopy water
- if(ixVegHyd/=integerMissing)then
- ! cross-derivative terms w.r.t. system temperatures (kg m-2 K-1)
- if(ixCasNrg/=integerMissing) aJac(ixOffDiag(ixVegHyd,ixCasNrg),ixCasNrg) = -dCanopyEvaporation_dTCanair*dt ! ixCasNrg: CORRECT
- if(ixVegNrg/=integerMissing) aJac(ixOffDiag(ixVegHyd,ixVegNrg),ixVegNrg) = -dCanopyEvaporation_dTCanopy*dt + dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy ! ixVegNrg: CORRECT
- if(ixTopNrg/=integerMissing) aJac(ixOffDiag(ixVegHyd,ixTopNrg),ixTopNrg) = -dCanopyEvaporation_dTGround*dt ! ixTopNrg: CORRECT
- ! cross-derivative terms w.r.t. canopy water (kg-1 m2)
- if(ixTopHyd/=integerMissing) aJac(ixOffDiag(ixTopHyd,ixVegHyd),ixVegHyd) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarFracLiqVeg*scalarCanopyLiqDeriv)/iden_water
- ! cross-derivative terms w.r.t. canopy liquid water (J m-1 kg-1)
- ! NOTE: dIce/dLiq = (1 - scalarFracLiqVeg); dIce*LH_fus/canopyDepth = J m-3; dLiq = kg m-2
- if(ixVegNrg/=integerMissing) aJac(ixOffDiag(ixVegNrg,ixVegHyd),ixVegHyd) = (dt/canopyDepth) *(-dCanopyNetFlux_dCanWat) - (1._dp - scalarFracLiqVeg)*LH_fus/canopyDepth ! dF/dLiq
- if(ixTopNrg/=integerMissing) aJac(ixOffDiag(ixTopNrg,ixVegHyd),ixVegHyd) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanWat)
- endif
-
- ! cross-derivative terms between surface hydrology and the temperature of the vegetation canopy (K-1)
- if(ixVegNrg/=integerMissing)then
- if(ixTopHyd/=integerMissing) aJac(ixOffDiag(ixTopHyd,ixVegNrg),ixVegNrg) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarCanopyLiqDeriv*dCanLiq_dTcanopy)/iden_water
- endif
-
- ! cross-derivative terms w.r.t. the temperature of the canopy air space (J m-3 K-1)
- if(ixCasNrg/=integerMissing)then
- if(ixVegNrg/=integerMissing) aJac(ixOffDiag(ixCasNrg,ixVegNrg),ixVegNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanopyTemp)
- if(ixTopNrg/=integerMissing) aJac(ixOffDiag(ixCasNrg,ixTopNrg),ixTopNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dGroundTemp)
- endif
-
- ! cross-derivative terms w.r.t. the temperature of the vegetation canopy (J m-3 K-1)
- if(ixVegNrg/=integerMissing)then
- if(ixCasNrg/=integerMissing) aJac(ixOffDiag(ixVegNrg,ixCasNrg),ixCasNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanairTemp)
- if(ixTopNrg/=integerMissing) aJac(ixOffDiag(ixVegNrg,ixTopNrg),ixTopNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dGroundTemp)
- endif
-
- ! cross-derivative terms w.r.t. the temperature of the surface (J m-3 K-1)
- if(ixTopNrg/=integerMissing)then
- if(ixCasNrg/=integerMissing) aJac(ixOffDiag(ixTopNrg,ixCasNrg),ixCasNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanairTemp)
- if(ixVegNrg/=integerMissing) aJac(ixOffDiag(ixTopNrg,ixVegNrg),ixVegNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanopyTemp)
- endif
-
- endif ! if there is a need to compute energy fluxes within vegetation
-
- ! -----
- ! * energy fluxes for the snow+soil domain...
- ! -------------------------------------------
- if(nSnowSoilNrg>0)then
- do iLayer=1,nLayers ! loop through all layers in the snow+soil domain
-
- ! check if the state is in the subset
- if(ixSnowSoilNrg(iLayer)==integerMissing) cycle
-
- ! - define index within the state subset and the full state vector
- jState = ixSnowSoilNrg(iLayer) ! index within the state subset
-
- ! - diagonal elements
- aJac(ixDiag,jState) = (dt/mLayerDepth(iLayer))*(-dNrgFlux_dTempBelow(iLayer-1) + dNrgFlux_dTempAbove(iLayer)) + dMat(jState)
-
- ! - lower-diagonal elements
- if(iLayer > 1)then
- if(ixSnowSoilNrg(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSnowSoilNrg(iLayer-1),jState),jState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dTempBelow(iLayer-1) )
- endif
-
- ! - upper diagonal elements
- if(iLayer < nLayers)then
- if(ixSnowSoilNrg(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSnowSoilNrg(iLayer+1),jState),jState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dTempAbove(iLayer ) )
+ ! data types
+ USE nrtype
+
+ ! derived types to define the data structures
+ USE data_types,only:&
+ var_ilength, & ! data vector with variable length dimension (i4b)
+ var_dlength ! data vector with variable length dimension (rkind)
+
+ ! named variables for structure elements
+ USE var_lookup,only:iLookPROG ! named variables for structure elements
+ USE var_lookup,only:iLookDIAG ! named variables for structure elements
+ USE var_lookup,only:iLookINDEX ! named variables for structure elements
+ USE var_lookup,only:iLookDERIV ! named variables for structure elements
+
+ ! access the global print flag
+ USE globalData,only:globalPrintFlag
+
+ ! access missing values
+ USE globalData,only:integerMissing ! missing integer
+ USE globalData,only:realMissing ! missing real number
+
+ ! domain types
+ USE globalData,only:iname_veg ! named variables for vegetation
+ USE globalData,only:iname_snow ! named variables for snow
+ USE globalData,only:iname_soil ! named variables for soil
+
+ ! 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 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
+
+ ! access named variables to describe the form and structure of the matrices used in the numerical solver
+ USE globalData,only: ku ! number of super-diagonal bands
+ USE globalData,only: kl ! number of sub-diagonal bands
+ USE globalData,only: ixDiag ! index for the diagonal band
+ USE globalData,only: nBands ! length of the leading dimension of the band diagonal matrix
+ USE globalData,only: ixFullMatrix ! named variable for the full Jacobian matrix
+ USE globalData,only: ixBandMatrix ! named variable for the band diagonal matrix
+ USE globalData,only: iJac1 ! first layer of the Jacobian to print
+ USE globalData,only: iJac2 ! last layer of the Jacobian to print
+
+ ! constants
+ USE multiconst,only:&
+ LH_fus, & ! latent heat of fusion (J kg-1)
+ iden_water ! intrinsic density of liquid water (kg m-3)
+
+ implicit none
+ ! define constants
+ real(rkind),parameter :: verySmall=tiny(1.0_rkind) ! a very small number
+ integer(i4b),parameter :: ixBandOffset=kl+ku+1 ! offset in the band Jacobian matrix
+
+ private
+ public::computJacob
+ contains
+
+ ! **********************************************************************************************************
+ ! public subroutine computJacob: compute the Jacobian matrix
+ ! **********************************************************************************************************
+ subroutine computJacob(&
+ ! input: model control
+ dt, & ! intent(in): length of the time step (seconds)
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation
+ computeBaseflow, & ! intent(in): flag to indicate if we need to compute baseflow
+ ixMatrix, & ! intent(in): form of the Jacobian matrix
+ ! input: data structures
+ indx_data, & ! intent(in): index data
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ diag_data, & ! intent(in): model diagnostic variables for a local HRU
+ deriv_data, & ! intent(in): derivatives in model fluxes w.r.t. relevant state variables
+ dBaseflow_dMatric, & ! intent(in): derivative in baseflow w.r.t. matric head (s-1)
+ ! input-output: Jacobian and its diagonal
+ dMat, & ! intent(inout): diagonal of the Jacobian matrix
+ aJac, & ! intent(out): Jacobian matrix
+ ! output: error control
+ err,message) ! intent(out): error code and error message
+ ! -----------------------------------------------------------------------------------------------------------------
+ implicit none
+ ! input: model control
+ real(rkind),intent(in) :: dt ! length of the time step (seconds)
+ integer(i4b),intent(in) :: nSnow ! number of snow layers
+ integer(i4b),intent(in) :: nSoil ! number of soil layers
+ integer(i4b),intent(in) :: nLayers ! total number of layers in the snow+soil domain
+ logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if computing fluxes over vegetation
+ logical(lgt),intent(in) :: computeBaseflow ! flag to indicate if computing baseflow
+ integer(i4b),intent(in) :: ixMatrix ! form of the Jacobian matrix
+ ! input: data structures
+ 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
+ type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU
+ type(var_dlength),intent(in) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
+ real(rkind),intent(in) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1)
+ ! input-output: Jacobian and its diagonal
+ real(rkind),intent(inout) :: dMat(:) ! diagonal of the Jacobian matrix
+ real(rkind),intent(out) :: aJac(:,:) ! Jacobian matrix
+ ! output variables
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! --------------------------------------------------------------
+ ! * local variables
+ ! --------------------------------------------------------------
+ ! indices of model state variables
+ integer(i4b) :: jState ! index of state within the state subset
+ integer(i4b) :: qState ! index of cross-derivative state variable for baseflow
+ integer(i4b) :: nrgState ! energy state variable
+ integer(i4b) :: watState ! hydrology state variable
+ integer(i4b) :: nState ! number of state variables
+ ! indices of model layers
+ integer(i4b) :: iLayer ! index of model layer
+ integer(i4b) :: jLayer ! index of model layer within the full state vector (hydrology)
+ integer(i4b) :: pLayer ! indices of soil layers (used for the baseflow derivatives)
+ ! conversion factors
+ real(rkind) :: convLiq2tot ! factor to convert liquid water derivative to total water derivative
+ ! --------------------------------------------------------------
+ ! associate variables from data structures
+ associate(&
+ ! indices of model state variables
+ ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable
+ ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable
+ ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass)
+ ixTopNrg => indx_data%var(iLookINDEX%ixTopNrg)%dat(1) ,& ! intent(in): [i4b] index of upper-most energy state in the snow+soil subdomain
+ ixTopHyd => indx_data%var(iLookINDEX%ixTopHyd)%dat(1) ,& ! intent(in): [i4b] index of upper-most hydrology state in the snow+soil subdomain
+ ixAqWat => indx_data%var(iLookINDEX%ixAqWat)%dat(1) ,& ! intent(in): [i4b] index of water storage in the aquifer
+ ! vectors of indices for specfic state types within specific sub-domains IN THE FULL STATE VECTOR
+ 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
+ ! vector of energy indices for the snow and soil domains
+ ! NOTE: states not in the subset are equal to integerMissing
+ ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow+soil domain
+ ixSnowOnlyNrg => indx_data%var(iLookINDEX%ixSnowOnlyNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow domain
+ ixSoilOnlyNrg => indx_data%var(iLookINDEX%ixSoilOnlyNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the soil domain
+ ! vector of hydrology indices for the snow and soil domains
+ ! NOTE: states not in the subset are equal to integerMissing
+ ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow+soil domain
+ ixSnowOnlyHyd => indx_data%var(iLookINDEX%ixSnowOnlyHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow domain
+ ixSoilOnlyHyd => indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the soil domain
+ ! number of state variables of a specific type
+ nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain
+ nSnowOnlyNrg => indx_data%var(iLookINDEX%nSnowOnlyNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow domain
+ nSoilOnlyNrg => indx_data%var(iLookINDEX%nSoilOnlyNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the soil domain
+ nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow+soil domain
+ nSnowOnlyHyd => indx_data%var(iLookINDEX%nSnowOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow domain
+ nSoilOnlyHyd => indx_data%var(iLookINDEX%nSoilOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the soil domain
+ ! type and index of model control volume
+ ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] index of the type of hydrology states in snow+soil domain
+ ixDomainType => indx_data%var(iLookINDEX%ixDomainType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the domain (iname_veg, iname_snow, iname_soil)
+ ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of the control volume for specific model domains
+ ! mapping between states and model layers
+ ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] list of indices in the full state vector that are in the state subset
+ ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset in each element of the full state vector
+ ! derivatives in net vegetation energy fluxes w.r.t. relevant state variables
+ dCanairNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. canopy air temperature
+ dCanairNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. canopy temperature
+ dCanairNetFlux_dGroundTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dGroundTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. ground temperature
+ dCanopyNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. canopy air temperature
+ dCanopyNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. canopy temperature
+ dCanopyNetFlux_dGroundTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dGroundTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. ground temperature
+ dCanopyNetFlux_dCanWat => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in net canopy fluxes w.r.t. canopy total water content
+ dGroundNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net ground flux w.r.t. canopy air temperature
+ dGroundNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net ground flux w.r.t. canopy temperature
+ dGroundNetFlux_dCanWat => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in net ground fluxes w.r.t. canopy total water content
+ ! derivatives in evaporative fluxes w.r.t. relevant state variables
+ dCanopyEvaporation_dTCanair => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTCanair )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy air temperature
+ dCanopyEvaporation_dTCanopy => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTCanopy )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy temperature
+ dCanopyEvaporation_dTGround => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTGround )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. ground temperature
+ dCanopyEvaporation_dCanWat => deriv_data%var(iLookDERIV%dCanopyEvaporation_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy total water content
+ dGroundEvaporation_dTCanair => deriv_data%var(iLookDERIV%dGroundEvaporation_dTCanair )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy air temperature
+ dGroundEvaporation_dTCanopy => deriv_data%var(iLookDERIV%dGroundEvaporation_dTCanopy )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy temperature
+ dGroundEvaporation_dTGround => deriv_data%var(iLookDERIV%dGroundEvaporation_dTGround )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. ground temperature
+ dGroundEvaporation_dCanWat => deriv_data%var(iLookDERIV%dGroundEvaporation_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy total water content
+ ! derivatives in canopy water w.r.t canopy temperature
+ dCanLiq_dTcanopy => deriv_data%var(iLookDERIV%dCanLiq_dTcanopy )%dat(1) ,& ! intent(in): [dp] derivative of canopy liquid storage w.r.t. temperature
+ dTheta_dTkCanopy => deriv_data%var(iLookDERIV%dTheta_dTkCanopy )%dat(1) ,& ! intent(in): [dp] derivative of volumetric liquid water content w.r.t. temperature
+ ! derivatives in canopy liquid fluxes w.r.t. canopy water
+ scalarCanopyLiqDeriv => deriv_data%var(iLookDERIV%scalarCanopyLiqDeriv )%dat(1) ,& ! intent(in): [dp] derivative in (throughfall + drainage) w.r.t. canopy liquid water
+ ! derivatives in energy fluxes at the interface of snow+soil layers w.r.t. temperature in layers above and below
+ dNrgFlux_dTempAbove => deriv_data%var(iLookDERIV%dNrgFlux_dTempAbove )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. temperature in the layer above
+ dNrgFlux_dTempBelow => deriv_data%var(iLookDERIV%dNrgFlux_dTempBelow )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. temperature in the layer below
+ ! derivatives in energy fluxes at the interface of snow+soil layers w.r.t. water state in layers above and below
+ dNrgFlux_dWatAbove => deriv_data%var(iLookDERIV%dNrgFlux_dWatAbove )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. water state in the layer above
+ dNrgFlux_dWatBelow => deriv_data%var(iLookDERIV%dNrgFlux_dWatBelow )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. water state in the layer below
+ ! derivatives in soil transpiration w.r.t. canopy state variables
+ mLayerdTrans_dTCanair => deriv_data%var(iLookDERIV%mLayerdTrans_dTCanair )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. canopy air temperature
+ mLayerdTrans_dTCanopy => deriv_data%var(iLookDERIV%mLayerdTrans_dTCanopy )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. canopy temperature
+ mLayerdTrans_dTGround => deriv_data%var(iLookDERIV%mLayerdTrans_dTGround )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. ground temperature
+ mLayerdTrans_dCanWat => deriv_data%var(iLookDERIV%mLayerdTrans_dCanWat )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. canopy total water
+ ! derivatives in aquifer transpiration w.r.t. canopy state variables
+ dAquiferTrans_dTCanair => deriv_data%var(iLookDERIV%dAquiferTrans_dTCanair )%dat(1) ,& ! intent(in): derivatives in the aquifer transpiration flux w.r.t. canopy air temperature
+ dAquiferTrans_dTCanopy => deriv_data%var(iLookDERIV%dAquiferTrans_dTCanopy )%dat(1) ,& ! intent(in): derivatives in the aquifer transpiration flux w.r.t. canopy temperature
+ dAquiferTrans_dTGround => deriv_data%var(iLookDERIV%dAquiferTrans_dTGround )%dat(1) ,& ! intent(in): derivatives in the aquifer transpiration flux w.r.t. ground temperature
+ dAquiferTrans_dCanWat => deriv_data%var(iLookDERIV%dAquiferTrans_dCanWat )%dat(1) ,& ! intent(in): derivatives in the aquifer transpiration flux w.r.t. canopy total water
+ ! derivative in liquid water fluxes at the interface of snow layers w.r.t. volumetric liquid water content in the layer above
+ iLayerLiqFluxSnowDeriv => deriv_data%var(iLookDERIV%iLayerLiqFluxSnowDeriv )%dat ,& ! intent(in): [dp(:)] derivative in vertical liquid water flux at layer interfaces
+ ! derivative in liquid water fluxes for the soil domain w.r.t hydrology state variables
+ dVolTot_dPsi0 => deriv_data%var(iLookDERIV%dVolTot_dPsi0 )%dat ,& ! intent(in): [dp(:)] derivative in total water content w.r.t. total water matric potential
+ dCompress_dPsi => deriv_data%var(iLookDERIV%dCompress_dPsi )%dat ,& ! intent(in): [dp(:)] derivative in compressibility w.r.t matric head
+ dq_dHydStateAbove => deriv_data%var(iLookDERIV%dq_dHydStateAbove )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer above
+ dq_dHydStateBelow => deriv_data%var(iLookDERIV%dq_dHydStateBelow )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer below
+ dq_dHydStateLayerSurfVec => deriv_data%var(iLookDERIV%dq_dHydStateLayerSurfVec )%dat ,& ! intent(in): [dp(:)] change in the flux in soil surface interface w.r.t. state variables in layers
+ ! derivative in baseflow flux w.r.t. aquifer storage
+ dBaseflow_dAquifer => deriv_data%var(iLookDERIV%dBaseflow_dAquifer )%dat(1) ,& ! intent(in): [dp(:)] erivative in baseflow flux w.r.t. aquifer storage (s-1)
+ ! derivative in liquid water fluxes for the soil domain w.r.t energy state variables
+ dq_dNrgStateAbove => deriv_data%var(iLookDERIV%dq_dNrgStateAbove )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer above
+ dq_dNrgStateBelow => deriv_data%var(iLookDERIV%dq_dNrgStateBelow )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer below
+ dq_dNrgStateLayerSurfVec => deriv_data%var(iLookDERIV%dq_dNrgStateLayerSurfVec )%dat ,& ! intent(in): [dp(:)] change in the flux in soil surface interface w.r.t. state variables in layers
+ ! derivative in liquid water fluxes for the soil and snow domain w.r.t temperature
+ mLayerdTheta_dTk => deriv_data%var(iLookDERIV%mLayerdTheta_dTk )%dat ,& ! intent(in): [dp(:)] derivative of volumetric liquid water content w.r.t. temperature
+ ! derivative in bulk heat capacity w.r.t. relevant state variables
+ dVolHtCapBulk_dPsi0 => deriv_data%var(iLookDERIV%dVolHtCapBulk_dPsi0 )%dat ,& ! intent(in): [dp(:)] derivative in bulk heat capacity w.r.t. matric potential
+ dVolHtCapBulk_dTheta => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTheta )%dat ,& ! intent(in): [dp(:)] derivative in bulk heat capacity w.r.t. volumetric water content
+ dVolHtCapBulk_dCanWat => deriv_data%var(iLookDERIV%dVolHtCapBulk_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in bulk heat capacity w.r.t. volumetric water content
+ dVolHtCapBulk_dTk => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTk )%dat ,& ! intent(in): [dp(:)] derivative in bulk heat capacity w.r.t. temperature
+ dVolHtCapBulk_dTkCanopy => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTkCanopy )%dat(1) ,& ! intent(in): [dp] derivative in bulk heat capacity w.r.t. temperature
+ ! derivatives in time
+ mLayerdTemp_dt => deriv_data%var(iLookDERIV%mLayerdTemp_dt )%dat ,& ! intent(in): [dp(:)] timestep change in layer temperature
+ scalarCanopydTemp_dt => deriv_data%var(iLookDERIV%scalarCanopydTemp_dt )%dat(1) ,& ! intent(in): [dp ] timestep change in canopy temperature
+ ! diagnostic variables
+ scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg)%dat(1) ,& ! intent(in): [dp] fraction of liquid water on vegetation (-)
+ scalarBulkVolHeatCapVeg => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg)%dat(1) ,& ! intent(in): [dp] bulk volumetric heat capacity of vegetation (J m-3 K-1)
+ mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow)%dat ,& ! intent(in): [dp(:)] fraction of liquid water in each snow layer (-)
+ mLayerVolHtCapBulk => diag_data%var(iLookDIAG%mLayerVolHtCapBulk)%dat ,& ! intent(in): [dp(:)] bulk volumetric heat capacity in each snow and soil layer (J m-3 K-1)
+ scalarSoilControl => diag_data%var(iLookDIAG%scalarSoilControl)%dat(1) ,& ! intent(in): [dp] soil control on infiltration, zero or one
+ ! canopy and layer depth
+ canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp ] canopy depth (m)
+ mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat & ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
+ ) ! making association with data in structures
+ ! --------------------------------------------------------------
+ ! initialize error control
+ err=0; message='computJacob/'
+
+ ! *********************************************************************************************************************************************************
+ ! *********************************************************************************************************************************************************
+ ! * PART 0: PRELIMINARIES (INITIALIZE JACOBIAN AND COMPUTE TIME-VARIABLE DIAGONAL TERMS)
+ ! *********************************************************************************************************************************************************
+ ! *********************************************************************************************************************************************************
+
+ ! get the number of state variables
+ nState = size(dMat)
+
+ ! initialize the Jacobian
+ ! NOTE: this needs to be done every time, since Jacobian matrix is modified in the solver
+ aJac(:,:) = 0._rkind ! analytical Jacobian matrix
+
+ ! compute terms in the Jacobian for vegetation (excluding fluxes)
+ ! NOTE: energy for vegetation is computed *within* the iteration loop as it includes phase change
+ if(ixVegNrg/=integerMissing)then
+ dMat(ixVegNrg) = scalarBulkVolHeatCapVeg + LH_fus*iden_water*dTheta_dTkCanopy &
+ + dVolHtCapBulk_dTkCanopy * scalarCanopydTemp_dt
endif
-
- end do ! (looping through energy states in the snow+soil domain)
- endif ! (if the subset includes energy state variables in the snow+soil domain)
-
- ! -----
- ! * liquid water fluxes for the snow domain...
- ! --------------------------------------------
- if(nSnowOnlyHyd>0)then
- do iLayer=1,nSnow ! loop through layers in the snow domain
-
- ! - check that the snow layer is desired
- if(ixSnowOnlyHyd(iLayer)==integerMissing) cycle
-
- ! - define state indices for the current layer
- watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset
-
- ! compute factor to convert liquid water derivative to total water derivative
- select case( ixHydType(iLayer) )
- case(iname_watLayer); convLiq2tot = mLayerFracLiqSnow(iLayer)
- case default; convLiq2tot = 1._dp
- end select
-
- ! - diagonal elements
- aJac(ixDiag,watState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot + dMat(watState)
-
- ! - lower-diagonal elements
- if(iLayer > 1)then
- if(ixSnowOnlyHyd(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSnowOnlyHyd(iLayer-1),watState),watState) = 0._dp ! sub-diagonal: no dependence on other layers
- endif
-
- ! - upper diagonal elements
- if(iLayer < nSnow)then
- if(ixSnowOnlyHyd(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSnowOnlyHyd(iLayer+1),watState),watState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot ! dVol(below)/dLiq(above) -- (-)
- endif
-
- ! - compute cross-derivative terms for energy
- ! NOTE: increase in volumetric liquid water content balanced by a decrease in volumetric ice content
- if(nSnowOnlyNrg>0)then
-
- ! (define the energy state)
- nrgState = ixSnowOnlyNrg(iLayer) ! index within the full state vector
- if(nrgstate/=integerMissing)then ! (energy state for the current layer is within the state subset)
-
- ! (cross-derivative terms for the current layer)
- aJac(ixOffDiag(nrgState,watState),watState) = -(1._dp - mLayerFracLiqSnow(iLayer))*LH_fus*iden_water ! (dF/dLiq)
- aJac(ixOffDiag(watState,nrgState),nrgState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! (dVol/dT)
-
- ! (cross-derivative terms for the layer below)
- if(iLayer < nSnow)then
- aJac(ixOffDiag(ixSnowOnlyHyd(iLayer+1),nrgState),nrgState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! dVol(below)/dT(above) -- K-1
- endif ! (if there is a water state in the layer below the current layer in the given state subset)
-
- endif ! (if the energy state for the current layer is within the state subset)
- endif ! (if state variables exist for energy in snow+soil layers)
-
- end do ! (looping through liquid water states in the snow domain)
- endif ! (if the subset includes hydrology state variables in the snow domain)
-
- ! -----
- ! * liquid water fluxes for the soil domain...
- ! --------------------------------------------
- if(nSoilOnlyHyd>0)then
- do iLayer=1,nSoil
-
- ! - check that the soil layer is desired
- if(ixSoilOnlyHyd(iLayer)==integerMissing) cycle
-
- ! - define state indices
- watState = ixSoilOnlyHyd(iLayer) ! hydrology state index within the state subset
-
- ! - define indices of the soil layers
- jLayer = iLayer+nSnow ! index of layer in the snow+soil vector
-
- ! - compute the diagonal elements
- ! all terms *excluding* baseflow
- aJac(ixDiag,watState) = (dt/mLayerDepth(jLayer))*(-dq_dHydStateBelow(iLayer-1) + dq_dHydStateAbove(iLayer)) + dMat(watState)
-
- ! - compute the lower-diagonal elements
- if(iLayer > 1)then
- if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(iLayer-1),watState),watState) = (dt/mLayerDepth(jLayer-1))*( dq_dHydStateBelow(iLayer-1))
- endif
-
- ! - compute the upper-diagonal elements
- if(iLayer<nSoil)then
- if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(iLayer+1),watState),watState) = (dt/mLayerDepth(jLayer+1))*(-dq_dHydStateAbove(iLayer))
- endif
-
- end do ! (looping through hydrology states in the soil domain)
- endif ! (if the subset includes hydrology state variables in the soil domain)
-
- ! -----
- ! * liquid water fluxes for the aquifer...
- ! ----------------------------------------
- if(ixAqWat/=integerMissing) aJac(ixDiag,ixAqWat) = -dBaseflow_dAquifer*dt + dMat(ixAqWat)
-
- ! -----
- ! * derivative in liquid water fluxes w.r.t. temperature for the soil domain...
- ! -----------------------------------------------------------------------------
- if(nSoilOnlyHyd>0 .and. nSoilOnlyNrg>0)then
- do iLayer=1,nSoilOnlyHyd
-
- ! - check that the soil layer is desired
- if(ixSoilOnlyHyd(iLayer)==integerMissing) cycle
-
- ! - define index of hydrology state variable within the state subset
- watState = ixSoilOnlyHyd(iLayer)
-
- ! - define indices of the soil layers
- jLayer = iLayer+nSnow ! index of layer in the snow+soil vector
-
- ! - define the energy state variable
- nrgState = ixNrgLayer(jLayer) ! index within the full state vector
-
- ! only compute derivatives if the energy state for the current layer is within the state subset
- if(nrgstate/=integerMissing)then
-
- ! - compute the Jacobian for the layer itself
- aJac(ixOffDiag(watState,nrgState),nrgState) = (dt/mLayerDepth(jLayer))*(-dq_dNrgStateBelow(iLayer-1) + dq_dNrgStateAbove(iLayer)) ! dVol/dT (K-1) -- flux depends on ice impedance
-
- ! - include derivatives w.r.t. ground evaporation
- if(nSnow==0 .and. iLayer==1)then ! upper-most soil layer
- if(computeVegFlux)then
- aJac(ixOffDiag(watState,ixVegHyd),ixVegHyd) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dCanWat/iden_water) ! dVol/dLiq (kg m-2)-1
- aJac(ixOffDiag(watState,ixCasNrg),ixCasNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanair/iden_water) ! dVol/dT (K-1)
- aJac(ixOffDiag(watState,ixVegNrg),ixVegNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanopy/iden_water) ! dVol/dT (K-1)
- endif
- aJac(ixOffDiag(watState,ixTopNrg),ixTopNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTGround/iden_water) + aJac(ixOffDiag(watState,ixTopNrg),ixTopNrg) ! dVol/dT (K-1)
+
+ ! compute additional terms for the Jacobian for the snow-soil domain (excluding fluxes)
+ ! NOTE: energy for snow+soil is computed *within* the iteration loop as it includes phase change
+ do iLayer=1,nLayers
+ if(ixSnowSoilNrg(iLayer)/=integerMissing)then
+ dMat(ixSnowSoilNrg(iLayer)) = mLayerVolHtCapBulk(iLayer) + LH_fus*iden_water*mLayerdTheta_dTk(iLayer) &
+ + dVolHtCapBulk_dTk(iLayer) * mLayerdTemp_dt(iLayer)
endif
-
- ! melt-freeze: compute derivative in energy with respect to mass
- if(mLayerdTheta_dTk(jLayer) > verySmall)then ! ice is present
- aJac(ixOffDiag(nrgState,watState),watState) = -dVolTot_dPsi0(iLayer)*LH_fus*iden_water ! dNrg/dMat (J m-3 m-1) -- dMat changes volumetric water, and hence ice content
- else
- aJac(ixOffDiag(nrgState,watState),watState) = 0._dp
+ end do
+
+ ! compute additional terms for the Jacobian for the soil domain (excluding fluxes)
+ do iLayer=1,nSoil
+ if(ixSoilOnlyHyd(iLayer)/=integerMissing)then
+ dMat(ixSoilOnlyHyd(iLayer)) = dVolTot_dPsi0(iLayer) + dCompress_dPsi(iLayer)
endif
-
- ! - compute lower diagonal elements
- if(iLayer>1)then
- if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(iLayer-1),nrgState),nrgState) = (dt/mLayerDepth(jLayer-1))*( dq_dNrgStateBelow(iLayer-1)) ! K-1
- endif
-
- ! compute upper-diagonal elements
- if(iLayer<nSoil)then
- if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(iLayer+1),nrgState),nrgState) = (dt/mLayerDepth(jLayer+1))*(-dq_dNrgStateAbove(iLayer)) ! K-1
- endif
-
- endif ! (if the energy state for the current layer is within the state subset)
-
- end do ! (looping through soil layers)
- endif ! (if there are state variables for both water and energy in the soil domain)
-
- if(globalPrintFlag)then
- print*, '** banded analytical Jacobian:'
- write(*,'(a4,1x,100(i17,1x))') 'xCol', (iLayer, iLayer=min(iJac1,nState),min(iJac2,nState))
- do iLayer=kl+1,nBands
- write(*,'(i4,1x,100(e17.10,1x))') iLayer, (aJac(iLayer,jLayer),jLayer=min(iJac1,nState),min(iJac2,nState))
- end do
- end if
- !print*, 'PAUSE: banded analytical Jacobian'; read(*,*)
-
- ! *********************************************************************************************************************************************************
- ! *********************************************************************************************************************************************************
- ! * PART 2: FULL MATRIX
- ! *********************************************************************************************************************************************************
- ! *********************************************************************************************************************************************************
- case(ixFullMatrix)
-
- ! check
- if(size(aJac,1)/=size(dMat) .or. size(aJac,2)/=size(dMat))then
- message=trim(message)//'unexpected shape of the Jacobian matrix: expect aJac(nState,nState)'
- err=20; return
- end if
-
- ! -----
- ! * energy and liquid fluxes over vegetation...
- ! ---------------------------------------------
- if(computeVegFlux)then ! (derivatives only defined when vegetation protrudes over the surface)
-
- ! * liquid water fluxes for vegetation canopy (-)
- if(ixVegHyd/=integerMissing) aJac(ixVegHyd,ixVegHyd) = -scalarFracLiqVeg*(dCanopyEvaporation_dCanWat - scalarCanopyLiqDeriv)*dt + 1._dp
-
- ! * cross-derivative terms for canopy water
- if(ixVegHyd/=integerMissing)then
- ! cross-derivative terms w.r.t. system temperatures (kg m-2 K-1)
- if(ixCasNrg/=integerMissing) aJac(ixVegHyd,ixCasNrg) = -dCanopyEvaporation_dTCanair*dt
- if(ixVegNrg/=integerMissing) aJac(ixVegHyd,ixVegNrg) = -dCanopyEvaporation_dTCanopy*dt + dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy
- if(ixTopNrg/=integerMissing) aJac(ixVegHyd,ixTopNrg) = -dCanopyEvaporation_dTGround*dt
- ! cross-derivative terms w.r.t. canopy water (kg-1 m2)
- if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegHyd) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarFracLiqVeg*scalarCanopyLiqDeriv)/iden_water
- ! cross-derivative terms w.r.t. canopy liquid water (J m-1 kg-1)
- ! NOTE: dIce/dLiq = (1 - scalarFracLiqVeg); dIce*LH_fus/canopyDepth = J m-3; dLiq = kg m-2
- if(ixVegNrg/=integerMissing) aJac(ixVegNrg,ixVegHyd) = (dt/canopyDepth) *(-dCanopyNetFlux_dCanWat) - (1._dp - scalarFracLiqVeg)*LH_fus/canopyDepth ! dF/dLiq
- if(ixTopNrg/=integerMissing) aJac(ixTopNrg,ixVegHyd) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanWat)
- endif
-
- ! cross-derivative terms w.r.t. canopy temperature (K-1)
- if(ixVegNrg/=integerMissing)then
- if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegNrg) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarCanopyLiqDeriv*dCanLiq_dTcanopy)/iden_water
- endif
-
- ! energy fluxes with the canopy air space (J m-3 K-1)
- if(ixCasNrg/=integerMissing)then
- aJac(ixCasNrg,ixCasNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanairTemp) + dMat(ixCasNrg)
- if(ixVegNrg/=integerMissing) aJac(ixCasNrg,ixVegNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanopyTemp)
- if(ixTopNrg/=integerMissing) aJac(ixCasNrg,ixTopNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dGroundTemp)
- endif
-
- ! energy fluxes with the vegetation canopy (J m-3 K-1)
- if(ixVegNrg/=integerMissing)then
- if(ixCasNrg/=integerMissing) aJac(ixVegNrg,ixCasNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanairTemp)
- aJac(ixVegNrg,ixVegNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanopyTemp) + dMat(ixVegNrg)
- if(ixTopNrg/=integerMissing) aJac(ixVegNrg,ixTopNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dGroundTemp)
- endif
-
- ! energy fluxes with the surface (J m-3 K-1)
- if(ixTopNrg/=integerMissing)then
- if(ixCasNrg/=integerMissing) aJac(ixTopNrg,ixCasNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanairTemp)
- if(ixVegNrg/=integerMissing) aJac(ixTopNrg,ixVegNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanopyTemp)
- endif
-
- endif ! if there is a need to compute energy fluxes within vegetation
-
- ! -----
- ! * energy fluxes for the snow+soil domain...
- ! -------------------------------------------
- if(nSnowSoilNrg>0)then
- do iLayer=1,nLayers ! loop through all layers in the snow+soil domain
-
- ! check if the state is in the subset
- if(ixSnowSoilNrg(iLayer)==integerMissing) cycle
-
- ! - define index within the state subset and the full state vector
- jState = ixSnowSoilNrg(iLayer) ! index within the state subset
-
- ! - diagonal elements
- aJac(jState,jState) = (dt/mLayerDepth(iLayer))*(-dNrgFlux_dTempBelow(iLayer-1) + dNrgFlux_dTempAbove(iLayer)) + dMat(jState)
-
- ! - lower-diagonal elements
- if(iLayer > 1)then
- if(ixSnowSoilNrg(iLayer-1)/=integerMissing) aJac(ixSnowSoilNrg(iLayer-1),jState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dTempBelow(iLayer-1) )
- endif
-
- ! - upper diagonal elements
- if(iLayer < nLayers)then
- if(ixSnowSoilNrg(iLayer+1)/=integerMissing) aJac(ixSnowSoilNrg(iLayer+1),jState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dTempAbove(iLayer ) )
- endif
-
- end do ! (looping through energy states in the snow+soil domain)
- endif ! (if the subset includes energy state variables in the snow+soil domain)
-
- ! -----
- ! * liquid water fluxes for the snow domain...
- ! --------------------------------------------
- if(nSnowOnlyHyd>0)then
- do iLayer=1,nSnow ! loop through layers in the snow domain
-
- ! - check that the snow layer is desired
- if(ixSnowOnlyHyd(iLayer)==integerMissing) cycle
-
- ! - define state indices for the current layer
- watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset
-
- ! compute factor to convert liquid water derivative to total water derivative
- select case( ixHydType(iLayer) )
- case(iname_watLayer); convLiq2tot = mLayerFracLiqSnow(iLayer)
- case default; convLiq2tot = 1._dp
- end select
-
- ! - diagonal elements
- aJac(watState,watState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot + dMat(watState)
-
- ! - lower-diagonal elements
- if(iLayer > 1)then
- if(ixSnowOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSnowOnlyHyd(iLayer-1),watState) = 0._dp ! sub-diagonal: no dependence on other layers
- endif
-
- ! - upper diagonal elements
- if(iLayer < nSnow)then
- if(ixSnowOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSnowOnlyHyd(iLayer+1),watState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot ! dVol(below)/dLiq(above) -- (-)
- endif
-
- ! - compute cross-derivative terms for energy
- ! NOTE: increase in volumetric liquid water content balanced by a decrease in volumetric ice content
- if(nSnowOnlyNrg>0)then
-
- ! (define the energy state)
- nrgState = ixSnowOnlyNrg(iLayer) ! index within the full state vector
- if(nrgstate/=integerMissing)then ! (energy state for the current layer is within the state subset)
-
- ! (cross-derivative terms for the current layer)
- aJac(nrgState,watState) = -(1._dp - mLayerFracLiqSnow(iLayer))*LH_fus*iden_water ! (dF/dLiq)
- aJac(watState,nrgState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! (dVol/dT)
-
- ! (cross-derivative terms for the layer below)
- if(iLayer < nSnow)then
- aJac(ixSnowOnlyHyd(iLayer+1),nrgState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! dVol(below)/dT(above) -- K-1
- endif ! (if there is a water state in the layer below the current layer in the given state subset)
-
- endif ! (if the energy state for the current layer is within the state subset)
- endif ! (if state variables exist for energy in snow+soil layers)
-
- end do ! (looping through liquid water states in the snow domain)
- endif ! (if the subset includes hydrology state variables in the snow domain)
-
- ! -----
- ! * liquid water fluxes for the soil domain...
- ! --------------------------------------------
- if(nSoilOnlyHyd>0)then
-
- do iLayer=1,nSoil
-
- ! - check that the soil layer is desired
- if(ixSoilOnlyHyd(iLayer)==integerMissing) cycle
-
- ! - define state indices
- watState = ixSoilOnlyHyd(iLayer) ! hydrology state index within the state subset
-
- ! - define indices of the soil layers
- jLayer = iLayer+nSnow ! index of layer in the snow+soil vector
-
- ! - compute the diagonal elements
- ! all terms *excluding* baseflow
- aJac(watState,watState) = (dt/mLayerDepth(jLayer))*(-dq_dHydStateBelow(iLayer-1) + dq_dHydStateAbove(iLayer)) + dMat(watState)
-
- ! - compute the lower-diagonal elements
- if(iLayer > 1)then
- if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer-1),watState) = (dt/mLayerDepth(jLayer-1))*( dq_dHydStateBelow(iLayer-1))
- endif
-
- ! - compute the upper-diagonal elements
- if(iLayer<nSoil)then
- if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer+1),watState) = (dt/mLayerDepth(jLayer+1))*(-dq_dHydStateAbove(iLayer))
- endif
-
- ! - include terms for baseflow
- if(computeBaseflow .and. nSoilOnlyHyd==nSoil)then
- do pLayer=1,nSoil
- qState = ixSoilOnlyHyd(pLayer) ! hydrology state index within the state subset
- aJac(watState,qState) = aJac(watState,qState) + (dt/mLayerDepth(jLayer))*dBaseflow_dMatric(iLayer,pLayer)
- end do
- endif
-
- end do ! (looping through hydrology states in the soil domain)
- endif ! (if the subset includes hydrology state variables in the soil domain)
-
- ! -----
- ! * liquid water fluxes for the aquifer...
- ! ----------------------------------------
- if(ixAqWat/=integerMissing) aJac(ixAqWat,ixAqWat) = -dBaseflow_dAquifer*dt + dMat(ixAqWat)
-
- ! -----
- ! * derivative in liquid water fluxes w.r.t. temperature for the soil domain...
- ! -----------------------------------------------------------------------------
- if(nSoilOnlyHyd>0 .and. nSoilOnlyNrg>0)then
- do iLayer=1,nSoilOnlyHyd
-
- ! - check that the soil layer is desired
- if(ixSoilOnlyHyd(iLayer)==integerMissing) cycle
-
- ! - define index of hydrology state variable within the state subset
- watState = ixSoilOnlyHyd(iLayer)
-
- ! - define indices of the soil layers
- jLayer = iLayer+nSnow ! index of layer in the snow+soil vector
-
- ! - define the energy state variable
- nrgState = ixNrgLayer(jLayer) ! index within the full state vector
-
- ! only compute derivatives if the energy state for the current layer is within the state subset
- if(nrgstate/=integerMissing)then
-
- ! - compute the Jacobian for the layer itself
- aJac(watState,nrgState) = (dt/mLayerDepth(jLayer))*(-dq_dNrgStateBelow(iLayer-1) + dq_dNrgStateAbove(iLayer)) ! dVol/dT (K-1) -- flux depends on ice impedance
-
- ! - include derivatives w.r.t. ground evaporation
- if(nSnow==0 .and. iLayer==1)then ! upper-most soil layer
- if(computeVegFlux)then
- aJac(watState,ixVegHyd) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dCanWat/iden_water) ! dVol/dLiq (kg m-2)-1
- aJac(watState,ixCasNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanair/iden_water) ! dVol/dT (K-1)
- aJac(watState,ixVegNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanopy/iden_water) ! dVol/dT (K-1)
+ end do
+
+ ! define the form of the matrix
+ select case(ixMatrix)
+
+ ! *********************************************************************************************************************************************************
+ ! *********************************************************************************************************************************************************
+ ! * PART 1: BAND MATRIX
+ ! *********************************************************************************************************************************************************
+ ! *********************************************************************************************************************************************************
+ case(ixBandMatrix)
+
+ ! check
+ if(size(aJac,1)/=nBands .or. size(aJac,2)/=size(dMat))then
+ message=trim(message)//'unexpected shape of the Jacobian matrix: expect aJac(nBands,nState)'
+ err=20; return
endif
- aJac(watState,ixTopNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTGround/iden_water) + aJac(watState,ixTopNrg) ! dVol/dT (K-1)
- endif
-
- ! melt-freeze: compute derivative in energy with respect to mass
- if(mLayerdTheta_dTk(jLayer) > verySmall)then ! ice is present
- aJac(nrgState,watState) = -dVolTot_dPsi0(iLayer)*LH_fus*iden_water ! dNrg/dMat (J m-3 m-1) -- dMat changes volumetric water, and hence ice content
- else
- aJac(nrgState,watState) = 0._dp
- endif
-
- ! - compute lower diagonal elements
- if(iLayer>1)then
- if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer-1),nrgState) = (dt/mLayerDepth(jLayer-1))*( dq_dNrgStateBelow(iLayer-1)) ! K-1
- endif
-
- ! compute upper-diagonal elements
- if(iLayer<nSoil)then
- if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer+1),nrgState) = (dt/mLayerDepth(jLayer+1))*(-dq_dNrgStateAbove(iLayer)) ! K-1
- endif
-
- endif ! (if the energy state for the current layer is within the state subset)
-
- end do ! (looping through soil layers)
- endif ! (if there are state variables for both water and energy in the soil domain)
-
- ! print the Jacobian
- if(globalPrintFlag)then
- print*, '** analytical Jacobian (full):'
- write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=min(iJac1,nState),min(iJac2,nState))
- do iLayer=min(iJac1,nState),min(iJac2,nState)
- write(*,'(i4,1x,100(e12.5,1x))') iLayer, aJac(min(iJac1,nState):min(iJac2,nState),iLayer)
- end do
- end if
-
- ! ***
- ! check
- case default; err=20; message=trim(message)//'unable to identify option for the type of matrix'; return
-
- end select ! type of matrix
-
- ! end association to variables in the data structures
- end associate
-
- end subroutine computJacob
-
-
- ! **********************************************************************************************************
- ! private function: get the off-diagonal index in the band-diagonal matrix
- ! **********************************************************************************************************
- function ixOffDiag(jState,iState)
- implicit none
- integer(i4b),intent(in) :: jState ! off-diagonal state
- integer(i4b),intent(in) :: iState ! diagonal state
- integer(i4b) :: ixOffDiag ! off-diagonal index in gthe band-diagonal matrix
- ixOffDiag = ixBandOffset + jState - iState
- end function ixOffDiag
-
-end module computJacob_module
+
+ ! -----
+ ! * energy and liquid fluxes over vegetation...
+ ! ---------------------------------------------
+ if(computeVegFlux)then ! (derivatives only defined when vegetation protrudes over the surface)
+
+ ! * energy fluxes with the canopy water
+ if(ixVegHyd/=integerMissing)then
+
+ ! * cross-derivative terms w.r.t. system temperatures (kg m-2 K-1)
+ if(ixCasNrg/=integerMissing) aJac(ixOffDiag(ixVegHyd,ixCasNrg),ixCasNrg) = -dCanopyEvaporation_dTCanair*dt
+ ! dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy is the derivative in throughfall and canopy drainage with canopy temperature
+ if(ixVegNrg/=integerMissing) aJac(ixOffDiag(ixVegHyd,ixVegNrg),ixVegNrg) = -dCanopyEvaporation_dTCanopy*dt + dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy
+ ! * liquid water fluxes for vegetation canopy (-), dt*scalarFracLiqVeg*scalarCanopyLiqDeriv is the derivative in throughfall and canopy drainage with canopy water
+ aJac(ixDiag, ixVegHyd) = -scalarFracLiqVeg*(dCanopyEvaporation_dCanWat - scalarCanopyLiqDeriv)*dt + 1._rkind
+ if(ixTopNrg/=integerMissing) aJac(ixOffDiag(ixVegHyd,ixTopNrg),ixTopNrg) = -dCanopyEvaporation_dTGround*dt
+
+ ! * cross-derivative terms w.r.t. canopy water (kg-1 m2)
+ if(ixTopHyd/=integerMissing) aJac(ixOffDiag(ixTopHyd,ixVegHyd),ixVegHyd) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarFracLiqVeg*scalarCanopyLiqDeriv)/iden_water
+
+ ! * cross-derivative terms w.r.t. canopy liquid water (J m-1 kg-1)
+ ! NOTE: dIce/dLiq = (1 - scalarFracLiqVeg); dIce*LH_fus/canopyDepth = J m-3; dLiq = kg m-2
+ if(ixVegNrg/=integerMissing) aJac(ixOffDiag(ixVegNrg,ixVegHyd),ixVegHyd) = (-1._rkind + scalarFracLiqVeg)*LH_fus/canopyDepth &
+ + dVolHtCapBulk_dCanWat * scalarCanopydTemp_dt - (dt/canopyDepth) * dCanopyNetFlux_dCanWat ! dF/dLiq
+ if(ixTopNrg/=integerMissing) aJac(ixOffDiag(ixTopNrg,ixVegHyd),ixVegHyd) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanWat)
+ endif
+
+ ! * cross-derivative terms w.r.t. canopy temperature (K-1)
+ if(ixVegNrg/=integerMissing)then
+ if(ixTopHyd/=integerMissing) aJac(ixOffDiag(ixTopHyd,ixVegNrg),ixVegNrg) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarCanopyLiqDeriv*dCanLiq_dTcanopy)/iden_water
+ endif
+
+ ! * energy fluxes with the canopy air space (J m-3 K-1)
+ if(ixCasNrg/=integerMissing)then
+ aJac(ixDiag, ixCasNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanairTemp) + dMat(ixCasNrg)
+ if(ixVegNrg/=integerMissing) aJac(ixOffDiag(ixCasNrg,ixVegNrg),ixVegNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanopyTemp)
+ if(ixTopNrg/=integerMissing) aJac(ixOffDiag(ixCasNrg,ixTopNrg),ixTopNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dGroundTemp)
+ endif
+
+ ! * energy fluxes with the vegetation canopy (J m-3 K-1)
+ if(ixVegNrg/=integerMissing)then
+ if(ixCasNrg/=integerMissing) aJac(ixOffDiag(ixVegNrg,ixCasNrg),ixCasNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanairTemp)
+ aJac(ixDiag, ixVegNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanopyTemp) + dMat(ixVegNrg)
+ if(ixTopNrg/=integerMissing) aJac(ixOffDiag(ixVegNrg,ixTopNrg),ixTopNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dGroundTemp)
+ endif
+
+ ! * energy fluxes with the surface (J m-3 K-1)
+ if(ixTopNrg/=integerMissing)then
+ if(ixCasNrg/=integerMissing) aJac(ixOffDiag(ixTopNrg,ixCasNrg),ixCasNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanairTemp)
+ if(ixVegNrg/=integerMissing) aJac(ixOffDiag(ixTopNrg,ixVegNrg),ixVegNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanopyTemp)
+ endif
+
+ endif ! if there is a need to compute energy fluxes within vegetation
+
+ ! -----
+ ! * energy fluxes for the snow+soil domain...
+ ! -------------------------------------------
+ if(nSnowSoilNrg>0)then
+ do iLayer=1,nLayers ! loop through all layers in the snow+soil domain
+
+ ! check if the state is in the subset
+ if(ixSnowSoilNrg(iLayer)==integerMissing) cycle
+
+ ! - define index within the state subset and the full state vector
+ jState = ixSnowSoilNrg(iLayer) ! index within the state subset
+
+ ! - diagonal elements
+ aJac(ixDiag,jState) = (dt/mLayerDepth(iLayer))*(-dNrgFlux_dTempBelow(iLayer-1) + dNrgFlux_dTempAbove(iLayer)) + dMat(jState)
+
+ ! - lower-diagonal elements
+ if(iLayer>1)then
+ if(ixSnowSoilNrg(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSnowSoilNrg(iLayer-1),jState),jState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dTempBelow(iLayer-1) )
+ endif
+
+ ! - upper diagonal elements
+ if(iLayer<nLayers)then
+ if(ixSnowSoilNrg(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSnowSoilNrg(iLayer+1),jState),jState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dTempAbove(iLayer ) )
+ endif
+
+ end do ! (looping through energy states in the snow+soil domain)
+ endif ! (if the subset includes energy state variables in the snow+soil domain)
+
+ ! -----
+ ! * liquid water fluxes for the snow domain...
+ ! --------------------------------------------
+ if(nSnowOnlyHyd>0)then
+ do iLayer=1,nSnow ! loop through layers in the snow domain
+
+ ! - check that the snow layer is desired
+ if(ixSnowOnlyHyd(iLayer)==integerMissing) cycle
+
+ ! - define state indices for the current layer
+ watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset
+
+ ! compute factor to convert liquid water derivative to total water derivative
+ select case( ixHydType(iLayer) )
+ case(iname_watLayer); convLiq2tot = mLayerFracLiqSnow(iLayer)
+ case default; convLiq2tot = 1._rkind
+ end select
+
+ ! - diagonal elements
+ aJac(ixDiag,watState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot + dMat(watState)
+
+ ! - lower-diagonal elements
+ if(iLayer>1)then
+ if(ixSnowOnlyHyd(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSnowOnlyHyd(iLayer-1),watState),watState) = 0._rkind ! sub-diagonal: no dependence on other layers
+ endif
+
+ ! - upper diagonal elements
+ if(iLayer<nSnow)then
+ if(ixSnowOnlyHyd(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSnowOnlyHyd(iLayer+1),watState),watState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot ! dVol(below)/dLiq(above) -- (-)
+ endif
+
+ end do ! (looping through liquid water states in the snow domain)
+ endif ! (if the subset includes hydrology state variables in the snow domain)
+
+ ! -----
+ ! * cross derivatives in the snow domain...
+ ! ----------------------------------------
+ if(nSnowOnlyHyd>0 .and. nSnowOnlyNrg>0)then
+ do iLayer=1,nSnow ! loop through layers in the snow domain
+
+ ! - check that the snow layer is desired
+ if(ixSnowOnlyNrg(iLayer)==integerMissing) cycle
+
+ ! (define the energy state)
+ nrgState = ixSnowOnlyNrg(iLayer) ! index within the full state vector
+
+ ! - define state indices for the current layer
+ watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset
+
+ if(watstate/=integerMissing)then ! (energy state for the current layer is within the state subset)
+
+ ! - include derivatives of energy fluxes w.r.t water fluxes for current layer
+ aJac(ixOffDiag(nrgState,watState),watState) = (-1._rkind + mLayerFracLiqSnow(iLayer))*LH_fus*iden_water &
+ + dVolHtCapBulk_dTheta(iLayer) * mLayerdTemp_dt(iLayer) &
+ + (dt/mLayerDepth(iLayer))*(-dNrgFlux_dWatBelow(iLayer-1) + dNrgFlux_dWatAbove(iLayer)) ! (dF/dLiq)
+
+ ! - include derivatives of water fluxes w.r.t energy fluxes for current layer
+ aJac(ixOffDiag(watState,nrgState),nrgState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! (dVol/dT)
+
+ ! (cross-derivative terms for the layer below)
+ if(iLayer < nSnow)then
+ aJac(ixOffDiag(ixSnowOnlyHyd(iLayer+1),nrgState),nrgState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! dVol(below)/dT(above) -- K-1
+ endif ! (if there is a water state in the layer below the current layer in the given state subset)
+
+ ! - include derivatives of heat capacity w.r.t water fluxes for surrounding layers starting with layer above
+ if(iLayer>1)then
+ if(ixSnowOnlyNrg(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSnowOnlyNrg(iLayer-1),watState),watState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dWatBelow(iLayer-1) )
+ endif
+
+ ! (cross-derivative terms for the layer below)
+ if(iLayer<nSnow)then
+ if(ixSnowOnlyNrg(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSnowOnlyNrg(iLayer+1),watState),watState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dWatAbove(iLayer ) )
+ elseif(iLayer==nSnow .and. nSoilOnlyNrg>0)then !bottom snow layer and there is soil below
+ if(ixSoilOnlyNrg(1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyNrg(1),watState),watState) = (dt/mLayerDepth(nSnow+1))*(-dNrgFlux_dWatAbove(nSnow) )
+ endif
+
+ endif ! (if the energy state for the current layer is within the state subset)
+
+ end do ! (looping through snow layers)
+ endif ! (if there are state variables for both water and energy in the snow domain)
+
+ ! -----
+ ! * liquid water fluxes for the soil domain...
+ ! --------------------------------------------
+ if(nSoilOnlyHyd>0)then
+ do iLayer=1,nSoil
+
+ ! - check that the soil layer is desired
+ if(ixSoilOnlyHyd(iLayer)==integerMissing) cycle
+
+ ! - define state indices
+ watState = ixSoilOnlyHyd(iLayer) ! hydrology state index within the state subset
+
+ ! - define indices of the soil layers
+ jLayer = iLayer+nSnow ! index of layer in the snow+soil vector
+
+ ! - compute the diagonal elements
+ ! all terms *excluding* baseflow
+ aJac(ixDiag,watState) = (dt/mLayerDepth(jLayer))*(-dq_dHydStateBelow(iLayer-1) + dq_dHydStateAbove(iLayer)) + dMat(watState)
+
+ ! - compute the lower-diagonal elements
+ if(iLayer>1)then
+ if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(iLayer-1),watState),watState) = (dt/mLayerDepth(jLayer-1))*( dq_dHydStateBelow(iLayer-1))
+ endif
+
+ ! - compute the upper-diagonal elements
+ if(iLayer<nSoil)then
+ if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(iLayer+1),watState),watState) = (dt/mLayerDepth(jLayer+1))*(-dq_dHydStateAbove(iLayer))
+ endif
+
+ ! - include terms for baseflow WHY NOT USE
+ !if(computeBaseflow .and. nSoilOnlyHyd==nSoil)then
+ ! do pLayer=1,nSoil
+ ! qState = ixSoilOnlyHyd(pLayer) ! hydrology state index within the state subset
+ ! aJac(ixOffDiag(watState,qState),qState) = (dt/mLayerDepth(jLayer))*dBaseflow_dMatric(iLayer,pLayer) + aJac(ixOffDiag(watState,qState),qState)
+ ! end do
+ !endif
+
+ ! - include terms for surface infiltration below surface
+ if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(1),watState),watState) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(iLayer) + aJac(ixOffDiag(ixSoilOnlyHyd(1),watState),watState)
+
+ end do ! (looping through hydrology states in the soil domain)
+
+ ! - include terms for surface infiltration above surface
+ if(nSnowOnlyHyd>0 .and. ixSnowOnlyHyd(nSnow)/=integerMissing)then
+ if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(1),ixSnowOnlyHyd(nSnow)),ixSnowOnlyHyd(nSnow)) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(0)
+ elseif(computeVegFlux .and. ixVegHyd/=integerMissing)then !ixTopHyd = ixSoilOnlyHyd(1)
+ if(ixTopHyd/=integerMissing) aJac(ixOffDiag(ixTopHyd,ixVegHyd),ixVegHyd) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(0) + aJac(ixOffDiag(ixTopHyd,ixVegHyd),ixVegHyd)
+ endif
+
+ endif ! (if the subset includes hydrology state variables in the soil domain)
+
+ ! -----
+ ! * liquid water fluxes for the aquifer...
+ ! ----------------------------------------
+ if(ixAqWat/=integerMissing) then
+ aJac(ixDiag,ixAqWat) = -dBaseflow_dAquifer*dt + dMat(ixAqWat)
+ aJac(ixOffDiag(ixAqWat,ixSoilOnlyNrg(nSoil)),ixSoilOnlyNrg(nSoil)) = -dq_dNrgStateAbove(nSoil)*dt ! dAquiferRecharge_dTk = d_iLayerLiqFluxSoil(nSoil)_dTk
+ aJac(ixOffDiag(ixAqWat,ixSoilOnlyHyd(nSoil)),ixSoilOnlyHyd(nSoil)) = -dq_dHydStateAbove(nSoil)*dt ! dAquiferRecharge_dWat = d_iLayerLiqFluxSoil(nSoil)_dWat
+ ! - include derivatives of energy and water w.r.t soil transpiration (dependent on canopy transpiration)
+ if(computeVegFlux)then
+ aJac(ixOffDiag(ixAqWat,ixCasNrg),ixCasNrg) = -dAquiferTrans_dTCanair*dt ! dVol/dT (K-1)
+ aJac(ixOffDiag(ixAqWat,ixVegNrg),ixVegNrg) = -dAquiferTrans_dTCanopy*dt ! dVol/dT (K-1)
+ aJac(ixOffDiag(ixAqWat,ixVegHyd),ixVegHyd) = -dAquiferTrans_dCanWat*dt ! dVol/dLiq (kg m-2)-1
+ aJac(ixOffDiag(ixAqWat,ixTopNrg),ixTopNrg) = -dAquiferTrans_dTGround*dt ! dVol/dT (K-1)
+ endif
+ endif
+
+ ! -----
+ ! * cross derivatives in the soil domain...
+ ! -----------------------------------------------------------------------------
+ if(nSoilOnlyHyd>0 .and. nSoilOnlyNrg>0)then
+ do iLayer=1,nSoilOnlyNrg
+
+ ! - check that the soil layer is desired
+ if(ixSoilOnlyNrg(iLayer)==integerMissing) cycle
+
+ ! - define indices of the soil layers
+ jLayer = iLayer+nSnow ! index of layer in the snow+soil vector
+
+ ! - define the energy state variable
+ nrgState = ixNrgLayer(jLayer) ! index within the full state vector
+
+ ! - define index of hydrology state variable within the state subset
+ watState = ixSoilOnlyHyd(iLayer)
+
+ ! only compute derivatives if the water state for the current layer is within the state subset
+ if(watstate/=integerMissing)then
+
+ ! - include derivatives in liquid water fluxes w.r.t. temperature for current layer
+ aJac(ixOffDiag(watState,nrgState),nrgState) = (dt/mLayerDepth(jLayer))*(-dq_dNrgStateBelow(iLayer-1) + dq_dNrgStateAbove(iLayer)) ! dVol/dT (K-1) -- flux depends on ice impedance
+
+ ! - compute lower diagonal elements
+ if(iLayer>1)then
+ if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(iLayer-1),nrgState),nrgState) = (dt/mLayerDepth(jLayer-1))*( dq_dNrgStateBelow(iLayer-1)) ! K-1
+ endif
+
+ ! compute upper-diagonal elements
+ if(iLayer<nSoil)then
+ if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(iLayer+1),nrgState),nrgState) = (dt/mLayerDepth(jLayer+1))*(-dq_dNrgStateAbove(iLayer)) ! K-1
+ endif
+
+ ! - include derivatives w.r.t. ground evaporation
+ if(nSnow==0 .and. iLayer==1)then ! upper-most soil layer
+ if(computeVegFlux)then
+ aJac(ixOffDiag(watState,ixCasNrg),ixCasNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanair/iden_water) ! dVol/dT (K-1)
+ aJac(ixOffDiag(watState,ixVegNrg),ixVegNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanopy/iden_water) ! dVol/dT (K-1)
+ aJac(ixOffDiag(watState,ixVegHyd),ixVegHyd) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dCanWat/iden_water) ! dVol/dLiq (kg m-2)-1
+ endif
+ aJac(ixOffDiag(watState,ixTopNrg),ixTopNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTGround/iden_water) + aJac(ixOffDiag(watState,ixTopNrg),ixTopNrg) ! dVol/dT (K-1)
+ endif
+
+ ! - include derivatives of energy and water w.r.t soil transpiration (dependent on canopy transpiration)
+ if(computeVegFlux)then
+ aJac(ixOffDiag(watState,ixCasNrg),ixCasNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTCanair(iLayer)) + aJac(ixOffDiag(watState,ixCasNrg),ixCasNrg) ! dVol/dT (K-1)
+ aJac(ixOffDiag(watState,ixVegNrg),ixVegNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTCanopy(iLayer)) + aJac(ixOffDiag(watState,ixVegNrg),ixVegNrg) ! dVol/dT (K-1)
+ aJac(ixOffDiag(watState,ixVegHyd),ixVegHyd) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dCanWat(iLayer)) + aJac(ixOffDiag(watState,ixVegHyd),ixVegHyd) ! dVol/dLiq (kg m-2)-1
+ aJac(ixOffDiag(watState,ixTopNrg),ixTopNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTGround(iLayer)) + aJac(ixOffDiag(watState,ixTopNrg),ixTopNrg) ! dVol/dT (K-1)
+ endif
+
+ ! - include derivatives in energy fluxes w.r.t. with respect to water for current layer
+ aJac(ixOffDiag(nrgState,watState),watState) = dVolHtCapBulk_dPsi0(iLayer) * mLayerdTemp_dt(jLayer) &
+ + (dt/mLayerDepth(jLayer))*(-dNrgFlux_dWatBelow(jLayer-1) + dNrgFlux_dWatAbove(jLayer))
+ if(mLayerdTheta_dTk(jLayer) > verySmall)then ! ice is present
+ aJac(ixOffDiag(nrgState,watState),watState) = -LH_fus*iden_water * dVolTot_dPsi0(iLayer) + aJac(ixOffDiag(nrgState,watState),watState) ! dNrg/dMat (J m-3 m-1) -- dMat changes volumetric water, and hence ice content
+ endif
+
+ ! - include derivatives of heat capacity w.r.t water fluxes for surrounding layers starting with layer above
+ if(iLayer>1)then
+ if(ixSoilOnlyNrg(iLayer-1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyNrg(iLayer-1),watState),watState) = (dt/mLayerDepth(jLayer-1))*( dNrgFlux_dWatBelow(jLayer-1) )
+ elseif(iLayer==1 .and. nSnowOnlyNrg>0)then !top soil layer and there is snow above
+ if(ixSnowOnlyNrg(nSnow)/=integerMissing) aJac(ixOffDiag(ixSnowOnlyNrg(nSnow),watState),watState) = (dt/mLayerDepth(nSnow))*( dNrgFlux_dWatBelow(nSnow) )
+ endif
+
+ ! (cross-derivative terms for the layer below)
+ if(iLayer<nSoil)then
+ if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyNrg(iLayer+1),watState),watState) = (dt/mLayerDepth(jLayer+1))*(-dNrgFlux_dWatAbove(jLayer ) )
+ endif
+
+ endif ! (if the water state for the current layer is within the state subset)
+
+ ! - include terms for surface infiltration below surface
+ if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(1),nrgState),nrgState) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(iLayer) + aJac(ixOffDiag(ixSoilOnlyHyd(1),nrgState),nrgState)
+
+ end do ! (looping through soil layers)
+
+ ! - include terms for surface infiltration above surface
+ if(nSnowOnlyHyd>0 .and. ixSnowOnlyNrg(nSnow)/=integerMissing)then
+ if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixOffDiag(ixSoilOnlyHyd(1),ixSnowOnlyNrg(nSnow)),ixSnowOnlyNrg(nSnow)) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(0)
+ elseif(computeVegFlux .and. ixVegNrg/=integerMissing) then !ixTopHyd = ixSoilOnlyHyd(1)
+ if(ixTopHyd/=integerMissing) aJac(ixOffDiag(ixTopHyd,ixVegNrg),ixVegNrg) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(0) + aJac(ixOffDiag(ixTopHyd,ixVegNrg),ixVegNrg)
+ endif
+
+ endif ! (if there are state variables for both water and energy in the soil domain)
+
+ if(globalPrintFlag)then
+ print*, '** banded analytical Jacobian:'
+ write(*,'(a4,1x,100(i17,1x))') 'xCol', (iLayer, iLayer=min(iJac1,nState),min(iJac2,nState))
+ do iLayer=kl+1,nBands
+ write(*,'(i4,1x,100(e17.10,1x))') iLayer, (aJac(iLayer,jLayer),jLayer=min(iJac1,nState),min(iJac2,nState))
+ end do
+ endif
+ !print*, 'PAUSE: banded analytical Jacobian'; read(*,*)
+
+ ! *********************************************************************************************************************************************************
+ ! *********************************************************************************************************************************************************
+ ! * PART 2: FULL MATRIX
+ ! *********************************************************************************************************************************************************
+ ! *********************************************************************************************************************************************************
+ case(ixFullMatrix)
+
+ ! check
+ if(size(aJac,1)/=size(dMat) .or. size(aJac,2)/=size(dMat))then
+ message=trim(message)//'unexpected shape of the Jacobian matrix: expect aJac(nState,nState)'
+ err=20; return
+ endif
+
+ ! -----
+ ! * energy and liquid fluxes over vegetation...
+ ! ---------------------------------------------
+ if(computeVegFlux)then ! (derivatives only defined when vegetation protrudes over the surface)
+
+ ! * energy fluxes with the canopy water
+ if(ixVegHyd/=integerMissing)then
+
+ ! cross-derivative terms w.r.t. system temperatures (kg m-2 K-1)
+ if(ixCasNrg/=integerMissing) aJac(ixVegHyd,ixCasNrg) = -dCanopyEvaporation_dTCanair*dt
+ ! dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy is the derivative in throughfall and canopy drainage with canopy temperature
+ if(ixVegNrg/=integerMissing) aJac(ixVegHyd,ixVegNrg) = -dCanopyEvaporation_dTCanopy*dt + dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy
+ ! * liquid water fluxes for vegetation canopy (-)
+ aJac(ixVegHyd,ixVegHyd) = -scalarFracLiqVeg*(dCanopyEvaporation_dCanWat - scalarCanopyLiqDeriv)*dt + 1._rkind
+ if(ixTopNrg/=integerMissing) aJac(ixVegHyd,ixTopNrg) = -dCanopyEvaporation_dTGround*dt
+
+ ! cross-derivative terms w.r.t. canopy water (kg-1 m2)
+ if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegHyd) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarFracLiqVeg*scalarCanopyLiqDeriv)/iden_water
+
+ ! cross-derivative terms w.r.t. canopy liquid water (J m-1 kg-1)
+ ! NOTE: dIce/dLiq = (1 - scalarFracLiqVeg); dIce*LH_fus/canopyDepth = J m-3; dLiq = kg m-2
+ if(ixVegNrg/=integerMissing) aJac(ixVegNrg,ixVegHyd) = (-1._rkind + scalarFracLiqVeg)*LH_fus/canopyDepth &
+ + dVolHtCapBulk_dCanWat * scalarCanopydTemp_dt - (dt/canopyDepth) * dCanopyNetFlux_dCanWat ! dF/dLiq
+ if(ixTopNrg/=integerMissing) aJac(ixTopNrg,ixVegHyd) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanWat)
+ endif
+
+ ! cross-derivative terms w.r.t. canopy temperature (K-1)
+ if(ixVegNrg/=integerMissing)then
+ if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegNrg) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarCanopyLiqDeriv*dCanLiq_dTcanopy)/iden_water
+ endif
+
+ ! energy fluxes with the canopy air space (J m-3 K-1)
+ if(ixCasNrg/=integerMissing)then
+ aJac(ixCasNrg,ixCasNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanairTemp) + dMat(ixCasNrg)
+ if(ixVegNrg/=integerMissing) aJac(ixCasNrg,ixVegNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanopyTemp)
+ if(ixTopNrg/=integerMissing) aJac(ixCasNrg,ixTopNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dGroundTemp)
+ endif
+
+ ! energy fluxes with the vegetation canopy (J m-3 K-1)
+ if(ixVegNrg/=integerMissing)then
+ if(ixCasNrg/=integerMissing) aJac(ixVegNrg,ixCasNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanairTemp)
+ aJac(ixVegNrg,ixVegNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanopyTemp) + dMat(ixVegNrg)
+ if(ixTopNrg/=integerMissing) aJac(ixVegNrg,ixTopNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dGroundTemp)
+ endif
+
+ ! energy fluxes with the surface (J m-3 K-1)
+ if(ixTopNrg/=integerMissing)then
+ if(ixCasNrg/=integerMissing) aJac(ixTopNrg,ixCasNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanairTemp)
+ if(ixVegNrg/=integerMissing) aJac(ixTopNrg,ixVegNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanopyTemp)
+ endif
+
+ endif ! if there is a need to compute energy fluxes within vegetation
+
+ ! -----
+ ! * energy fluxes for the snow+soil domain...
+ ! -------------------------------------------
+ if(nSnowSoilNrg>0)then
+ do iLayer=1,nLayers ! loop through all layers in the snow+soil domain
+
+ ! check if the state is in the subset
+ if(ixSnowSoilNrg(iLayer)==integerMissing) cycle
+
+ ! - define index within the state subset and the full state vector
+ jState = ixSnowSoilNrg(iLayer) ! index within the state subset
+
+ ! - diagonal elements
+ aJac(jState,jState) = (dt/mLayerDepth(iLayer))*(-dNrgFlux_dTempBelow(iLayer-1) + dNrgFlux_dTempAbove(iLayer)) + dMat(jState)
+
+ ! - lower-diagonal elements
+ if(iLayer>1)then
+ if(ixSnowSoilNrg(iLayer-1)/=integerMissing) aJac(ixSnowSoilNrg(iLayer-1),jState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dTempBelow(iLayer-1) )
+ endif
+
+ ! - upper diagonal elements
+ if(iLayer<nLayers)then
+ if(ixSnowSoilNrg(iLayer+1)/=integerMissing) aJac(ixSnowSoilNrg(iLayer+1),jState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dTempAbove(iLayer ) )
+ endif
+
+ end do ! (looping through energy states in the snow+soil domain)
+ endif ! (if the subset includes energy state variables in the snow+soil domain)
+
+ ! -----
+ ! * liquid water fluxes for the snow domain...
+ ! --------------------------------------------
+ if(nSnowOnlyHyd>0)then
+ do iLayer=1,nSnow ! loop through layers in the snow domain
+
+ ! - check that the snow layer is desired
+ if(ixSnowOnlyHyd(iLayer)==integerMissing) cycle
+
+ ! - define state indices for the current layer
+ watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset
+
+ ! compute factor to convert liquid water derivative to total water derivative
+ select case( ixHydType(iLayer) )
+ case(iname_watLayer); convLiq2tot = mLayerFracLiqSnow(iLayer)
+ case default; convLiq2tot = 1._rkind
+ end select
+
+ ! - diagonal elements
+ aJac(watState,watState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot + dMat(watState)
+
+ ! - lower-diagonal elements
+ if(iLayer>1)then
+ if(ixSnowOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSnowOnlyHyd(iLayer-1),watState) = 0._rkind ! sub-diagonal: no dependence on other layers
+ endif
+
+ ! - upper diagonal elements
+ if(iLayer<nSnow)then
+ if(ixSnowOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSnowOnlyHyd(iLayer+1),watState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot ! dVol(below)/dLiq(above) -- (-)
+ endif
+
+ end do ! (looping through liquid water states in the snow domain)
+ endif ! (if the subset includes hydrology state variables in the snow domain)
+
+ ! -----
+ ! * cross derivatives in the snow domain...
+ ! ----------------------------------------
+ if(nSnowOnlyHyd>0 .and. nSnowOnlyNrg>0)then
+ do iLayer=1,nSnow ! loop through layers in the snow domain
+
+ ! - check that the snow layer is desired
+ if(ixSnowOnlyNrg(iLayer)==integerMissing) cycle
+
+ ! (define the energy state)
+ nrgState = ixSnowOnlyNrg(iLayer) ! index within the full state vector
+
+ ! - define state indices for the current layer
+ watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset
+
+ if(watstate/=integerMissing)then ! (energy state for the current layer is within the state subset)
+
+ ! - include derivatives of energy fluxes w.r.t water fluxes for current layer
+ aJac(nrgState,watState) = (-1._rkind + mLayerFracLiqSnow(iLayer))*LH_fus*iden_water &
+ + dVolHtCapBulk_dTheta(iLayer) * mLayerdTemp_dt(iLayer) &
+ + (dt/mLayerDepth(iLayer))*(-dNrgFlux_dWatBelow(iLayer-1) + dNrgFlux_dWatAbove(iLayer)) ! (dF/dLiq)
+
+ ! - include derivatives of water fluxes w.r.t energy fluxes for current layer
+ aJac(watState,nrgState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! (dVol/dT)
+
+ ! (cross-derivative terms for the layer below)
+ if(iLayer < nSnow)then
+ aJac(ixSnowOnlyHyd(iLayer+1),nrgState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! dVol(below)/dT(above) -- K-1
+ endif ! (if there is a water state in the layer below the current layer in the given state subset)
+
+ ! - include derivatives of heat capacity w.r.t water fluxes for surrounding layers starting with layer above
+ if(iLayer>1)then
+ if(ixSnowOnlyNrg(iLayer-1)/=integerMissing) aJac(ixSnowOnlyNrg(iLayer-1),watState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dWatBelow(iLayer-1) )
+ endif
+
+ ! (cross-derivative terms for the layer below)
+ if(iLayer<nSnow)then
+ if(ixSnowOnlyNrg(iLayer+1)/=integerMissing) aJac(ixSnowOnlyNrg(iLayer+1),watState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dWatAbove(iLayer ) )
+ elseif(iLayer==nSnow .and. nSoilOnlyNrg>0)then !bottom snow layer and there is soil below
+ if(ixSoilOnlyNrg(1)/=integerMissing) aJac(ixSoilOnlyNrg(1),watState) = (dt/mLayerDepth(nSnow+1))*(-dNrgFlux_dWatAbove(nSnow) )
+ endif
+
+ endif ! (if the energy state for the current layer is within the state subset)
+
+ end do ! (looping through snow layers)
+ endif ! (if there are state variables for both water and energy in the snow domain)
+
+ ! -----
+ ! * liquid water fluxes for the soil domain...
+ ! --------------------------------------------
+ if(nSoilOnlyHyd>0)then
+ do iLayer=1,nSoil
+
+ ! - check that the soil layer is desired
+ if(ixSoilOnlyHyd(iLayer)==integerMissing) cycle
+
+ ! - define state indices
+ watState = ixSoilOnlyHyd(iLayer) ! hydrology state index within the state subset
+
+ ! - define indices of the soil layers
+ jLayer = iLayer+nSnow ! index of layer in the snow+soil vector
+
+ ! - compute the diagonal elements
+ ! all terms *excluding* baseflow
+ aJac(watState,watState) = (dt/mLayerDepth(jLayer))*(-dq_dHydStateBelow(iLayer-1) + dq_dHydStateAbove(iLayer)) + dMat(watState)
+
+ ! - compute the lower-diagonal elements
+ if(iLayer>1)then
+ if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer-1),watState) = (dt/mLayerDepth(jLayer-1))*( dq_dHydStateBelow(iLayer-1))
+ endif
+
+ ! - compute the upper-diagonal elements
+ if(iLayer<nSoil)then
+ if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer+1),watState) = (dt/mLayerDepth(jLayer+1))*(-dq_dHydStateAbove(iLayer))
+ endif
+
+ ! - include terms for baseflow
+ if(computeBaseflow .and. nSoilOnlyHyd==nSoil)then
+ do pLayer=1,nSoil
+ qState = ixSoilOnlyHyd(pLayer) ! hydrology state index within the state subset
+ aJac(watState,qState) = (dt/mLayerDepth(jLayer))*dBaseflow_dMatric(iLayer,pLayer) + aJac(watState,qState)
+ end do
+ endif
+
+ ! - include terms for surface infiltration below surface
+ if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),watState) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(iLayer) + aJac(ixSoilOnlyHyd(1),watState)
+
+ end do ! (looping through hydrology states in the soil domain)
+
+ ! - include terms for surface infiltration above surface
+ if(nSnowOnlyHyd>0 .and. ixSnowOnlyHyd(nSnow)/=integerMissing)then
+ if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),ixSnowOnlyHyd(nSnow)) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(0)
+ elseif(computeVegFlux .and. ixVegHyd/=integerMissing)then !ixTopHyd = ixSoilOnlyHyd(1)
+ if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegHyd) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(0) + aJac(ixTopHyd,ixVegHyd)
+ endif
+
+ endif ! (if the subset includes hydrology state variables in the soil domain)
+
+ ! -----
+ ! * liquid water fluxes for the aquifer...
+ ! ----------------------------------------
+ if(ixAqWat/=integerMissing) then
+ aJac(ixAqWat,ixAqWat) = -dBaseflow_dAquifer*dt + dMat(ixAqWat)
+ aJac(ixAqWat,ixSoilOnlyNrg(nSoil)) = -dq_dNrgStateAbove(nSoil)*dt ! dAquiferRecharge_dTk = d_iLayerLiqFluxSoil(nSoil)_dTk
+ aJac(ixAqWat,ixSoilOnlyHyd(nSoil)) = -dq_dHydStateAbove(nSoil)*dt ! dAquiferRecharge_dWat = d_iLayerLiqFluxSoil(nSoil)_dWat
+ ! - include derivatives of energy and water w.r.t soil transpiration (dependent on canopy transpiration)
+ if(computeVegFlux)then
+ aJac(ixAqWat,ixCasNrg) = -dAquiferTrans_dTCanair*dt ! dVol/dT (K-1)
+ aJac(ixAqWat,ixVegNrg) = -dAquiferTrans_dTCanopy*dt ! dVol/dT (K-1)
+ aJac(ixAqWat,ixVegHyd) = -dAquiferTrans_dCanWat*dt ! dVol/dLiq (kg m-2)-1
+ aJac(ixAqWat,ixTopNrg) = -dAquiferTrans_dTGround*dt ! dVol/dT (K-1)
+ endif
+ endif
+
+ ! -----
+ ! * cross derivatives in the soil domain...
+ ! -----------------------------------------------------------------------------
+ if(nSoilOnlyHyd>0 .and. nSoilOnlyNrg>0)then
+ do iLayer=1,nSoilOnlyNrg
+
+ ! - check that the soil layer is desired
+ if(ixSoilOnlyNrg(iLayer)==integerMissing) cycle
+
+ ! - define indices of the soil layers
+ jLayer = iLayer+nSnow ! index of layer in the snow+soil vector
+
+ ! - define the energy state variable
+ nrgState = ixNrgLayer(jLayer) ! index within the full state vector
+
+ ! - define index of hydrology state variable within the state subset
+ watState = ixSoilOnlyHyd(iLayer)
+
+ ! only compute derivatives if the water state for the current layer is within the state subset
+ if(watstate/=integerMissing)then
+
+ ! - include derivatives in liquid water fluxes w.r.t. temperature for current layer
+ aJac(watState,nrgState) = (dt/mLayerDepth(jLayer))*(-dq_dNrgStateBelow(iLayer-1) + dq_dNrgStateAbove(iLayer)) ! dVol/dT (K-1) -- flux depends on ice impedance
+
+ ! - compute lower diagonal elements
+ if(iLayer>1)then
+ if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer-1),nrgState) = (dt/mLayerDepth(jLayer-1))*( dq_dNrgStateBelow(iLayer-1)) ! K-1
+ endif
+
+ ! compute upper-diagonal elements
+ if(iLayer<nSoil)then
+ if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer+1),nrgState) = (dt/mLayerDepth(jLayer+1))*(-dq_dNrgStateAbove(iLayer)) ! K-1
+ endif
+
+ ! - include derivatives w.r.t. ground evaporation
+ if(nSnow==0 .and. iLayer==1)then ! upper-most soil layer
+ if(computeVegFlux)then
+ aJac(watState,ixCasNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanair/iden_water) ! dVol/dT (K-1)
+ aJac(watState,ixVegNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanopy/iden_water) ! dVol/dT (K-1)
+ aJac(watState,ixVegHyd) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dCanWat/iden_water) ! dVol/dLiq (kg m-2)-1
+ endif
+ aJac(watState,ixTopNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTGround/iden_water) + aJac(watState,ixTopNrg) ! dVol/dT (K-1)
+ endif
+
+ ! - include derivatives of energy and water w.r.t soil transpiration (dependent on canopy transpiration)
+ if(computeVegFlux)then
+ aJac(watState,ixCasNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTCanair(iLayer)) + aJac(watState,ixCasNrg) ! dVol/dT (K-1)
+ aJac(watState,ixVegNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTCanopy(iLayer)) + aJac(watState,ixVegNrg) ! dVol/dT (K-1)
+ aJac(watState,ixVegHyd) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dCanWat(iLayer)) + aJac(watState,ixVegHyd) ! dVol/dLiq (kg m-2)-1
+ aJac(watState,ixTopNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTGround(iLayer)) + aJac(watState,ixTopNrg) ! dVol/dT (K-1)
+ endif
+
+ ! - include derivatives in energy fluxes w.r.t. with respect to water for current layer
+ aJac(nrgState,watState) = dVolHtCapBulk_dPsi0(iLayer) * mLayerdTemp_dt(jLayer) &
+ + (dt/mLayerDepth(jLayer))*(-dNrgFlux_dWatBelow(jLayer-1) + dNrgFlux_dWatAbove(jLayer))
+ if(mLayerdTheta_dTk(jLayer) > verySmall)then ! ice is present
+ aJac(nrgState,watState) = -LH_fus*iden_water * dVolTot_dPsi0(iLayer) + aJac(nrgState,watState) ! dNrg/dMat (J m-3 m-1) -- dMat changes volumetric water, and hence ice content
+ endif
+
+ ! - include derivatives of heat capacity w.r.t water fluxes for surrounding layers starting with layer above
+ if(iLayer>1)then
+ if(ixSoilOnlyNrg(iLayer-1)/=integerMissing) aJac(ixSoilOnlyNrg(iLayer-1),watState) = (dt/mLayerDepth(jLayer-1))*( dNrgFlux_dWatBelow(jLayer-1) )
+ elseif(iLayer==1 .and. nSnowOnlyNrg>0)then !top soil layer and there is snow above
+ if(ixSnowOnlyNrg(nSnow)/=integerMissing) aJac(ixSnowOnlyNrg(nSnow),watState) = (dt/mLayerDepth(nSnow))*( dNrgFlux_dWatBelow(nSnow) )
+ endif
+
+ ! (cross-derivative terms for the layer below)
+ if(iLayer<nSoil)then
+ if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyNrg(iLayer+1),watState) = (dt/mLayerDepth(jLayer+1))*(-dNrgFlux_dWatAbove(jLayer ) )
+ endif
+
+ endif ! (if the water state for the current layer is within the state subset)
+
+ ! - include terms for surface infiltration below surface
+ if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),nrgState) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(iLayer) + aJac(ixSoilOnlyHyd(1),nrgState)
+
+ end do ! (looping through soil layers)
+
+ ! - include terms for surface infiltration above surface
+ if(nSnowOnlyHyd>0 .and. ixSnowOnlyNrg(nSnow)/=integerMissing)then
+ if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),ixSnowOnlyNrg(nSnow)) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(0)
+ elseif(computeVegFlux .and. ixVegNrg/=integerMissing) then !ixTopHyd = ixSoilOnlyHyd(1)
+ if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegNrg) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(0) + aJac(ixTopHyd,ixVegNrg)
+ endif
+
+ endif ! (if there are state variables for both water and energy in the soil domain)
+
+ ! print the Jacobian
+ if(globalPrintFlag)then
+ print*, '** analytical Jacobian (full):'
+ write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=min(iJac1,nState),min(iJac2,nState))
+ do iLayer=min(iJac1,nState),min(iJac2,nState)
+ write(*,'(i4,1x,100(e12.5,1x))') iLayer, aJac(min(iJac1,nState):min(iJac2,nState),iLayer)
+ end do
+ endif
+
+ !print*, '** analytical Jacobian (full):'
+ !write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=1,size(aJac,2))
+ !do iLayer=1,size(aJac,1)
+ ! write(*,'(i4,1x,100(e12.5,1x))') iLayer, aJac(iLayer,1:size(aJac,2))
+ !end do
+ !print *, '--------------------------------------------------------------'
+
+ ! ***
+ ! check
+ case default; err=20; message=trim(message)//'unable to identify option for the type of matrix'; return
+
+ end select ! type of matrix
+
+ if(any(isNan(aJac)))then
+ print *, '******************************* WE FOUND NAN IN JACOBIAN ************************************'
+ stop 1
+ message=trim(message)//'we found NaN'
+ err=20; return
+ endif
+
+ ! end association to variables in the data structures
+ end associate
+
+ end subroutine computJacob
+
+
+ ! **********************************************************************************************************
+ ! private function: get the off-diagonal index in the band-diagonal matrix
+ ! **********************************************************************************************************
+ function ixOffDiag(jState,iState)
+ implicit none
+ integer(i4b),intent(in) :: jState ! off-diagonal state
+ integer(i4b),intent(in) :: iState ! diagonal state
+ integer(i4b) :: ixOffDiag ! off-diagonal index in gthe band-diagonal matrix
+ ixOffDiag = ixBandOffset + jState - iState
+ end function ixOffDiag
+
+ end module computJacob_module
+
\ No newline at end of file
diff --git a/build/source/engine/conv_funcs.f90 b/build/source/engine/conv_funcs.f90
index b4bc8f27042cb414bb75555ed7683f2c5dcd690a..b0481d58c21aaf092e20fdec6c678b8ae0b724e8 100755
--- a/build/source/engine/conv_funcs.f90
+++ b/build/source/engine/conv_funcs.f90
@@ -70,26 +70,33 @@ end function vapPress
! NOTE: temperature units are degC !!!!
! ***************************************************************************************************************
subroutine satVapPress(TC, SVP, dSVP_dT)
-IMPLICIT NONE
-! input
-real(dp), intent(in) :: TC ! temperature (C)
-! output
-real(dp), intent(out) :: SVP ! saturation vapor pressure (Pa)
-real(dp), intent(out) :: dSVP_dT ! d(SVP)/dT
-! local
-real(dp), parameter :: X1 = 17.27_dp
-real(dp), parameter :: X2 = 237.30_dp
-! local (use to test derivative calculations)
-real(dp),parameter :: dx = 1.e-8_dp ! finite difference increment
-logical(lgt),parameter :: testDeriv=.false. ! flag to test the derivative
-!---------------------------------------------------------------------------------------------------
-! Units note : Pa = N m-2 = kg m-1 s-2
-! SATVPFRZ= 610.8 ! Saturation water vapour pressure at 273.16K (Pa)
-
-SVP = SATVPFRZ * EXP( (X1*TC)/(X2 + TC) ) ! Saturated Vapour Press (Pa)
-dSVP_dT = SVP * (X1/(X2 + TC) - X1*TC/(X2 + TC)**2._dp)
-if(testDeriv) print*, 'dSVP_dT check... ', SVP, dSVP_dT, (SATVPRESS(TC+dx) - SVP)/dx
+ IMPLICIT NONE
+ ! input
+ real(rkind), intent(in) :: TC ! temperature (C)
+ ! output
+ real(rkind), intent(out) :: SVP ! saturation vapor pressure (Pa)
+ real(rkind), intent(out) :: dSVP_dT ! d(SVP)/dT
+ ! local
+ real(rkind), parameter :: X1 = 17.27_rkind
+ real(rkind), parameter :: X2 = 237.30_rkind
+ ! local (use to test derivative calculations)
+ real(rkind),parameter :: dx = 1.e-8_rkind ! finite difference increment
+ logical(lgt),parameter :: testDeriv=.false. ! flag to test the derivative
+ !---------------------------------------------------------------------------------------------------
+ ! Units note : Pa = N m-2 = kg m-1 s-2
+ ! SATVPFRZ= 610.8 ! Saturation water vapour pressure at 273.16K (Pa)
+
+ if(X2 + TC < 0)then
+ !print*, "error, canopy temperature is very low, satVapPress=Inf" !will fail as SVP=inf
+ SVP = 0._rkind
+ dSVP_dT = 0._rkind
+ else
+ SVP = SATVPFRZ * EXP( (X1*TC)/(X2 + TC) ) ! Saturated Vapour Press (Pa)
+ dSVP_dT = SVP * (X1/(X2 + TC) - X1*TC/(X2 + TC)**2._rkind)
+ end if
+ if(testDeriv) print*, 'dSVP_dT check... ', SVP, dSVP_dT, (SATVPRESS(TC+dx) - SVP)/dx
END SUBROUTINE satVapPress
+
! ***************************************************************************************************************
diff --git a/build/source/engine/eval8summa.f90 b/build/source/engine/eval8summa.f90
index 23816632f6913c5e5649e39ee214f0aa69c13262..2c198cfe959f6d4c6b687e75d8f797d0959c976c 100755
--- a/build/source/engine/eval8summa.f90
+++ b/build/source/engine/eval8summa.f90
@@ -1,5 +1,5 @@
! SUMMA - Structure for Unifying Multiple Modeling Alternatives
-! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
+! Copyright (C) 2014-2015 NCAR/RAL
!
! This file is part of SUMMA
!
@@ -20,390 +20,395 @@
module eval8summa_module
-! data types
-USE nrtype
-
-! access missing values
-USE globalData,only:integerMissing ! missing integer
-USE globalData,only:realMissing ! missing double precision number
-USE globalData,only:quadMissing ! missing quadruple precision number
-
-! access the global print flag
-USE globalData,only:globalPrintFlag
-
-! define access to state variables to print
-USE globalData,only: iJac1 ! first layer of the Jacobian to print
-USE globalData,only: iJac2 ! last layer of the Jacobian to print
-
-! domain types
-USE globalData,only:iname_veg ! named variables for vegetation
-USE globalData,only:iname_snow ! named variables for snow
-USE globalData,only:iname_soil ! named variables for soil
-
-! 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 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
-
-! constants
-USE multiconst,only:&
- Tfreeze, & ! temperature at freezing (K)
- LH_fus, & ! latent heat of fusion (J kg-1)
- LH_vap, & ! latent heat of vaporization (J kg-1)
- LH_sub, & ! latent heat of sublimation (J kg-1)
- Cp_air, & ! specific heat of air (J kg-1 K-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 (dp)
- var_ilength, & ! data vector with variable length dimension (i4b)
- var_dlength, & ! data vector with variable length dimension (dp)
- zLookup, &
- model_options ! defines the model decisions
-
-! indices that define elements of the data structures
-USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure
-USE var_lookup,only:iLookPARAM ! named variables for structure elements
-USE var_lookup,only:iLookPROG ! named variables for structure elements
-USE var_lookup,only:iLookINDEX ! named variables for structure elements
-USE var_lookup,only:iLookDIAG ! named variables for structure elements
-USE var_lookup,only:iLookFLUX ! named variables for structure elements
-USE var_lookup,only:iLookDERIV ! named variables for structure elements
-
-! look-up values for the choice of groundwater representation (local-column, or single-basin)
-USE mDecisions_module,only: &
- localColumn, & ! separate groundwater representation in each local soil column
- singleBasin ! single groundwater store over the entire basin
-
-! look-up values for the choice of groundwater parameterization
-USE mDecisions_module,only: &
- qbaseTopmodel, & ! TOPMODEL-ish baseflow parameterization
- bigBucket, & ! a big bucket (lumped aquifer model)
- noExplicit ! no explicit groundwater parameterization
-
-! look-up values for the form of Richards' equation
-USE mDecisions_module,only: &
- moisture, & ! moisture-based form of Richards' equation
- mixdform ! mixed form of Richards' equation
-
-implicit none
-private
-public::eval8summa
-
-contains
-
-! **********************************************************************************************************
-! public subroutine eval8summa: compute the residual vector and the Jacobian matrix
-! **********************************************************************************************************
-subroutine eval8summa(&
- ! input: model control
- dt, & ! intent(in): length of the time step (seconds)
- nSnow, & ! intent(in): number of snow layers
- nSoil, & ! intent(in): number of soil layers
- nLayers, & ! intent(in): total number of layers
- nState, & ! intent(in): total number of state variables
- firstSubStep, & ! intent(in): flag to indicate if we are processing the first sub-step
- firstFluxCall, & ! intent(inout): flag to indicate if we are processing the first flux call
- firstSplitOper, & ! intent(in): flag to indicate if we are processing the first flux call in a splitting operation
- computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation
- scalarSolution, & ! intent(in): flag to indicate the scalar solution
- ! input: state vectors
- stateVecTrial, & ! intent(in): model state vector
- fScale, & ! intent(in): function scaling vector
- sMul, & ! intent(in): state vector multiplier (used in the residual calculations)
- ! input: data structures
- model_decisions, & ! intent(in): model decisions
- lookup_data, & ! intent(in): lookup tables
- type_data, & ! intent(in): type of vegetation and soil
- attr_data, & ! intent(in): spatial attributes
- mpar_data, & ! intent(in): model parameters
- forc_data, & ! intent(in): model forcing data
- bvar_data, & ! intent(in): average model variables for the entire basin
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- ! input-output: data structures
- indx_data, & ! intent(inout): index data
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- flux_data, & ! intent(inout): model fluxes for a local HRU
- deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables
- ! input-output: baseflow
- ixSaturation, & ! intent(inout): index of the lowest saturated layer (NOTE: only computed on the first iteration)
- dBaseflow_dMatric, & ! intent(out): derivative in baseflow w.r.t. matric head (s-1)
- ! output: flux and residual vectors
- feasible, & ! intent(out): flag to denote the feasibility of the solution
- fluxVec, & ! intent(out): flux vector
- resSink, & ! intent(out): additional (sink) terms on the RHS of the state equation
- resVec, & ! intent(out): residual vector
- fEval, & ! intent(out): function evaluation
- err,message) ! intent(out): error control
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! provide access to subroutines
- USE getVectorz_module, only:varExtract ! extract variables from the state vector
- USE updateVars_module, only:updateVars ! update prognostic variables
- USE t2enthalpy_module, only:t2enthalpy ! compute enthalpy
- USE computFlux_module, only:soilCmpres ! compute soil compression
- USE computFlux_module, only:computFlux ! compute fluxes given a state vector
- USE computResid_module,only:computResid ! compute residuals given a state vector
+ ! data types
+ USE nrtype
+
+ ! access missing values
+ USE globalData,only:integerMissing ! missing integer
+ USE globalData,only:realMissing ! missing double precision number
+ USE globalData,only:quadMissing ! missing quadruple precision number
+
+ ! access the global print flag
+ USE globalData,only:globalPrintFlag
+
+ ! define access to state variables to print
+ USE globalData,only: iJac1 ! first layer of the Jacobian to print
+ USE globalData,only: iJac2 ! last layer of the Jacobian to print
+
+ ! domain types
+ USE globalData,only:iname_veg ! named variables for vegetation
+ USE globalData,only:iname_snow ! named variables for snow
+ USE globalData,only:iname_soil ! named variables for soil
+
+ ! 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 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
+
+ ! constants
+ USE multiconst,only:&
+ Tfreeze, & ! temperature at freezing (K)
+ LH_fus, & ! latent heat of fusion (J kg-1)
+ LH_vap, & ! latent heat of vaporization (J kg-1)
+ LH_sub, & ! latent heat of sublimation (J kg-1)
+ Cp_air, & ! specific heat of air (J kg-1 K-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)
+ zLookup, & ! data vector with variable length dimension (rkind)
+ model_options ! defines the model decisions
+
+ ! indices that define elements of the data structures
+ USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure
+ USE var_lookup,only:iLookPARAM ! named variables for structure elements
+ USE var_lookup,only:iLookPROG ! named variables for structure elements
+ USE var_lookup,only:iLookINDEX ! named variables for structure elements
+ USE var_lookup,only:iLookDIAG ! named variables for structure elements
+ USE var_lookup,only:iLookFLUX ! named variables for structure elements
+ USE var_lookup,only:iLookDERIV ! named variables for structure elements
+
+ ! look-up values for the choice of groundwater representation (local-column, or single-basin)
+ USE mDecisions_module,only: &
+ localColumn, & ! separate groundwater representation in each local soil column
+ singleBasin ! single groundwater store over the entire basin
+
+ ! look-up values for the choice of groundwater parameterization
+ USE mDecisions_module,only: &
+ qbaseTopmodel, & ! TOPMODEL-ish baseflow parameterization
+ bigBucket, & ! a big bucket (lumped aquifer model)
+ noExplicit ! no explicit groundwater parameterization
+
+ ! look-up values for the form of Richards' equation
+ USE mDecisions_module,only: &
+ moisture, & ! moisture-based form of Richards' equation
+ mixdform ! mixed form of Richards' equation
+
implicit none
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! input: model control
- real(dp),intent(in) :: dt ! length of the time step (seconds)
- integer(i4b),intent(in) :: nSnow ! number of snow layers
- integer(i4b),intent(in) :: nSoil ! number of soil layers
- integer(i4b),intent(in) :: nLayers ! total number of layers
- integer(i4b),intent(in) :: nState ! total number of state variables
- logical(lgt),intent(in) :: firstSubStep ! flag to indicate if we are processing the first sub-step
- logical(lgt),intent(inout) :: firstFluxCall ! flag to indicate if we are processing the first flux call
- logical(lgt),intent(in) :: firstSplitOper ! flag to indicate if we are processing the first flux call in a splitting operation
- logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if computing fluxes over vegetation
- logical(lgt),intent(in) :: scalarSolution ! flag to denote if implementing the scalar solution
- ! input: state vectors
- real(dp),intent(in) :: stateVecTrial(:) ! model state vector
- real(dp),intent(in) :: fScale(:) ! function scaling vector
- real(qp),intent(in) :: sMul(:) ! NOTE: qp ! state vector multiplier (used in the residual calculations)
- ! input: data structures
- type(model_options),intent(in) :: model_decisions(:) ! model decisions
- type(zLookup), intent(in) :: lookup_data ! lookup tables
- type(var_i), intent(in) :: type_data ! type of vegetation and soil
- type(var_d), intent(in) :: attr_data ! spatial attributes
- type(var_dlength), intent(in) :: mpar_data ! model parameters
- type(var_d), intent(in) :: forc_data ! model forcing data
- type(var_dlength), intent(in) :: bvar_data ! model variables for the local basin
- type(var_dlength), intent(in) :: prog_data ! prognostic variables for a local HRU
- ! output: data structures
- type(var_ilength),intent(inout) :: indx_data ! indices defining model states and layers
- type(var_dlength),intent(inout) :: diag_data ! diagnostic variables for a local HRU
- type(var_dlength),intent(inout) :: flux_data ! model fluxes for a local HRU
- type(var_dlength),intent(inout) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
- ! input-output: baseflow
- integer(i4b),intent(inout) :: ixSaturation ! index of the lowest saturated layer (NOTE: only computed on the first iteration)
- real(dp),intent(out) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1)
- ! output: flux and residual vectors
- logical(lgt),intent(out) :: feasible ! flag to denote the feasibility of the solution
- real(dp),intent(out) :: fluxVec(:) ! flux vector
- real(dp),intent(out) :: resSink(:) ! sink terms on the RHS of the flux equation
- real(qp),intent(out) :: resVec(:) ! NOTE: qp ! residual vector
- real(dp),intent(out) :: fEval ! function evaluation
- ! output: error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! local variables
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! state variables
- real(dp) :: scalarCanairTempTrial ! trial value for temperature of the canopy air space (K)
- real(dp) :: scalarCanopyTempTrial ! trial value for temperature of the vegetation canopy (K)
- real(dp) :: scalarCanopyWatTrial ! trial value for liquid water storage in the canopy (kg m-2)
- real(dp),dimension(nLayers) :: mLayerTempTrial ! trial value for temperature of layers in the snow and soil domains (K)
- real(dp),dimension(nLayers) :: mLayerVolFracWatTrial ! trial value for volumetric fraction of total water (-)
- real(dp),dimension(nSoil) :: mLayerMatricHeadTrial ! trial value for total water matric potential (m)
- real(dp),dimension(nSoil) :: mLayerMatricHeadLiqTrial ! trial value for liquid water matric potential (m)
- real(dp) :: scalarAquiferStorageTrial ! trial value of storage of water in the aquifer (m)
- ! diagnostic variables
- logical(lgt),parameter :: needEnthalpy=.true. ! flag to compute enthalpy
- real(dp) :: scalarCanopyLiqTrial ! trial value for mass of liquid water on the vegetation canopy (kg m-2)
- real(dp) :: scalarCanopyIceTrial ! trial value for mass of ice on the vegetation canopy (kg m-2)
- real(dp),dimension(nLayers) :: mLayerVolFracLiqTrial ! trial value for volumetric fraction of liquid water (-)
- real(dp),dimension(nLayers) :: mLayerVolFracIceTrial ! trial value for volumetric fraction of ice (-)
- ! other local variables
- integer(i4b) :: iLayer ! index of model layer in the snow+soil domain
- integer(i4b) :: jState(1) ! index of model state for the scalar solution within the soil domain
- integer(i4b) :: ixBeg,ixEnd ! index of indices for the soil compression routine
- integer(i4b),parameter :: ixVegVolume=1 ! index of the desired vegetation control volumne (currently only one veg layer)
- real(dp) :: xMin,xMax ! minimum and maximum values for water content
- real(dp) :: scalarCanopyHydTrial ! trial value for mass of water on the vegetation canopy (kg m-2)
- real(dp),parameter :: canopyTempMax=500._dp ! expected maximum value for the canopy temperature (K)
- real(dp),dimension(nLayers) :: mLayerVolFracHydTrial ! trial value for volumetric fraction of water (-), general vector merged from Wat and Liq
- real(dp),dimension(nState) :: rVecScaled ! scaled residual vector
- character(LEN=256) :: cmessage ! error message of downwind routine
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! association to variables in the data structures
- ! --------------------------------------------------------------------------------------------------------------------------------
- associate(&
- ! model decisions
- ixRichards => model_decisions(iLookDECISIONS%f_Richards)%iDecision ,& ! intent(in): [i4b] index of the form of Richards' equation
- ! snow parameters
- snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) ,& ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1)
- ! soil parameters
- theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat ,& ! intent(in): [dp(:)] soil porosity (-)
- specificStorage => mpar_data%var(iLookPARAM%specificStorage)%dat(1) ,& ! intent(in): [dp] specific storage coefficient (m-1)
- ! canopy and layer depth
- canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp ] canopy depth (m)
- mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
- ! model state variables
- scalarSfcMeltPond => prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1) ,& ! intent(in): [dp] ponded water caused by melt of the "snow without a layer" (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 (-)
- mLayerMatricHeadLiq => diag_data%var(iLookDIAG%mLayerMatricHeadLiq)%dat ,& ! intent(in): [dp(:)] liquid water matric potential (m)
- ! model diagnostic variables
- scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg)%dat(1) ,& ! intent(in): [dp] fraction of liquid water on vegetation (-)
- mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow)%dat ,& ! intent(in): [dp(:)] fraction of liquid water in each snow layer (-)
- ! enthalpy
- scalarCanairEnthalpy => diag_data%var(iLookDIAG%scalarCanairEnthalpy)%dat(1) ,& ! intent(out): [dp] enthalpy of the canopy air space (J m-3)
- scalarCanopyEnthalpy => diag_data%var(iLookDIAG%scalarCanopyEnthalpy)%dat(1) ,& ! intent(out): [dp] enthalpy of the vegetation canopy (J m-3)
- mLayerEnthalpy => diag_data%var(iLookDIAG%mLayerEnthalpy)%dat ,& ! intent(out): [dp(:)] enthalpy of the snow+soil layers (J m-3)
- ! soil compression
- scalarSoilCompress => diag_data%var(iLookDIAG%scalarSoilCompress)%dat(1) ,& ! intent(in): [dp] total change in storage associated with compression of the soil matrix (kg m-2)
- mLayerCompress => diag_data%var(iLookDIAG%mLayerCompress)%dat ,& ! intent(in): [dp(:)] change in storage associated with compression of the soil matrix (-)
- ! derivatives
- dVolTot_dPsi0 => deriv_data%var(iLookDERIV%dVolTot_dPsi0)%dat ,& ! intent(in): [dp(:)] derivative in total water content w.r.t. total water matric potential
- dCompress_dPsi => deriv_data%var(iLookDERIV%dCompress_dPsi)%dat ,& ! intent(in): [dp(:)] derivative in compressibility w.r.t. matric head (m-1)
- ! mapping
- ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] mapping of full state vector to the state subset
- ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of control volume for different domains (veg, snow, soil)
- ! indices
- ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable (nrg)
- ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable (nrg)
- ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass)
- ixSnowOnlyNrg => indx_data%var(iLookINDEX%ixSnowOnlyNrg)%dat ,& ! intent(in): [i4b(:)] indices for energy states in the snow subdomain
- ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the snow+soil subdomain
- ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (iname_nrgLayer...)
- ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] index of the hydrology states in the canopy domain
- ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] index of the type of hydrology states in snow+soil domain
- layerType => indx_data%var(iLookINDEX%layerType)%dat & ! intent(in): [i4b(:)] layer type (iname_soil or iname_snow)
- ) ! association to variables in the data structures
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message="eval8summa/"
-
- ! check the feasibility of the solution
- feasible=.true.
-
- ! check that the canopy air space temperature is reasonable
- if(ixCasNrg/=integerMissing)then
+ private
+ public::eval8summa
+
+ contains
+
+
+ ! **********************************************************************************************************
+ ! public subroutine eval8summa: compute the residual vector and the Jacobian matrix
+ ! **********************************************************************************************************
+ subroutine eval8summa(&
+ ! input: model control
+ dt, & ! intent(in): length of the time step (seconds)
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ nState, & ! intent(in): total number of state variables
+ firstSubStep, & ! intent(in): flag to indicate if we are processing the first sub-step
+ firstFluxCall, & ! intent(inout): flag to indicate if we are processing the first flux call
+ firstSplitOper, & ! intent(in): flag to indicate if we are processing the first flux call in a splitting operation
+ computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation
+ scalarSolution, & ! intent(in): flag to indicate the scalar solution
+ ! input: state vectors
+ stateVecTrial, & ! intent(in): model state vector
+ fScale, & ! intent(in): function scaling vector
+ sMul, & ! intent(in): state vector multiplier (used in the residual calculations)
+ ! input: data structures
+ model_decisions, & ! intent(in): model decisions
+ lookup_data, & ! intent(in): lookup tables
+ type_data, & ! intent(in): type of vegetation and soil
+ attr_data, & ! intent(in): spatial attributes
+ mpar_data, & ! intent(in): model parameters
+ forc_data, & ! intent(in): model forcing data
+ bvar_data, & ! intent(in): average model variables for the entire basin
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ ! input-output: data structures
+ indx_data, & ! intent(inout): index data
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+ deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables
+ ! input-output: baseflow
+ ixSaturation, & ! intent(inout): index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ dBaseflow_dMatric, & ! intent(out): derivative in baseflow w.r.t. matric head (s-1)
+ ! output: flux and residual vectors
+ feasible, & ! intent(out): flag to denote the feasibility of the solution
+ fluxVec, & ! intent(out): flux vector
+ resSink, & ! intent(out): additional (sink) terms on the RHS of the state equation
+ resVec, & ! intent(out): residual vector
+ fEval, & ! intent(out): function evaluation
+ err,message) ! intent(out): error control
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! provide access to subroutines
+ USE getVectorz_module, only:varExtract ! extract variables from the state vector
+ USE updateVars_module, only:updateVars ! update prognostic variables
+ USE t2enthalpy_module, only:t2enthalpy ! compute enthalpy
+ USE computFlux_module, only:soilCmpres ! compute soil compression, use non-sundials version because sundials version needs mLayerMatricHeadPrime
+ USE computFlux_module, only:computFlux ! compute fluxes given a state vector
+ USE computResid_module,only:computResid ! compute residuals given a state vector
+ implicit none
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! input: model control
+ real(rkind),intent(in) :: dt ! length of the time step (seconds)
+ integer(i4b),intent(in) :: nSnow ! number of snow layers
+ integer(i4b),intent(in) :: nSoil ! number of soil layers
+ integer(i4b),intent(in) :: nLayers ! total number of layers
+ integer(i4b),intent(in) :: nState ! total number of state variables
+ logical(lgt),intent(in) :: firstSubStep ! flag to indicate if we are processing the first sub-step
+ logical(lgt),intent(inout) :: firstFluxCall ! flag to indicate if we are processing the first flux call
+ logical(lgt),intent(in) :: firstSplitOper ! flag to indicate if we are processing the first flux call in a splitting operation
+ logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if computing fluxes over vegetation
+ logical(lgt),intent(in) :: scalarSolution ! flag to denote if implementing the scalar solution
+ ! input: state vectors
+ real(rkind),intent(in) :: stateVecTrial(:) ! model state vector
+ real(rkind),intent(in) :: fScale(:) ! function scaling vector
+ real(rkind),intent(in) :: sMul(:) ! NOTE: qp ! state vector multiplier (used in the residual calculations)
+ ! input: data structures
+ type(model_options),intent(in) :: model_decisions(:) ! model decisions
+ type(zLookup), intent(in) :: lookup_data ! lookup tables
+ type(var_i), intent(in) :: type_data ! type of vegetation and soil
+ type(var_d), intent(in) :: attr_data ! spatial attributes
+ type(var_dlength), intent(in) :: mpar_data ! model parameters
+ type(var_d), intent(in) :: forc_data ! model forcing data
+ type(var_dlength), intent(in) :: bvar_data ! model variables for the local basin
+ type(var_dlength), intent(in) :: prog_data ! prognostic variables for a local HRU
+ ! output: data structures
+ type(var_ilength),intent(inout) :: indx_data ! indices defining model states and layers
+ type(var_dlength),intent(inout) :: diag_data ! diagnostic variables for a local HRU
+ type(var_dlength),intent(inout) :: flux_data ! model fluxes for a local HRU
+ type(var_dlength),intent(inout) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
+ ! input-output: baseflow
+ integer(i4b),intent(inout) :: ixSaturation ! index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ real(rkind),intent(out) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1)
+ ! output: flux and residual vectors
+ logical(lgt),intent(out) :: feasible ! flag to denote the feasibility of the solution
+ real(rkind),intent(out) :: fluxVec(:) ! flux vector
+ real(rkind),intent(out) :: resSink(:) ! sink terms on the RHS of the flux equation
+ real(rkind),intent(out) :: resVec(:) ! NOTE: qp ! residual vector
+ real(rkind),intent(out) :: fEval ! function evaluation
+ ! output: error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! local variables
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! state variables
+ real(rkind) :: scalarCanairTempTrial ! trial value for temperature of the canopy air space (K)
+ real(rkind) :: scalarCanopyTempTrial ! trial value for temperature of the vegetation canopy (K)
+ real(rkind) :: scalarCanopyWatTrial ! trial value for liquid water storage in the canopy (kg m-2)
+ real(rkind),dimension(nLayers) :: mLayerTempTrial ! trial value for temperature of layers in the snow and soil domains (K)
+ real(rkind),dimension(nLayers) :: mLayerVolFracWatTrial ! trial value for volumetric fraction of total water (-)
+ real(rkind),dimension(nSoil) :: mLayerMatricHeadTrial ! trial value for total water matric potential (m)
+ real(rkind),dimension(nSoil) :: mLayerMatricHeadLiqTrial ! trial value for liquid water matric potential (m)
+ real(rkind) :: scalarAquiferStorageTrial ! trial value of storage of water in the aquifer (m)
+ ! diagnostic variables
+ logical(lgt),parameter :: needEnthalpy=.true. ! flag to compute enthalpy
+ real(rkind) :: scalarCanopyLiqTrial ! trial value for mass of liquid water on the vegetation canopy (kg m-2)
+ real(rkind) :: scalarCanopyIceTrial ! trial value for mass of ice on the vegetation canopy (kg m-2)
+ real(rkind),dimension(nLayers) :: mLayerVolFracLiqTrial ! trial value for volumetric fraction of liquid water (-)
+ real(rkind),dimension(nLayers) :: mLayerVolFracIceTrial ! trial value for volumetric fraction of ice (-)
+ ! other local variables
+ integer(i4b) :: iLayer ! index of model layer in the snow+soil domain
+ integer(i4b) :: jState(1) ! index of model state for the scalar solution within the soil domain
+ integer(i4b) :: ixBeg,ixEnd ! index of indices for the soil compression routine
+ integer(i4b),parameter :: ixVegVolume=1 ! index of the desired vegetation control volumne (currently only one veg layer)
+ real(rkind) :: xMin,xMax ! minimum and maximum values for water content
+ real(rkind) :: scalarCanopyHydTrial ! trial value for mass of water on the vegetation canopy (kg m-2)
+ real(rkind),parameter :: canopyTempMax=500._rkind ! expected maximum value for the canopy temperature (K)
+ real(rkind),dimension(nLayers) :: mLayerVolFracHydTrial ! trial value for volumetric fraction of water (-), general vector merged from Wat and Liq
+ real(rkind),dimension(nState) :: rVecScaled ! scaled residual vector
+ character(LEN=256) :: cmessage ! error message of downwind routine
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! association to variables in the data structures
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ associate(&
+ ! model decisions
+ ixRichards => model_decisions(iLookDECISIONS%f_Richards)%iDecision ,& ! intent(in): [i4b] index of the form of Richards' equation
+ ! snow parameters
+ snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) ,& ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1)
+ ! soil parameters
+ 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(:)] residual volumetric water content (-)
+ specificStorage => mpar_data%var(iLookPARAM%specificStorage)%dat(1) ,& ! intent(in): [dp] specific storage coefficient (m-1)
+ ! canopy and layer depth
+ canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp ] canopy depth (m)
+ mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
+ ! model state variables from a previous solution
+ mLayerMatricHead => prog_data%var(iLookPROG%mLayerMatricHead)%dat ,& ! intent(in): [dp(:)] total water matric potential (m)
+ ! model diagnostic variables from a previous solution
+ scalarSfcMeltPond => prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1) ,& ! intent(in): [dp] ponded water caused by melt of the "snow without a layer" (kg m-2)
+ mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(in): [dp(:)] volumetric fraction of ice (-)
+ ! enthalpy
+ scalarCanairEnthalpy => diag_data%var(iLookDIAG%scalarCanairEnthalpy)%dat(1) ,& ! intent(out): [dp] enthalpy of the canopy air space (J m-3)
+ scalarCanopyEnthalpy => diag_data%var(iLookDIAG%scalarCanopyEnthalpy)%dat(1) ,& ! intent(out): [dp] enthalpy of the vegetation canopy (J m-3)
+ mLayerEnthalpy => diag_data%var(iLookDIAG%mLayerEnthalpy)%dat ,& ! intent(out): [dp(:)] enthalpy of the snow+soil layers (J m-3)
+ ! soil compression
+ scalarSoilCompress => diag_data%var(iLookDIAG%scalarSoilCompress)%dat(1) ,& ! intent(in): [dp] total change in storage associated with compression of the soil matrix (kg m-2)
+ mLayerCompress => diag_data%var(iLookDIAG%mLayerCompress)%dat ,& ! intent(in): [dp(:)] change in storage associated with compression of the soil matrix (-)
+ ! derivatives
+ dVolTot_dPsi0 => deriv_data%var(iLookDERIV%dVolTot_dPsi0)%dat ,& ! intent(in): [dp(:)] derivative in total water content w.r.t. total water matric potential
+ dCompress_dPsi => deriv_data%var(iLookDERIV%dCompress_dPsi)%dat ,& ! intent(in): [dp(:)] derivative in compressibility w.r.t. matric head (m-1)
+ ! mapping
+ ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] mapping of full state vector to the state subset
+ ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of control volume for different domains (veg, snow, soil)
+ ! indices
+ ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable (nrg)
+ ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable (nrg)
+ ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass)
+ ixSnowOnlyNrg => indx_data%var(iLookINDEX%ixSnowOnlyNrg)%dat ,& ! intent(in): [i4b(:)] indices for energy states in the snow subdomain
+ ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the snow+soil subdomain
+ ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (iname_nrgLayer...)
+ ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] index of the hydrology states in the canopy domain
+ ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] index of the type of hydrology states in snow+soil domain
+ layerType => indx_data%var(iLookINDEX%layerType)%dat & ! intent(in): [i4b(:)] layer type (iname_soil or iname_snow)
+ ) ! association to variables in the data structures
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message="eval8summa/"
+
+ ! check the feasibility of the solution
+ feasible=.true.
+
+ ! check that the canopy air space temperature is reasonable
+ if(ixCasNrg/=integerMissing)then
if(stateVecTrial(ixCasNrg) > canopyTempMax) feasible=.false.
- endif
-
- ! check that the canopy air space temperature is reasonable
- if(ixVegNrg/=integerMissing)then
+ if(stateVecTrial(ixCasNrg) > canopyTempMax) message=trim(message)//'canopy air space temp high,'
+ if(.not.feasible) write(*,'(a,1x,L1,1x,10(f20.10,1x))') 'feasible, max, stateVecTrial( ixCasNrg )', feasible, canopyTempMax, stateVecTrial(ixCasNrg)
+ endif
+
+ ! check that the canopy air space temperature is reasonable
+ if(ixVegNrg/=integerMissing)then
if(stateVecTrial(ixVegNrg) > canopyTempMax) feasible=.false.
- endif
-
- ! check canopy liquid water is not negative
- if(ixVegHyd/=integerMissing)then
- if(stateVecTrial(ixVegHyd) < 0._dp) feasible=.false.
- end if
-
- ! check snow temperature is below freezing
- if(count(ixSnowOnlyNrg/=integerMissing)>0)then
+ if(stateVecTrial(ixVegNrg) > canopyTempMax) message=trim(message)//'canopy temp high,'
+ if(.not.feasible) write(*,'(a,1x,L1,1x,10(f20.10,1x))') 'feasible, max, stateVecTrial( ixVegNrg )', feasible, canopyTempMax, stateVecTrial(ixVegNrg)
+ endif
+
+ ! check canopy liquid water is not negative
+ if(ixVegHyd/=integerMissing)then
+ if(stateVecTrial(ixVegHyd) < 0._rkind) feasible=.false.
+ if(stateVecTrial(ixVegHyd) < 0._rkind) message=trim(message)//'canopy water negative,'
+ if(.not.feasible) write(*,'(a,1x,L1,1x,10(f20.10,1x))') 'feasible, min, stateVecTrial( ixVegHyd )', feasible, 0._rkind, stateVecTrial(ixVegHyd)
+
+ end if
+
+ ! check snow temperature is below freezing
+ if(count(ixSnowOnlyNrg/=integerMissing)>0)then
if(any(stateVecTrial( pack(ixSnowOnlyNrg,ixSnowOnlyNrg/=integerMissing) ) > Tfreeze)) feasible=.false.
- endif
-
- ! loop through non-missing hydrology state variables in the snow+soil domain
- do concurrent (iLayer=1:nLayers,ixSnowSoilHyd(iLayer)/=integerMissing)
-
+ if(any(stateVecTrial( pack(ixSnowOnlyNrg,ixSnowOnlyNrg/=integerMissing) ) > Tfreeze)) message=trim(message)//'snow temp above freezing,'
+ do iLayer=1,nSnow
+ if(.not.feasible) write(*,'(a,1x,i4,1x,L1,1x,10(f20.10,1x))') 'iLayer, feasible, max, stateVecTrial( ixSnowOnlyNrg(iLayer) )', iLayer, feasible, Tfreeze, stateVecTrial( ixSnowOnlyNrg(iLayer) )
+ enddo
+ endif
+
+ ! loop through non-missing hydrology state variables in the snow+soil domain
+ do concurrent (iLayer=1:nLayers,ixSnowSoilHyd(iLayer)/=integerMissing)
+
! check the minimum and maximum water constraints
if(ixHydType(iLayer)==iname_watLayer .or. ixHydType(iLayer)==iname_liqLayer)then
-
- ! --> minimum
- if (layerType(iLayer) == iname_soil) then
- xMin = theta_sat(iLayer-nSnow)
- else
- xMin = 0._dp
- endif
-
- ! --> maximum
- select case( layerType(iLayer) )
- case(iname_snow); xMax = merge(iden_ice, 1._dp - mLayerVolFracIce(iLayer), ixHydType(iLayer)==iname_watLayer)
- case(iname_soil); xMax = merge(theta_sat(iLayer-nSnow), theta_sat(iLayer-nSnow) - mLayerVolFracIce(iLayer), ixHydType(iLayer)==iname_watLayer)
- end select
-
- ! --> check
- if(stateVecTrial( ixSnowSoilHyd(iLayer) ) < xMin .or. stateVecTrial( ixSnowSoilHyd(iLayer) ) > xMax) feasible=.false.
- !if(.not.feasible) write(*,'(a,1x,i4,1x,L1,1x,10(f20.10,1x))') 'iLayer, feasible, stateVecTrial( ixSnowSoilHyd(iLayer) ), xMin, xMax = ', iLayer, feasible, stateVecTrial( ixSnowSoilHyd(iLayer) ), xMin, xMax
-
+
+ ! --> minimum
+ if (layerType(iLayer) == iname_soil) then
+ xMin = theta_res(iLayer-nSnow)
+ else
+ xMin = 0._rkind
+ endif
+
+ ! --> maximum
+ select case( layerType(iLayer) )
+ case(iname_snow); xMax = merge(iden_ice, 1._rkind - mLayerVolFracIce(iLayer), ixHydType(iLayer)==iname_watLayer)
+ case(iname_soil); xMax = merge(theta_sat(iLayer-nSnow), theta_sat(iLayer-nSnow) - mLayerVolFracIce(iLayer), ixHydType(iLayer)==iname_watLayer)
+ end select
+
+ ! --> check
+ if(stateVecTrial( ixSnowSoilHyd(iLayer) ) < xMin .or. stateVecTrial( ixSnowSoilHyd(iLayer) ) > xMax) feasible=.false.
+ if(stateVecTrial( ixSnowSoilHyd(iLayer) ) < xMin .or. stateVecTrial( ixSnowSoilHyd(iLayer) ) > xMax) message=trim(message)//'layer water outside bounds,'
+ if(.not.feasible) write(*,'(a,1x,i4,1x,L1,1x,10(f20.10,1x))') 'iLayer, feasible, stateVecTrial( ixSnowSoilHyd(iLayer) ), xMin, xMax = ', iLayer, feasible, stateVecTrial( ixSnowSoilHyd(iLayer) ), xMin, xMax
+
endif ! if water states
-
- end do ! loop through non-missing hydrology state variables in the snow+soil domain
-
- ! early return for non-feasible solutions
- if(.not.feasible)then
+
+ end do ! loop through non-missing hydrology state variables in the snow+soil domain
+
+ ! early return for non-feasible solutions
+ if(.not.feasible)then
fluxVec(:) = realMissing
resVec(:) = quadMissing
fEval = realMissing
- return
- end if
-
- ! get the start and end indices for the soil compression calculations
- if(scalarSolution)then
+ message=trim(message)//'non-feasible'
+ err=20; return
+ endif
+
+ ! get the start and end indices for the soil compression calculations
+ if(scalarSolution)then
jState = pack(ixControlVolume, ixMapFull2Subset/=integerMissing)
ixBeg = jState(1)
ixEnd = jState(1)
- else
+ else
ixBeg = 1
ixEnd = nSoil
- endif
-
- ! extract variables from the model state vector
- call varExtract(&
- ! input
- stateVecTrial, & ! intent(in): model state vector (mixed units)
- diag_data, & ! intent(in): model diagnostic variables for a local HRU
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- indx_data, & ! intent(in): indices defining model states and layers
- ! output: variables for the vegetation canopy
- scalarCanairTempTrial, & ! intent(out): trial value of canopy air temperature (K)
- scalarCanopyTempTrial, & ! intent(out): trial value of canopy temperature (K)
- scalarCanopyWatTrial, & ! intent(out): trial value of canopy total water (kg m-2)
- scalarCanopyLiqTrial, & ! intent(out): trial value of canopy liquid water (kg m-2)
- scalarCanopyIceTrial, & ! intent(out): trial value of canopy ice content (kg m-2)
- ! output: variables for the snow-soil domain
- mLayerTempTrial, & ! intent(out): trial vector of layer temperature (K)
- mLayerVolFracWatTrial, & ! intent(out): trial vector of volumetric total water content (-)
- mLayerVolFracLiqTrial, & ! intent(out): trial vector of volumetric liquid water content (-)
- mLayerVolFracIceTrial, & ! intent(out): trial vector of volumetric ice water content (-)
- mLayerMatricHeadTrial, & ! intent(out): trial vector of total water matric potential (m)
- mLayerMatricHeadLiqTrial, & ! intent(out): trial vector of liquid water matric potential (m)
- ! output: variables for the aquifer
- scalarAquiferStorageTrial,& ! intent(out): trial value of storage of water in the aquifer (m)
- ! output: error control
- err,cmessage) ! intent(out): error control
- if(err/=0)then
- message=trim(message)//trim(cmessage)
- print*, message
- return
- end if ! (check for errors)
-
- ! update diagnostic variables
- call updateVars(&
- ! input
- .false., & ! 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
- 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)
- ! 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)
- ! output: error control
- err,cmessage) ! intent(out): error control
- if(err/=0)then
- message=trim(message)//trim(cmessage)
- print*, message
- return
- end if ! (check for errors)
-
- ! compute enthalpy (J m-3)
- if(needEnthalpy)then
+ endif
+
+ ! extract variables from the model state vector
+ call varExtract(&
+ ! input
+ stateVecTrial, & ! intent(in): model state vector (mixed units)
+ diag_data, & ! intent(in): model diagnostic variables for a local HRU
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ indx_data, & ! intent(in): indices defining model states and layers
+ ! output: variables for the vegetation canopy
+ scalarCanairTempTrial, & ! intent(out): trial value of canopy air temperature (K)
+ scalarCanopyTempTrial, & ! intent(out): trial value of canopy temperature (K)
+ scalarCanopyWatTrial, & ! intent(out): trial value of canopy total water (kg m-2)
+ scalarCanopyLiqTrial, & ! intent(out): trial value of canopy liquid water (kg m-2)
+ scalarCanopyIceTrial, & ! intent(out): trial value of canopy ice content (kg m-2)
+ ! output: variables for the snow-soil domain
+ mLayerTempTrial, & ! intent(out): trial vector of layer temperature (K)
+ mLayerVolFracWatTrial, & ! intent(out): trial vector of volumetric total water content (-)
+ mLayerVolFracLiqTrial, & ! intent(out): trial vector of volumetric liquid water content (-)
+ mLayerVolFracIceTrial, & ! intent(out): trial vector of volumetric ice water content (-)
+ mLayerMatricHeadTrial, & ! intent(out): trial vector of total water matric potential (m)
+ mLayerMatricHeadLiqTrial, & ! intent(out): trial vector of liquid water matric potential (m)
+ ! output: variables for the aquifer
+ scalarAquiferStorageTrial,& ! intent(out): trial value of storage of water in the aquifer (m)
+ ! output: error control
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors)
+
+ ! update diagnostic variables and derivatives
+ call updateVars(&
+ ! input
+ .false., & ! intent(in): logical flag to adjust temperature to account for the energy used in melt+freeze
+ lookup_data, & ! intent(in): lookup tables for a local HRU
+ 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
+ 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)
+ ! 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)
+ ! output: error control
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors)
+
+ ! compute enthalpy (J m-3)
+ if(needEnthalpy)then
call t2enthalpy(&
! input: data structures
diag_data, & ! intent(in): model diagnostic variables for a local HRU
@@ -426,126 +431,114 @@ subroutine eval8summa(&
mLayerEnthalpy, & ! intent(out): enthalpy of each snow+soil layer (J m-3)
! output: error control
err,cmessage) ! intent(out): error control
- if(err/=0)then
- message=trim(message)//trim(cmessage)
- print*, message
- return
- endif
- endif ! if computing enthalpy
-
-
- ! print the states in the canopy domain
- !print*, 'dt = ', dt
- !write(*,'(a,1x,10(f20.10,1x))') 'scalarCanopyTempTrial = ', scalarCanopyTempTrial
- !write(*,'(a,1x,10(f20.10,1x))') 'scalarCanopyWatTrial = ', scalarCanopyWatTrial
- !write(*,'(a,1x,10(f20.10,1x))') 'scalarCanopyLiqTrial = ', scalarCanopyLiqTrial
- !write(*,'(a,1x,10(f20.10,1x))') 'scalarCanopyIceTrial = ', scalarCanopyIceTrial
-
- ! print the states in the snow+soil domain
- !write(*,'(a,1x,10(f20.10,1x))') 'mLayerTempTrial = ', mLayerTempTrial(iJac1:min(nLayers,iJac2))
- !write(*,'(a,1x,10(f20.10,1x))') 'mLayerVolFracWatTrial = ', mLayerVolFracWatTrial(iJac1:min(nLayers,iJac2))
- !write(*,'(a,1x,10(f20.10,1x))') 'mLayerVolFracLiqTrial = ', mLayerVolFracLiqTrial(iJac1:min(nLayers,iJac2))
- !write(*,'(a,1x,10(f20.10,1x))') 'mLayerVolFracIceTrial = ', mLayerVolFracIceTrial(iJac1:min(nLayers,iJac2))
- !write(*,'(a,1x,10(f20.10,1x))') 'mLayerMatricHeadTrial = ', mLayerMatricHeadTrial(iJac1:min(nSoil,iJac2))
- !write(*,'(a,1x,10(f20.10,1x))') 'mLayerMatricHeadLiqTrial = ', mLayerMatricHeadLiqTrial(iJac1:min(nSoil,iJac2))
-
- ! print the water content
- if(globalPrintFlag)then
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+ endif ! if computing enthalpy
+
+ ! print the states in the canopy domain
+ !print*, 'dt = ', dt
+ !write(*,'(a,1x,10(f20.10,1x))') 'scalarCanopyTempTrial = ', scalarCanopyTempTrial
+ !write(*,'(a,1x,10(f20.10,1x))') 'scalarCanopyWatTrial = ', scalarCanopyWatTrial
+ !write(*,'(a,1x,10(f20.10,1x))') 'scalarCanopyLiqTrial = ', scalarCanopyLiqTrial
+ !write(*,'(a,1x,10(f20.10,1x))') 'scalarCanopyIceTrial = ', scalarCanopyIceTrial
+
+ ! print the states in the snow+soil domain
+ !write(*,'(a,1x,10(f20.10,1x))') 'mLayerTempTrial = ', mLayerTempTrial(iJac1:min(nLayers,iJac2))
+ !write(*,'(a,1x,10(f20.10,1x))') 'mLayerVolFracWatTrial = ', mLayerVolFracWatTrial(iJac1:min(nLayers,iJac2))
+ !write(*,'(a,1x,10(f20.10,1x))') 'mLayerVolFracLiqTrial = ', mLayerVolFracLiqTrial(iJac1:min(nLayers,iJac2))
+ !write(*,'(a,1x,10(f20.10,1x))') 'mLayerVolFracIceTrial = ', mLayerVolFracIceTrial(iJac1:min(nLayers,iJac2))
+ !write(*,'(a,1x,10(f20.10,1x))') 'mLayerMatricHeadTrial = ', mLayerMatricHeadTrial(iJac1:min(nSoil,iJac2))
+ !write(*,'(a,1x,10(f20.10,1x))') 'mLayerMatricHeadLiqTrial = ', mLayerMatricHeadLiqTrial(iJac1:min(nSoil,iJac2))
+
+ ! print the water content
+ if(globalPrintFlag)then
if(iJac1<nSnow) write(*,'(a,10(f16.10,1x))') 'mLayerVolFracWatTrial = ', mLayerVolFracWatTrial(iJac1:min(iJac2,nSnow))
if(iJac1<nSnow) write(*,'(a,10(f16.10,1x))') 'mLayerVolFracLiqTrial = ', mLayerVolFracLiqTrial(iJac1:min(iJac2,nSnow))
if(iJac1<nSnow) write(*,'(a,10(f16.10,1x))') 'mLayerVolFracIceTrial = ', mLayerVolFracIceTrial(iJac1:min(iJac2,nSnow))
- endif
-
- ! save the number of flux calls per time step
- indx_data%var(iLookINDEX%numberFluxCalc)%dat(1) = indx_data%var(iLookINDEX%numberFluxCalc)%dat(1) + 1
-
- ! compute the fluxes for a given state vector
- call computFlux(&
- ! input-output: model control
- nSnow, & ! intent(in): number of snow layers
- nSoil, & ! intent(in): number of soil layers
- nLayers, & ! intent(in): total number of layers
- firstSubStep, & ! intent(in): flag to indicate if we are processing the first sub-step
- firstFluxCall, & ! intent(inout): flag to denote the first flux call
- firstSplitOper, & ! intent(in): flag to indicate if we are processing the first flux call in a splitting operation
- computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation
- scalarSolution, & ! intent(in): flag to indicate the scalar solution
- .true., & ! intent(in): balance longwave
- scalarSfcMeltPond/dt, & ! intent(in): drainage from the surface melt pond (kg m-2 s-1)
- ! input: state variables
- scalarCanairTempTrial, & ! intent(in): trial value for the temperature of the canopy air space (K)
- scalarCanopyTempTrial, & ! intent(in): trial value for the temperature of the vegetation canopy (K)
- mLayerTempTrial, & ! intent(in): trial value for the temperature of each snow and soil layer (K)
- mLayerMatricHeadLiqTrial, & ! intent(in): trial value for the liquid water matric potential in each soil layer (m)
- mLayerMatricHeadTrial, & ! intent(in): trial vector of total water matric potential (m)
- scalarAquiferStorageTrial, & ! intent(in): trial value of storage of water in the aquifer (m)
- ! input: diagnostic variables defining the liquid water and ice content
- scalarCanopyLiqTrial, & ! intent(in): trial value for the liquid water on the vegetation canopy (kg m-2)
- scalarCanopyIceTrial, & ! intent(in): trial value for the ice on the vegetation canopy (kg m-2)
- mLayerVolFracLiqTrial, & ! intent(in): trial value for the volumetric liquid water content in each snow and soil layer (-)
- mLayerVolFracIceTrial, & ! intent(in): trial value for the volumetric ice in each snow and soil layer (-)
- ! input: data structures
- model_decisions, & ! intent(in): model decisions
- type_data, & ! intent(in): type of vegetation and soil
- attr_data, & ! intent(in): spatial attributes
- mpar_data, & ! intent(in): model parameters
- forc_data, & ! intent(in): model forcing data
- bvar_data, & ! intent(in): average model variables for the entire basin
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- indx_data, & ! intent(in): index data
- ! input-output: data structures
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- flux_data, & ! intent(inout): model fluxes for a local HRU
- deriv_data, & ! intent(out): derivatives in model fluxes w.r.t. relevant state variables
- ! input-output: flux vector and baseflow derivatives
- ixSaturation, & ! intent(inout): index of the lowest saturated layer (NOTE: only computed on the first iteration)
- dBaseflow_dMatric, & ! intent(out): derivative in baseflow w.r.t. matric head (s-1)
- fluxVec, & ! intent(out): flux vector (mixed units)
- ! output: error control
- err,cmessage) ! intent(out): error code and error message
- if(err/=0)then
- message=trim(message)//trim(cmessage)
- print*, message
- return
- end if ! (check for errors)
-
- ! compute soil compressibility (-) and its derivative w.r.t. matric head (m)
- ! NOTE: we already extracted trial matrix head and volumetric liquid water as part of the flux calculations
- call soilCmpres(&
- ! input:
- ixRichards, & ! intent(in): choice of option for Richards' equation
- ixBeg,ixEnd, & ! intent(in): start and end indices defining desired layers
- mLayerMatricHeadLiq(1:nSoil), & ! intent(in): matric head at the start of the time step (m)
- mLayerMatricHeadLiqTrial(1:nSoil), & ! intent(in): trial value of matric head (m)
- mLayerVolFracLiqTrial(nSnow+1:nLayers), & ! intent(in): trial value for the volumetric liquid water content in each soil layer (-)
- mLayerVolFracIceTrial(nSnow+1:nLayers), & ! intent(in): trial value for the volumetric ice content in each soil layer (-)
- specificStorage, & ! intent(in): specific storage coefficient (m-1)
- theta_sat, & ! intent(in): soil porosity (-)
- ! output:
- mLayerCompress, & ! intent(inout): compressibility of the soil matrix (-)
- dCompress_dPsi, & ! intent(inout): derivative in compressibility w.r.t. matric head (m-1)
- err,cmessage) ! intent(out): error code and error message
- if(err/=0)then
- message=trim(message)//trim(cmessage)
- print*, message
- return
- end if ! (check for errors)
-
- ! compute the total change in storage associated with compression of the soil matrix (kg m-2)
- scalarSoilCompress = sum(mLayerCompress(1:nSoil)*mLayerDepth(nSnow+1:nLayers))*iden_water
-
- ! vegetation domain: get the correct water states (total water, or liquid water, depending on the state type)
- if(computeVegFlux)then
+ endif
+
+ ! save the number of flux calls per time step
+ indx_data%var(iLookINDEX%numberFluxCalc)%dat(1) = indx_data%var(iLookINDEX%numberFluxCalc)%dat(1) + 1
+
+ ! compute the fluxes for a given state vector
+ call computFlux(&
+ ! input-output: model control
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ firstSubStep, & ! intent(in): flag to indicate if we are processing the first sub-step
+ firstFluxCall, & ! intent(inout): flag to denote the first flux call
+ firstSplitOper, & ! intent(in): flag to indicate if we are processing the first flux call in a splitting operation
+ computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation
+ scalarSolution, & ! intent(in): flag to indicate the scalar solution
+ .false., & ! intent(in): not inside Sundials solver loop
+ scalarSfcMeltPond/dt, & ! intent(in): drainage from the surface melt pond (kg m-2 s-1)
+ ! input: state variables
+ scalarCanairTempTrial, & ! intent(in): trial value for the temperature of the canopy air space (K)
+ scalarCanopyTempTrial, & ! intent(in): trial value for the temperature of the vegetation canopy (K)
+ mLayerTempTrial, & ! intent(in): trial value for the temperature of each snow and soil layer (K)
+ mLayerMatricHeadLiqTrial, & ! intent(in): trial value for the liquid water matric potential in each soil layer (m)
+ mLayerMatricHeadTrial, & ! intent(in): trial vector of total water matric potential (m)
+ scalarAquiferStorageTrial, & ! intent(in): trial value of storage of water in the aquifer (m)
+ ! input: diagnostic variables defining the liquid water and ice content
+ scalarCanopyLiqTrial, & ! intent(in): trial value for the liquid water on the vegetation canopy (kg m-2)
+ scalarCanopyIceTrial, & ! intent(in): trial value for the ice on the vegetation canopy (kg m-2)
+ mLayerVolFracLiqTrial, & ! intent(in): trial value for the volumetric liquid water content in each snow and soil layer (-)
+ mLayerVolFracIceTrial, & ! intent(in): trial value for the volumetric ice in each snow and soil layer (-)
+ ! input: data structures
+ model_decisions, & ! intent(in): model decisions
+ type_data, & ! intent(in): type of vegetation and soil
+ attr_data, & ! intent(in): spatial attributes
+ mpar_data, & ! intent(in): model parameters
+ forc_data, & ! intent(in): model forcing data
+ bvar_data, & ! intent(in): average model variables for the entire basin
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ indx_data, & ! intent(in): index data
+ ! input-output: data structures
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+ deriv_data, & ! intent(out): derivatives in model fluxes w.r.t. relevant state variables
+ ! input-output: flux vector and baseflow derivatives
+ ixSaturation, & ! intent(inout): index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ dBaseflow_dMatric, & ! intent(out): derivative in baseflow w.r.t. matric head (s-1)
+ fluxVec, & ! intent(out): flux vector (mixed units)
+ ! output: error control
+ err,cmessage) ! intent(out): error code and error message
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors)
+
+ ! compute soil compressibility (-) and its derivative w.r.t. matric head (m)
+ ! NOTE: we already extracted trial matrix head and volumetric liquid water as part of the flux calculations
+ ! use non-sundials version because sundials version needs mLayerMatricHeadPrime
+ call soilCmpres(&
+ ! input:
+ ixRichards, & ! intent(in): choice of option for Richards' equation
+ ixBeg,ixEnd, & ! intent(in): start and end indices defining desired layers
+ mLayerMatricHead(1:nSoil), & ! intent(in): matric head at the start of the time step (m)
+ mLayerMatricHeadTrial(1:nSoil), & ! intent(in): trial value of matric head (m)
+ mLayerVolFracLiqTrial(nSnow+1:nLayers), & ! intent(in): trial value for the volumetric liquid water content in each soil layer (-)
+ mLayerVolFracIceTrial(nSnow+1:nLayers), & ! intent(in): trial value for the volumetric ice content in each soil layer (-)
+ specificStorage, & ! intent(in): specific storage coefficient (m-1)
+ theta_sat, & ! intent(in): soil porosity (-)
+ ! output:
+ mLayerCompress, & ! intent(inout): compressibility of the soil matrix (-)
+ dCompress_dPsi, & ! intent(inout): derivative in compressibility w.r.t. matric head (m-1)
+ err,cmessage) ! intent(out): error code and error message
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors)
+
+ ! compute the total change in storage associated with compression of the soil matrix (kg m-2)
+ scalarSoilCompress = sum(mLayerCompress(1:nSoil)*mLayerDepth(nSnow+1:nLayers))*iden_water
+
+ ! vegetation domain: get the correct water states (total water, or liquid water, depending on the state type)
+ if(computeVegFlux)then
scalarCanopyHydTrial = merge(scalarCanopyWatTrial, scalarCanopyLiqTrial, (ixStateType( ixHydCanopy(ixVegVolume) )==iname_watCanopy) )
- else
+ else
scalarCanopyHydTrial = realMissing
- endif
-
- ! snow+soil domain: get the correct water states (total water, or liquid water, depending on the state type)
- mLayerVolFracHydTrial = merge(mLayerVolFracWatTrial, mLayerVolFracLiqTrial, (ixHydType==iname_watLayer .or. ixHydType==iname_matLayer) )
-
- ! compute the residual vector
- call computResid(&
+ endif
+
+ ! snow+soil domain: get the correct water states (total water, or liquid water, depending on the state type)
+ mLayerVolFracHydTrial = merge(mLayerVolFracWatTrial, mLayerVolFracLiqTrial, (ixHydType==iname_watLayer .or. ixHydType==iname_matLayer) )
+
+ ! compute the residual vector
+ call computResid(&
! input: model control
dt, & ! intent(in): length of the time step (seconds)
nSnow, & ! intent(in): number of snow layers
@@ -573,18 +566,16 @@ subroutine eval8summa(&
resSink, & ! intent(out): additional (sink) terms on the RHS of the state equation
resVec, & ! intent(out): residual vector
err,cmessage) ! intent(out): error control
- if(err/=0)then
- message=trim(message)//trim(cmessage)
- print*, message
- return
- end if ! (check for errors)
-
- ! compute the function evaluation
- rVecScaled = fScale(:)*real(resVec(:), dp) ! scale the residual vector (NOTE: residual vector is in quadruple precision)
- fEval = 0.5_dp*dot_product(rVecScaled,rVecScaled)
-
- ! end association with the information in the data structures
- end associate
-
-end subroutine eval8summa
-end module eval8summa_module
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors)
+
+ ! compute the function evaluation
+ rVecScaled = fScale(:)*real(resVec(:), rkind) ! scale the residual vector (NOTE: residual vector is in quadruple precision)
+ fEval = 0.5_rkind*dot_product(rVecScaled,rVecScaled)
+
+ ! end association with the information in the data structures
+ end associate
+
+ end subroutine eval8summa
+
+ end module eval8summa_module
+
\ No newline at end of file
diff --git a/build/source/engine/updateVars.f90 b/build/source/engine/updateVars.f90
index 5b91649e31812cd3fe7dc6afa89f988d803c112c..eadc29a8319848345fcb20381d1103b488f20d91 100755
--- a/build/source/engine/updateVars.f90
+++ b/build/source/engine/updateVars.f90
@@ -1,5 +1,5 @@
! SUMMA - Structure for Unifying Multiple Modeling Alternatives
-! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
+! Copyright (C) 2014-2015 NCAR/RAL
!
! This file is part of SUMMA
!
@@ -20,716 +20,740 @@
module updateVars_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 (dp)
- var_ilength, & ! data vector with variable length dimension (i4b)
- var_dlength, & ! data vector with variable length dimension (dp)
- zLookup
-
-! 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 updatState_module,only:updateSnow ! update snow states
-USE updatState_module,only:updateSoil ! 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_utils_module,only:liquidHead ! compute the liquid water matric potential
-
-! IEEE check
-USE, intrinsic :: ieee_arithmetic ! check values (NaN, etc.)
-
-implicit none
-private
-public::updateVars
-
-contains
-
- ! **********************************************************************************************************
- ! public subroutine updateVars: compute diagnostic variables
- ! **********************************************************************************************************
- subroutine updateVars(&
- ! input
- 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
- 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)
- ! 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)
- ! output: error control
- err,message) ! intent(out): error control
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------------------------------------------------------
- implicit none
- ! input
- 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
- 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(dp),intent(inout) :: scalarCanopyTempTrial ! trial value of canopy temperature (K)
- real(dp),intent(inout) :: scalarCanopyWatTrial ! trial value of canopy total water (kg m-2)
- real(dp),intent(inout) :: scalarCanopyLiqTrial ! trial value of canopy liquid water (kg m-2)
- real(dp),intent(inout) :: scalarCanopyIceTrial ! trial value of canopy ice content (kg m-2)
- ! output: variables for the snow-soil domain
- real(dp),intent(inout) :: mLayerTempTrial(:) ! trial vector of layer temperature (K)
- real(dp),intent(inout) :: mLayerVolFracWatTrial(:) ! trial vector of volumetric total water content (-)
- real(dp),intent(inout) :: mLayerVolFracLiqTrial(:) ! trial vector of volumetric liquid water content (-)
- real(dp),intent(inout) :: mLayerVolFracIceTrial(:) ! trial vector of volumetric ice water content (-)
- real(dp),intent(inout) :: mLayerMatricHeadTrial(:) ! trial vector of total water matric potential (m)
- real(dp),intent(inout) :: mLayerMatricHeadLiqTrial(:) ! trial vector of 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(dp) :: scalarVolFracLiq ! volumetric fraction of liquid water (-)
- real(dp) :: scalarVolFracIce ! volumetric fraction of ice (-)
- real(dp) :: Tcrit ! critical soil temperature below which ice exists (K)
- real(dp) :: xTemp ! temporary temperature (K)
- real(dp) :: effSat ! effective saturation (-)
- real(dp) :: 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(dp) :: meltNrg ! energy for melt+freeze (J m-3)
- real(dp) :: residual ! residual in the energy equation (J m-3)
- real(dp) :: derivative ! derivative in the energy equation (J m-3 K-1)
- real(dp) :: 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(dp),parameter :: nrgConvTol=1.e-4_dp ! convergence tolerance for energy (J m-3)
- real(dp),parameter :: tempConvTol=1.e-6_dp ! convergence tolerance for temperature (K)
- real(dp) :: critDiff ! temperature difference from critical (K)
- real(dp) :: tempMin ! minimum bracket for temperature (K)
- real(dp) :: tempMax ! maximum bracket for temperature (K)
- logical(lgt) :: bFlag ! flag to denote that iteration increment was constrained using bi-section
- real(dp),parameter :: epsT=1.e-7_dp ! 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='updateVars/'
-
- ! 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
-
- ! 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
- case(iname_snow); mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer)*iden_ice/iden_water
- case(iname_soil); mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer) ! no volume expansion
- 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._dp,1._dp,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
- mLayerMatricHeadTrial(ixControlIndex) = matricHead(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
- !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))
- ! --> 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))
- 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._dp
- tempMax = xTemp + 10._dp
-
- ! 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_dp*(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._dp,1._dp,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._dp
- case(iname_snow, iname_soil); mLayerdTheta_dTk(iLayer) = 0._dp
- 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 volumetric fraction of liquid water and ice
- call updateSnow(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)
- scalarVolFracLiq, & ! intent(out) : trial volumetric fraction of liquid water (-)
- scalarVolFracIce, & ! intent(out) : trial volumetric fraction if 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
- scalarCanopyIceTrial = (1._dp - scalarFracLiqVeg)*scalarCanopyWatTrial
-
- ! *** snow layers
- case(iname_snow)
-
- ! compute volumetric fraction of liquid water and ice
- call updateSnow(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)
- mLayerVolFracLiqTrial(iLayer), & ! intent(out) : trial volumetric fraction of liquid water (-)
- mLayerVolFracIceTrial(iLayer), & ! intent(out) : trial volumetric fraction if ice (-)
- 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 updateSoil(xTemp, & ! intent(in) : temperature (K)
- mLayerMatricHeadTrial(ixControlIndex), & ! intent(in) : total water matric potential (m)
- vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),theta_sat(ixControlIndex),theta_res(ixControlIndex),vGn_m(ixControlIndex), & ! intent(in) : soil parameters
- 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 (-)
- 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...
- ! ------------------------
-
- ! check the need to adjust temperature
- if(do_adjustTemp)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._dp)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_dp ! 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._dp)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 (include tolerance)
- if(mLayerVolFracLiqTrial(iLayer) > theta_sat(ixControlIndex)+epsT )then
- message=trim(message)//'liquid water greater than porosity'
- print*,'---------------'
- print*,'porosity(theta_sat)=', theta_sat(ixControlIndex)
- print*,'liq water =',mLayerVolFracLiqTrial(iLayer)
- print*,'layer =',iLayer
- print*,'---------------'
- 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._dp ! exact correspondence (psiLiq=psi0)
- dPsiLiq_dTemp(ixControlIndex) = 0._dp ! 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 liquidHead(&
- ! input
- mLayerMatricHeadTrial(ixControlIndex) ,& ! intent(in) : total water matric potential (m)
- 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)
- ! output
- mLayerMatricHeadLiqTrial(ixControlIndex) ,& ! intent(out): liquid water matric potential (m)
- 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 updateVars
-
-
- ! **********************************************************************************************************
- ! 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(dp),intent(in) :: meltNrg ! energy for melt+freeze (J m-3)
- real(dp),intent(in) :: heatCap ! volumetric heat capacity (J m-3 K-1)
- real(dp),intent(in) :: tempInit ! initial temperature (K)
- real(dp),intent(in) :: volFracIceInit ! initial volumetric fraction of ice (-)
- ! input-output: trial values
- real(dp),intent(inout) :: xTemp ! trial value for temperature
- real(dp),intent(in) :: dLiq_dT ! derivative in liquid water content w.r.t. temperature (K-1)
- real(dp),intent(in) :: volFracIceTrial ! trial value for the volumetric fraction of ice (-)
- ! output: residual and derivative
- real(dp),intent(out) :: residual ! residual (J m-3)
- real(dp),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
-
- ! check validity of residual ...
- ! informational only: if nan, the sim will start to error out from calling routine
- if( ieee_is_nan(residual) )then
- print*, '--------'
- print*, 'ERROR: residual is not valid in xTempSolve'
- print*, 'heatCap', heatCap
- print*, 'xTemp', xTemp
- print*, 'tempInit', tempInit
- print*, 'meltNrg', meltNrg
- print*, 'volFracIceTrial', volFracIceTrial
- print*, 'volFracIceInit', volFracIceInit
- print*, 'dLiq_dT', dLiq_dT
- print*, '--------'
- endif
-
- end subroutine xTempSolve
-
-end module updateVars_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)
+ 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)
+ zLookup, & ! data vector with variable length dimension (rkind)
+ 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 the routines to calculate enthalpy
+ USE t2enthalpy_module,only:t2enthalpy
+
+ ! provide access to routines to update states
+ USE updatState_module,only:updateSnow ! update snow states
+ USE updatState_module,only:updateSoil ! 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_utils_module,only:liquidHead ! compute the liquid water matric potential
+
+ ! IEEE check
+ USE, intrinsic :: ieee_arithmetic ! check values (NaN, etc.)
+
+ implicit none
+ private
+ public::updateVars
+
+ contains
+
+ ! **********************************************************************************************************
+ ! public subroutine updateVars: compute diagnostic variables and derivatives
+ ! **********************************************************************************************************
+ subroutine updateVars(&
+ ! input
+ do_adjustTemp, & ! intent(in): logical flag to adjust temperature to account for the energy used in melt+freeze
+ lookup_data, & ! intent(in): lookup tables for a local HRU
+ 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
+ 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)
+ ! 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)
+ ! output: error control
+ err,message) ! intent(out): error control
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ implicit none
+ ! input
+ logical(lgt) ,intent(in) :: do_adjustTemp ! flag to adjust temperature to account for the energy used in melt+freeze
+ type(zLookup), intent(in) :: lookup_data ! lookup tables for a local HRU
+ 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
+ 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)
+ ! 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)
+ ! 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) :: 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
+ ! 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
+ mLayerdTemp_dt => deriv_data%var(iLookDERIV%mLayerdTemp_dt )%dat ,& ! intent(out): [dp(:)] timestep change in layer temperature
+ scalarCanopydTemp_dt => deriv_data%var(iLookDERIV%scalarCanopydTemp_dt )%dat(1) & ! intent(out): [dp ] timestep change in canopy temperature
+ ) ! association with variables in the data structures
+
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------------------------------------------------------
+
+ ! initialize error control
+ err=0; message='updateVars/'
+
+ ! 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
+
+ ! 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
+ case(iname_snow); mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer)*iden_ice/iden_water
+ case(iname_soil); mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer) ! no volume expansion
+ 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
+ mLayerMatricHeadTrial(ixControlIndex) = matricHead(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))
+ !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))
+ ! --> 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))
+ 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 temperature time derivatives
+ select case(ixDomainType)
+ case(iname_veg); scalarCanopydTemp_dt = xTemp - scalarCanopyTemp
+ case(iname_snow, iname_soil); mLayerdTemp_dt(iLayer) = xTemp - mLayerTemp(iLayer)
+ end select
+
+ ! 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
+ select case(ixDomainType)
+ case(iname_veg)
+ fLiq = fracLiquid(xTemp,snowfrz_scale)
+ dVolHtCapBulk_dCanWat = ( -Cp_ice*( fLiq-1._rkind ) + Cp_water*fLiq )/canopyDepth !this is iden_water/(iden_water*canopyDepth)
+ case(iname_snow)
+ 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
+ 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
+ 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
+ ! dVolHtCapBulk_dPsi0 uses the derivative in water retention curve above critical temp w.r.t.state variable dVolTot_dPsi0
+ 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)
+ end select
+
+ ! 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)
+ dVolHtCapBulk_dTkCanopy = iden_water * (-Cp_ice + Cp_water) * dTheta_dTkCanopy !same as snow but there is no derivative in air
+ case(iname_snow)
+ mLayerdTheta_dTk(iLayer) = dFracLiq_dTk(xTemp,snowfrz_scale)*mLayerVolFracWatTrial(iLayer)
+ dVolHtCapBulk_dTk(iLayer) = ( iden_water * (-Cp_ice + Cp_water) + iden_air * (iden_water/iden_ice - 1._rkind) * Cp_air ) * mLayerdTheta_dTk(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))
+ dVolHtCapBulk_dTk(iLayer) = (-iden_ice * Cp_ice + iden_water * Cp_water) * mLayerdTheta_dTk(iLayer)
+ 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; dVolHtCapBulk_dTkCanopy = 0._rkind
+ case(iname_snow, iname_soil); mLayerdTheta_dTk(iLayer) = 0._rkind; dVolHtCapBulk_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 volumetric fraction of liquid water and ice
+ call updateSnow(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)
+ scalarVolFracLiq, & ! intent(out) : trial volumetric fraction of liquid water (-)
+ scalarVolFracIce, & ! intent(out) : trial volumetric fraction if 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
+ scalarCanopyIceTrial = (1._rkind - scalarFracLiqVeg)*scalarCanopyWatTrial !(kg m-2), scalarVolFracIce* iden_ice *canopyDepth
+
+ ! *** snow layers
+ case(iname_snow)
+
+ ! compute volumetric fraction of liquid water and ice
+ call updateSnow(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)
+ mLayerVolFracLiqTrial(iLayer), & ! intent(out) : trial volumetric fraction of liquid water (-)
+ mLayerVolFracIceTrial(iLayer), & ! intent(out) : trial volumetric fraction if ice (-)
+ 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 updateSoil(xTemp, & ! intent(in) : temperature (K)
+ mLayerMatricHeadTrial(ixControlIndex), & ! intent(in) : total water matric potential (m)
+ vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),theta_sat(ixControlIndex),theta_res(ixControlIndex),vGn_m(ixControlIndex), & ! intent(in) : soil parameters
+ 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 (-)
+ 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...
+ ! ------------------------
+
+ ! check the need to adjust temperature
+ if(do_adjustTemp)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 (include tolerance)
+ if(mLayerVolFracLiqTrial(iLayer) > theta_sat(ixControlIndex)+epsT )then
+ message=trim(message)//'liquid water greater than porosity'
+ print*,'---------------'
+ print*,'porosity(theta_sat)=', theta_sat(ixControlIndex)
+ print*,'liq water =',mLayerVolFracLiqTrial(iLayer)
+ print*,'layer =',iLayer
+ print*,'---------------'
+ 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 liquidHead(&
+ ! input
+ mLayerMatricHeadTrial(ixControlIndex) ,& ! intent(in) : total water matric potential (m)
+ 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)
+ ! output
+ mLayerMatricHeadLiqTrial(ixControlIndex) ,& ! intent(out): liquid water matric potential (m)
+ 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 updateVars
+
+
+ ! **********************************************************************************************************
+ ! 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 updateVars_module
+
\ No newline at end of file
diff --git a/build/source/engine/varSubstep.f90 b/build/source/engine/varSubstep.f90
index 030edb69dcbd785e7d6ccbd9814ece6ce2ab63ca..d5dd079ff7e53d72f1c554b0d6fe8a4aa5c185ba 100755
--- a/build/source/engine/varSubstep.f90
+++ b/build/source/engine/varSubstep.f90
@@ -403,7 +403,7 @@ subroutine varSubstep(&
! update prognostic variables
call updateProg(dtSubstep,nSnow,nSoil,nLayers,doAdjustTemp,computeVegFlux,untappedMelt,stateVecTrial,checkMassBalance, & ! input: model control
- mpar_data,indx_data,flux_temp,prog_data,diag_data,deriv_data, & ! input-output: data structures
+ lookup_data,mpar_data,indx_data,flux_temp,prog_data,diag_data,deriv_data, & ! input-output: data structures
waterBalanceError,nrgFluxModified,tooMuchMelt,err,cmessage) ! output: flags and error control
if(err/=0)then
message=trim(message)//trim(cmessage)
@@ -551,7 +551,7 @@ subroutine varSubstep(&
! private subroutine updateProg: update prognostic variables
! **********************************************************************************************************
subroutine updateProg(dt,nSnow,nSoil,nLayers,doAdjustTemp,computeVegFlux,untappedMelt,stateVecTrial,checkMassBalance, & ! input: model control
- mpar_data,indx_data,flux_data,prog_data,diag_data,deriv_data, & ! input-output: data structures
+ lookup_data,mpar_data,indx_data,flux_data,prog_data,diag_data,deriv_data, & ! input-output: data structures
waterBalanceError,nrgFluxModified,tooMuchMelt,err,message) ! output: flags and error control
USE getVectorz_module,only:varExtract ! extract variables from the state vector
USE updateVars_module,only:updateVars ! update prognostic variables
@@ -567,6 +567,7 @@ subroutine varSubstep(&
real(dp) ,intent(in) :: stateVecTrial(:) ! trial state vector (mixed units)
logical(lgt) ,intent(in) :: checkMassBalance ! flag to check the mass balance
! data structures
+ type(zLookup), intent(in) :: lookup_data ! lookup tables
type(var_dlength),intent(in) :: mpar_data ! model parameters
type(var_ilength),intent(in) :: indx_data ! indices for a local HRU
type(var_dlength),intent(inout) :: flux_data ! model fluxes for a local HRU
@@ -713,6 +714,7 @@ subroutine varSubstep(&
call updateVars(&
! input
doAdjustTemp, & ! intent(in): logical flag to adjust temperature to account for the energy used in melt+freeze
+ lookup_data, & ! intent(in): lookup tables for a local HRU
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
diff --git a/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp2_G1-1_timestep.nc b/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp2_G1-1_timestep.nc
index 47a90093471d0726dd9294c067e9bec7fc075edd..7137052004bd3af9cc32bd35ec5952c4530a6f18 100644
Binary files a/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp2_G1-1_timestep.nc and b/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp2_G1-1_timestep.nc differ
diff --git a/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt b/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt
index 25874f090a496bf205b34945ebc0f70159fdebef..a4bcab1640d74292f06bf1205217ac7d3547830e 100644
--- a/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt
+++ b/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt
@@ -1 +1 @@
- Run start time on system: ccyy=2022 - mm=09 - dd=09 - hh=20 - mi=06 - ss=06.185
+ Run start time on system: ccyy=2022 - mm=09 - dd=09 - hh=20 - mi=18 - ss=37.358