CH-Earth · martynpclark · Apr 19, 2022
diff --git a/build/source/driver/summa_setup.f90 b/build/source/driver/summa_setup.f90
@@ -255,12 +255,6 @@ subroutine summa_paramSetup(summa1_struc, err, message)
  ! loop through GRUs
  do iGRU=1,nGRU
 
-  ! calculate the fraction of runoff in future time steps
-  call fracFuture(bparStruct%gru(iGRU)%var,    &  ! vector of basin-average model parameters
-                  bvarStruct%gru(iGRU),        &  ! data structure of basin-average variables
-                  err,cmessage)                   ! error control
-  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
-
   ! loop through local HRUs
   do iHRU=1,gru_struc(iGRU)%hruCount
 
@@ -315,6 +309,12 @@ subroutine summa_paramSetup(summa1_struc, err, message)
   end do
   end associate
 
+  ! calculate the fraction of runoff in future time steps
+  call fracFuture(bparStruct%gru(iGRU)%var,    &  ! vector of basin-average model parameters
+                  bvarStruct%gru(iGRU),        &  ! data structure of basin-average variables
+                  err,cmessage)                   ! error control
+  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+
  end do ! GRU
 
  ! identify the end of the initialization

diff --git a/build/source/dshare/get_ixname.f90 b/build/source/dshare/get_ixname.f90
@@ -859,6 +859,7 @@ function get_ixbpar(varName)
   case('basin__aquiferScaleFactor'); get_ixbpar = iLookBPAR%basin__aquiferScaleFactor ! scaling factor for aquifer storage in the big bucket (m)
   case('basin__aquiferBaseflowExp'); get_ixbpar = iLookBPAR%basin__aquiferBaseflowExp ! baseflow exponent for the big bucket (-)
   ! sub-grid routing
+  case('routingVelocity'          ); get_ixbpar = iLookBPAR%routingVelocity           ! velocity parameter in Gamma distribution used for sub-grid routing (-)
   case('routingGammaShape'        ); get_ixbpar = iLookBPAR%routingGammaShape         ! shape parameter in Gamma distribution used for sub-grid routing (-)
   case('routingGammaScale'        ); get_ixbpar = iLookBPAR%routingGammaScale         ! scale parameter in Gamma distribution used for sub-grid routing (s)
   ! get to here if cannot find the variable

