GwfSfrCommon.f90

Go to the documentation of this file.

  !> @brief Kinematic routing equation residual
    !!
    !! Method to calculate the kinematic-wave routing
    !! residual.
    !!
  !<
  function kinematic_residual(qa, qb, qc, qd, &
                              aa, ab, ac, ad, &
                              qsrc, length, weight, delt, &
                              courant)
    use kindmodule, only: dp
    use constantsmodule, only: dzero, dhalf, done
    ! -- return variable
    real(dp) :: kinematic_residual !< kinematic-wave residual
    ! -- dummy variables
    real(dp), intent(in) :: qa !< upstream flow at previous time
    real(dp), intent(in) :: qb !< downstream flow at previous time
    real(dp), intent(in) :: qc !< upstream flow
    real(dp), intent(in) :: qd !< downstream flow
    real(dp), intent(in) :: aa !< upstream area at previous time
    real(dp), intent(in) :: ab !< downstream area at previous time
    real(dp), intent(in) :: ac !< upstream area
    real(dp), intent(in) :: ad !< downstream area
    real(dp), intent(in) :: qsrc !< lateral flow term (L3T-1L-1)
    real(dp), intent(in) :: length !< reach length
    real(dp), intent(in) :: weight !< temporal weight
    real(dp), intent(in) :: delt !< time step length
    real(dp), intent(in) :: courant !< courant number
    ! --local variables
    real(dp) :: f11
    real(dp) :: f12
 
    if (courant > done) then
      f11 = (qd - qc) / length
      f12 = (ac - aa) / delt
    else
      f11 = (weight * (qd - qc) + (done - weight) * (qb - qa)) / length
      f12 = dhalf * ((ad - ab) + (ac - aa)) / delt
    end if
    kinematic_residual = f11 + f12 - qsrc
 
  end function kinematic_residual
 
  !> @brief Kinematic routing equation storage term
    !!
    !! Method to calculate the kinematic-wave routing
    !! storage term.
    !!
  !<
  function kinematic_storage(aa, ab, ac, ad, length, delt, courant)
    use kindmodule, only: dp
    use constantsmodule, only: dhalf
 
    real(dp) :: kinematic_storage !< kinematic-wave storage change
 
    real(dp), intent(in) :: aa !< upstream area at previous time
    real(dp), intent(in) :: ab !< downstream area at previous time
    real(dp), intent(in) :: ac !< upstream area
    real(dp), intent(in) :: ad !< downstream area
    real(dp), intent(in) :: length !< reach length
    real(dp), intent(in) :: delt !< time step length
    real(dp), intent(in) :: courant !< courant number
 
    if (courant > done) then
      kinematic_storage = (aa - ac) * length / delt
    else
      kinematic_storage = dhalf * ((ac + ad) - (aa + ab)) * length / delt
    end if
 
  end function kinematic_storage