diff --git a/build/source/engine/groundwatr.f90 b/build/source/engine/groundwatr.f90
index 0e16b27ae0796c48894cac41884a0f2571da5921..6c9c91744ed1a6f2e045ca27f8908dc27b5d2ae0 100755
--- a/build/source/engine/groundwatr.f90
+++ b/build/source/engine/groundwatr.f90
@@ -28,8 +28,8 @@ USE multiconst,only:iden_water ! density of water (kg m-3)
! derived types to define the data structures
USE data_types,only:&
- var_d, & ! data vector (dp)
- var_dlength ! data vector with variable length dimension (dp)
+ var_d, & ! data vector (rkind)
+ var_dlength ! data vector with variable length dimension (rkind)
! named variables defining elements in the data structures
USE var_lookup,only:iLookATTR ! named variables for structure elements
@@ -40,496 +40,484 @@ USE var_lookup,only:iLookPARAM ! named variables for structure elements
! 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
+ qbaseTopmodel, & ! TOPMODEL-ish baseflow parameterization
+ bigBucket, & ! a big bucket (lumped aquifer model)
+ noExplicit ! no explicit groundwater parameterization
! privacy
implicit none
! constant parameters
-real(dp),parameter :: valueMissing=-9999._dp ! missing value parameter
-real(dp),parameter :: verySmall=epsilon(1.0_dp) ! a very small number (used to avoid divide by zero)
-real(dp),parameter :: dx=1.e-8_dp ! finite difference increment
+real(rkind),parameter :: valueMissing=-9999._rkind ! missing value parameter
+real(rkind),parameter :: verySmall=epsilon(1.0_rkind) ! a very small number (used to avoid divide by zero)
+real(rkind),parameter :: dx=1.e-8_rkind ! finite difference increment
private
public::groundwatr
contains
- ! ************************************************************************************************
- ! public subroutine groundwatr: compute the groundwater sink term in Richards' equation
- ! ************************************************************************************************
- !
- ! Method
- ! ------
- !
- ! Here we assume that water avaialble for shallow groundwater flow includes is all water above
- ! "field capacity" below the depth zCrit, where zCrit is defined as the lowest point in the soil
- ! profile where the volumetric liquid water content is less than field capacity.
- !
- ! We further assume that transmssivity (m2 s-1) for each layer is defined asuming that the water
- ! available for saturated flow is located at the bottom of the soil profile. Specifically:
- ! trTotal(iLayer) = tran0*(zActive(iLayer)/soilDepth)**zScale_TOPMODEL
- ! trSoil(iLayer) = trTotal(iLayer) - trTotal(iLayer+1)
- ! where zActive(iLayer) is the effective water table thickness for all layers up to and including
- ! the current layer (working from the bottom to the top).
- !
- ! The outflow from each layer is then (m3 s-1)
- ! mLayerOutflow(iLayer) = trSoil(iLayer)*tan_slope*contourLength
- ! where contourLength is the width of a hillslope (m) parallel to a stream
- !
- ! ************************************************************************************************
- subroutine groundwatr(&
-
- ! input: model control
- nSnow, & ! intent(in): number of snow layers
- nSoil, & ! intent(in): number of soil layers
- nLayers, & ! intent(in): total number of layers
- getSatDepth, & ! intent(in): logical flag to compute index of the lowest saturated layer
-
- ! input: state and diagnostic variables
- mLayerdTheta_dPsi, & ! intent(in): derivative in the soil water characteristic w.r.t. matric head in each layer (m-1)
- mLayerMatricHeadLiq, & ! intent(in): liquid water matric potential (m)
- mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water (-)
- mLayerVolFracIce, & ! intent(in): volumetric fraction of ice (-)
-
- ! input/output: data structures
- attr_data, & ! intent(in): spatial attributes
- mpar_data, & ! intent(in): model parameters
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- diag_data, & ! intent(in): model diagnostic variables for a local HRU
- flux_data, & ! intent(inout): model fluxes for a local HRU
-
- ! output: baseflow
- ixSaturation, & ! intent(inout) index of lowest saturated layer (NOTE: only computed on the first iteration)
- mLayerBaseflow, & ! intent(out): baseflow from each soil layer (m s-1)
- dBaseflow_dMatric, & ! intent(out): derivative in baseflow w.r.t. matric head (s-1)
-
- ! output: error control
- err,message) ! intent(out): error control
- ! ---------------------------------------------------------------------------------------
- ! utility modules
- USE soil_utils_module,only:volFracLiq ! compute volumetric fraction of liquid water as a function of matric head
- USE soil_utils_module,only:hydCond_psi ! compute hydraulic conductivity as a function of matric head
- implicit none
- ! ---------------------------------------------------------------------------------------
- ! * dummy variables
- ! ---------------------------------------------------------------------------------------
- ! input: model control
- integer(i4b),intent(in) :: nSnow ! number of snow layers
- integer(i4b),intent(in) :: nSoil ! number of soil layers
- integer(i4b),intent(in) :: nLayers ! total number of layers
- logical(lgt),intent(in) :: getSatDepth ! logical flag to compute index of the lowest saturated layer
- ! input: state and diagnostic variables
- real(dp),intent(in) :: mLayerdTheta_dPsi(:) ! derivative in the soil water characteristic w.r.t. matric head in each layer (m-1)
- real(dp),intent(in) :: mLayerMatricHeadLiq(:) ! matric head in each layer at the current iteration (m)
- real(dp),intent(in) :: mLayerVolFracLiq(:) ! volumetric fraction of liquid water (-)
- real(dp),intent(in) :: mLayerVolFracIce(:) ! volumetric fraction of ice (-)
- ! input/output: data structures
- type(var_d),intent(in) :: attr_data ! spatial attributes
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- 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(inout) :: flux_data ! model fluxes for a local HRU
- ! output: baseflow
- integer(i4b),intent(inout) :: ixSaturation ! index of lowest saturated layer (NOTE: only computed on the first iteration)
- real(dp),intent(out) :: mLayerBaseflow(:) ! baseflow from each soil layer (m s-1)
- real(dp),intent(out) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1)
- ! output: error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! ---------------------------------------------------------------------------------------
- ! * local variables
- ! ---------------------------------------------------------------------------------------
- ! general local variables
- integer(i4b) :: iLayer ! index of soil layer
- real(dp),dimension(nSoil,nSoil) :: dBaseflow_dVolLiq ! derivative in the baseflow flux w.r.t. volumetric liquid water content (m s-1)
- ! local variables to compute the numerical Jacobian
- logical(lgt),parameter :: doNumericalJacobian=.false. ! flag to compute the numerical Jacobian
- real(dp),dimension(nSoil) :: mLayerMatricHeadPerturbed ! perturbed matric head (m)
- real(dp),dimension(nSoil) :: mLayerVolFracLiqPerturbed ! perturbed volumetric fraction of liquid water (-)
- real(dp),dimension(nSoil) :: mLayerBaseflowPerturbed ! perturbed baseflow (m s-1)
- real(dp),dimension(nSoil,nSoil) :: nJac ! numerical Jacobian (s-1)
- ! ***************************************************************************************
- ! ***************************************************************************************
- ! initialize error control
- err=0; message='groundwatr/'
- ! ---------------------------------------------------------------------------------------
- ! ---------------------------------------------------------------------------------------
- ! associate variables in data structures
- associate(&
-
- ! input: baseflow parameters
- fieldCapacity => mpar_data%var(iLookPARAM%fieldCapacity)%dat(1), & ! intent(in): [dp] field capacity (-)
- 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 (-)
-
- ! input: van Genuchten soil parametrers
- vGn_alpha => mpar_data%var(iLookPARAM%vGn_alpha)%dat, & ! intent(in): [dp] van Genutchen "alpha" parameter (m-1)
- vGn_n => mpar_data%var(iLookPARAM%vGn_n)%dat, & ! intent(in): [dp] van Genutchen "n" parameter (-)
- vGn_m => diag_data%var(iLookDIAG%scalarVGn_m)%dat, & ! intent(in): [dp] van Genutchen "m" parameter (-)
-
- ! output: diagnostic variables
- scalarExfiltration => flux_data%var(iLookFLUX%scalarExfiltration)%dat(1), & ! intent(out):[dp] exfiltration from the soil profile (m s-1)
- mLayerColumnOutflow => flux_data%var(iLookFLUX%mLayerColumnOutflow)%dat & ! intent(out):[dp(:)] column outflow from each soil layer (m3 s-1)
-
- ) ! end association to variables in data structures
-
- ! ************************************************************************************************
- ! (1) compute the "active" portion of the soil profile
- ! ************************************************************************************************
-
- ! get index of the lowest saturated layer
- if(getSatDepth)then ! NOTE: only compute for the first flux call
- ixSaturation = nSoil+1 ! unsaturated profile when ixSaturation>nSoil
- do iLayer=nSoil,1,-1 ! start at the lowest soil layer and work upwards to the top layer
- if(mLayerVolFracLiq(iLayer) > fieldCapacity)then; ixSaturation = iLayer ! index of saturated layer -- keeps getting over-written as move upwards
- else; exit; end if ! (only consider saturated layer at the bottom of the soil profile)
- end do ! (looping through soil layers)
- end if
-
- ! check for an early return (no layers are "active")
- if(ixSaturation > nSoil)then
- scalarExfiltration = 0._dp ! exfiltration from the soil profile (m s-1)
- mLayerColumnOutflow(:) = 0._dp ! column outflow from each soil layer (m3 s-1)
- mLayerBaseflow(:) = 0._dp ! baseflow from each soil layer (m s-1)
- dBaseflow_dMatric(:,:) = 0._dp ! derivative in baseflow w.r.t. matric head (s-1)
- return
- end if ! if some layers are saturated
-
- ! ************************************************************************************************
- ! (2) compute the baseflow flux and its derivative w.r.t volumetric liquid water content
- ! ************************************************************************************************
-
- ! use private subroutine to compute baseflow (for multiple calls for numerical Jacobian)
- call computeBaseflow(&
- ! input: control and state variables
- nSnow, & ! intent(in): number of snow layers
- nSoil, & ! intent(in): number of soil layers
- nLayers, & ! intent(in): total number of layers
- .true., & ! intent(in): .true. if analytical derivatives are desired
- ixSaturation, & ! intent(in): index of upper-most "saturated" layer
- mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water in each soil layer (-)
- mLayerVolFracIce, & ! intent(in): volumetric fraction of ice in each soil layer (-)
- ! input/output: data structures
- attr_data, & ! intent(in): spatial attributes
- mpar_data, & ! intent(in): model parameters
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- flux_data, & ! intent(inout): model fluxes for a local HRU
- ! output: fluxes and derivatives
- mLayerBaseflow, & ! intent(out): baseflow flux in each soil layer (m s-1)
- dBaseflow_dVolLiq) ! intent(out): derivative in baseflow w.r.t. volumetric liquid water content (s-1)
-
- ! use the chain rule to compute the baseflow derivative w.r.t. matric head (s-1)
- do iLayer=1,nSoil
- dBaseflow_dMatric(1:iLayer,iLayer) = dBaseflow_dVolLiq(1:iLayer,iLayer)*mLayerdTheta_dPsi(iLayer)
- if(iLayer<nSoil) dBaseflow_dMatric(iLayer+1:nSoil,iLayer) = 0._dp
- end do
-
- ! ************************************************************************************************
- ! (4) compute the numerical Jacobian
- ! ************************************************************************************************
- if(doNumericalJacobian)then
-
- ! first, print the analytical Jacobian
- print*, '** analytical Jacobian:'
- write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=1,nSoil)
- do iLayer=1,nSoil; write(*,'(i4,1x,100(e12.5,1x))') iLayer, dBaseflow_dMatric(1:nSoil,iLayer); end do
-
- ! get a copy of the state vector to perturb
- mLayerMatricHeadPerturbed(:) = mLayerMatricHeadLiq(:)
- mLayerVolFracLiqPerturbed(:) = mLayerVolFracLiq(:)
-
- ! loop through the state variables
- do iLayer=1,nSoil
-
- ! perturb state vector
- mLayerMatricHeadPerturbed(iLayer) = mLayerMatricHeadPerturbed(iLayer) + dx
-
- ! compute the columetruc liquid water content
- mLayerVolFracLiqPerturbed(iLayer) = volFracLiq(mLayerMatricHeadPerturbed(iLayer),vGn_alpha(iLayer),theta_res(iLayer),theta_sat(iLayer),vGn_n(iLayer),vGn_m(iLayer))
-
- ! compute baseflow flux
- ! NOTE: This is an optional second call to computeBaseflow that is invoked when computing numerical derivatives.
- ! Since the purpose here is to compute the numerical derivatives, we do not need to compute analytical derivatives also.
- ! Hence, analytical derivatives are not desired
- call computeBaseflow(&
- ! input: control and state variables
- nSnow, & ! intent(in): number of snow layers
- nSoil, & ! intent(in): number of soil layers
- nLayers, & ! intent(in): total number of layers
- .false., & ! intent(in): .true. if analytical derivatives are desired
- ixSaturation, & ! intent(in): index of upper-most "saturated" layer
- mLayerVolFracLiqPerturbed, & ! intent(in): volumetric fraction of liquid water in each soil layer (-)
- mLayerVolFracIce, & ! intent(in): volumetric fraction of ice in each soil layer (-)
- ! input/output: data structures
- attr_data, & ! intent(in): spatial attributes
- mpar_data, & ! intent(in): model parameters
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- flux_data, & ! intent(inout): model fluxes for a local HRU
- ! output: fluxes and derivatives
- mLayerBaseflowPerturbed, & ! intent(out): baseflow flux in each soil layer (m s-1)
- dBaseflow_dVolLiq) ! intent(out): ** NOT USED ** derivative in baseflow w.r.t. volumetric liquid water content (s-1)
-
- ! compute the numerical Jacobian
- nJac(:,iLayer) = (mLayerBaseflowPerturbed(:) - mLayerBaseflow(:))/dx ! compute the Jacobian
-
- ! set the state back to the input value
- mLayerMatricHeadPerturbed(iLayer) = mLayerMatricHeadLiq(iLayer)
- mLayerVolFracLiqPerturbed(iLayer) = mLayerVolFracLiq(iLayer)
-
- end do ! looping through state variables
-
- ! print the numerical Jacobian
- print*, '** numerical Jacobian:'
- write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=1,nSoil)
- do iLayer=1,nSoil; write(*,'(i4,1x,100(e12.5,1x))') iLayer, nJac(1:nSoil,iLayer); end do
- !pause 'testing Jacobian'
-
- end if ! if desire to compute the Jacobian
-
- ! end association to variables in data structures
- end associate
-
- end subroutine groundwatr
-
-
- ! ***********************************************************************************************************************
- ! * private subroutine computeBaseflow: private subroutine so can be used to test the numerical jacobian
- ! ***********************************************************************************************************************
- subroutine computeBaseflow(&
- ! input: control and state variables
- nSnow, & ! intent(in): number of snow layers
- nSoil, & ! intent(in): number of soil layers
- nLayers, & ! intent(in): total number of layers
- derivDesired, & ! intent(in): .true. if derivatives are desired
- ixSaturation, & ! intent(in): index of upper-most "saturated" layer
- mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water in each soil layer (-)
- mLayerVolFracIce, & ! intent(in): volumetric fraction of ice in each soil layer (-)
- ! input/output: data structures
- attr_data, & ! intent(in): spatial attributes
- mpar_data, & ! intent(in): model parameters
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- flux_data, & ! intent(inout): model fluxes for a local HRU
- ! output: fluxes and derivatives
- mLayerBaseflow, & ! intent(out): baseflow flux in each soil layer (m s-1)
- dBaseflow_dVolLiq) ! intent(out): derivative in baseflow w.r.t. volumetric liquid water content (s-1)
- implicit none
- ! ---------------------------------------------------------------------------------------
- ! * dummy variables
- ! ---------------------------------------------------------------------------------------
- ! input: control and state variables
- integer(i4b),intent(in) :: nSnow ! number of snow layers
- integer(i4b),intent(in) :: nSoil ! number of soil layers
- integer(i4b),intent(in) :: nLayers ! total number of layers
- logical(lgt),intent(in) :: derivDesired ! .true. if derivatives are desired
- integer(i4b),intent(in) :: ixSaturation ! index of upper-most "saturated" layer
- real(dp),intent(in) :: mLayerVolFracLiq(:) ! volumetric fraction of liquid water (-)
- real(dp),intent(in) :: mLayerVolFracIce(:) ! volumetric fraction of ice (-)
- ! input/output: data structures
- type(var_d),intent(in) :: attr_data ! spatial attributes
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- type(var_dlength),intent(in) :: prog_data ! prognostic variables for a local HRU
- type(var_dlength),intent(inout) :: flux_data ! model fluxes for a local HRU
- ! output: baseflow
- real(dp),intent(out) :: mLayerBaseflow(:) ! baseflow from each soil layer (m s-1)
- real(dp),intent(out) :: dBaseflow_dVolLiq(:,:) ! derivative in baseflow w.r.t. matric head (s-1)
- ! ---------------------------------------------------------------------------------------
- ! * local variables
- ! ---------------------------------------------------------------------------------------
- ! general local variables
- integer(i4b) :: iLayer,jLayer ! index of model layer
- ! local variables for the exfiltration
- real(dp) :: totalColumnInflow ! total column inflow (m s-1)
- real(dp) :: totalColumnOutflow ! total column outflow (m s-1)
- real(dp) :: availStorage ! available storage (m)
- real(dp),parameter :: xMinEval=0.002_dp ! minimum value to evaluate the exfiltration function (m)
- real(dp),parameter :: xCenter=0.001_dp ! center of the exfiltration function (m)
- real(dp),parameter :: xWidth=0.0001_dp ! width of the exfiltration function (m)
- real(dp) :: expF,logF ! logistic smoothing function (-)
- ! local variables for the lateral flux among soil columns
- real(dp) :: activePorosity ! "active" porosity associated with storage above a threshold (-)
- real(dp) :: drainableWater ! drainable water in eaxch layer (m)
- real(dp) :: tran0 ! maximum transmissivity (m2 s-1)
- real(dp),dimension(nSoil) :: zActive ! water table thickness associated with storage below and including the given layer (m)
- real(dp),dimension(nSoil) :: trTotal ! total transmissivity associated with total water table depth zActive (m2 s-1)
- real(dp),dimension(nSoil) :: trSoil ! transmissivity of water in a given layer (m2 s-1)
- ! local variables for the derivatives
- real(dp) :: qbTotal ! total baseflow (m s-1)
- real(dp) :: length2area ! ratio of hillslope width to hillslope area (m m-2)
- real(dp),dimension(nSoil) :: depth2capacity ! ratio of layer depth to total subsurface storage capacity (-)
- real(dp),dimension(nSoil) :: dXdS ! change in dimensionless flux w.r.t. change in dimensionless storage (-)
- real(dp),dimension(nSoil) :: dLogFunc_dLiq ! derivative in the logistic function w.r.t. volumetric liquid water content (-)
- real(dp),dimension(nSoil) :: dExfiltrate_dVolLiq ! derivative in exfiltration w.r.t. volumetric liquid water content (-)
- ! local variables for testing (debugging)
- logical(lgt),parameter :: printFlag=.false. ! flag for printing (debugging)
- real(dp) :: xDepth,xTran,xFlow ! temporary variables (depth, transmissivity, flow)
- ! ---------------------------------------------------------------------------------------
- ! * association to data in structures
- ! ---------------------------------------------------------------------------------------
- associate(&
-
- ! input: coordinate variables
- soilDepth => prog_data%var(iLookPROG%iLayerHeight)%dat(nLayers), & ! intent(in): [dp] total soil depth (m)
- mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat(nSnow+1:nLayers),& ! intent(in): [dp(:)] depth of each soil layer (m)
-
- ! input: diagnostic variables
- surfaceHydCond => flux_data%var(iLookFLUX%mLayerSatHydCondMP)%dat(1), & ! intent(in): [dp] saturated hydraulic conductivity at the surface (m s-1)
- mLayerColumnInflow => flux_data%var(iLookFLUX%mLayerColumnInflow)%dat, & ! intent(in): [dp(:)] inflow into each soil layer (m3/s)
-
- ! input: local attributes
- HRUarea => attr_data%var(iLookATTR%HRUarea), & ! intent(in): [dp] HRU area (m2)
- tan_slope => attr_data%var(iLookATTR%tan_slope), & ! intent(in): [dp] tan water table slope, taken as tan local ground surface slope (-)
- contourLength => attr_data%var(iLookATTR%contourLength), & ! intent(in): [dp] length of contour at downslope edge of HRU (m)
-
- ! input: baseflow parameters
- zScale_TOPMODEL => mpar_data%var(iLookPARAM%zScale_TOPMODEL)%dat(1), & ! intent(in): [dp] TOPMODEL exponent (-)
- kAnisotropic => mpar_data%var(iLookPARAM%kAnisotropic)%dat(1), & ! intent(in): [dp] anisotropy factor for lateral hydraulic conductivity (-
- fieldCapacity => mpar_data%var(iLookPARAM%fieldCapacity)%dat(1), & ! intent(in): [dp] field capacity (-)
- theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat, & ! intent(in): [dp(:)] soil porosity (-)
-
- ! output: diagnostic variables
- scalarExfiltration => flux_data%var(iLookFLUX%scalarExfiltration)%dat(1), & ! intent(out):[dp] exfiltration from the soil profile (m s-1)
- mLayerColumnOutflow => flux_data%var(iLookFLUX%mLayerColumnOutflow)%dat & ! intent(out):[dp(:)] column outflow from each soil layer (m3 s-1)
-
- ) ! end association to variables in data structures
- ! ***********************************************************************************************************************
- ! ***********************************************************************************************************************
- ! start routine here
-
- ! ***********************************************************************************************************************
- ! (1) compute the baseflow flux in each soil layer
- ! ***********************************************************************************************************************
-
- ! compute the maximum transmissivity
- ! NOTE: this can be done as a pre-processing step
- tran0 = kAnisotropic*surfaceHydCond*soilDepth/zScale_TOPMODEL ! maximum transmissivity (m2 s-1)
-
- ! compute the water table thickness (m) and transmissivity in each layer (m2 s-1)
- do iLayer=nSoil,ixSaturation,-1 ! loop through "active" soil layers, from lowest to highest
- ! define drainable water in each layer (m)
- activePorosity = theta_sat(iLayer) - fieldCapacity ! "active" porosity (-)
- drainableWater = mLayerDepth(iLayer)*(max(0._dp,mLayerVolFracLiq(iLayer) - fieldCapacity))/activePorosity
- ! compute layer transmissivity
- if(iLayer==nSoil)then
- zActive(iLayer) = drainableWater ! water table thickness associated with storage in a given layer (m)
- trTotal(iLayer) = tran0*(zActive(iLayer)/soilDepth)**zScale_TOPMODEL ! total transmissivity for total depth zActive (m2 s-1)
- trSoil(iLayer) = trTotal(iLayer) ! transmissivity of water in a given layer (m2 s-1)
- else
- zActive(iLayer) = zActive(iLayer+1) + drainableWater
- trTotal(iLayer) = tran0*(zActive(iLayer)/soilDepth)**zScale_TOPMODEL
- trSoil(iLayer) = trTotal(iLayer) - trTotal(iLayer+1)
- end if
- !write(*,'(a,1x,i4,1x,10(f20.15,1x))') 'iLayer, mLayerMatricHeadLiq(iLayer), mLayerVolFracLiq(iLayer), zActive(iLayer), trTotal(iLayer), trSoil(iLayer) = ', &
- ! iLayer, mLayerMatricHeadLiq(iLayer), mLayerVolFracLiq(iLayer), zActive(iLayer), trTotal(iLayer), trSoil(iLayer)
- end do ! looping through soil layers
-
- ! set un-used portions of the vectors to zero
- if(ixSaturation>1)then
- zActive(1:ixSaturation-1) = 0._dp
- trTotal(1:ixSaturation-1) = 0._dp
- trSoil(1:ixSaturation-1) = 0._dp
- end if
-
- ! compute the outflow from each layer (m3 s-1)
- mLayerColumnOutflow(1:nSoil) = trSoil(1:nSoil)*tan_slope*contourLength
-
- ! compute total column inflow and total column outflow (m s-1)
- totalColumnInflow = sum(mLayerColumnInflow(1:nSoil))/HRUarea
- totalColumnOutflow = sum(mLayerColumnOutflow(1:nSoil))/HRUarea
-
- ! compute the available storage (m)
- availStorage = sum(mLayerDepth(1:nSoil)*(theta_sat - (mLayerVolFracLiq(1:nSoil)+mLayerVolFracIce(1:nSoil))) )
-
- ! compute the smoothing function (-)
- if(availStorage < xMinEval)then
- ! (compute the logistic function)
- expF = exp((availStorage - xCenter)/xWidth)
- logF = 1._dp / (1._dp + expF)
- ! (compute the derivative in the logistic function w.r.t. volumetric liquid water content in each soil layer)
- dLogFunc_dLiq(1:nSoil) = mLayerDepth(1:nSoil)*(expF/xWidth)/(1._dp + expF)**2._dp
- else
- logF = 0._dp
- dLogFunc_dLiq(:) = 0._dp
- end if
-
- ! compute the exfiltartion (m s-1)
- if(totalColumnInflow > totalColumnOutflow .and. logF > tiny(1._dp))then
- scalarExfiltration = logF*(totalColumnInflow - totalColumnOutflow) ! m s-1
- !write(*,'(a,1x,10(f30.20,1x))') 'scalarExfiltration = ', scalarExfiltration
- else
- scalarExfiltration = 0._dp
- end if
-
- ! check
- !write(*,'(a,1x,10(f30.20,1x))') 'zActive(1), soilDepth, availStorage, logF, scalarExfiltration = ', &
- ! zActive(1), soilDepth, availStorage, logF, scalarExfiltration
- !if(scalarExfiltration > tiny(1.0_dp)) pause 'exfiltrating'
-
- ! compute the baseflow in each layer (m s-1)
- mLayerBaseflow(1:nSoil) = (mLayerColumnOutflow(1:nSoil) - mLayerColumnInflow(1:nSoil))/HRUarea
-
- ! compute the total baseflow
- qbTotal = sum(mLayerBaseflow)
-
- ! add exfiltration to the baseflow flux at the top layer
- mLayerBaseflow(1) = mLayerBaseflow(1) + scalarExfiltration
- mLayerColumnOutflow(1) = mLayerColumnOutflow(1) + scalarExfiltration*HRUarea
-
- ! test
- if(printFlag)then
- xDepth = sum(mLayerDepth(ixSaturation:nSoil)*(mLayerVolFracLiq(ixSaturation:nSoil) - fieldCapacity))/sum(theta_sat(ixSaturation:nSoil) - fieldCapacity) ! "effective" water table thickness (m)
- xTran = tran0*(xDepth/soilDepth)**zScale_TOPMODEL ! transmissivity for the entire aquifer (m2 s-1)
- xFlow = xTran*tan_slope*contourLength/HRUarea ! total column outflow (m s-1)
- print*, 'ixSaturation = ', ixSaturation
- write(*,'(a,1x,5(f30.20,1x))') 'tran0, soilDepth = ', tran0, soilDepth
- write(*,'(a,1x,5(f30.20,1x))') 'surfaceHydCond, zScale_TOPMODEL = ', surfaceHydCond, zScale_TOPMODEL
- write(*,'(a,1x,5(f30.20,1x))') 'xDepth, zActive(ixSaturation) = ', xDepth, zActive(ixSaturation)
- write(*,'(a,1x,5(f30.20,1x))') 'xTran, trTotal(ixSaturation) = ', xTran, trTotal(ixSaturation)
- write(*,'(a,1x,5(f30.20,1x))') 'xFlow, totalColumnOutflow = ', xFlow, sum(mLayerColumnOutflow(:))/HRUarea
- !pause 'check groundwater'
- end if
-
- ! ***********************************************************************************************************************
- ! (2) compute the derivative in the baseflow flux w.r.t. volumetric liquid water content (m s-1)
- ! ***********************************************************************************************************************
-
- ! initialize the derivative matrix
- dBaseflow_dVolLiq(:,:) = 0._dp
-
- ! check if derivatives are actually required
- if(.not.derivDesired) return
-
- ! compute ratio of hillslope width to hillslope area (m m-2)
- length2area = tan_slope*contourLength/HRUarea
-
- ! compute the ratio of layer depth to maximum water holding capacity (-)
- depth2capacity(1:nSoil) = mLayerDepth(1:nSoil)/sum( (theta_sat(1:nSoil) - fieldCapacity)*mLayerDepth(1:nSoil) )
-
- ! compute the change in dimensionless flux w.r.t. change in dimensionless storage (-)
- dXdS(1:nSoil) = zScale_TOPMODEL*(zActive(1:nSoil)/SoilDepth)**(zScale_TOPMODEL - 1._dp)
-
- ! loop through soil layers
- do iLayer=1,nSoil
- ! compute diagonal terms (s-1)
- dBaseflow_dVolLiq(iLayer,iLayer) = tran0*dXdS(iLayer)*depth2capacity(iLayer)*length2area
- ! compute off-diagonal terms
- do jLayer=iLayer+1,nSoil ! (only dependent on layers below)
- dBaseflow_dVolLiq(iLayer,jLayer) = tran0*(dXdS(iLayer) - dXdS(iLayer+1))*depth2capacity(jLayer)*length2area
- end do ! looping through soil layers
- end do ! looping through soil layers
-
- ! compute the derivative in the exfiltration flux w.r.t. volumetric liquid water content (m s-1)
- if(qbTotal < 0._dp)then
- do iLayer=1,nSoil
- dExfiltrate_dVolLiq(iLayer) = dBaseflow_dVolLiq(iLayer,iLayer)*logF + dLogFunc_dLiq(iLayer)*qbTotal
- end do ! looping through soil layers
- dBaseflow_dVolLiq(1,1:nSoil) = dBaseflow_dVolLiq(1,1:nSoil) - dExfiltrate_dVolLiq(1:nSoil)
- end if
-
- ! end association to data in structures
- end associate
-
- end subroutine computeBaseflow
-
+! ************************************************************************************************
+! public subroutine groundwatr: compute the groundwater sink term in Richards' equation
+! ************************************************************************************************
+!
+! Method
+! ------
+!
+! Here we assume that water avaialble for shallow groundwater flow includes is all water above
+! "field capacity" below the depth zCrit, where zCrit is defined as the lowest point in the soil
+! profile where the volumetric liquid water content is less than field capacity.
+!
+! We further assume that transmssivity (m2 s-1) for each layer is defined asuming that the water
+! available for saturated flow is located at the bottom of the soil profile. Specifically:
+! trTotal(iLayer) = tran0*(zActive(iLayer)/soilDepth)**zScale_TOPMODEL
+! trSoil(iLayer) = trTotal(iLayer) - trTotal(iLayer+1)
+! where zActive(iLayer) is the effective water table thickness for all layers up to and including
+! the current layer (working from the bottom to the top).
+!
+! The outflow from each layer is then (m3 s-1)
+! mLayerOutflow(iLayer) = trSoil(iLayer)*tan_slope*contourLength
+! where contourLength is the width of a hillslope (m) parallel to a stream
+!
+! ************************************************************************************************
+subroutine groundwatr(&
+
+ ! input: model control
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ getSatDepth, & ! intent(in): logical flag to compute index of the lowest saturated layer
+
+ ! input: state and diagnostic variables
+ mLayerdTheta_dPsi, & ! intent(in): derivative in the soil water characteristic w.r.t. matric head in each layer (m-1)
+ mLayerMatricHeadLiq, & ! intent(in): liquid water matric potential (m)
+ mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water (-)
+ mLayerVolFracIce, & ! intent(in): volumetric fraction of ice (-)
+
+ ! input/output: data structures
+ attr_data, & ! intent(in): spatial attributes
+ mpar_data, & ! intent(in): model parameters
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ diag_data, & ! intent(in): model diagnostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+
+ ! output: baseflow
+ ixSaturation, & ! intent(inout) index of lowest saturated layer (NOTE: only computed on the first iteration)
+ mLayerBaseflow, & ! intent(out): baseflow from each soil layer (m s-1)
+ dBaseflow_dMatric, & ! intent(out): derivative in baseflow w.r.t. matric head (s-1)
+
+ ! output: error control
+ err,message) ! intent(out): error control
+ ! ---------------------------------------------------------------------------------------
+ ! utility modules
+ USE soil_utils_module,only:volFracLiq ! compute volumetric fraction of liquid water as a function of matric head
+ USE soil_utils_module,only:hydCond_psi ! compute hydraulic conductivity as a function of matric head
+ implicit none
+ ! ---------------------------------------------------------------------------------------
+ ! * dummy variables
+ ! ---------------------------------------------------------------------------------------
+ ! input: model control
+ integer(i4b),intent(in) :: nSnow ! number of snow layers
+ integer(i4b),intent(in) :: nSoil ! number of soil layers
+ integer(i4b),intent(in) :: nLayers ! total number of layers
+ logical(lgt),intent(in) :: getSatDepth ! logical flag to compute index of the lowest saturated layer
+ ! input: state and diagnostic variables
+ real(rkind),intent(in) :: mLayerdTheta_dPsi(:) ! derivative in the soil water characteristic w.r.t. matric head in each layer (m-1)
+ real(rkind),intent(in) :: mLayerMatricHeadLiq(:) ! matric head in each layer at the current iteration (m)
+ real(rkind),intent(in) :: mLayerVolFracLiq(:) ! volumetric fraction of liquid water (-)
+ real(rkind),intent(in) :: mLayerVolFracIce(:) ! volumetric fraction of ice (-)
+ ! input/output: data structures
+ type(var_d),intent(in) :: attr_data ! spatial attributes
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ 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(inout) :: flux_data ! model fluxes for a local HRU
+ ! output: baseflow
+ integer(i4b),intent(inout) :: ixSaturation ! index of lowest saturated layer (NOTE: only computed on the first iteration)
+ real(rkind),intent(out) :: mLayerBaseflow(:) ! baseflow from each soil layer (m s-1)
+ real(rkind),intent(out) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1)
+ ! output: error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! ---------------------------------------------------------------------------------------
+ ! * local variables
+ ! ---------------------------------------------------------------------------------------
+ ! general local variables
+ integer(i4b) :: iLayer ! index of soil layer
+ real(rkind),dimension(nSoil,nSoil) :: dBaseflow_dVolLiq ! derivative in the baseflow flux w.r.t. volumetric liquid water content (m s-1)
+ ! local variables to compute the numerical Jacobian
+ logical(lgt),parameter :: doNumericalJacobian=.false. ! flag to compute the numerical Jacobian
+ real(rkind),dimension(nSoil) :: mLayerMatricHeadPerturbed ! perturbed matric head (m)
+ real(rkind),dimension(nSoil) :: mLayerVolFracLiqPerturbed ! perturbed volumetric fraction of liquid water (-)
+ real(rkind),dimension(nSoil) :: mLayerBaseflowPerturbed ! perturbed baseflow (m s-1)
+ real(rkind),dimension(nSoil,nSoil) :: nJac ! numerical Jacobian (s-1)
+ ! ***************************************************************************************
+ ! ***************************************************************************************
+ ! initialize error control
+ err=0; message='groundwatr/'
+ ! ---------------------------------------------------------------------------------------
+ ! ---------------------------------------------------------------------------------------
+ ! associate variables in data structures
+ associate(&
+
+ ! input: baseflow parameters
+ fieldCapacity => mpar_data%var(iLookPARAM%fieldCapacity)%dat(1), & ! intent(in): [dp] field capacity (-)
+ 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 (-)
+
+ ! input: van Genuchten soil parametrers
+ vGn_alpha => mpar_data%var(iLookPARAM%vGn_alpha)%dat, & ! intent(in): [dp] van Genutchen "alpha" parameter (m-1)
+ vGn_n => mpar_data%var(iLookPARAM%vGn_n)%dat, & ! intent(in): [dp] van Genutchen "n" parameter (-)
+ vGn_m => diag_data%var(iLookDIAG%scalarVGn_m)%dat, & ! intent(in): [dp] van Genutchen "m" parameter (-)
+
+ ! output: diagnostic variables
+ scalarExfiltration => flux_data%var(iLookFLUX%scalarExfiltration)%dat(1), & ! intent(out):[dp] exfiltration from the soil profile (m s-1)
+ mLayerColumnOutflow => flux_data%var(iLookFLUX%mLayerColumnOutflow)%dat & ! intent(out):[dp(:)] column outflow from each soil layer (m3 s-1)
+
+ ) ! end association to variables in data structures
+
+ ! ************************************************************************************************
+ ! (1) compute the "active" portion of the soil profile
+ ! ************************************************************************************************
+
+ ! get index of the lowest saturated layer
+ if(getSatDepth)then ! NOTE: only compute for the first flux call
+ ixSaturation = nSoil+1 ! unsaturated profile when ixSaturation>nSoil
+ do iLayer=nSoil,1,-1 ! start at the lowest soil layer and work upwards to the top layer
+ if(mLayerVolFracLiq(iLayer) > fieldCapacity)then; ixSaturation = iLayer ! index of saturated layer -- keeps getting over-written as move upwards
+ else; exit; end if ! (only consider saturated layer at the bottom of the soil profile)
+ end do ! (looping through soil layers)
+ end if
+
+ ! check for an early return (no layers are "active")
+ if(ixSaturation > nSoil)then
+ scalarExfiltration = 0._rkind ! exfiltration from the soil profile (m s-1)
+ mLayerColumnOutflow(:) = 0._rkind ! column outflow from each soil layer (m3 s-1)
+ mLayerBaseflow(:) = 0._rkind ! baseflow from each soil layer (m s-1)
+ dBaseflow_dMatric(:,:) = 0._rkind ! derivative in baseflow w.r.t. matric head (s-1)
+ return
+ end if ! if some layers are saturated
+
+ ! ************************************************************************************************
+ ! (2) compute the baseflow flux and its derivative w.r.t volumetric liquid water content
+ ! ************************************************************************************************
+
+ ! use private subroutine to compute baseflow (for multiple calls for numerical Jacobian)
+ call computeBaseflow(&
+ ! input: control and state variables
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ .true., & ! intent(in): .true. if analytical derivatives are desired
+ ixSaturation, & ! intent(in): index of upper-most "saturated" layer
+ mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water in each soil layer (-)
+ mLayerVolFracIce, & ! intent(in): volumetric fraction of ice in each soil layer (-)
+ ! input/output: data structures
+ attr_data, & ! intent(in): spatial attributes
+ mpar_data, & ! intent(in): model parameters
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+ ! output: fluxes and derivatives
+ mLayerBaseflow, & ! intent(out): baseflow flux in each soil layer (m s-1)
+ dBaseflow_dVolLiq) ! intent(out): derivative in baseflow w.r.t. volumetric liquid water content (s-1)
+
+ ! use the chain rule to compute the baseflow derivative w.r.t. matric head (s-1)
+ do iLayer=1,nSoil
+ dBaseflow_dMatric(1:iLayer,iLayer) = dBaseflow_dVolLiq(1:iLayer,iLayer)*mLayerdTheta_dPsi(iLayer)
+ if(iLayer<nSoil) dBaseflow_dMatric(iLayer+1:nSoil,iLayer) = 0._rkind
+ end do
+
+ ! ************************************************************************************************
+ ! (4) compute the numerical Jacobian
+ ! ************************************************************************************************
+ if(doNumericalJacobian)then
+
+ ! first, print the analytical Jacobian
+ print*, '** analytical Jacobian:'
+ write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=1,nSoil)
+ do iLayer=1,nSoil; write(*,'(i4,1x,100(e12.5,1x))') iLayer, dBaseflow_dMatric(1:nSoil,iLayer); end do
+
+ ! get a copy of the state vector to perturb
+ mLayerMatricHeadPerturbed(:) = mLayerMatricHeadLiq(:)
+ mLayerVolFracLiqPerturbed(:) = mLayerVolFracLiq(:)
+
+ ! loop through the state variables
+ do iLayer=1,nSoil
+
+ ! perturb state vector
+ mLayerMatricHeadPerturbed(iLayer) = mLayerMatricHeadPerturbed(iLayer) + dx
+
+ ! compute the columetruc liquid water content
+ mLayerVolFracLiqPerturbed(iLayer) = volFracLiq(mLayerMatricHeadPerturbed(iLayer),vGn_alpha(iLayer),theta_res(iLayer),theta_sat(iLayer),vGn_n(iLayer),vGn_m(iLayer))
+
+ ! compute baseflow flux
+ ! NOTE: This is an optional second call to computeBaseflow that is invoked when computing numerical derivatives.
+ ! Since the purpose here is to compute the numerical derivatives, we do not need to compute analytical derivatives also.
+ ! Hence, analytical derivatives are not desired
+ call computeBaseflow(&
+ ! input: control and state variables
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ .false., & ! intent(in): .true. if analytical derivatives are desired
+ ixSaturation, & ! intent(in): index of upper-most "saturated" layer
+ mLayerVolFracLiqPerturbed, & ! intent(in): volumetric fraction of liquid water in each soil layer (-)
+ mLayerVolFracIce, & ! intent(in): volumetric fraction of ice in each soil layer (-)
+ ! input/output: data structures
+ attr_data, & ! intent(in): spatial attributes
+ mpar_data, & ! intent(in): model parameters
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+ ! output: fluxes and derivatives
+ mLayerBaseflowPerturbed, & ! intent(out): baseflow flux in each soil layer (m s-1)
+ dBaseflow_dVolLiq) ! intent(out): ** NOT USED ** derivative in baseflow w.r.t. volumetric liquid water content (s-1)
+
+ ! compute the numerical Jacobian
+ nJac(:,iLayer) = (mLayerBaseflowPerturbed(:) - mLayerBaseflow(:))/dx ! compute the Jacobian
+
+ ! set the state back to the input value
+ mLayerMatricHeadPerturbed(iLayer) = mLayerMatricHeadLiq(iLayer)
+ mLayerVolFracLiqPerturbed(iLayer) = mLayerVolFracLiq(iLayer)
+
+ end do ! looping through state variables
+
+ ! print the numerical Jacobian
+ print*, '** numerical Jacobian:'
+ write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=1,nSoil)
+ do iLayer=1,nSoil; write(*,'(i4,1x,100(e12.5,1x))') iLayer, nJac(1:nSoil,iLayer); end do
+ !pause 'testing Jacobian'
+
+ end if ! if desire to compute the Jacobian
+
+ ! end association to variables in data structures
+ end associate
+
+end subroutine groundwatr
+
+
+! ***********************************************************************************************************************
+! * private subroutine computeBaseflow: private subroutine so can be used to test the numerical jacobian
+! ***********************************************************************************************************************
+subroutine computeBaseflow(&
+ ! input: control and state variables
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ derivDesired, & ! intent(in): .true. if derivatives are desired
+ ixSaturation, & ! intent(in): index of upper-most "saturated" layer
+ mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water in each soil layer (-)
+ mLayerVolFracIce, & ! intent(in): volumetric fraction of ice in each soil layer (-)
+ ! input/output: data structures
+ attr_data, & ! intent(in): spatial attributes
+ mpar_data, & ! intent(in): model parameters
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+ ! output: fluxes and derivatives
+ mLayerBaseflow, & ! intent(out): baseflow flux in each soil layer (m s-1)
+ dBaseflow_dVolLiq) ! intent(out): derivative in baseflow w.r.t. volumetric liquid water content (s-1)
+ implicit none
+ ! ---------------------------------------------------------------------------------------
+ ! * dummy variables
+ ! ---------------------------------------------------------------------------------------
+ ! input: control and state variables
+ integer(i4b),intent(in) :: nSnow ! number of snow layers
+ integer(i4b),intent(in) :: nSoil ! number of soil layers
+ integer(i4b),intent(in) :: nLayers ! total number of layers
+ logical(lgt),intent(in) :: derivDesired ! .true. if derivatives are desired
+ integer(i4b),intent(in) :: ixSaturation ! index of upper-most "saturated" layer
+ real(rkind),intent(in) :: mLayerVolFracLiq(:) ! volumetric fraction of liquid water (-)
+ real(rkind),intent(in) :: mLayerVolFracIce(:) ! volumetric fraction of ice (-)
+ ! input/output: data structures
+ type(var_d),intent(in) :: attr_data ! spatial attributes
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ type(var_dlength),intent(in) :: prog_data ! prognostic variables for a local HRU
+ type(var_dlength),intent(inout) :: flux_data ! model fluxes for a local HRU
+ ! output: baseflow
+ real(rkind),intent(out) :: mLayerBaseflow(:) ! baseflow from each soil layer (m s-1)
+ real(rkind),intent(out) :: dBaseflow_dVolLiq(:,:) ! derivative in baseflow w.r.t. matric head (s-1)
+ ! ---------------------------------------------------------------------------------------
+ ! * local variables
+ ! ---------------------------------------------------------------------------------------
+ ! general local variables
+ integer(i4b) :: iLayer,jLayer ! index of model layer
+ ! local variables for the exfiltration
+ real(rkind) :: totalColumnInflow ! total column inflow (m s-1)
+ real(rkind) :: totalColumnOutflow ! total column outflow (m s-1)
+ real(rkind) :: availStorage ! available storage (m)
+ real(rkind),parameter :: xMinEval=0.002_rkind ! minimum value to evaluate the exfiltration function (m)
+ real(rkind),parameter :: xCenter=0.001_rkind ! center of the exfiltration function (m)
+ real(rkind),parameter :: xWidth=0.0001_rkind ! width of the exfiltration function (m)
+ real(rkind) :: expF,logF ! logistic smoothing function (-)
+ ! local variables for the lateral flux among soil columns
+ real(rkind) :: activePorosity ! "active" porosity associated with storage above a threshold (-)
+ real(rkind) :: drainableWater ! drainable water in eaxch layer (m)
+ real(rkind) :: tran0 ! maximum transmissivity (m2 s-1)
+ real(rkind),dimension(nSoil) :: zActive ! water table thickness associated with storage below and including the given layer (m)
+ real(rkind),dimension(nSoil) :: trTotal ! total transmissivity associated with total water table depth zActive (m2 s-1)
+ real(rkind),dimension(nSoil) :: trSoil ! transmissivity of water in a given layer (m2 s-1)
+ ! local variables for the derivatives
+ real(rkind) :: qbTotal ! total baseflow (m s-1)
+ real(rkind) :: length2area ! ratio of hillslope width to hillslope area (m m-2)
+ real(rkind),dimension(nSoil) :: depth2capacity ! ratio of layer depth to total subsurface storage capacity (-)
+ real(rkind),dimension(nSoil) :: dXdS ! change in dimensionless flux w.r.t. change in dimensionless storage (-)
+ real(rkind),dimension(nSoil) :: dLogFunc_dLiq ! derivative in the logistic function w.r.t. volumetric liquid water content (-)
+ real(rkind),dimension(nSoil) :: dExfiltrate_dVolLiq ! derivative in exfiltration w.r.t. volumetric liquid water content (-)
+ ! local variables for testing (debugging)
+ logical(lgt),parameter :: printFlag=.false. ! flag for printing (debugging)
+ real(rkind) :: xDepth,xTran,xFlow ! temporary variables (depth, transmissivity, flow)
+ ! ---------------------------------------------------------------------------------------
+ ! * association to data in structures
+ ! ---------------------------------------------------------------------------------------
+ associate(&
+
+ ! input: coordinate variables
+ soilDepth => prog_data%var(iLookPROG%iLayerHeight)%dat(nLayers), & ! intent(in): [dp] total soil depth (m)
+ mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat(nSnow+1:nLayers),& ! intent(in): [dp(:)] depth of each soil layer (m)
+
+ ! input: diagnostic variables
+ surfaceHydCond => flux_data%var(iLookFLUX%mLayerSatHydCondMP)%dat(1), & ! intent(in): [dp] saturated hydraulic conductivity at the surface (m s-1)
+ mLayerColumnInflow => flux_data%var(iLookFLUX%mLayerColumnInflow)%dat, & ! intent(in): [dp(:)] inflow into each soil layer (m3/s)
+
+ ! input: local attributes
+ HRUarea => attr_data%var(iLookATTR%HRUarea), & ! intent(in): [dp] HRU area (m2)
+ tan_slope => attr_data%var(iLookATTR%tan_slope), & ! intent(in): [dp] tan water table slope, taken as tan local ground surface slope (-)
+ contourLength => attr_data%var(iLookATTR%contourLength), & ! intent(in): [dp] length of contour at downslope edge of HRU (m)
+
+ ! input: baseflow parameters
+ zScale_TOPMODEL => mpar_data%var(iLookPARAM%zScale_TOPMODEL)%dat(1), & ! intent(in): [dp] TOPMODEL exponent (-)
+ kAnisotropic => mpar_data%var(iLookPARAM%kAnisotropic)%dat(1), & ! intent(in): [dp] anisotropy factor for lateral hydraulic conductivity (-
+ fieldCapacity => mpar_data%var(iLookPARAM%fieldCapacity)%dat(1), & ! intent(in): [dp] field capacity (-)
+ theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat, & ! intent(in): [dp(:)] soil porosity (-)
+
+ ! output: diagnostic variables
+ scalarExfiltration => flux_data%var(iLookFLUX%scalarExfiltration)%dat(1), & ! intent(out):[dp] exfiltration from the soil profile (m s-1)
+ mLayerColumnOutflow => flux_data%var(iLookFLUX%mLayerColumnOutflow)%dat & ! intent(out):[dp(:)] column outflow from each soil layer (m3 s-1)
+
+ ) ! end association to variables in data structures
+
+ ! ***********************************************************************************************************************
+ ! (1) compute the baseflow flux in each soil layer
+ ! ***********************************************************************************************************************
+
+ ! compute the maximum transmissivity
+ ! NOTE: this can be done as a pre-processing step
+ tran0 = kAnisotropic*surfaceHydCond*soilDepth/zScale_TOPMODEL ! maximum transmissivity (m2 s-1)
+
+ ! compute the water table thickness (m) and transmissivity in each layer (m2 s-1)
+ do iLayer=nSoil,ixSaturation,-1 ! loop through "active" soil layers, from lowest to highest
+ ! define drainable water in each layer (m)
+ activePorosity = theta_sat(iLayer) - fieldCapacity ! "active" porosity (-)
+ drainableWater = mLayerDepth(iLayer)*(max(0._rkind,mLayerVolFracLiq(iLayer) - fieldCapacity))/activePorosity
+ ! compute layer transmissivity
+ if(iLayer==nSoil)then
+ zActive(iLayer) = drainableWater ! water table thickness associated with storage in a given layer (m)
+ trTotal(iLayer) = tran0*(zActive(iLayer)/soilDepth)**zScale_TOPMODEL ! total transmissivity for total depth zActive (m2 s-1)
+ trSoil(iLayer) = trTotal(iLayer) ! transmissivity of water in a given layer (m2 s-1)
+ else
+ zActive(iLayer) = zActive(iLayer+1) + drainableWater
+ trTotal(iLayer) = tran0*(zActive(iLayer)/soilDepth)**zScale_TOPMODEL
+ trSoil(iLayer) = trTotal(iLayer) - trTotal(iLayer+1)
+ end if
+ end do ! looping through soil layers
+
+ ! set un-used portions of the vectors to zero
+ if(ixSaturation>1)then
+ zActive(1:ixSaturation-1) = 0._rkind
+ trTotal(1:ixSaturation-1) = 0._rkind
+ trSoil(1:ixSaturation-1) = 0._rkind
+ end if
+
+ ! compute the outflow from each layer (m3 s-1)
+ mLayerColumnOutflow(1:nSoil) = trSoil(1:nSoil)*tan_slope*contourLength
+
+ ! compute total column inflow and total column outflow (m s-1)
+ totalColumnInflow = sum(mLayerColumnInflow(1:nSoil))/HRUarea
+ totalColumnOutflow = sum(mLayerColumnOutflow(1:nSoil))/HRUarea
+
+ ! compute the available storage (m)
+ availStorage = sum(mLayerDepth(1:nSoil)*(theta_sat - (mLayerVolFracLiq(1:nSoil)+mLayerVolFracIce(1:nSoil))) )
+
+ ! compute the smoothing function (-)
+ if(availStorage < xMinEval)then
+ ! (compute the logistic function)
+ expF = exp((availStorage - xCenter)/xWidth)
+ logF = 1._rkind / (1._rkind + expF)
+ ! (compute the derivative in the logistic function w.r.t. volumetric liquid water content in each soil layer)
+ dLogFunc_dLiq(1:nSoil) = mLayerDepth(1:nSoil)*(expF/xWidth)/(1._rkind + expF)**2._rkind
+ else
+ logF = 0._rkind
+ dLogFunc_dLiq(:) = 0._rkind
+ end if
+
+ ! compute the exfiltartion (m s-1)
+ if(totalColumnInflow > totalColumnOutflow .and. logF > tiny(1._rkind))then
+ scalarExfiltration = logF*(totalColumnInflow - totalColumnOutflow) ! m s-1
+ else
+ scalarExfiltration = 0._rkind
+ end if
+
+ ! compute the baseflow in each layer (m s-1)
+ mLayerBaseflow(1:nSoil) = (mLayerColumnOutflow(1:nSoil) - mLayerColumnInflow(1:nSoil))/HRUarea
+
+ ! compute the total baseflow
+ qbTotal = sum(mLayerBaseflow)
+
+ ! add exfiltration to the baseflow flux at the top layer
+ mLayerBaseflow(1) = mLayerBaseflow(1) + scalarExfiltration
+ mLayerColumnOutflow(1) = mLayerColumnOutflow(1) + scalarExfiltration*HRUarea
+
+ ! test
+ if(printFlag)then
+ xDepth = sum(mLayerDepth(ixSaturation:nSoil)*(mLayerVolFracLiq(ixSaturation:nSoil) - fieldCapacity))/sum(theta_sat(ixSaturation:nSoil) - fieldCapacity) ! "effective" water table thickness (m)
+ xTran = tran0*(xDepth/soilDepth)**zScale_TOPMODEL ! transmissivity for the entire aquifer (m2 s-1)
+ xFlow = xTran*tan_slope*contourLength/HRUarea ! total column outflow (m s-1)
+ print*, 'ixSaturation = ', ixSaturation
+ write(*,'(a,1x,5(f30.20,1x))') 'tran0, soilDepth = ', tran0, soilDepth
+ write(*,'(a,1x,5(f30.20,1x))') 'surfaceHydCond, zScale_TOPMODEL = ', surfaceHydCond, zScale_TOPMODEL
+ write(*,'(a,1x,5(f30.20,1x))') 'xDepth, zActive(ixSaturation) = ', xDepth, zActive(ixSaturation)
+ write(*,'(a,1x,5(f30.20,1x))') 'xTran, trTotal(ixSaturation) = ', xTran, trTotal(ixSaturation)
+ write(*,'(a,1x,5(f30.20,1x))') 'xFlow, totalColumnOutflow = ', xFlow, sum(mLayerColumnOutflow(:))/HRUarea
+ !pause 'check groundwater'
+ end if
+
+ ! ***********************************************************************************************************************
+ ! (2) compute the derivative in the baseflow flux w.r.t. volumetric liquid water content (m s-1)
+ ! ***********************************************************************************************************************
+
+ ! initialize the derivative matrix
+ dBaseflow_dVolLiq(:,:) = 0._rkind
+
+ ! check if derivatives are actually required
+ if(.not.derivDesired) return
+
+ ! compute ratio of hillslope width to hillslope area (m m-2)
+ length2area = tan_slope*contourLength/HRUarea
+
+ ! compute the ratio of layer depth to maximum water holding capacity (-)
+ depth2capacity(1:nSoil) = mLayerDepth(1:nSoil)/sum( (theta_sat(1:nSoil) - fieldCapacity)*mLayerDepth(1:nSoil) )
+
+ ! compute the change in dimensionless flux w.r.t. change in dimensionless storage (-)
+ dXdS(1:nSoil) = zScale_TOPMODEL*(zActive(1:nSoil)/SoilDepth)**(zScale_TOPMODEL - 1._rkind)
+
+ ! loop through soil layers
+ do iLayer=1,nSoil
+ ! compute diagonal terms (s-1)
+ dBaseflow_dVolLiq(iLayer,iLayer) = tran0*dXdS(iLayer)*depth2capacity(iLayer)*length2area
+ ! compute off-diagonal terms
+ do jLayer=iLayer+1,nSoil ! (only dependent on layers below)
+ dBaseflow_dVolLiq(iLayer,jLayer) = tran0*(dXdS(iLayer) - dXdS(iLayer+1))*depth2capacity(jLayer)*length2area
+ end do ! looping through soil layers
+ end do ! looping through soil layers
+
+ ! compute the derivative in the exfiltration flux w.r.t. volumetric liquid water content (m s-1)
+ if(qbTotal < 0._rkind)then
+ do iLayer=1,nSoil
+ dExfiltrate_dVolLiq(iLayer) = dBaseflow_dVolLiq(iLayer,iLayer)*logF + dLogFunc_dLiq(iLayer)*qbTotal
+ end do ! looping through soil layers
+ dBaseflow_dVolLiq(1,1:nSoil) = dBaseflow_dVolLiq(1,1:nSoil) - dExfiltrate_dVolLiq(1:nSoil)
+ end if
+
+ ! end association to data in structures
+ end associate
+
+end subroutine computeBaseflow
end module groundwatr_module
diff --git a/build/source/engine/snowLiqFlx.f90 b/build/source/engine/snowLiqFlx.f90
index 48566b9ccfd840778a30f1254914beb9b4adaf47..e31612523722a9a7622570594ea15b933c0c73ac 100644
--- a/build/source/engine/snowLiqFlx.f90
+++ b/build/source/engine/snowLiqFlx.f90
@@ -20,92 +20,92 @@
module snowLiqFlx_module
- ! access modules
- USE nrtype ! numerical recipes data types
- USE multiconst,only:iden_ice,iden_water ! intrinsic density of ice and water (kg m-3)
-
- ! access missing values
- USE globalData,only:integerMissing ! missing integer
- USE globalData,only:realMissing ! missing real number
-
- ! named variables
- USE var_lookup,only:iLookINDEX ! named variables for structure elements
- USE var_lookup,only:iLookPARAM ! 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
-
- ! data types
- USE data_types,only:var_d ! x%var(:) (rkind)
- USE data_types,only:var_dlength ! x%var(:)%dat (rkind)
- USE data_types,only:var_ilength ! x%var(:)%dat (i4b)
-
- ! privacy
+! access modules
+USE nrtype ! numerical recipes data types
+USE multiconst,only:iden_ice,iden_water ! intrinsic density of ice and water (kg m-3)
+
+! access missing values
+USE globalData,only:integerMissing ! missing integer
+USE globalData,only:realMissing ! missing real number
+
+! named variables
+USE var_lookup,only:iLookINDEX ! named variables for structure elements
+USE var_lookup,only:iLookPARAM ! 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
+
+! data types
+USE data_types,only:var_d ! x%var(:) (rkind)
+USE data_types,only:var_dlength ! x%var(:)%dat (rkind)
+USE data_types,only:var_ilength ! x%var(:)%dat (i4b)
+
+! privacy
+implicit none
+private
+public::snowLiqFlx
+public::snowLiqFlxSundials
+contains
+
+
+! ************************************************************************************************
+! public subroutine snowLiqFlx: compute liquid water flux through the snowpack
+! ************************************************************************************************
+subroutine snowLiqFlx(&
+ ! input: model control
+ nSnow, & ! intent(in): number of snow layers
+ firstFluxCall, & ! intent(in): the first flux call
+ scalarSolution, & ! intent(in): flag to indicate the scalar solution
+ ! input: forcing for the snow domain
+ scalarThroughfallRain, & ! intent(in): rain that reaches the snow surface without ever touching vegetation (kg m-2 s-1)
+ scalarCanopyLiqDrainage, & ! intent(in): liquid drainage from the vegetation canopy (kg m-2 s-1)
+ ! input: model state vector
+ mLayerVolFracLiqTrial, & ! intent(in): trial value of volumetric fraction of liquid water at the current iteration (-)
+ ! input-output: data structures
+ indx_data, & ! intent(in): model indices
+ mpar_data, & ! intent(in): model parameters
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ ! output: fluxes and derivatives
+ iLayerLiqFluxSnow, & ! intent(inout): vertical liquid water flux at layer interfaces (m s-1)
+ iLayerLiqFluxSnowDeriv, & ! intent(inout): derivative in vertical liquid water flux at layer interfaces (m s-1)
+ ! output: error control
+ err,message) ! intent(out): error control
implicit none
- private
- public::snowLiqFlx
- public::snowLiqFlxSundials
- contains
-
-
- ! ************************************************************************************************
- ! public subroutine snowLiqFlx: compute liquid water flux through the snowpack
- ! ************************************************************************************************
- subroutine snowLiqFlx(&
- ! input: model control
- nSnow, & ! intent(in): number of snow layers
- firstFluxCall, & ! intent(in): the first flux call
- scalarSolution, & ! intent(in): flag to indicate the scalar solution
- ! input: forcing for the snow domain
- scalarThroughfallRain, & ! intent(in): rain that reaches the snow surface without ever touching vegetation (kg m-2 s-1)
- scalarCanopyLiqDrainage, & ! intent(in): liquid drainage from the vegetation canopy (kg m-2 s-1)
- ! input: model state vector
- mLayerVolFracLiqTrial, & ! intent(in): trial value of volumetric fraction of liquid water at the current iteration (-)
- ! input-output: data structures
- indx_data, & ! intent(in): model indices
- mpar_data, & ! intent(in): model parameters
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- ! output: fluxes and derivatives
- iLayerLiqFluxSnow, & ! intent(inout): vertical liquid water flux at layer interfaces (m s-1)
- iLayerLiqFluxSnowDeriv, & ! intent(inout): derivative in vertical liquid water flux at layer interfaces (m s-1)
- ! output: error control
- err,message) ! intent(out): error control
- implicit none
- ! input: model control
- integer(i4b),intent(in) :: nSnow ! number of snow layers
- logical(lgt),intent(in) :: firstFluxCall ! the first flux call
- logical(lgt),intent(in) :: scalarSolution ! flag to denote if implementing the scalar solution
- ! input: forcing for the snow domain
- real(rkind),intent(in) :: scalarThroughfallRain ! computed throughfall rate (kg m-2 s-1)
- real(rkind),intent(in) :: scalarCanopyLiqDrainage ! computed drainage of liquid water (kg m-2 s-1)
- ! input: model state vector
- real(rkind),intent(in) :: mLayerVolFracLiqTrial(:) ! trial value of volumetric fraction of liquid water at the current iteration (-)
- ! input-output: data structures
- type(var_ilength),intent(in) :: indx_data ! model indices
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- 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
- ! output: fluxes and derivatives
- real(rkind),intent(inout) :: iLayerLiqFluxSnow(0:) ! vertical liquid water flux at layer interfaces (m s-1)
- real(rkind),intent(inout) :: iLayerLiqFluxSnowDeriv(0:) ! derivative in vertical liquid water flux at layer interfaces (m s-1)
- ! output: error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! ------------------------------------------------------------------------------------------------------------------------------------------
- ! local variables
- integer(i4b) :: i ! search index for scalar solution
- integer(i4b) :: iLayer ! layer index
- integer(i4b) :: ixTop ! top layer in subroutine call
- integer(i4b) :: ixBot ! bottom layer in subroutine call
- real(rkind) :: multResid ! multiplier for the residual water content (-)
- real(rkind),parameter :: residThrs=550._rkind ! ice density threshold to reduce residual liquid water content (kg m-3)
- real(rkind),parameter :: residScal=10._rkind ! scaling factor for residual liquid water content reduction factor (kg m-3)
- real(rkind),parameter :: maxVolIceContent=0.7_rkind ! maximum volumetric ice content to store water (-)
- real(rkind) :: availCap ! available storage capacity [0,1] (-)
- real(rkind) :: relSaturn ! relative saturation [0,1] (-)
- ! ------------------------------------------------------------------------------------------------------------------------------------------
- ! make association of local variables with information in the data structures
- associate(&
+ ! input: model control
+ integer(i4b),intent(in) :: nSnow ! number of snow layers
+ logical(lgt),intent(in) :: firstFluxCall ! the first flux call
+ logical(lgt),intent(in) :: scalarSolution ! flag to denote if implementing the scalar solution
+ ! input: forcing for the snow domain
+ real(rkind),intent(in) :: scalarThroughfallRain ! computed throughfall rate (kg m-2 s-1)
+ real(rkind),intent(in) :: scalarCanopyLiqDrainage ! computed drainage of liquid water (kg m-2 s-1)
+ ! input: model state vector
+ real(rkind),intent(in) :: mLayerVolFracLiqTrial(:) ! trial value of volumetric fraction of liquid water at the current iteration (-)
+ ! input-output: data structures
+ type(var_ilength),intent(in) :: indx_data ! model indices
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ 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
+ ! output: fluxes and derivatives
+ real(rkind),intent(inout) :: iLayerLiqFluxSnow(0:) ! vertical liquid water flux at layer interfaces (m s-1)
+ real(rkind),intent(inout) :: iLayerLiqFluxSnowDeriv(0:) ! derivative in vertical liquid water flux at layer interfaces (m s-1)
+ ! output: error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! ------------------------------------------------------------------------------------------------------------------------------------------
+ ! local variables
+ integer(i4b) :: i ! search index for scalar solution
+ integer(i4b) :: iLayer ! layer index
+ integer(i4b) :: ixTop ! top layer in subroutine call
+ integer(i4b) :: ixBot ! bottom layer in subroutine call
+ real(rkind) :: multResid ! multiplier for the residual water content (-)
+ real(rkind),parameter :: residThrs=550._rkind ! ice density threshold to reduce residual liquid water content (kg m-3)
+ real(rkind),parameter :: residScal=10._rkind ! scaling factor for residual liquid water content reduction factor (kg m-3)
+ real(rkind),parameter :: maxVolIceContent=0.7_rkind ! maximum volumetric ice content to store water (-)
+ real(rkind) :: availCap ! available storage capacity [0,1] (-)
+ real(rkind) :: relSaturn ! relative saturation [0,1] (-)
+ ! ------------------------------------------------------------------------------------------------------------------------------------------
+ ! make association of local variables with information in the data structures
+ associate(&
! input: layer indices
ixLayerState => indx_data%var(iLookINDEX%ixLayerState)%dat, & ! intent(in): list of indices for all model layers
ixSnowOnlyHyd => indx_data%var(iLookINDEX%ixSnowOnlyHyd)%dat, & ! intent(in): index in the state subset for hydrology state variables in the snow domain
@@ -117,151 +117,140 @@ module snowLiqFlx_module
! input/output: diagnostic variables -- only computed for the first iteration
mLayerPoreSpace => diag_data%var(iLookDIAG%mLayerPoreSpace)%dat, & ! intent(inout): pore space in each snow layer (-)
mLayerThetaResid => diag_data%var(iLookDIAG%mLayerThetaResid)%dat & ! intent(inout): esidual volumetric liquid water content in each snow layer (-)
- ) ! association of local variables with information in the data structures
- ! ------------------------------------------------------------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message='snowLiqFlx/'
-
- ! check that the input vectors match nSnow
- if(size(mLayerVolFracLiqTrial)/=nSnow .or. size(mLayerVolFracIce)/=nSnow .or. &
- size(iLayerLiqFluxSnow)/=nSnow+1 .or. size(iLayerLiqFluxSnowDeriv)/=nSnow+1) then
- err=20; message=trim(message)//'size mismatch of input/output vectors'; return
- end if
-
- ! check the meltwater exponent is >=1
- if(mw_exp<1._rkind)then; err=20; message=trim(message)//'meltwater exponent < 1'; return; end if
-
- ! get the indices for the snow+soil layers
- ixTop = integerMissing
- if(scalarSolution)then
- ! WARNING: Previously this was implemented as:
- ! ixLayerDesired = pack(ixLayerState, ixSnowOnlyHyd/=integerMissing)
- ! ixTop = ixLayerDesired(1)
- ! ixBot = ixLayerDesired(1)
- ! This implementation can result in a segfault when using JRDN layering.
- ! The segfault occurs when trying to access `mw_exp` in:
- ! iLayerLiqFluxSnow(iLayer) = k_snow*relSaturn**mw_exp
- ! Debugging found that the `pack` statement caused `mw_exp` to no longer be accessible.
- ! We have not been able to determine the underlying reason for this segfault.
- do i=1,size(ixSnowOnlyHyd)
- if(ixSnowOnlyHyd(i) /= integerMissing)then
- ixTop=ixLayerState(i)
- ixBot=ixTop
- exit ! break out of loop once found
- endif
- end do
- if(ixTop == integerMissing)then
- err=20; message=trim(message)//'Unable to identify snow layer for scalar solution!'; return
+ ) ! association of local variables with information in the data structures
+ ! ------------------------------------------------------------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message='snowLiqFlx/'
+
+ ! check that the input vectors match nSnow
+ if(size(mLayerVolFracLiqTrial)/=nSnow .or. size(mLayerVolFracIce)/=nSnow .or. &
+ size(iLayerLiqFluxSnow)/=nSnow+1 .or. size(iLayerLiqFluxSnowDeriv)/=nSnow+1) then
+ err=20; message=trim(message)//'size mismatch of input/output vectors'; return
end if
- else
- ixTop = 1
- ixBot = nSnow
- endif
-
- ! define the liquid flux at the upper boundary (m s-1)
- iLayerLiqFluxSnow(0) = (scalarThroughfallRain + scalarCanopyLiqDrainage)/iden_water
- iLayerLiqFluxSnowDeriv(0) = 0._rkind !computed inside computJacobSundials_module
-
- ! compute properties fixed over the time step
- if(firstFluxCall)then
- ! loop through snow layers
- do iLayer=1,nSnow
- ! compute the reduction in liquid water holding capacity at high snow density (-)
- multResid = 1._rkind / ( 1._rkind + exp( (mLayerVolFracIce(iLayer)*iden_ice - residThrs) / residScal) )
- ! compute the pore space (-)
- mLayerPoreSpace(iLayer) = 1._rkind - mLayerVolFracIce(iLayer)
- ! compute the residual volumetric liquid water content (-)
- mLayerThetaResid(iLayer) = Fcapil*mLayerPoreSpace(iLayer) * multResid
- end do ! (looping through snow layers)
- end if ! (if the first flux call)
-
- ! compute fluxes
- do iLayer=ixTop,ixBot ! (loop through snow layers)
- ! check that flow occurs
- if(mLayerVolFracLiqTrial(iLayer) > mLayerThetaResid(iLayer))then
- ! compute the relative saturation (-)
- availCap = mLayerPoreSpace(iLayer) - mLayerThetaResid(iLayer) ! available capacity
- relSaturn = (mLayerVolFracLiqTrial(iLayer) - mLayerThetaResid(iLayer)) / availCap ! relative saturation
- iLayerLiqFluxSnow(iLayer) = k_snow*relSaturn**mw_exp
- iLayerLiqFluxSnowDeriv(iLayer) = ( (k_snow*mw_exp)/availCap ) * relSaturn**(mw_exp - 1._rkind)
- if(mLayerVolFracIce(iLayer) > maxVolIceContent)then ! NOTE: use start-of-step ice content, to avoid convergence problems
- ! ** allow liquid water to pass through under very high ice density
- iLayerLiqFluxSnow(iLayer) = iLayerLiqFluxSnow(iLayer) + iLayerLiqFluxSnow(iLayer-1) !NOTE: derivative may need to be updated in future.
- end if
- else ! flow does not occur
- iLayerLiqFluxSnow(iLayer) = 0._rkind
- iLayerLiqFluxSnowDeriv(iLayer) = 0._rkind
- endif ! storage above residual content
- end do ! loop through snow layers
-
- ! end association of local variables with information in the data structures
- end associate
-
- end subroutine snowLiqFlx
-
-
-
- ! ************************************************************************************************
- ! public subroutine snowLiqFlxSundials: compute liquid water flux through the snowpack
- ! ************************************************************************************************
- subroutine snowLiqFlxSundials(&
- ! input: model control
- nSnow, & ! intent(in): number of snow layers
- firstFluxCall, & ! intent(in): the first flux call
- scalarSolution, & ! intent(in): flag to indicate the scalar solution
- ! input: forcing for the snow domain
- scalarThroughfallRain, & ! intent(in): rain that reaches the snow surface without ever touching vegetation (kg m-2 s-1)
- scalarCanopyLiqDrainage, & ! intent(in): liquid drainage from the vegetation canopy (kg m-2 s-1)
- ! input: model state vector
- mLayerVolFracIce, & ! intent(in)
- mLayerVolFracLiqTrial, & ! intent(in): trial value of volumetric fraction of liquid water at the current iteration (-)
- ! input-output: data structures
- indx_data, & ! intent(in): model indices
- mpar_data, & ! intent(in): model parameters
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- ! output: fluxes and derivatives
- iLayerLiqFluxSnow, & ! intent(inout): vertical liquid water flux at layer interfaces (m s-1)
- iLayerLiqFluxSnowDeriv, & ! intent(inout): derivative in vertical liquid water flux at layer interfaces (m s-1)
- ! output: error control
- err,message) ! intent(out): error control
- implicit none
- ! input: model control
- integer(i4b),intent(in) :: nSnow ! number of snow layers
- logical(lgt),intent(in) :: firstFluxCall ! the first flux call
- logical(lgt),intent(in) :: scalarSolution ! flag to denote if implementing the scalar solution
- ! input: forcing for the snow domain
- real(rkind),intent(in) :: scalarThroughfallRain ! computed throughfall rate (kg m-2 s-1)
- real(rkind),intent(in) :: scalarCanopyLiqDrainage ! computed drainage of liquid water (kg m-2 s-1)
- ! input: model state vector
- real(rkind),intent(in) :: mLayerVolFracIce(:)
- real(rkind),intent(in) :: mLayerVolFracLiqTrial(:) ! trial value of volumetric fraction of liquid water at the current iteration (-)
- ! input-output: data structures
- type(var_ilength),intent(in) :: indx_data ! model indices
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- 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
- ! output: fluxes and derivatives
- real(rkind),intent(inout) :: iLayerLiqFluxSnow(0:) ! vertical liquid water flux at layer interfaces (m s-1)
- real(rkind),intent(inout) :: iLayerLiqFluxSnowDeriv(0:) ! derivative in vertical liquid water flux at layer interfaces (m s-1)
- ! output: error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! ------------------------------------------------------------------------------------------------------------------------------------------
- ! local variables
- integer(i4b) :: i ! search index for scalar solution
- integer(i4b) :: iLayer ! layer index
- integer(i4b) :: ixTop ! top layer in subroutine call
- integer(i4b) :: ixBot ! bottom layer in subroutine call
- real(rkind) :: multResid ! multiplier for the residual water content (-)
- real(rkind),parameter :: residThrs=550._rkind ! ice density threshold to reduce residual liquid water content (kg m-3)
- real(rkind),parameter :: residScal=10._rkind ! scaling factor for residual liquid water content reduction factor (kg m-3)
- real(rkind),parameter :: maxVolIceContent=0.7_rkind ! maximum volumetric ice content to store water (-)
- real(rkind) :: availCap ! available storage capacity [0,1] (-)
- real(rkind) :: relSaturn ! relative saturation [0,1] (-)
- ! ------------------------------------------------------------------------------------------------------------------------------------------
- ! make association of local variables with information in the data structures
- associate(&
+
+ ! check the meltwater exponent is >=1
+ if(mw_exp<1._rkind)then; err=20; message=trim(message)//'meltwater exponent < 1'; return; end if
+
+ ! get the indices for the snow+soil layers
+ ixTop = integerMissing
+ if(scalarSolution)then
+ do i=1,size(ixSnowOnlyHyd)
+ if(ixSnowOnlyHyd(i) /= integerMissing)then
+ ixTop=ixLayerState(i)
+ ixBot=ixTop
+ exit ! break out of loop once found
+ endif
+ end do
+ if(ixTop == integerMissing)then
+ err=20; message=trim(message)//'Unable to identify snow layer for scalar solution!'; return
+ end if
+ else
+ ixTop = 1
+ ixBot = nSnow
+ endif
+
+ ! define the liquid flux at the upper boundary (m s-1)
+ iLayerLiqFluxSnow(0) = (scalarThroughfallRain + scalarCanopyLiqDrainage)/iden_water
+ iLayerLiqFluxSnowDeriv(0) = 0._rkind !computed inside computJacobSundials_module
+
+ ! compute properties fixed over the time step
+ if(firstFluxCall)then
+ ! loop through snow layers
+ do iLayer=1,nSnow
+ ! compute the reduction in liquid water holding capacity at high snow density (-)
+ multResid = 1._rkind / ( 1._rkind + exp( (mLayerVolFracIce(iLayer)*iden_ice - residThrs) / residScal) )
+ ! compute the pore space (-)
+ mLayerPoreSpace(iLayer) = 1._rkind - mLayerVolFracIce(iLayer)
+ ! compute the residual volumetric liquid water content (-)
+ mLayerThetaResid(iLayer) = Fcapil*mLayerPoreSpace(iLayer) * multResid
+ end do ! (looping through snow layers)
+ end if ! (if the first flux call)
+
+ ! compute fluxes
+ do iLayer=ixTop,ixBot ! (loop through snow layers)
+ ! check that flow occurs
+ if(mLayerVolFracLiqTrial(iLayer) > mLayerThetaResid(iLayer))then
+ ! compute the relative saturation (-)
+ availCap = mLayerPoreSpace(iLayer) - mLayerThetaResid(iLayer) ! available capacity
+ relSaturn = (mLayerVolFracLiqTrial(iLayer) - mLayerThetaResid(iLayer)) / availCap ! relative saturation
+ iLayerLiqFluxSnow(iLayer) = k_snow*relSaturn**mw_exp
+ iLayerLiqFluxSnowDeriv(iLayer) = ( (k_snow*mw_exp)/availCap ) * relSaturn**(mw_exp - 1._rkind)
+ if(mLayerVolFracIce(iLayer) > maxVolIceContent)then ! NOTE: use start-of-step ice content, to avoid convergence problems
+ ! ** allow liquid water to pass through under very high ice density
+ iLayerLiqFluxSnow(iLayer) = iLayerLiqFluxSnow(iLayer) + iLayerLiqFluxSnow(iLayer-1) !NOTE: derivative may need to be updated in future.
+ end if
+ else ! flow does not occur
+ iLayerLiqFluxSnow(iLayer) = 0._rkind
+ iLayerLiqFluxSnowDeriv(iLayer) = 0._rkind
+ endif ! storage above residual content
+ end do ! loop through snow layers
+
+ ! end association of local variables with information in the data structures
+ end associate
+
+end subroutine snowLiqFlx
+
+! ************************************************************************************************
+! public subroutine snowLiqFlxSundials: compute liquid water flux through the snowpack
+! ************************************************************************************************
+subroutine snowLiqFlxSundials(&
+ ! input: model control
+ nSnow, & ! intent(in): number of snow layers
+ firstFluxCall, & ! intent(in): the first flux call
+ scalarSolution, & ! intent(in): flag to indicate the scalar solution
+ ! input: forcing for the snow domain
+ scalarThroughfallRain, & ! intent(in): rain that reaches the snow surface without ever touching vegetation (kg m-2 s-1)
+ scalarCanopyLiqDrainage, & ! intent(in): liquid drainage from the vegetation canopy (kg m-2 s-1)
+ ! input: model state vector
+ mLayerVolFracIce, & ! intent(in)
+ mLayerVolFracLiqTrial, & ! intent(in): trial value of volumetric fraction of liquid water at the current iteration (-)
+ ! input-output: data structures
+ indx_data, & ! intent(in): model indices
+ mpar_data, & ! intent(in): model parameters
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ ! output: fluxes and derivatives
+ iLayerLiqFluxSnow, & ! intent(inout): vertical liquid water flux at layer interfaces (m s-1)
+ iLayerLiqFluxSnowDeriv, & ! intent(inout): derivative in vertical liquid water flux at layer interfaces (m s-1)
+ ! output: error control
+ err,message) ! intent(out): error control
+ implicit none
+ ! input: model control
+ integer(i4b),intent(in) :: nSnow ! number of snow layers
+ logical(lgt),intent(in) :: firstFluxCall ! the first flux call
+ logical(lgt),intent(in) :: scalarSolution ! flag to denote if implementing the scalar solution
+ ! input: forcing for the snow domain
+ real(rkind),intent(in) :: scalarThroughfallRain ! computed throughfall rate (kg m-2 s-1)
+ real(rkind),intent(in) :: scalarCanopyLiqDrainage ! computed drainage of liquid water (kg m-2 s-1)
+ ! input: model state vector
+ real(rkind),intent(in) :: mLayerVolFracIce(:)
+ real(rkind),intent(in) :: mLayerVolFracLiqTrial(:) ! trial value of volumetric fraction of liquid water at the current iteration (-)
+ ! input-output: data structures
+ type(var_ilength),intent(in) :: indx_data ! model indices
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ 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
+ ! output: fluxes and derivatives
+ real(rkind),intent(inout) :: iLayerLiqFluxSnow(0:) ! vertical liquid water flux at layer interfaces (m s-1)
+ real(rkind),intent(inout) :: iLayerLiqFluxSnowDeriv(0:) ! derivative in vertical liquid water flux at layer interfaces (m s-1)
+ ! output: error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! ------------------------------------------------------------------------------------------------------------------------------------------
+ ! local variables
+ integer(i4b) :: i ! search index for scalar solution
+ integer(i4b) :: iLayer ! layer index
+ integer(i4b) :: ixTop ! top layer in subroutine call
+ integer(i4b) :: ixBot ! bottom layer in subroutine call
+ real(rkind) :: multResid ! multiplier for the residual water content (-)
+ real(rkind),parameter :: residThrs=550._rkind ! ice density threshold to reduce residual liquid water content (kg m-3)
+ real(rkind),parameter :: residScal=10._rkind ! scaling factor for residual liquid water content reduction factor (kg m-3)
+ real(rkind),parameter :: maxVolIceContent=0.7_rkind ! maximum volumetric ice content to store water (-)
+ real(rkind) :: availCap ! available storage capacity [0,1] (-)
+ real(rkind) :: relSaturn ! relative saturation [0,1] (-)
+ ! ------------------------------------------------------------------------------------------------------------------------------------------
+ ! make association of local variables with information in the data structures
+ associate(&
! input: layer indices
ixLayerState => indx_data%var(iLookINDEX%ixLayerState)%dat, & ! intent(in): list of indices for all model layers
ixSnowOnlyHyd => indx_data%var(iLookINDEX%ixSnowOnlyHyd)%dat, & ! intent(in): index in the state subset for hydrology state variables in the snow domain
@@ -272,88 +261,75 @@ module snowLiqFlx_module
! input/output: diagnostic variables -- only computed for the first iteration
mLayerPoreSpace => diag_data%var(iLookDIAG%mLayerPoreSpace)%dat, & ! intent(inout): pore space in each snow layer (-)
mLayerThetaResid => diag_data%var(iLookDIAG%mLayerThetaResid)%dat & ! intent(inout): esidual volumetric liquid water content in each snow layer (-)
- ) ! association of local variables with information in the data structures
- ! ------------------------------------------------------------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message='snowLiqFlxSundials/'
-
- ! check that the input vectors match nSnow
- if(size(mLayerVolFracLiqTrial)/=nSnow .or. size(mLayerVolFracIce)/=nSnow .or. &
- size(iLayerLiqFluxSnow)/=nSnow+1 .or. size(iLayerLiqFluxSnowDeriv)/=nSnow+1) then
- err=20; message=trim(message)//'size mismatch of input/output vectors'; return
- end if
-
- ! check the meltwater exponent is >=1
- if(mw_exp<1._rkind)then; err=20; message=trim(message)//'meltwater exponent < 1'; return; end if
-
- ! get the indices for the snow+soil layers
- ixTop = integerMissing
- if(scalarSolution)then
- ! WARNING: Previously this was implemented as:
- ! ixLayerDesired = pack(ixLayerState, ixSnowOnlyHyd/=integerMissing)
- ! ixTop = ixLayerDesired(1)
- ! ixBot = ixLayerDesired(1)
- ! This implementation can result in a segfault when using JRDN layering.
- ! The segfault occurs when trying to access `mw_exp` in:
- ! iLayerLiqFluxSnow(iLayer) = k_snow*relSaturn**mw_exp
- ! Debugging found that the `pack` statement caused `mw_exp` to no longer be accessible.
- ! We have not been able to determine the underlying reason for this segfault.
- do i=1,size(ixSnowOnlyHyd)
- if(ixSnowOnlyHyd(i) /= integerMissing)then
- ixTop=ixLayerState(i)
- ixBot=ixTop
- exit ! break out of loop once found
- endif
- end do
- if(ixTop == integerMissing)then
- err=20; message=trim(message)//'Unable to identify snow layer for scalar solution!'; return
+ ) ! association of local variables with information in the data structures
+ ! ------------------------------------------------------------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message='snowLiqFlxSundials/'
+
+ ! check that the input vectors match nSnow
+ if(size(mLayerVolFracLiqTrial)/=nSnow .or. size(mLayerVolFracIce)/=nSnow .or. &
+ size(iLayerLiqFluxSnow)/=nSnow+1 .or. size(iLayerLiqFluxSnowDeriv)/=nSnow+1) then
+ err=20; message=trim(message)//'size mismatch of input/output vectors'; return
end if
- else
- ixTop = 1
- ixBot = nSnow
- endif
-
- ! define the liquid flux at the upper boundary (m s-1)
- iLayerLiqFluxSnow(0) = (scalarThroughfallRain + scalarCanopyLiqDrainage)/iden_water
- iLayerLiqFluxSnowDeriv(0) = 0._rkind
-
- ! compute properties fixed over the time step
- if(firstFluxCall)then
- ! loop through snow layers
- do iLayer=1,nSnow
- ! compute the reduction in liquid water holding capacity at high snow density (-)
- multResid = 1._rkind / ( 1._rkind + exp( (mLayerVolFracIce(iLayer)*iden_ice - residThrs) / residScal) )
- ! compute the pore space (-)
- mLayerPoreSpace(iLayer) = 1._rkind - mLayerVolFracIce(iLayer)
- ! compute the residual volumetric liquid water content (-)
- mLayerThetaResid(iLayer) = Fcapil*mLayerPoreSpace(iLayer) * multResid
- end do ! (looping through snow layers)
- end if ! (if the first flux call)
-
- ! compute fluxes
- do iLayer=ixTop,ixBot ! (loop through snow layers)
- ! check that flow occurs
- if(mLayerVolFracLiqTrial(iLayer) > mLayerThetaResid(iLayer))then
- ! compute the relative saturation (-)
- availCap = mLayerPoreSpace(iLayer) - mLayerThetaResid(iLayer) ! available capacity
- relSaturn = (mLayerVolFracLiqTrial(iLayer) - mLayerThetaResid(iLayer)) / availCap ! relative saturation
- iLayerLiqFluxSnow(iLayer) = k_snow*relSaturn**mw_exp
- iLayerLiqFluxSnowDeriv(iLayer) = ( (k_snow*mw_exp)/availCap ) * relSaturn**(mw_exp - 1._rkind)
- if(mLayerVolFracIce(iLayer) > maxVolIceContent)then ! NOTE: use start-of-step ice content, to avoid convergence problems
- ! ** allow liquid water to pass through under very high ice density
- iLayerLiqFluxSnow(iLayer) = iLayerLiqFluxSnow(iLayer) + iLayerLiqFluxSnow(iLayer-1) !NOTE: derivative may need to be updated in future.
- end if
- else ! flow does not occur
- iLayerLiqFluxSnow(iLayer) = 0._rkind
- iLayerLiqFluxSnowDeriv(iLayer) = 0._rkind
- endif ! storage above residual content
- end do ! loop through snow layers
-
- ! end association of local variables with information in the data structures
- end associate
-
- end subroutine snowLiqFlxSundials
-
-
- end module snowLiqFlx_module
-
\ No newline at end of file
+
+ ! check the meltwater exponent is >=1
+ if(mw_exp<1._rkind)then; err=20; message=trim(message)//'meltwater exponent < 1'; return; end if
+
+ ! get the indices for the snow+soil layers
+ ixTop = integerMissing
+ if(scalarSolution)then
+ do i=1,size(ixSnowOnlyHyd)
+ if(ixSnowOnlyHyd(i) /= integerMissing)then
+ ixTop=ixLayerState(i)
+ ixBot=ixTop
+ exit ! break out of loop once found
+ endif
+ end do
+ if(ixTop == integerMissing)then
+ err=20; message=trim(message)//'Unable to identify snow layer for scalar solution!'; return
+ end if
+ else
+ ixTop = 1
+ ixBot = nSnow
+ endif
+
+ ! define the liquid flux at the upper boundary (m s-1)
+ iLayerLiqFluxSnow(0) = (scalarThroughfallRain + scalarCanopyLiqDrainage)/iden_water
+ iLayerLiqFluxSnowDeriv(0) = 0._rkind
+
+ ! compute properties fixed over the time step
+ if(firstFluxCall)then
+ ! loop through snow layers
+ do iLayer=1,nSnow
+ ! compute the reduction in liquid water holding capacity at high snow density (-)
+ multResid = 1._rkind / ( 1._rkind + exp( (mLayerVolFracIce(iLayer)*iden_ice - residThrs) / residScal) )
+ ! compute the pore space (-)
+ mLayerPoreSpace(iLayer) = 1._rkind - mLayerVolFracIce(iLayer)
+ ! compute the residual volumetric liquid water content (-)
+ mLayerThetaResid(iLayer) = Fcapil*mLayerPoreSpace(iLayer) * multResid
+ end do ! (looping through snow layers)
+ end if ! (if the first flux call)
+
+ ! compute fluxes
+ do iLayer=ixTop,ixBot ! (loop through snow layers)
+ ! check that flow occurs
+ if(mLayerVolFracLiqTrial(iLayer) > mLayerThetaResid(iLayer))then
+ ! compute the relative saturation (-)
+ availCap = mLayerPoreSpace(iLayer) - mLayerThetaResid(iLayer) ! available capacity
+ relSaturn = (mLayerVolFracLiqTrial(iLayer) - mLayerThetaResid(iLayer)) / availCap ! relative saturation
+ iLayerLiqFluxSnow(iLayer) = k_snow*relSaturn**mw_exp
+ iLayerLiqFluxSnowDeriv(iLayer) = ( (k_snow*mw_exp)/availCap ) * relSaturn**(mw_exp - 1._rkind)
+ if(mLayerVolFracIce(iLayer) > maxVolIceContent)then ! NOTE: use start-of-step ice content, to avoid convergence problems
+ ! ** allow liquid water to pass through under very high ice density
+ iLayerLiqFluxSnow(iLayer) = iLayerLiqFluxSnow(iLayer) + iLayerLiqFluxSnow(iLayer-1) !NOTE: derivative may need to be updated in future.
+ end if
+ else ! flow does not occur
+ iLayerLiqFluxSnow(iLayer) = 0._rkind
+ iLayerLiqFluxSnowDeriv(iLayer) = 0._rkind
+ endif ! storage above residual content
+ end do ! loop through snow layers
+ ! end association of local variables with information in the data structures
+ end associate
+
+end subroutine snowLiqFlxSundials
+end module snowLiqFlx_module