diff --git a/build/source/dshare/popMetadat.f90 b/build/source/dshare/popMetadat.f90
@@ -291,6 +291,7 @@ subroutine popMetadat(err,message)
  bpar_meta(iLookBPAR%basin__aquiferHydCond)           = var_info('basin__aquiferHydCond'    , 'hydraulic conductivity of the aquifer'                          , 'm s-1', get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  bpar_meta(iLookBPAR%basin__aquiferScaleFactor)       = var_info('basin__aquiferScaleFactor', 'scaling factor for aquifer storage in the big bucket'           , 'm'    , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  bpar_meta(iLookBPAR%basin__aquiferBaseflowExp)       = var_info('basin__aquiferBaseflowExp', 'baseflow exponent for the big bucket'                           , '-'    , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
+ bpar_meta(iLookBPAR%routingVelocity)                 = var_info('routingVelocity'          , 'constant velocity parameter used for sub-grid routing'          , 'm s-1', get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  bpar_meta(iLookBPAR%routingGammaShape)               = var_info('routingGammaShape'        , 'shape parameter in Gamma distribution used for sub-grid routing', '-'    , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  bpar_meta(iLookBPAR%routingGammaScale)               = var_info('routingGammaScale'        , 'scale parameter in Gamma distribution used for sub-grid routing', 's'    , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
 

diff --git a/build/source/dshare/var_lookup.f90 b/build/source/dshare/var_lookup.f90
@@ -688,6 +688,7 @@ MODULE var_lookup
   integer(i4b)    :: basin__aquiferScaleFactor  = integerMissing ! scaling factor for aquifer storage in the big bucket (m)
   integer(i4b)    :: basin__aquiferBaseflowExp  = integerMissing ! baseflow exponent for the big bucket (-)
   ! within-grid routing
+  integer(i4b)    :: routingVelocity            = integerMissing ! constant velocity parameter used for sub-grid routing (-)
   integer(i4b)    :: routingGammaShape          = integerMissing ! shape parameter in Gamma distribution used for sub-grid routing (-)
   integer(i4b)    :: routingGammaScale          = integerMissing ! scale parameter in Gamma distribution used for sub-grid routing (s)
  endtype iLook_bpar
@@ -839,7 +840,7 @@ MODULE var_lookup
                                                                          51, 52, 53, 54, 55, 56, 57, 58, 59, 60)
 
  ! named variables: basin-average parameters
- type(iLook_bpar),    public,parameter :: iLookBPAR     =ilook_bpar    (  1,  2,  3,  4,  5)
+ type(iLook_bpar),    public,parameter :: iLookBPAR     =ilook_bpar    (  1,  2,  3,  4,  5,  6)
 
  ! named variables: basin-average variables
  type(iLook_bvar),    public,parameter :: iLookBVAR     =ilook_bvar    (  1,  2,  3,  4,  5,  6,  7,  8,  9, 10,&

diff --git a/build/source/engine/mDecisions.f90 b/build/source/engine/mDecisions.f90
@@ -136,7 +136,8 @@ module mDecisions_module
 integer(i4b),parameter,public :: localColumn          = 291    ! separate groundwater representation in each local soil column
 integer(i4b),parameter,public :: singleBasin          = 292    ! single groundwater store over the entire basin
 ! look-up values for the choice of sub-grid routing method
-integer(i4b),parameter,public :: timeDelay            = 301    ! time-delay histogram
+integer(i4b),parameter,public :: constVelocity        = 300    ! time-delay histogram based on constant velocity
+integer(i4b),parameter,public :: timeDelay            = 301    ! time-delay histogram computed using a gamma distribution with shape and scale parameters
 integer(i4b),parameter,public :: qInstant             = 302    ! instantaneous routing
 ! look-up values for the choice of new snow density method
 integer(i4b),parameter,public :: constDens            = 311    ! Constant new snow density
@@ -603,7 +604,8 @@ subroutine mDecisions(err,message)
 
  ! choice of routing method
  select case(trim(model_decisions(iLookDECISIONS%subRouting)%cDecision))
-  case('timeDlay'); model_decisions(iLookDECISIONS%subRouting)%iDecision = timeDelay           ! time-delay histogram
+  case('constVel'); model_decisions(iLookDECISIONS%subRouting)%iDecision = constVelocity       ! time-delay histogram computed based on constant velocity
+  case('timeDlay'); model_decisions(iLookDECISIONS%subRouting)%iDecision = timeDelay           ! time-delay histogram computed using a gamma distribution with shape and scale parameters
   case('qInstant'); model_decisions(iLookDECISIONS%subRouting)%iDecision = qInstant            ! instantaneous routing
   case default
    err=10; message=trim(message)//"unknown option for sub-grid routing [option="//trim(model_decisions(iLookDECISIONS%subRouting)%cDecision)//"]"; return

diff --git a/build/source/engine/qTimeDelay.f90 b/build/source/engine/qTimeDelay.f90
@@ -25,8 +25,9 @@ module qTimeDelay_module
 
 ! look-up values for the sub-grid routing method
 USE mDecisions_module,only:      &
- timeDelay,&  ! time-delay histogram
- qInstant     ! instantaneous routing
+ constVelocity,&  ! time-delay histogram computed based on constant velocity
+ timeDelay,&      ! time-delay histogram computed using a gamma distribution with shape and scale parameter
+ qInstant         ! instantaneous routing
 
 implicit none
 private
@@ -74,7 +75,7 @@ subroutine qOverland(&
    averageRoutedRunoff = averageInstantRunoff
 
   ! ** time delay histogram
-  case(timeDelay)
+  case(timeDelay,constVelocity)
    ! identify number of points in the time-delay histogram
    nTDH = size(qFuture)
    ! place a fraction of runoff in future steps

diff --git a/build/source/engine/read_pinit.f90 b/build/source/engine/read_pinit.f90
@@ -21,10 +21,10 @@
 module read_pinit_module
 USE nrtype
 ! check for when model decisions are undefined
-USE mDecisions_module,only: unDefined
+USE mDecisions_module,only: unDefined, constVelocity
 USE globalData,only:model_decisions
 USE globalData,only:realMissing
-USE var_lookup,only:iLookDECISIONS,iLookPARAM
+USE var_lookup,only:iLookDECISIONS,iLookPARAM,iLookBPAR
 implicit none
 private
 public::read_pinit
@@ -34,7 +34,7 @@ module read_pinit_module
  ! ************************************************************************************************
  ! public subroutine read_pinit: read default model parameter values and constraints
  ! ************************************************************************************************
- subroutine read_pinit(filenm,isLocal,mpar_meta,parFallback,err,message)
+ subroutine read_pinit(filenm,isHRU,mpar_meta,parFallback,err,message)
  ! used to read metadata on the forcing data file
  USE summaFileManager,only:SETTINGS_PATH   ! path for input parameter and other configuration files
  USE ascii_util_module,only:file_open      ! open ascii file
@@ -46,7 +46,7 @@ subroutine read_pinit(filenm,isLocal,mpar_meta,parFallback,err,message)
  implicit none
  ! define input
  character(*),intent(in)                :: filenm         ! name of file containing default values and constraints of model parameters
- logical(lgt),intent(in)                :: isLocal        ! .true. if the file describes local column parameters
+ logical(lgt),intent(in)                :: isHRU          ! .true. if the file describes HRU parameters
  type(var_info),intent(in)              :: mpar_meta(:)   ! metadata for model parameters
  ! define output
  type(par_info),intent(out)             :: parFallback(:) ! default values and constraints of model parameters
@@ -111,7 +111,7 @@ subroutine read_pinit(filenm,isLocal,mpar_meta,parFallback,err,message)
   read(temp,trim(ffmt),iostat=err) varname, dLim, parTemp%default_val, dLim, parTemp%lower_limit, dLim, parTemp%upper_limit
   if (err/=0) then; err=30; message=trim(message)//"errorReadLine"; return; end if
   ! (identify the index of the variable in the data structure)
-  if(isLocal)then
+  if(isHRU)then
    ivar = get_ixParam(trim(varname))
   else
    ivar = get_ixBpar(trim(varname))
@@ -129,24 +129,32 @@ subroutine read_pinit(filenm,isLocal,mpar_meta,parFallback,err,message)
    err=40; message=trim(message)//"variable in parameter file not present in data structure [var="//trim(varname)//"]"; return
   end if
  end do  ! (looping through lines in the file)
+ ! populate HRU parameters that were not included in the original control files
+ if(isHRU)then
+  if(model_decisions(iLookDECISIONS%cIntercept)%iDecision == unDefined)then
+   parFallback(iLookPARAM%canopyWettingFactor)%default_val = 1._rkind             ! maximum wetted fraction of the canopy (-)
+   parFallback(iLookPARAM%canopyWettingExp)%default_val    = 0.666666667_rkind    ! exponent in canopy wetting function (-)
+  end if
+ ! populate GRU parameters that were not included in the original control files
+ else
+  if(parFallback(iLookBPAR%routingVelocity)%default_val < 0.99_rkind*realMissing)then
+   parFallback(iLookBPAR%routingVelocity)%default_val = 1._dp
+   if(model_decisions(iLookDECISIONS%subRouting)%iDecision == constVelocity)then
+    print*, trim(message)//'WARNING: '// &
+                           'We noticed that you have implemented the constant velocity parameterization for sub-grid routing '// &
+                           'but you do not have the "routingVelocity" parameter in your basinPARAM file. We set the constant '// &
+                           'velocity parameter to 1.0. We consider this to be an immense act of philanthropy. You are welcome.'
+   endif  ! if using the constant velocity routing
+  endif  ! if the parameter is not populated
+ endif  ! if dealing with GRU parameters
  ! check we have populated all variables
  ! NOTE: ultimately need a need a parameter dictionary to ensure that the parameters used are populated
- if(.not.backwardsCompatible)then  ! if we add new variables in future versions of the code, then some may be missing in the input file
-  if(any(parFallback(:)%default_val < 0.99_rkind*realMissing))then
-   do ivar=1,size(parFallback)
-    if(parFallback(ivar)%default_val < 0.99_rkind*realMissing)then
-     err=40; message=trim(message)//"variableNonexistent[var="//trim(mpar_meta(ivar)%varname)//"]"; return
-    end if
-   end do
-  end if
- ! populate parameters that were not included in the original control files
- else ! (need backwards compatibility)
-  if(isLocal)then
-   if(model_decisions(iLookDECISIONS%cIntercept)%iDecision == unDefined)then
-    parFallback(iLookPARAM%canopyWettingFactor)%default_val = 1._rkind             ! maximum wetted fraction of the canopy (-)
-    parFallback(iLookPARAM%canopyWettingExp)%default_val    = 0.666666667_rkind    ! exponent in canopy wetting function (-)
+ if(any(parFallback(:)%default_val < 0.99_rkind*realMissing))then
+  do ivar=1,size(parFallback)
+   if(parFallback(ivar)%default_val < 0.99_rkind*realMissing)then
+    err=40; message=trim(message)//"variableNonexistent[var="//trim(mpar_meta(ivar)%varname)//"]"; return
    end if
-  end if
+  end do
  end if
  ! close file unit
  close(unt)

diff --git a/build/source/engine/var_derive.f90 b/build/source/engine/var_derive.f90
@@ -59,8 +59,9 @@ module var_derive_module
 
 ! look-up values for the sub-grid routing method
 USE mDecisions_module,only:      &
- timeDelay,&  ! time-delay histogram
- qInstant     ! instantaneous routing
+ constVelocity,&  ! time-delay histogram computed based on constant velocity
+ timeDelay,&      ! time-delay histogram computed using a gamma distribution with shape and scale parameter
+ qInstant         ! instantaneous routing
 
 ! privacy
 implicit none
@@ -384,15 +385,16 @@ subroutine fracFuture(bpar_data,bvar_data,err,message)
 
  implicit none
  ! input variables
- real(rkind),intent(in)             :: bpar_data(:)           ! vector of basin-average model parameters
+ real(rkind),intent(inout)          :: bpar_data(:)           ! vector of basin-average model parameters
  ! output variables
- type(var_dlength),intent(inout) :: bvar_data              ! data structure of basin-average model variables
- integer(i4b),intent(out)        :: err                    ! error code
- character(*),intent(out)        :: message                ! error message
+ type(var_dlength),intent(inout)    :: bvar_data              ! data structure of basin-average model variables
+ integer(i4b),intent(out)           :: err                    ! error code
+ character(*),intent(out)           :: message                ! error message
  ! internal
  real(rkind)                        :: dt                     ! data time step (s)
- integer(i4b)                    :: nTDH                   ! number of points in the time-delay histogram
- integer(i4b)                    :: iFuture                ! index in time delay histogram
+ real(rkind)                        :: chanLength             ! GRU channel length (m)
+ integer(i4b)                       :: nTDH                   ! number of points in the time-delay histogram
+ integer(i4b)                       :: iFuture                ! index in time delay histogram
  real(rkind)                        :: tFuture                ! future time (end of step)
  real(rkind)                        :: pSave                  ! cumulative probability at the start of the step
  real(rkind)                        :: cumProb                ! cumulative probability at the end of the step
@@ -404,8 +406,10 @@ subroutine fracFuture(bpar_data,bvar_data,err,message)
  ! associate variables in data structure
  associate(&
  ixRouting         => model_decisions(iLookDECISIONS%subRouting)%iDecision, & ! index for routing method
+ routingVelocity   => bpar_data(iLookBPAR%routingVelocity),                 & ! velocity parameter in Gamma distribution used for sub-grid routing (-)
  routingGammaShape => bpar_data(iLookBPAR%routingGammaShape),               & ! shape parameter in Gamma distribution used for sub-grid routing (-)
  routingGammaScale => bpar_data(iLookBPAR%routingGammaScale),               & ! scale parameter in Gamma distribution used for sub-grid routing (s)
+ totalArea         => bvar_data%var(iLookBVAR%basin__TotalArea)%dat(1),     & ! total basin area (m2)
  runoffFuture      => bvar_data%var(iLookBVAR%routingRunoffFuture)%dat,     & ! runoff in future time steps (m s-1)
  fractionFuture    => bvar_data%var(iLookBVAR%routingFractionFuture)%dat    & ! fraction of runoff in future time steps (-)
  ) ! end associate
@@ -414,6 +418,12 @@ subroutine fracFuture(bpar_data,bvar_data,err,message)
  ! define time step
  dt =  data_step ! length of the data step (s)
 
+ ! check if we need to compute the scale parameter for the Gamma distribution
+ if(ixRouting == constVelocity)then
+  chanLength        = sqrt(totalArea/PI_D)        ! channel length (m)
+  routingGammaScale = chanLength/routingVelocity  ! scale parameter (s)
+ endif  ! if need to compute the shape and scale parameters
+
  ! identify number of points in the time-delay runoff variable (should be allocated match nTimeDelay)
  nTDH = size(runoffFuture)
 
@@ -429,7 +439,7 @@ subroutine fracFuture(bpar_data,bvar_data,err,message)
    fractionFuture(2:nTDH) = 0._rkind
 
   ! ** time delay histogram
-  case(timeDelay)
+  case(timeDelay,constVelocity)
    ! initialize
    pSave   = 0._rkind ! cumulative probability at the start of the step
    if(routingGammaShape <= 0._rkind .or. routingGammaScale <= 0._rkind)then