diff --git a/build/source/dshare/globalData.f90 b/build/source/dshare/globalData.f90
index 205788e200250c50aa5686176dc792977edccae1..1e1cbf74b546f14a96753ee129038803b75079e6 100755
--- a/build/source/dshare/globalData.f90
+++ b/build/source/dshare/globalData.f90
@@ -315,7 +315,7 @@ MODULE globalData
! define fixed dimensions
integer(i4b),parameter,public :: nBand=2 ! number of spectral bands
- integer(i4b),parameter,public :: nTimeDelay=2000 ! number of hours in the time delay histogram (default: ~1 season = 24*365/4)
+ integer(i4b),parameter,public :: nTimeDelay=3000 ! number of hours in the time delay histogram (default: ~1 season = 24*365/4)
character(len=1024),public,save :: fname ! temporary filename
diff --git a/build/source/engine/computFlux.f90 b/build/source/engine/computFlux.f90
index c2be017999a60896dd330974dc1be1471e726344..9d443f370006e5729d674ef5d31a5e2230ba3df1 100755
--- a/build/source/engine/computFlux.f90
+++ b/build/source/engine/computFlux.f90
@@ -620,6 +620,10 @@ subroutine computFlux(&
! check the need to compute liquid water fluxes through snow
if(nSnowOnlyHyd>0)then
+ ! print*, "scalarThroughfallRain = ", scalarThroughfallRain
+ ! print*, "scalarCanopyLiqDrainage = ", scalarCanopyLiqDrainage
+ ! print*, "mLayerVolFracLiqTrial(1) =", mLayerVolFracLiqTrial(1)
+
! compute liquid fluxes through snow
call snowLiqFlx(&
! input: model control
@@ -655,6 +659,7 @@ subroutine computFlux(&
! compute drainage from the soil zone (needed for mass balance checks)
scalarSnowDrainage = iLayerLiqFluxSnow(nSnow)
+ ! print*, "scalarSnowDrainage = ", scalarSnowDrainage
! save bottom layer of snow derivatives
above_soilLiqFluxDeriv = iLayerLiqFluxSnowDeriv(nSnow) ! derivative in vertical liquid water flux at bottom snow layer interface
diff --git a/build/source/engine/coupled_em.f90 b/build/source/engine/coupled_em.f90
index 76f05fe29ad340c2d84126a4295cbb1d767b8141..c93d332cdbc9587a2f36164d51e9b218a66da601 100755
--- a/build/source/engine/coupled_em.f90
+++ b/build/source/engine/coupled_em.f90
@@ -253,6 +253,14 @@ subroutine coupled_em(&
! initialize error control
err=0; message="coupled_em/"
+ ! print*, "scalarCanopyWat"
+ ! print*, "mLayerVolFracWat"
+ ! print*, "mLayerMatricHead"
+ ! print*, ""
+ ! print*, ""
+ ! print*, ""
+ ! print*, ""
+
! check that the decision is supported
if(model_decisions(iLookDECISIONS%groundwatr)%iDecision==bigBucket .and. &
model_decisions(iLookDECISIONS%spatial_gw)%iDecision/=localColumn)then
@@ -1355,22 +1363,32 @@ subroutine coupled_em(&
newSWE = prog_data%var(iLookPROG%scalarSWE)%dat(1)
delSWE = newSWE - (oldSWE - sfcMeltPond)
massBalance = delSWE - (effSnowfall + effRainfall + averageSnowSublimation - averageSnowDrainage*iden_water)*data_step
- if(abs(massBalance) > 1.d-6)then
- print*, 'nSnow = ', nSnow
- print*, 'nSub = ', nSub
- write(*,'(a,1x,f20.10)') 'data_step = ', data_step
- write(*,'(a,1x,f20.10)') 'oldSWE = ', oldSWE
- write(*,'(a,1x,f20.10)') 'newSWE = ', newSWE
- write(*,'(a,1x,f20.10)') 'delSWE = ', delSWE
- write(*,'(a,1x,f20.10)') 'effRainfall = ', effRainfall*data_step
- write(*,'(a,1x,f20.10)') 'effSnowfall = ', effSnowfall*data_step
- write(*,'(a,1x,f20.10)') 'sublimation = ', averageSnowSublimation*data_step
- write(*,'(a,1x,f20.10)') 'snwDrainage = ', averageSnowDrainage*iden_water*data_step
- write(*,'(a,1x,f20.10)') 'sfcMeltPond = ', sfcMeltPond
- write(*,'(a,1x,f20.10)') 'massBalance = ', massBalance
- message=trim(message)//'SWE does not balance'
- print*,message
- err=20; return
+ ! print*, "effSnowfall = ", effSnowfall
+ ! print*, "effRainfall = ", effRainfall
+ ! print*, "averageSnowSublimation = ", averageSnowSublimation
+ ! print*, "averageSnowDrainage = ", averageSnowDrainage
+ ! print*, "iden_water = ", iden_water
+ ! print*, "newSWE = ", newSWE
+ ! print*, "delSWE = ", delSWE
+ ! print*, "massBalance = ", massBalance
+
+
+ if(abs(massBalance) > absConvTol_liquid*iden_water*10._dp)then
+ print*, 'nSnow = ', nSnow
+ print*, 'nSub = ', nSub
+ write(*,'(a,1x,f20.10)') 'data_step = ', data_step
+ write(*,'(a,1x,f20.10)') 'oldSWE = ', oldSWE
+ write(*,'(a,1x,f20.10)') 'newSWE = ', newSWE
+ write(*,'(a,1x,f20.10)') 'delSWE = ', delSWE
+ write(*,'(a,1x,f20.10)') 'effRainfall = ', effRainfall*data_step
+ write(*,'(a,1x,f20.10)') 'effSnowfall = ', effSnowfall*data_step
+ write(*,'(a,1x,f20.10)') 'sublimation = ', averageSnowSublimation*data_step
+ write(*,'(a,1x,f20.10)') 'snwDrainage = ', averageSnowDrainage*iden_water*data_step
+ write(*,'(a,1x,f20.10)') 'sfcMeltPond = ', sfcMeltPond
+ write(*,'(a,1x,f20.10)') 'massBalance = ', massBalance
+ message=trim(message)//'SWE does not balance'
+ print*,message
+ err=20; return
endif ! if failed mass balance check
endif ! if snow layers exist
diff --git a/build/source/engine/opSplittin.f90 b/build/source/engine/opSplittin.f90
index c22a32ecf3b0ac97a5c2aeba998fa71a05f4d736..90393b6254f160b20dc66a900cb2067c56b8eeff 100755
--- a/build/source/engine/opSplittin.f90
+++ b/build/source/engine/opSplittin.f90
@@ -365,6 +365,18 @@ subroutine opSplittin(&
! ---------------------------------------------------------------------------------------
! initialize error control
err=0; message="opSplittin/"
+ ! print*, "BEFORE******"
+ ! print*, "scalarCanairTemp = ", scalarCanairTemp
+ ! print*, "scalarCanopyTemp = ", scalarCanopyTemp
+ ! print*, "scalarCanopyIce = ", scalarCanopyIce
+ ! print*, "scalarCanopyLiq = ", scalarCanopyLiq
+ ! print*, "scalarCanopyWat = ", scalarCanopyWat
+ ! print*, "mLayerTemp = ", mLayerTemp(1)
+ ! print*, "mLayerVolFracIce = ", mLayerVolFracIce(1)
+ ! print*, "mLayerVolFracLiq = ", mLayerVolFracLiq(1)
+ ! print*, "mLayerVolFracWat = ", mLayerVolFracWat(1)
+ ! print*, "mLayerMatricHead = ", mLayerMatricHead(1)
+ ! print*, "mLayerMatricHeadLiq = ", mLayerMatricHeadLiq(1)
! we just solve the fully coupled problem by ida
select case(model_decisions(iLookDECISIONS%diffEqSolv)%iDecision)
@@ -919,9 +931,9 @@ subroutine opSplittin(&
if(err>0) return
endif ! (check for errors)
- ! print*, trim(message)//'after varSubstep: scalarSnowDrainage = ', flux_data%var(iLookFLUX%scalarSnowDrainage)%dat
- ! print*, trim(message)//'after varSubstep: iLayerLiqFluxSnow = ', flux_data%var(iLookFLUX%iLayerLiqFluxSnow)%dat
- ! print*, trim(message)//'after varSubstep: iLayerLiqFluxSoil = ', flux_data%var(iLookFLUX%iLayerLiqFluxSoil)%dat
+ ! print*, trim(message)//'after varSubstep: scalarSnowDrainage = ', flux_data%var(iLookFLUX%scalarSnowDrainage)%dat
+ ! print*, trim(message)//'after varSubstep: iLayerLiqFluxSnow = ', flux_data%var(iLookFLUX%iLayerLiqFluxSnow)%dat
+ ! print*, trim(message)//'after varSubstep: iLayerLiqFluxSoil = ', flux_data%var(iLookFLUX%iLayerLiqFluxSoil)%dat
! check
!if(ixSolution==scalar)then
@@ -1069,20 +1081,9 @@ subroutine opSplittin(&
! ==========================================================================================================================================
! success = exit the coupling loop
- ! terminate DO loop early if fullyCoupled returns a solution,
- ! so that the loop does not proceed to ixCoupling = stateTypeSplit
+
if(ixCoupling==fullyCoupled .and. .not. failure) exit coupling
-
- ! if we reach stateTypeSplit, terminating the DO loop here is cleaner
- ! than letting the loop complete, because in the latter case the coupling
- ! loop will end with ixCoupling = nCoupling+1 = 3 (a FORTRAN loop
- ! increments the index variable at the end of each iteration and stops
- ! the loop if the index > specified stop value). Variable ixCoupling is
- ! used for error reporting in coupled_em.f90 in the balance checks and
- ! we thus need to make sure ixCoupling is not incremented to be larger
- ! than nCoupling.
- if(ixCoupling==stateTypeSplit .and. .not. failure) exit coupling
-
+
end do coupling ! coupling method
! check that all state variables were updated
@@ -1105,9 +1106,24 @@ subroutine opSplittin(&
if(ixCoupling/=fullyCoupled .or. nSubsteps>1) dtMultiplier=0.5_dp
! compute the melt in each snow and soil layer
- if(nSnow>0) mLayerMeltFreeze( 1:nSnow ) = -(mLayerVolFracIce( 1:nSnow ) - mLayerVolFracIceInit( 1:nSnow ))*iden_ice
- mLayerMeltFreeze(nSnow+1:nLayers) = -(mLayerVolFracIce(nSnow+1:nLayers) - mLayerVolFracIceInit(nSnow+1:nLayers))*iden_water
+ if(nSnow>0) then
+ diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat(1:nSnow) = -( mLayerVolFracIce(1:nSnow) - mLayerVolFracIceInit(1:nSnow) ) * iden_ice
+ diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat(nSnow+1:nLayers) = -(mLayerVolFracIce(nSnow+1:nLayers) - mLayerVolFracIceInit(nSnow+1:nLayers))*iden_water
+ endif
+
+ ! print*, "After******"
+ ! print*, "scalarCanairTemp = ", scalarCanairTemp
+ ! print*, "scalarCanopyTemp = ", scalarCanopyTemp
+ ! print*, "scalarCanopyIce = ", scalarCanopyIce
+ ! print*, "scalarCanopyLiq = ", scalarCanopyLiq
+ ! print*, "scalarCanopyWat = ", scalarCanopyWat
+ ! print*, "mLayerTemp = ", mLayerTemp(1)
+ ! print*, "mLayerVolFracIce = ", mLayerVolFracIce(1)
+ ! print*, "mLayerVolFracLiq = ", mLayerVolFracLiq(1)
+ ! print*, "mLayerVolFracWat = ", mLayerVolFracWat(1)
+ ! print*, "mLayerMatricHead = ", mLayerMatricHead(1)
+ ! print*, "mLayerMatricHeadLiq = ", mLayerMatricHeadLiq(1)
! end associate statements
end associate globalVars
diff --git a/build/source/engine/snowLiqFlx.f90 b/build/source/engine/snowLiqFlx.f90
old mode 100755
new mode 100644
index 53b4fb29a127dbf68ccc95dffe65986000de8a88..6778561f876a97ebb584ba918c37c4373be3a9da
--- a/build/source/engine/snowLiqFlx.f90
+++ b/build/source/engine/snowLiqFlx.f90
@@ -35,14 +35,15 @@ 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(:) (dp)
-USE data_types,only:var_dlength ! x%var(:)%dat (dp)
+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
@@ -75,18 +76,18 @@ contains
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(dp),intent(in) :: scalarThroughfallRain ! computed throughfall rate (kg m-2 s-1)
- real(dp),intent(in) :: scalarCanopyLiqDrainage ! computed drainage of liquid water (kg m-2 s-1)
+ 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(dp),intent(in) :: mLayerVolFracLiqTrial(:) ! trial value of volumetric fraction of liquid water at the current iteration (-)
+ 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(dp),intent(inout) :: iLayerLiqFluxSnow(0:) ! vertical liquid water flux at layer interfaces (m s-1)
- real(dp),intent(inout) :: iLayerLiqFluxSnowDeriv(0:) ! derivative in vertical liquid water flux at layer interfaces (m s-1)
+ 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
@@ -96,12 +97,12 @@ contains
integer(i4b) :: iLayer ! layer index
integer(i4b) :: ixTop ! top layer in subroutine call
integer(i4b) :: ixBot ! bottom layer in subroutine call
- real(dp) :: multResid ! multiplier for the residual water content (-)
- real(dp),parameter :: residThrs=550._dp ! ice density threshold to reduce residual liquid water content (kg m-3)
- real(dp),parameter :: residScal=10._dp ! scaling factor for residual liquid water content reduction factor (kg m-3)
- real(dp),parameter :: maxVolIceContent=0.7_dp ! maximum volumetric ice content to store water (-)
- real(dp) :: availCap ! available storage capacity [0,1] (-)
- real(dp) :: relSaturn ! relative saturation [0,1] (-)
+ 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(&
@@ -128,7 +129,7 @@ contains
end if
! check the meltwater exponent is >=1
- if(mw_exp<1._dp)then; err=20; message=trim(message)//'meltwater exponent < 1'; return; end if
+ 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
@@ -159,16 +160,16 @@ contains
! define the liquid flux at the upper boundary (m s-1)
iLayerLiqFluxSnow(0) = (scalarThroughfallRain + scalarCanopyLiqDrainage)/iden_water
- iLayerLiqFluxSnowDeriv(0) = 0._dp
+ 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._dp / ( 1._dp + exp( (mLayerVolFracIce(iLayer)*iden_ice - residThrs) / residScal) )
+ multResid = 1._rkind / ( 1._rkind + exp( (mLayerVolFracIce(iLayer)*iden_ice - residThrs) / residScal) )
! compute the pore space (-)
- mLayerPoreSpace(iLayer) = 1._dp - mLayerVolFracIce(iLayer)
+ 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)
@@ -182,14 +183,14 @@ contains
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._dp)
+ 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._dp
- iLayerLiqFluxSnowDeriv(iLayer) = 0._dp
+ iLayerLiqFluxSnow(iLayer) = 0._rkind
+ iLayerLiqFluxSnowDeriv(iLayer) = 0._rkind
endif ! storage above residual content
end do ! loop through snow layers
@@ -199,4 +200,159 @@ contains
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.
+ 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
diff --git a/build/source/engine/snowLiqFlx_old.f90 b/build/source/engine/snowLiqFlx_old.f90
new file mode 100755
index 0000000000000000000000000000000000000000..53b4fb29a127dbf68ccc95dffe65986000de8a88
--- /dev/null
+++ b/build/source/engine/snowLiqFlx_old.f90
@@ -0,0 +1,202 @@
+! SUMMA - Structure for Unifying Multiple Modeling Alternatives
+! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
+!
+! This file is part of SUMMA
+!
+! For more information see: http://www.ral.ucar.edu/projects/summa
+!
+! This program is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+module 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(:) (dp)
+USE data_types,only:var_dlength ! x%var(:)%dat (dp)
+USE data_types,only:var_ilength ! x%var(:)%dat (i4b)
+
+! privacy
+implicit none
+private
+public::snowLiqFlx
+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(dp),intent(in) :: scalarThroughfallRain ! computed throughfall rate (kg m-2 s-1)
+ real(dp),intent(in) :: scalarCanopyLiqDrainage ! computed drainage of liquid water (kg m-2 s-1)
+ ! input: model state vector
+ real(dp),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(dp),intent(inout) :: iLayerLiqFluxSnow(0:) ! vertical liquid water flux at layer interfaces (m s-1)
+ real(dp),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(dp) :: multResid ! multiplier for the residual water content (-)
+ real(dp),parameter :: residThrs=550._dp ! ice density threshold to reduce residual liquid water content (kg m-3)
+ real(dp),parameter :: residScal=10._dp ! scaling factor for residual liquid water content reduction factor (kg m-3)
+ real(dp),parameter :: maxVolIceContent=0.7_dp ! maximum volumetric ice content to store water (-)
+ real(dp) :: availCap ! available storage capacity [0,1] (-)
+ real(dp) :: 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._dp)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._dp
+
+ ! 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._dp / ( 1._dp + exp( (mLayerVolFracIce(iLayer)*iden_ice - residThrs) / residScal) )
+ ! compute the pore space (-)
+ mLayerPoreSpace(iLayer) = 1._dp - 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._dp)
+ 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._dp
+ iLayerLiqFluxSnowDeriv(iLayer) = 0._dp
+ 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
+
+
+end module snowLiqFlx_module
diff --git a/build/source/engine/var_derive.f90 b/build/source/engine/var_derive.f90
index 0432a41cadbe4be71faeb68ac354ebc5a6ed7b36..7e3eed052a7330c6ae2d2989ea795b01d6aa90a6 100755
--- a/build/source/engine/var_derive.f90
+++ b/build/source/engine/var_derive.f90
@@ -458,6 +458,10 @@ contains
if(abs(1._dp - sumFrac) > tolerFrac)then
write(*,*) 'fraction of basin runoff histogram being accounted for by time delay vector is ', sumFrac
write(*,*) 'this is less than allowed by tolerFrac = ', tolerFrac
+ write(*,*) 'Solutions:'
+ write(*,*) ' (1) Check that the values of routingGammaShape and routingGammaScale are appropriate (and fix if necessary); or'
+ write(*,*) ' (2) Increase the hard coded parameter nTimeDelay in globalData.f90 (currently nTimeDelay is set to ', nTDH, ')'
+ write(*,*) ' -- note that nTimeDelay defines the number of time steps in the time delay histogram'
message=trim(message)//'not enough bins for the time delay histogram -- fix hard-coded parameter in globalData.f90'
err=20; return
end if
diff --git a/utils/laugh_tests/celia1990/verification_data/runinfo.txt b/utils/laugh_tests/celia1990/verification_data/runinfo.txt
index 2e294c3db47c5c2e65b97fb0ec50b4f6c9bf0658..4f7e767fdcbbaa6083334e085da26299ed886975 100644
--- a/utils/laugh_tests/celia1990/verification_data/runinfo.txt
+++ b/utils/laugh_tests/celia1990/verification_data/runinfo.txt
@@ -1 +1 @@
- Run start time on system: ccyy=2022 - mm=09 - dd=08 - hh=18 - mi=34 - ss=44.626
+ Run start time on system: ccyy=2022 - mm=09 - dd=08 - hh=18 - mi=40 - ss=09.307
diff --git a/utils/laugh_tests/colbeck1976/actors_out.txt b/utils/laugh_tests/colbeck1976/actors_out.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3fe835e32740d3f165c31c89b51a1a9ed761ca11
--- /dev/null
+++ b/utils/laugh_tests/colbeck1976/actors_out.txt
@@ -0,0 +1,349 @@
+SETTINGS FOR SUMMA_ACTOR
+Output Structure Size = 250
+Max GRUs Per Job = 250
+
+SETTINGS FOR JOB_ACTOR
+File Manager Path = /Summa-Actors/utils/laugh_tests/colbeck1976/settings/summa_fileManager_colbeck1976-exp1.txt
+output_csv = false
+Job Actor Initalized
+
+----------File_Access_Actor Started----------
+Initalizing Output Structure
+GRU Actor Has Started
+
+SETTINGS FOR HRU_ACTOR
+Print Output = true
+Print Output every 1 timesteps
+
+All Forcing Files Loaded
+1 - Timestep = 1
+1 - Timestep = 2
+1 - Timestep = 3
+1 - Timestep = 4
+1 - Timestep = 5
+1 - Timestep = 6
+1 - Timestep = 7
+1 - Timestep = 8
+1 - Timestep = 9
+1 - Timestep = 10
+1 - Timestep = 11
+ /Summa-Actors/utils/laugh_tests/colbeck1976/settings/summa_defineModelOutput.txt
+ HERE
+
+INFO: aspect not found in the input attribute file, continuing ...
+
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 0.0000000000000000
+ iden_water = 1000.0000000000000
+ newSWE = 300.60000000000002
+ delSWE = 0.60000000000002274
+ massBalance = 2.2759572004815709E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 4.9587017904807871E-017
+ iden_water = 1000.0000000000000
+ newSWE = 301.19999999999783
+ delSWE = 0.59999999999780584
+ massBalance = 7.8104189782379763E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 3.9772892257868010E-016
+ iden_water = 1000.0000000000000
+ newSWE = 301.79999999997375
+ delSWE = 0.59999999997592113
+ massBalance = -2.1516122217235534E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.3458248922576075E-015
+ iden_water = 1000.0000000000000
+ newSWE = 302.39999999989305
+ delSWE = 0.59999999991930508
+ massBalance = 5.4511950509095186E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 3.1983853384144842E-015
+ iden_water = 1000.0000000000000
+ newSWE = 302.99999999970117
+ delSWE = 0.59999999980811936
+ massBalance = 2.2537527399890678E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 6.2630410898092297E-015
+ iden_water = 1000.0000000000000
+ newSWE = 303.59999999932541
+ delSWE = 0.59999999962423090
+ massBalance = 1.3322676295501878E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.0850552546761924E-014
+ iden_water = 1000.0000000000000
+ newSWE = 304.19999999867446
+ delSWE = 0.59999999934905190
+ massBalance = 8.5043083686286991E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.7274818811323277E-014
+ iden_water = 1000.0000000000000
+ newSWE = 304.79999999763805
+ delSWE = 0.59999999896359668
+ massBalance = 8.5820239803524601E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 2.5852886556788660E-014
+ iden_water = 1000.0000000000000
+ newSWE = 305.39999999608659
+ delSWE = 0.59999999844853846
+ massBalance = -2.8843594179761567E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 3.6904959605757223E-014
+ iden_water = 1000.0000000000000
+ newSWE = 305.99999999387285
+ delSWE = 0.59999999778625579
+ massBalance = 5.5333515547317802E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 5.0754385602336477E-014
+ iden_water = 1000.0000000000000
+ newSWE = 306.59999999082720
+ delSWE = 0.59999999695435235
+ massBalance = -3.8458125573015423E-013
+1 - Timestep = 12
+1 - Timestep = 13
+1 - Timestep = 14
+1 - Timestep = 15
+1 - Timestep = 16
+1 - Timestep = 17
+1 - Timestep = 18
+1 - Timestep = 19
+1 - Timestep = 20
+1 - Timestep = 21
+1 - Timestep = 22
+1 - Timestep = 23
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 6.7727690522882357E-014
+ iden_water = 1000.0000000000000
+ newSWE = 307.19999998676326
+ delSWE = 0.59999999593605935
+ massBalance = -2.7933211299568939E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 8.8154616957988339E-014
+ iden_water = 1000.0000000000000
+ newSWE = 307.79999998147423
+ delSWE = 0.59999999471096999
+ massBalance = 2.4702462297909733E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.1236808114754007E-013
+ iden_water = 1000.0000000000000
+ newSWE = 308.39999997473200
+ delSWE = 0.59999999325776798
+ massBalance = -1.4710455076283324E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.4070416063307484E-013
+ iden_water = 1000.0000000000000
+ newSWE = 308.99999996628952
+ delSWE = 0.59999999155752448
+ massBalance = -2.2593038551121936E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.7350215945493872E-013
+ iden_water = 1000.0000000000000
+ newSWE = 309.59999995587907
+ delSWE = 0.59999998958954848
+ massBalance = -3.2196467714129540E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 2.1110465794128700E-013
+ iden_water = 1000.0000000000000
+ newSWE = 310.19999994321324
+ delSWE = 0.59999998733417215
+ massBalance = 4.5163872641751368E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 2.5385736837767799E-013
+ iden_water = 1000.0000000000000
+ newSWE = 310.79999992798213
+ delSWE = 0.59999998476888550
+ massBalance = 3.2762681456688370E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 3.0210926319812266E-013
+ iden_water = 1000.0000000000000
+ newSWE = 311.39999990985564
+ delSWE = 0.59999998187350911
+ massBalance = 6.4837024638109142E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 3.5621257403774553E-013
+ iden_water = 1000.0000000000000
+ newSWE = 311.99999988848282
+ delSWE = 0.59999997862718146
+ massBalance = -6.4170890823334048E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 4.1652276488339124E-013
+ iden_water = 1000.0000000000000
+ newSWE = 312.59999986349135
+ delSWE = 0.59999997500852942
+ massBalance = -1.0469403122215226E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 4.8339858993693743E-013
+ iden_water = 1000.0000000000000
+ newSWE = 313.19999983448764
+ delSWE = 0.59999997099629354
+ massBalance = 2.0894397323445446E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 5.5720191270652321E-013
+ iden_water = 1000.0000000000000
+ newSWE = 313.79999980105532
+ delSWE = 0.59999996656767962
+ massBalance = -2.0561330416057899E-013
+1 - Timestep = 24
+1 - Timestep = 25
+1 - Timestep = 26
+1 - Timestep = 27
+1 - Timestep = 28
+1 - Timestep = 29
+1 - Timestep = 30
+1 - Timestep = 31
+1 - Timestep = 32
+1 - Timestep = 33
+1 - Timestep = 34
+1 - Timestep = 35
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 6.3829798560848913E-013
+ iden_water = 1000.0000000000000
+ newSWE = 314.39999976275755
+ delSWE = 0.59999996170222403
+ massBalance = 1.0313971898767704E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 7.2705534112198951E-013
+ iden_water = 1000.0000000000000
+ newSWE = 314.99999971913411
+ delSWE = 0.59999995637656411
+ massBalance = -1.1546319456101628E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 8.2384578888474908E-013
+ iden_water = 1000.0000000000000
+ newSWE = 315.59999966970321
+ delSWE = 0.59999995056909938
+ massBalance = -1.5332179970073412E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 9.2904443946227717E-013
+ iden_water = 1000.0000000000000
+ newSWE = 316.19999961396070
+ delSWE = 0.59999994425749037
+ massBalance = 1.5665246877460959E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.0430297137078835E-012
+ iden_water = 1000.0000000000000
+ newSWE = 316.79999955137873
+ delSWE = 0.59999993741803337
+ massBalance = -1.8385293287792592E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.1661833323718264E-012
+ iden_water = 1000.0000000000000
+ newSWE = 317.39999948140763
+ delSWE = 0.59999993002890051
+ massBalance = -9.9587005308876542E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.2988901893932967E-012
+ iden_water = 1000.0000000000000
+ newSWE = 317.99999940347448
+ delSWE = 0.59999992206684283
+ massBalance = 2.5424107263916085E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.4415385941891714E-012
+ iden_water = 1000.0000000000000
+ newSWE = 318.59999931698201
+ delSWE = 0.59999991350753135
+ massBalance = -1.5309975509580909E-013
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.5945202341461499E-012
+ iden_water = 1000.0000000000000
+ newSWE = 319.19999922131086
+ delSWE = 0.59999990432885397
+ massBalance = 6.8056671409522096E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.7582301232143246E-012
+ iden_water = 1000.0000000000000
+ newSWE = 319.79999911581712
+ delSWE = 0.59999989450625435
+ massBalance = 6.1728400169158704E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 1.9330666721022761E-012
+ iden_water = 1000.0000000000000
+ newSWE = 320.39999899983303
+ delSWE = 0.59999988401591509
+ massBalance = -8.4598994476436928E-014
+ effSnowfall = 0.0000000000000000
+ effRainfall = 1.0000000000000000E-002
+ averageSnowSublimation = 0.0000000000000000
+ averageSnowDrainage = 2.1194312945133812E-012
+ iden_water = 1000.0000000000000
+ newSWE = 320.99999887266722
+ delSWE = 0.59999987283418932
+ massBalance = 6.6946448384896939E-014
+1 - Timestep = 36
+1 - Timestep = 37
+1 - Timestep = 38
+1 - Timestep = 39
+1 - Timestep = 40
+1 - Timestep = 41
+1 - Timestep = 42
+1 - Timestep = 43
+1 - Timestep = 44
+Error: RunPhysics - HRU = 1 - indxGRU = 1 - refGRU = 1 - Timestep = 44
+*** unexpected message [id: 10, name: user.scheduled-actor]: message(run_failure(), 11@38DDCD1CBB03D559767271D57A7B827F2487A98D#36110, 1, 20)
diff --git a/utils/laugh_tests/colbeck1976/run_test_summa.sh b/utils/laugh_tests/colbeck1976/run_test_summa.sh
index 4b56489e34e4b2eb0b717c8bb33a4acb9c2a0173..00d3c26faa980899982615ee7c639ca3dacb7708 100755
--- a/utils/laugh_tests/colbeck1976/run_test_summa.sh
+++ b/utils/laugh_tests/colbeck1976/run_test_summa.sh
@@ -1,5 +1,5 @@
#! /bin/bash
-/SUMMA/bin/summa_sundials.exe -g 1 1 -m /Summa-Actors/utils/laugh_tests/colbeck1976/settings/summa_fileManager_verify_colbeck1976-exp1.txt
+/SUMMA/bin/summa_sundials.exe -g 1 1 -m /Summa-Actors/utils/laugh_tests/colbeck1976/settings/summa_fileManager_verify_colbeck1979-exp1.txt
/SUMMA/bin/summa_sundials.exe -g 1 1 -m /Summa-Actors/utils/laugh_tests/colbeck1976/settings/summa_fileManager_verify_colbeck1976-exp2.txt
/SUMMA/bin/summa_sundials.exe -g 1 1 -m /Summa-Actors/utils/laugh_tests/colbeck1976/settings/summa_fileManager_verify_colbeck1976-exp3.txt
\ No newline at end of file
diff --git a/utils/laugh_tests/colbeck1976/summa_out.txt b/utils/laugh_tests/colbeck1976/summa_out.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp1_G1-1_timestep.nc b/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp1_G1-1_timestep.nc
index 9e6ae158fabce58a11b81c45cffc294b752b523c..772be806da8b0d04af8fed0ed001a661552e6606 100644
Binary files a/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp1_G1-1_timestep.nc and b/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp1_G1-1_timestep.nc differ
diff --git a/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp2_G1-1_timestep.nc b/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp2_G1-1_timestep.nc
index 27e99c3dc1da709bf1e3e307406e1804486f4feb..8752d3307f1284e7a2e875f0458895b4e388ab0d 100644
Binary files a/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp2_G1-1_timestep.nc and b/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp2_G1-1_timestep.nc differ
diff --git a/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp3_G1-1_timestep.nc b/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp3_G1-1_timestep.nc
index 96eef89531f15d3d7649199e9d5dceedb46a74ae..c563189156bf58e2071bf9d4c33adbee91bd6cda 100644
Binary files a/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp3_G1-1_timestep.nc and b/utils/laugh_tests/colbeck1976/verification_data/colbeck1976-exp3_G1-1_timestep.nc differ
diff --git a/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt b/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt
index 9ec6f2ad3e6b7c139f46cd38ed5365d905e58f52..4370b2e1c2243faa068074cccc4a2c33e4c1637d 100644
--- a/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt
+++ b/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt
@@ -1 +1 @@
- Run start time on system: ccyy=2022 - mm=09 - dd=08 - hh=15 - mi=27 - ss=10.981
+ Run start time on system: ccyy=2022 - mm=09 - dd=08 - hh=22 - mi=30 - ss=05.308
diff --git a/utils/laugh_tests/colbeck1976/verify_colbeck.py b/utils/laugh_tests/colbeck1976/verify_colbeck.py
index 0394fb583cad070f25de0435eda4477de119be88..e690e9ed1a9cfd933eaa0f8f922a99d1f177a6fe 100644
--- a/utils/laugh_tests/colbeck1976/verify_colbeck.py
+++ b/utils/laugh_tests/colbeck1976/verify_colbeck.py
@@ -80,13 +80,13 @@ mLayerDepth = "mLayerDepth"
output_variables = [scalarRainfall, scalarSnowfall, scalarRainPlusMelt, mLayerVolFracLiq, \
mLayerVolFracIce, iLayerNrgFlux, iLayerHeight, mLayerDepth]
-verified_data_path = Path("./verification_data/colbeck1976-exp1_G1-1_timestep.nc")
-data_to_compare_path = Path("./output/colbeck1976-exp1GRU1-1_timestep.nc")
-verify(verified_data_path, data_to_compare_path, output_variables, numHRU)
+# verified_data_path = Path("./verification_data/colbeck1976-exp1_G1-1_timestep.nc")
+# data_to_compare_path = Path("./output/colbeck1976-exp1GRU1-1_timestep.nc")
+# verify(verified_data_path, data_to_compare_path, output_variables, numHRU)
-verified_data_path = Path("./verification_data/colbeck1976-exp2_G1-1_timestep.nc")
-data_to_compare_path = Path("./output/colbeck1976-exp2GRU1-1_timestep.nc")
-verify(verified_data_path, data_to_compare_path, output_variables, numHRU)
+# verified_data_path = Path("./verification_data/colbeck1976-exp2_G1-1_timestep.nc")
+# data_to_compare_path = Path("./output/colbeck1976-exp2GRU1-1_timestep.nc")
+# verify(verified_data_path, data_to_compare_path, output_variables, numHRU)
verified_data_path = Path("./verification_data/colbeck1976-exp3_G1-1_timestep.nc")
data_to_compare_path = Path("./output/colbeck1976-exp3GRU1-1_timestep.nc")