var_derive.f90

! 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 var_derive_module

! data types
USE nrtype

! derived types to define the data structures
USE data_types,only:var_ilength    ! x%var(:)%dat (i4b)
USE data_types,only:var_dlength    ! x%var(:)%dat (dp)

! named variables for snow and soil
USE globalData,only:iname_snow     ! named variables for snow
USE globalData,only:iname_soil     ! named variables for soil

! named variables
USE globalData,only:data_step      ! time step of forcing data

! named variables
USE var_lookup,only:iLookPARAM,iLookINDEX,iLookPROG,iLookDIAG,iLookFLUX        ! HRU: named variables for structure elements
USE var_lookup,only:iLookBVAR,iLookBPAR                                        ! GRU: named variables for structure elements

! model decision structures
USE globalData,only:model_decisions        ! model decision structure
USE var_lookup,only:iLookDECISIONS         ! named variables for elements of the decision structure

! look-up values for the choice of the rooting profile
USE mDecisions_module,only: &
 powerLaw,                   & ! simple power-law rooting profile
 doubleExp                     ! the double exponential function of Xeng et al. (JHM 2001)

! look-up values for the choice of groundwater parameterization
USE mDecisions_module,only: &
 bigBucket,                  & ! a big bucket (lumped aquifer model)
 noExplicit                    ! no explicit groundwater parameterization

! look-up values for the choice of groundwater parameterization
USE mDecisions_module,only: &
 constant,                  & ! constant hydraulic conductivity with depth
 powerLaw_profile             ! power-law profile

! look-up values for the sub-grid routing method
USE mDecisions_module,only:      &
 timeDelay,&  ! time-delay histogram
 qInstant     ! instantaneous routing

