Solve Muskingum-Cunge-Todini equation. When parameters k and x are passed, C1, C2, and C3 are computed using the passed values. When x = 0, linear reservoir model. When x = 0.5 and k = dt translation only
References:
Todini, E. (2007) A mass conservative and water storage consistent variable parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659.
Ponce, V.M., (1994) Engineering Hydrology: Principles and Practices, Prentice Hall
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer(kind=short), | intent(in) | :: | dt |
time increment (s) |
||
real(kind=float), | intent(in) | :: | dx |
reach length (m) |
||
real(kind=float), | intent(in) | :: | Qin |
input discharge at current time (m3/s) |
||
real(kind=float), | intent(in) | :: | Pin |
input and output at previous time (m3/s) |
||
real(kind=float), | intent(in) | :: | Pout |
input and output at previous time (m3/s) |
||
real(kind=float), | intent(in) | :: | Qlat |
lateral flow at current and previous time (m3/s) |
||
real(kind=float), | intent(in) | :: | Plat |
lateral flow at current and previous time (m3/s) |
||
real(kind=float), | intent(in) | :: | B0 |
bottom width, = 0 for triangular section (m) |
||
real(kind=float), | intent(in) | :: | alpha |
angle formed by dykes over a horizontal plane (deg) |
||
real(kind=float), | intent(in) | :: | slope |
bed slope (m/m) |
||
real(kind=float), | intent(in) | :: | n |
manning roughness coefficient (s m^(-1/3)) |
||
real(kind=float), | intent(out) | :: | Qout |
output discharge at current time |
||
real(kind=float), | intent(out) | :: | depth | |||
real(kind=float), | intent(out) | :: | width | |||
real(kind=float), | intent(in), | optional | :: | k |
flood wave travel time used to compute C1, C2, and C3 (s) |
|
real(kind=float), | intent(in), | optional | :: | x |
used to compute C1, C2, and C3 (-) |
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
real(kind=float), | public | :: | B | ||||
real(kind=float), | public | :: | C1 | ||||
real(kind=float), | public | :: | C2 | ||||
real(kind=float), | public | :: | C3 | ||||
real(kind=float), | public | :: | C4 | ||||
real(kind=float), | public | :: | Cstar |
dimensionless Courant number at time t+deltat |
|||
real(kind=float), | public | :: | Dstar |
Cell Reynolds number at time t+deltat |
|||
real(kind=float), | public | :: | PCstar |
dimensionless Courant number at time t |
|||
real(kind=float), | public | :: | PDstar |
Cell Reynolds number at time t |
|||
real(kind=float), | public | :: | Pref |
reference discharge at time t |
|||
real(kind=float), | public | :: | Qoutp | ||||
real(kind=float), | public | :: | Qref |
reference discharge at time t+deltat |
|||
real(kind=float), | public | :: | aveWetArea | ||||
real(kind=float), | public | :: | betaf |
correction factor |
|||
real(kind=float), | public | :: | cel |
kinematic celerity (m/s) |
|||
integer, | public | :: | i | ||||
integer, | public | :: | maxiter | = | 20 | ||
real(kind=float), | public | :: | normal_depth |
normal flow depth (m) |
|||
real(kind=float), | public | :: | storage | ||||
real(kind=float), | public | :: | tol | = | 0.000001 |
SUBROUTINE MuskingumCungeTodini & ! ( dt, dx, Qin, Pin, Pout, Qlat, Plat, B0, alpha, slope, n, Qout, depth, width, k, x ) !arguments with intent in: INTEGER (KIND = short), INTENT (IN) :: dt !!time increment (s) REAL (KIND = float), INTENT (IN) :: dx !!reach length (m) REAL (KIND = float), INTENT (IN) :: Qin !!input discharge at current time (m3/s) REAL (KIND = float), INTENT (IN) :: Pin, Pout !!input and output at previous time (m3/s) REAL (KIND = float), INTENT (IN) :: Qlat, Plat !!lateral flow at current and previous time (m3/s) REAL (KIND = float), INTENT (IN) :: B0 !!bottom width, = 0 for triangular section (m) REAL (KIND = float), INTENT (IN) :: alpha !!angle formed by dykes over a horizontal plane (deg) REAL (KIND = float), INTENT (IN) :: slope !!bed slope (m/m) REAL (KIND = float), INTENT (IN) :: n !!manning roughness coefficient (s m^(-1/3)) REAL (KIND = float), OPTIONAL, INTENT(IN) :: k !!flood wave travel time used to compute C1, C2, and C3 (s) REAL (KIND = float), OPTIONAL, INTENT(IN) :: x !!used to compute C1, C2, and C3 (-) !arguments with intent out: REAL (KIND = float), INTENT (OUT) :: Qout !!output discharge at current time REAL (KIND = float), INTENT (OUT) :: depth !average water depth REAL (KIND = float), INTENT (OUT) :: width !wetted width !local declarations REAL (KIND = float) :: Cstar !!dimensionless Courant number at time t+deltat REAL (KIND = float) :: Dstar !!Cell Reynolds number at time t+deltat REAL (KIND = float) :: PCstar !!dimensionless Courant number at time t REAL (KIND = float) :: PDstar !!Cell Reynolds number at time t REAL (KIND = float) :: cel !!kinematic celerity (m/s) REAL (KIND = float) :: normal_depth !!normal flow depth (m) REAL (KIND = float) :: Qref !!reference discharge at time t+deltat REAL (KIND = float) :: Pref !!reference discharge at time t REAL (KIND = float) :: betaf !!correction factor REAL (KIND = float) :: B ! topwidth REAL (KIND = float) :: C1, C2, C3, C4 REAL (KIND = float) :: tol = 0.000001 !tolerance for checking convergence REAL (KIND = float) :: Qoutp !Qout of the previous iteration step REAL (KIND = float) :: storage !storage at time t + deltat REAL (KIND = float) :: aveWetArea !average wetted area in the reach INTEGER :: maxiter = 20 !maximum number of iterations INTEGER :: i !-------------------------------end of declaration----------------------------- IF (PRESENT (x) .AND. PRESENT (k) ) THEN C1 = (dt-2.*k*x) / (2.*k*(1.-x) + dt) C2 = (dt+2.*k*x) / (2.*k*(1.-x) + dt) C3 = (2.*k*(1.-x) - dt) / (2.*k*(1.-x) + dt) C4 = (2.*dt) / (2.*k*(1.-x) + dt) Qout = C1 * Qin + C2 * Pin + C3 * Pout + C4 * 0.5 * (Qlat + Plat) !Compute the storage storage = k * x * Qin + k * (1. - x) * Qout IF (storage > 0.) THEN !Compute the average wetted area aveWetArea = storage / dx !Compute average water depth depth = DepthFromArea (aveWetArea, B0, alpha) !Compute width width = TopWidth ( depth, B0, alpha ) ELSE depth = 0. width = 0. END IF RETURN !no need to go on END IF !compute reference discharge Pref = (Pin + Pout) / 2. IF ( (Qin > 0. .OR. Qlat > 0. ) .AND. Pref == 0.) THEN !prevent division by zero Pref = 0.001 END IF IF (Pref > 0.) THEN !compute normal flow depth normal_depth = NormalDepth ( Pref, B0, alpha, slope, n ) !compute topwidth B = TopWidth ( normal_depth, B0, alpha ) !compute celerity cel = Celerity ( normal_depth, B0, alpha, slope, n ) !compute beta betaf = Beta ( normal_depth, B0, alpha, slope, n ) !compute dimensionless Courant number PCstar = cel * dt / (dx * betaf) !Compute cell Reynolds number IF ( (B * slope * cel * betaf * dx) > 0.) THEN PDstar = Pref / (B * slope * cel * betaf * dx) ELSE PDstar = 0. END IF END IF !first guess estimate Qout = Pout + (Qin + Qlat - Pin) IF (Qout < 0.) THEN Qout = 0. END IF IF (Qin > 0. .OR. Qlat > 0.) THEN DO i = 1, maxiter Qoutp = Qout !save previous iteration value !compute reference discharge Qref = (Qin + Qout) / 2. !compute normal flow depth normal_depth = NormalDepth ( Qref, B0, alpha, slope, n ) !compute topwidth B = TopWidth ( normal_depth, B0, alpha ) !compute celerity cel = Celerity ( normal_depth, B0, alpha, slope, n ) !compute beta betaf = Beta ( normal_depth, B0, alpha, slope, n ) !compute dimensionless Courant number IF (dx * betaf > 0.) THEN Cstar = cel * dt / (dx * betaf) ELSE Cstar = 0. END IF !Compute cell Reynolds number IF ( (B * slope * cel * betaf * dx) > 0.) THEN Dstar = Qref / (B * slope * cel * betaf * dx) ELSE Dstar = 0. END IF !M-C coefficients C1 = (-1. + Cstar + Dstar) / (1. + Cstar + Dstar) C2 = (1. + PCstar - PDstar) / (1. + Cstar + Dstar) * Cstar / PCstar C3 = (1. - PCstar + PDstar) / (1. + Cstar + Dstar) * Cstar / PCstar C4 = 2. * Cstar / (1. + Cstar + Dstar) Qout = C1 * Qin + C2 * Pin + C3 * Pout + C4 * 0.5 * (Qlat + Plat) !IF (Qout < 0.) THEN ! Qout = 0. !END IF !check if other iteration is needed IF (ABS (Qout - Qoutp) / Qoutp <= tol) EXIT END DO !last check IF (Qout < 0.) then Qout = 0. END IF !Compute the storage storage = (1. - Dstar) * dt * Qin / (2. * Cstar) + & (1. + Dstar) * dt * Qout / (2. * Cstar) IF (storage > 0.) THEN !Compute the average wetted area aveWetArea = storage / dx !Compute average water depth depth = DepthFromArea (aveWetArea, B0, alpha) !Compute width width = TopWidth ( depth, B0, alpha ) ELSE depth = 0. width = 0. END IF ELSE !M-C coefficients C1 = 1. / 3. C2 = 1. / 3. C3 = 1. / 3. Qout = C1 * Qin + C2 * Pin + C3 * Pout depth = 0. END IF RETURN END SUBROUTINE MuskingumCungeTodini