diff --git a/build/source/driver/summa_setup.f90 b/build/source/driver/summa_setup.f90 index e406dfd4b..416807b76 100755 --- 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 index 8243d20a2..6b0f12e60 100755 --- 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 index 0a0eff43e..e47d5a360 100755 --- 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 index 73d7d93f7..1b88de275 100755 --- 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 index b658d7c0d..53597150f 100755 --- 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 index ba5a5aa16..2e6ee170c 100755 --- 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 index d0103d615..e779fdd91 100755 --- 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 index 88b5a5edd..1199421d7 100755 --- 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