! privacy
implicit none
private
public::calcHeight
public::rootDensty
public::satHydCond
public::fracFuture
public::v_shortcut
contains


 ! **********************************************************************************************************
 ! public subroutine calcHeight: compute snow height
 ! **********************************************************************************************************
 subroutine calcHeight(&
                       ! input/output: data structures
                       indx_data,   & ! intent(in): layer type
                       prog_data,   & ! intent(inout): model variables for a local HRU
                       ! output: error control
                       err,message)
 implicit none
 ! ----------------------------------------------------------------------------------
 ! dummy variables
 ! input/output: data structures
 type(var_ilength),intent(in)       :: indx_data      ! type of model layer
 type(var_dlength),intent(inout)    :: prog_data      ! model variables for a local HRU
 ! output: error control
 integer(i4b),intent(out)           :: err            ! error code
 character(*),intent(out)           :: message        ! error message
 ! local variables
 integer(i4b)                       :: iLayer         ! loop through layers
 integer(i4b)                       :: ixLower(1)     ! index of the lower bound
 ! ----------------------------------------------------------------------------------
 ! initialize error control
 err=0; message='calcHeight/'
 ! ----------------------------------------------------------------------------------
 ! associate variables in data structure
 associate(&
 ! associate the model index structures
 nLayers        => indx_data%var(iLookINDEX%nLayers)%dat(1),  &   ! total number of layers
 layerType      => indx_data%var(iLookINDEX%layerType)%dat,   &   ! layer type (iname_soil or iname_snow)
 ! associate the values in the model variable structures
 mLayerDepth    => prog_data%var(iLookPROG%mLayerDepth)%dat,  &   ! depth of the layer (m)
 mLayerHeight   => prog_data%var(iLookPROG%mLayerHeight)%dat, &   ! height of the layer mid-point (m)
 iLayerHeight   => prog_data%var(iLookPROG%iLayerHeight)%dat  &   ! height of the layer interface (m)
 ) ! end associate
 ! ----------------------------------------------------------------------------------

 ! initialize layer height as the top of the snowpack -- positive downward
 ixLower=lbound(iLayerHeight); if(ixLower(1) > 0)then; err=20; message=trim(message)//'unexpected lower bound for iLayerHeight'; return; endif
 iLayerHeight(0) = -sum(mLayerDepth, mask=layerType==iname_snow)

 ! loop through layers
 do iLayer=1,nLayers
  ! compute the height at the layer midpoint
  mLayerHeight(iLayer) = iLayerHeight(iLayer-1) + mLayerDepth(iLayer)/2._dp
  ! compute the height at layer interfaces
  iLayerHeight(iLayer) = iLayerHeight(iLayer-1) + mLayerDepth(iLayer)
 end do ! (looping through layers)

 !print*, 'layerType   = ',  layerType
 !print*, 'mLayerDepth = ',  mLayerDepth
 !print*, 'mLayerHeight = ', mLayerHeight
 !print*, 'iLayerHeight = ', iLayerHeight
 !print*, '************** '

 ! end association to variables in the data structure
 end associate

 end subroutine calcHeight


 ! **********************************************************************************************************
 ! public subroutine rootDensty: compute vertical distribution of root density
 ! **********************************************************************************************************
 subroutine rootDensty(mpar_data,indx_data,prog_data,diag_data,err,message)
 implicit none
 ! declare input variables
 type(var_dlength),intent(in)    :: mpar_data       ! data structure of model parameters for a local HRU
 type(var_ilength),intent(in)    :: indx_data       ! data structure of model indices for a local HRU
 type(var_dlength),intent(in)    :: prog_data       ! data structure of model prognostic (state) variables for a local HRU
 type(var_dlength),intent(inout) :: diag_data       ! data structure of model diagnostic variables for a local HRU
 ! declare output variables
 integer(i4b),intent(out)        :: err             ! error code
 character(*),intent(out)        :: message         ! error message
 ! declare local variables
 integer(i4b)                    :: iLayer          ! loop through layers
 real(dp)                        :: fracRootLower   ! fraction of the rooting depth at the lower interface
 real(dp)                        :: fracRootUpper   ! fraction of the rooting depth at the upper interface
 real(dp), parameter             :: rootTolerance = 0.05_dp ! tolerance for error in doubleExp rooting option
 real(dp)                        :: error           ! machine precision error in rooting distribution
 ! initialize error control
 err=0; message='rootDensty/'

 ! ----------------------------------------------------------------------------------
 ! associate variables in data structure
 associate(&
 ! associate the model decisions
 ixRootProfile         =>model_decisions(iLookDECISIONS%rootProfil)%iDecision,  & ! choice of the rooting profile
 ixGroundwater         =>model_decisions(iLookDECISIONS%groundwatr)%iDecision,  & ! choice of groundwater parameterization
 ! associate the values in the model parameter structures
 rootScaleFactor1      =>mpar_data%var(iLookPARAM%rootScaleFactor1)%dat(1),     & ! 1st scaling factor (m-1)
 rootScaleFactor2      =>mpar_data%var(iLookPARAM%rootScaleFactor2)%dat(1),     & ! 2nd scaling factor (m-1)
 rootingDepth          =>mpar_data%var(iLookPARAM%rootingDepth)%dat(1),         & ! rooting depth (m)
 rootDistExp           =>mpar_data%var(iLookPARAM%rootDistExp)%dat(1),          & ! root distribution exponent (-)
 ! associate the model index structures
 nSoil                 =>indx_data%var(iLookINDEX%nSoil)%dat(1),                & ! number of soil layers
 nSnow                 =>indx_data%var(iLookINDEX%nSnow)%dat(1),                & ! number of snow layers
 nLayers               =>indx_data%var(iLookINDEX%nLayers)%dat(1),              & ! total number of layers
 iLayerHeight          =>prog_data%var(iLookPROG%iLayerHeight)%dat,             & ! height of the layer interface (m)
 ! associate the values in the model variable structures
 scalarAquiferRootFrac =>diag_data%var(iLookDIAG%scalarAquiferRootFrac)%dat(1), & ! fraction of roots below the soil profile (in the aquifer)
 mLayerRootDensity     =>diag_data%var(iLookDIAG%mLayerRootDensity)%dat         & ! fraction of roots in each soil layer (-)
 ) ! end associate
 ! ----------------------------------------------------------------------------------

 !print*, 'nSnow   = ', nSnow
 !print*, 'nLayers = ', nLayers

 ! compute the fraction of roots in each soil layer
 do iLayer=nSnow+1,nLayers

  ! different options for the rooting profile
  select case(ixRootProfile)

   ! ** option 1: simple power-law profile
   case(powerLaw)
    if(iLayerHeight(iLayer-1)<rootingDepth)then
     ! compute the fraction of the rooting depth at the lower and upper interfaces
     if(iLayer==nSnow+1)then  ! height=0; avoid precision issues
      fracRootLower = 0._dp
     else
      fracRootLower = iLayerHeight(iLayer-1)/rootingDepth
     end if
     fracRootUpper = iLayerHeight(iLayer)/rootingDepth
     if(fracRootUpper>1._dp) fracRootUpper=1._dp
     ! compute the root density
     mLayerRootDensity(iLayer-nSnow) = fracRootUpper**rootDistExp - fracRootLower**rootDistExp
    else
     mLayerRootDensity(iLayer-nSnow) = 0._dp
    end if
    !write(*,'(a,10(f11.5,1x))') 'mLayerRootDensity(iLayer-nSnow), fracRootUpper, fracRootLower = ', &
    !                             mLayerRootDensity(iLayer-nSnow), fracRootUpper, fracRootLower

   ! ** option 2: double expoential profile of Zeng et al. (JHM 2001)
   case(doubleExp)
    ! compute the cumulative fraction of roots at the top and bottom of the layer
    fracRootLower = 1._dp - 0.5_dp*(exp(-iLayerHeight(iLayer-1)*rootScaleFactor1) + exp(-iLayerHeight(iLayer-1)*rootScaleFactor2) )
    fracRootUpper = 1._dp - 0.5_dp*(exp(-iLayerHeight(iLayer  )*rootScaleFactor1) + exp(-iLayerHeight(iLayer  )*rootScaleFactor2) )
    ! compute the root density
    mLayerRootDensity(iLayer-nSnow) = fracRootUpper - fracRootLower
    !write(*,'(a,10(f11.5,1x))') 'mLayerRootDensity(iLayer-nSnow), fracRootUpper, fracRootLower = ', &
    !                             mLayerRootDensity(iLayer-nSnow), fracRootUpper, fracRootLower

   ! ** check
   case default; err=20; message=trim(message)//'unable to identify option for rooting profile'; return

  end select

 end do  ! (looping thru layers)

 ! check that root density is within some reaosnable version of machine tolerance
 ! This is the case when root density is greater than 1. Can only happen with powerLaw option.
 error = sum(mLayerRootDensity) - 1._dp
 if (error > 2._dp*epsilon(rootingDepth)) then
  message=trim(message)//'problem with the root density calaculation'
  err=20; return
 else
  mLayerRootDensity = mLayerRootDensity - error/real(nSoil,kind(dp))
 end if

 ! compute fraction of roots in the aquifer
 if(sum(mLayerRootDensity) < 1._dp)then
  scalarAquiferRootFrac = 1._dp - sum(mLayerRootDensity)
 else
  scalarAquiferRootFrac = 0._dp
 end if

 ! check that roots in the aquifer are appropriate
 if ((ixGroundwater /= bigBucket).and.(scalarAquiferRootFrac > 2._dp*epsilon(rootingDepth)))then
  if(scalarAquiferRootFrac < rootTolerance) then
   mLayerRootDensity = mLayerRootDensity + scalarAquiferRootFrac/real(nSoil, kind(dp))
   scalarAquiferRootFrac = 0._dp
  else
   select case(ixRootProfile)
    case(powerLaw);  message=trim(message)//'roots in the aquifer only allowed for the big bucket gw parameterization: check that rooting depth < soil depth'
    case(doubleExp); message=trim(message)//'roots in the aquifer only allowed for the big bucket gw parameterization: increase soil depth to alow for exponential roots'
   end select
   err=10; return
  end if  ! if roots in the aquifer
 end if  ! if not the big bucket

 end associate

 end subroutine rootDensty


 ! **********************************************************************************************************
 ! public subroutine satHydCond: compute vertical profile of saturated hydraulic conductivity
 ! **********************************************************************************************************
 subroutine satHydCond(mpar_data,indx_data,prog_data,flux_data,err,message)
 implicit none
 ! declare input variables
 type(var_dlength),intent(in)    :: mpar_data           ! data structure of model parameters for a local HRU
 type(var_ilength),intent(in)    :: indx_data           ! data structure of model indices for a local HRU
 type(var_dlength),intent(in)    :: prog_data           ! data structure of model prognostic (state) variables for a local HRU
 type(var_dlength),intent(inout) :: flux_data           ! data structure of model fluxes for a local HRU
 ! declare output variables
 integer(i4b),intent(out)        :: err                 ! error code
 character(*),intent(out)        :: message             ! error message
 ! declare local variables
 integer(i4b)                    :: iLayer              ! loop through layers
 real(dp)                        :: ifcDepthScaleFactor ! depth scaling factor (layer interfaces)
 real(dp)                        :: midDepthScaleFactor ! depth scaling factor (layer midpoints)
 ! initialize error control
 err=0; message='satHydCond/'
 ! ----------------------------------------------------------------------------------
 ! associate variables in data structure
 associate(&
 ! associate the values in the parameter structures
 k_soil             => mpar_data%var(iLookPARAM%k_soil)%dat,            & ! saturated hydraulic conductivity at the compacted depth (m s-1)
 k_macropore        => mpar_data%var(iLookPARAM%k_macropore)%dat,       & ! saturated hydraulic conductivity at the compacted depth for macropores (m s-1)
 compactedDepth     => mpar_data%var(iLookPARAM%compactedDepth)%dat(1), & ! the depth at which k_soil reaches the compacted value given by CH78 (m)
 zScale_TOPMODEL    => mpar_data%var(iLookPARAM%zScale_TOPMODEL)%dat(1),& ! exponent for the TOPMODEL-ish baseflow parameterization (-)
 ! associate the model index structures
 nSnow              => indx_data%var(iLookINDEX%nSnow)%dat(1),          & ! number of snow layers
 nSoil              => indx_data%var(iLookINDEX%nSoil)%dat(1),          & ! number of soil layers
 nLayers            => indx_data%var(iLookINDEX%nLayers)%dat(1),        & ! total number of layers
 ! associate the coordinate variables
 mLayerHeight       => prog_data%var(iLookPROG%mLayerHeight)%dat,       & ! height at the mid-point of each layer (m)
 iLayerHeight       => prog_data%var(iLookPROG%iLayerHeight)%dat,       & ! height at the interface of each layer (m)
 ! associate the values in the model variable structures
 mLayerSatHydCondMP => flux_data%var(ilookFLUX%mlayersathydcondmp)%dat, & ! saturated hydraulic conductivity for macropores at the mid-point of each layer (m s-1)
 mLayerSatHydCond   => flux_data%var(ilookFLUX%mlayersathydcond)%dat,   & ! saturated hydraulic conductivity at the mid-point of each layer (m s-1)
 iLayerSatHydCond   => flux_data%var(ilookFLUX%ilayersathydcond)%dat    & ! saturated hydraulic conductivity at the interface of each layer (m s-1)
 ) ! end associate
 ! ----------------------------------------------------------------------------------

 ! loop through soil layers
 ! NOTE: could do constant profile with the power-law profile with exponent=1, but keep constant profile decision for clarity
 do iLayer=nSnow,nLayers
  select case(model_decisions(iLookDECISIONS%hc_profile)%iDecision)

   ! constant hydraulic conductivity with depth
   case(constant)
    ! - conductivity at layer interfaces
    !   --> NOTE: Do we need a weighted average based on layer depth for interior layers?
    if(iLayer==nSnow)then
     iLayerSatHydCond(iLayer-nSnow) = k_soil(1)
    else
     if(iLayer==nLayers)then
      iLayerSatHydCond(iLayer-nSnow) = k_soil(nSoil)
     else
      iLayerSatHydCond(iLayer-nSnow)   = 0.5_dp * (k_soil(iLayer-nSnow) + k_soil(iLayer+1-nSnow) )
     endif
     ! - conductivity at layer midpoints
     mLayerSatHydCond(iLayer-nSnow)   = k_soil(iLayer-nSnow)
     mLayerSatHydCondMP(iLayer-nSnow) = k_macropore(iLayer-nSnow)
    end if ! if iLayer>nSnow

   ! power-law profile
   case(powerLaw_profile)
    ! - conductivity at layer interfaces
    !   --> NOTE: Do we need a weighted average based on layer depth for interior layers?
    
    if(compactedDepth/iLayerHeight(nLayers) /= 1._dp) then    ! avoid divide by zero
     ifcDepthScaleFactor = ( (1._dp - iLayerHeight(iLayer)/iLayerHeight(nLayers))**(zScale_TOPMODEL - 1._dp) ) / &
                           ( (1._dp -       compactedDepth/iLayerHeight(nLayers))**(zScale_TOPMODEL - 1._dp) )
    else
     ifcDepthScaleFactor = 1.0_dp
    endif                           
    if(iLayer==nSnow)then
     iLayerSatHydCond(iLayer-nSnow) = k_soil(1) * ifcDepthScaleFactor
    else   ! if the mid-point of a layer
     if(iLayer==nLayers)then
      iLayerSatHydCond(iLayer-nSnow) = k_soil(nSoil) * ifcDepthScaleFactor
     else
      iLayerSatHydCond(iLayer-nSnow)   = 0.5_dp * (k_soil(iLayer-nSnow) + k_soil(iLayer+1-nSnow) ) * ifcDepthScaleFactor
     endif
     ! - conductivity at layer midpoints
     if(compactedDepth/iLayerHeight(nLayers) /= 1._dp) then    ! avoid divide by zero
      midDepthScaleFactor = ( (1._dp - mLayerHeight(iLayer)/iLayerHeight(nLayers))**(zScale_TOPMODEL - 1._dp) ) / &
                            ( (1._dp -       compactedDepth/iLayerHeight(nLayers))**(zScale_TOPMODEL - 1._dp) )
     else
      midDepthScaleFactor = 1.0_dp
     endif                            
     mLayerSatHydCond(iLayer-nSnow)   = k_soil(iLayer-nSnow)      * midDepthScaleFactor
     mLayerSatHydCondMP(iLayer-nSnow) = k_macropore(iLayer-nSnow) * midDepthScaleFactor
    end if  
 
    !print*, 'compactedDepth = ', compactedDepth
    !print*, 'k_macropore    = ', k_macropore
    !print*, 'mLayerHeight(iLayer) = ', mLayerHeight(iLayer)
    !print*, 'iLayerHeight(nLayers) = ', iLayerHeight(nLayers)
    !print*, 'iLayer, mLayerSatHydCondMP(iLayer-nSnow) = ', mLayerSatHydCondMP(iLayer-nSnow)

   ! error check (errors checked earlier also, so should not get here)
   case default
    message=trim(message)//"unknown hydraulic conductivity profile [option="//trim(model_decisions(iLookDECISIONS%hc_profile)%cDecision)//"]"
    err=10; return

  end select
  !if(iLayer > nSnow)& ! avoid layer 0
  ! write(*,'(a,1x,i4,1x,2(f11.5,1x,e20.10,1x))') 'satHydCond: ', iLayer, mLayerHeight(iLayer), mLayerSatHydCond(iLayer-nSnow), iLayerHeight(iLayer), iLayerSatHydCond(iLayer-nSnow)
 end do  ! looping through soil layers
 !print*, trim(model_decisions(iLookDECISIONS%hc_profile)%cDecision)
 !print*, 'k_soil, k_macropore, zScale_TOPMODEL = ', k_soil, k_macropore, zScale_TOPMODEL
 !pause ' in satHydCond'
 end associate

 end subroutine satHydCond


 ! **********************************************************************************************************
 ! public subroutine fracFuture: compute the fraction of runoff in future time steps
 ! **********************************************************************************************************
 subroutine fracFuture(bpar_data,bvar_data,err,message)
 ! external functions
 USE soil_utils_module,only:gammp                     ! compute the cumulative probabilty based on the Gamma distribution

 implicit none
 ! input variables
 real(dp),intent(in)             :: 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
 ! internal
 real(dp)                        :: 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(dp)                        :: aLambda                ! scale parameter in the Gamma distribution
 real(dp)                        :: tFuture                ! future time (end of step)
 real(dp)                        :: pSave                  ! cumulative probability at the start of the step
 real(dp)                        :: cumProb                ! cumulative probability at the end of the step
 real(dp)                        :: sumFrac                ! sum of runoff fractions in all steps
 real(dp),parameter              :: tolerFrac=0.01_dp      ! tolerance for missing fractional runoff by truncating histogram
 ! initialize error control
 err=0; message='fracFuture/'
 ! ----------------------------------------------------------------------------------
 ! associate variables in data structure
 associate(&
 ixRouting         => model_decisions(iLookDECISIONS%subRouting)%iDecision, & ! index for routing method
 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)
 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
 ! ----------------------------------------------------------------------------------

 ! define time step
 dt =  data_step ! length of the data step (s)

 ! identify number of points in the time-delay runoff variable (should be allocated match nTimeDelay)
 nTDH = size(runoffFuture)

 ! initialize runoffFuture (will be overwritten by initial conditions file values if present)
 runoffFuture(1:nTDH) = 0._dp

 ! select option for sub-grid routing
 select case(ixRouting)

  ! ** instantaneous routing
  case(qInstant)
   fractionFuture(1)      = 1._dp
   fractionFuture(2:nTDH) = 0._dp

  ! ** time delay histogram
  case(timeDelay)
   ! initialize
   pSave   = 0._dp ! cumulative probability at the start of the step
   aLambda = routingGammaShape / routingGammaScale
   if(routingGammaShape <= 0._dp .or. aLambda < 0._dp)then
    message=trim(message)//'bad arguments for the Gamma distribution'
    err=20; return
   end if
   ! loop through time steps and compute fraction of runoff in future steps
   do iFuture = 1,nTDH
    ! get weight for a given bin
    tFuture = real(iFuture, kind(dt))*dt                  ! future time (end of step)
    cumProb = gammp(routingGammaShape,aLambda*tFuture)    ! cumulative probability at the end of the step
    fractionFuture(iFuture) = max(0._dp, cumProb - pSave) ! fraction of runoff in the current step
    pSave   = cumProb                                     ! save the cumulative probability for use in the next step
    !write(*,'(a,1x,i4,1x,3(f20.10,1x))') trim(message), iFuture, tFuture, cumProb, fractionFuture(iFuture)
    ! set remaining bins to zero
    if(fractionFuture(iFuture) < tiny(dt))then
     fractionFuture(iFuture:nTDH) = 0._dp
     exit
    end if
   end do ! (looping through future time steps)

   ! check that we have enough bins
   sumFrac  = sum(fractionFuture)
   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
    message=trim(message)//'not enough bins for the time delay histogram -- fix hard-coded parameter in globalData.f90'
    err=20; return
   end if
   ! ensure the fraction sums to one
   fractionFuture = fractionFuture/sumFrac

  ! ** error checking
  case default; err=20; message=trim(message)//'cannot find option for sub-grid routing'; return

 end select ! (select option for sub-grid routing)

 end associate

 end subroutine fracFuture


 ! **********************************************************************************************************
 ! public subroutine v_shortcut: compute "short-cut" variables
 ! **********************************************************************************************************
 subroutine v_shortcut(mpar_data,diag_data,err,message)
 implicit none
 ! declare input variables
 type(var_dlength),intent(in)    :: mpar_data       ! data structure of model parameters for a local HRU
 type(var_dlength),intent(inout) :: diag_data       ! data structure of model variables for a local HRU
 ! declare output variables
 integer(i4b),intent(out)        :: err             ! error code
 character(*),intent(out)        :: message         ! error message
 ! initialize error control
 err=0; message='v_shortcut/'
 ! ----------------------------------------------------------------------------------
 ! associate variables in data structure
 associate(&
 ! associate values in the parameter structures
 vGn_n          =>mpar_data%var(iLookPARAM%vGn_n)%dat,                  & ! van Genutchen "n" parameter (-)
 vGn_m          =>diag_data%var(iLookDIAG%scalarVGn_m)%dat              & ! van Genutchen "m" parameter (-)
 ) ! end associate
 ! ----------------------------------------------------------------------------------

 ! compute the van Genutchen "m" parameter
 vGn_m = 1._dp - 1._dp/vGn_n
 end associate

 end subroutine v_shortcut


end module var_derive_module