diff --git a/build/source/engine/snowLiqFlx.f90 b/build/source/engine/snowLiqFlx.f90
index 6778561f876a97ebb584ba918c37c4373be3a9da..48566b9ccfd840778a30f1254914beb9b4adaf47 100644
--- a/build/source/engine/snowLiqFlx.f90
+++ b/build/source/engine/snowLiqFlx.f90
@@ -20,339 +20,340 @@
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
-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: 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
- ! input: snow properties and parameters
- mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow), & ! intent(in): volumetric ice content at the start of the time step (-)
- Fcapil => mpar_data%var(iLookPARAM%Fcapil)%dat(1), & ! intent(in): capillary retention as a fraction of the total pore volume (-)
- k_snow => mpar_data%var(iLookPARAM%k_snow)%dat(1), & ! intent(in): hydraulic conductivity of snow (m s-1), 0.0055 = approx. 20 m/hr, from UEB
- mw_exp => mpar_data%var(iLookPARAM%mw_exp)%dat(1), & ! intent(in): exponent for meltwater flow (-)
- ! 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
- 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 computJacDAE_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.
+ ! 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
+ ! 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
+ ! input: snow properties and parameters
+ mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow), & ! intent(in): volumetric ice content at the start of the time step (-)
+ Fcapil => mpar_data%var(iLookPARAM%Fcapil)%dat(1), & ! intent(in): capillary retention as a fraction of the total pore volume (-)
+ k_snow => mpar_data%var(iLookPARAM%k_snow)%dat(1), & ! intent(in): hydraulic conductivity of snow (m s-1), 0.0055 = approx. 20 m/hr, from UEB
+ mw_exp => mpar_data%var(iLookPARAM%mw_exp)%dat(1), & ! intent(in): exponent for meltwater flow (-)
+ ! 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
- 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
- ! input: snow properties and parameters
- Fcapil => mpar_data%var(iLookPARAM%Fcapil)%dat(1), & ! intent(in): capillary retention as a fraction of the total pore volume (-)
- k_snow => mpar_data%var(iLookPARAM%k_snow)%dat(1), & ! intent(in): hydraulic conductivity of snow (m s-1), 0.0055 = approx. 20 m/hr, from UEB
- mw_exp => mpar_data%var(iLookPARAM%mw_exp)%dat(1), & ! intent(in): exponent for meltwater flow (-)
- ! 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
- 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.
+
+ ! 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
+ 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
+ ! input: snow properties and parameters
+ Fcapil => mpar_data%var(iLookPARAM%Fcapil)%dat(1), & ! intent(in): capillary retention as a fraction of the total pore volume (-)
+ k_snow => mpar_data%var(iLookPARAM%k_snow)%dat(1), & ! intent(in): hydraulic conductivity of snow (m s-1), 0.0055 = approx. 20 m/hr, from UEB
+ mw_exp => mpar_data%var(iLookPARAM%mw_exp)%dat(1), & ! intent(in): exponent for meltwater flow (-)
+ ! 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
- 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
+
+ ! 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
+ 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