Muskingum hydrologic routing in Microsoft Excel (UDF)

MuskyRouting.bas

Function Muskingum(timeSeries As Range, _
                   K As Double, X As Double, _
                   Optional Reaches As Integer = 1) As Variant
'Muskingum routing in Excel
'Returns a column array the same length as timeSeries (use array formula or spill)
'K (travel time) is in same units as the time step of the input timeSeries
'X must be between 0 and 0.5
'Reaches will be automatically adjusted to avoid negative coefficients if K is < 1 / (2*(1-X)) or K > 1/2X
Dim Coeffs(1 To 5) As Double

If X > 0.5 Or X < 0 Then
    MsgBox "Muskingum K must be between 0 and 0.5"
    Exit Function
End If

Kmin = 1 / (2 * (1 - X))
Kmax = 1 / (2 * X)
Ksub = K / Reaches

If Ksub < Kmin Then
    If Reaches > 1 Then
        Do While Reaches > 1 And Ksub < Kmin
            Reaches = Reaches - 1
            Ksub = K / Reaches
        Loop
    Else
        MsgBox "Muskingum routing requires K > 1/(2(1-x)) = " & Kmin
        Exit Function
    End If
End If

Do While Ksub > Kmax
    Reaches = Reaches + 1
    Ksub = K / Reaches
Loop

Denom = 2 * Ksub * (1 - X) + 1
CC = (Denom - 2) / Denom

Coeffs(1) = (1 - Ksub * X) / Denom
Coeffs(2) = Coeffs(1) * CC + (1 + Ksub * X) / Denom
Coeffs(3) = Coeffs(2) * CC
Coeffs(4) = Coeffs(3) * CC
Coeffs(5) = 1 - WorksheetFunction.Sum(Coeffs) 'Ensure Coeffs Sum to 1

TS = timeSeries.Value
TSOut = TS

Do While Reaches > 0
    For i = 1 To UBound(TS)
        v = 0#
        For c = 1 To 5
            TSindex = WorksheetFunction.Max(1, i - c + 1)
            v = v + TS(TSindex, 1) * Coeffs(c)
        Next c
        TSOut(i, 1) = v
    Next i
    Reaches = Reaches - 1
    TS = TSOut
Loop

Muskingum = TSOut
End Function

References:

HEC OPROUT (1982)

HEC-5 User Manual (1998)

HEC OPROUT, updated link