MuskingumCungeTodini Subroutine

public subroutine MuskingumCungeTodini(dt, dx, Qin, Pin, Pout, Qlat, Plat, B0, alpha, slope, n, Qout, depth, width, k, x)

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

Arguments

Type IntentOptional 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 (-)


Variables

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

Source Code

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