RoutingModels.f90 Source File
Implement models for discharge routing
Source Code
!! Implement models for discharge routing !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.4 - 25th April 2025 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 03/Oct/2016 | Original code | ! | 1.1 | 05/May/2023 | bug corrected in RungeKutta third order | ! | 1.2 | 28/Nov/2023 | manage reservoir when level is high | ! | 1.3 | 12/Dec/2023 | ecological flow is defined at daily scale | ! | 1.4 | 24/Apr/2025 | added 'noleap' when computing doy in DiverFlow | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description ! This module includes models for discharge routing ! MODULE RoutingModels ! Modules used: USE DataTypeSizes, ONLY : & ! Imported Parameters: float, short USE Units, ONLY : & ! Imported parameters: degToRad USE Loglib, ONLY : & !Imported routines: Catch USE Chronos, ONLY : & !Imported definitions: DateTime , & !Imported routines: DayOfYear USE TableLib, ONLY : & !Imported definitions: Table, & !Imported routines: TableGetValue USE StringManipulation, ONLY : & !Imported routines: ToString USE Utilities, ONLY : & !Imported routines: LinearInterp USE Reservoirs, ONLY: & !Imported definition: Reservoir USE Diversions, ONLY : & !Imported definition: Diversion IMPLICIT NONE ! Global (i.e. public) Declarations: !Public routines PUBLIC :: MuskingumCungeTodini PUBLIC :: LevelPool PUBLIC :: DivertFlow !Local routines PRIVATE :: WettedArea PRIVATE :: WettedPerimeter PRIVATE :: TopWidth PRIVATE :: Discharge PRIVATE :: Velocity PRIVATE :: Celerity PRIVATE :: Beta PRIVATE :: NormalDepth PRIVATE :: DepthFromArea PRIVATE :: RungeKutta !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! 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 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 !============================================================================== !| Description: ! Solve Muskingum-Cunge equation. ! ! References: ! ! Ponce, V.M., (1994) Engineering Hydrology: Principles and Practices, ! Prentice Hall SUBROUTINE MuskingumCunge & ! ( 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) :: xx, kk !!used for computing C1, C2, and C3 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) :: B ! topwidth REAL (KIND = float) :: C1, C2, C3, C4 REAL (KIND = float) :: storage !storage at time t + deltat REAL (KIND = float) :: aveWetArea !average wetted area in the reach !-------------------------------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 Qref = (Qin + Pin + Pout) / 3. IF (Qref == 0. .AND. Qlat > 0.) THEN Qref = Qlat END IF IF (Qref > 0.) THEN !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 k and x kk = dx / cel xx = 0.5 * ( 1. - Qref / (B * cel * slope * dx) ) !compute weigths C1 = (dt-2.*kk*xx) / (2.*kk*(1.-xx)+dt) C2 = (dt+2.*kk*xx) / (2.*kk*(1.-xx)+dt) C3 = (2.*kk*(1.-xx)-dt) / (2.*kk*(1.-xx)+dt) C4 = 2.*dt / (2.*kk*(1.-xx)+dt) !compute Qout Qout = C1 * Qin + C2 * Pin + C3 * Pout + C4 * (Qlat+Plat) / 2. !Qout = C1 * Qin + C2 * Pin + C3 * Pout + C4 * Qlat !last check IF (Qout < 0.) then Qout = 0. END IF !compute the storage storage = kk *xx * Qin + kk * (1. - xx) * Qout !compute the stage 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 Qout = 0. depth = 0. END IF RETURN END SUBROUTINE MuskingumCunge !============================================================================== !| Description: ! Discharge routing through diversion. Discharge is split ! between river and diversion channel SUBROUTINE DivertFlow & ! ( time, div, Qin, Qout) !arguments with intent in: TYPE(DateTime), INTENT(IN) :: time REAL(KIND = float), INTENT(IN) :: Qin !!Input discharge at time t + dt (m3/s) !arguments with intent inout: TYPE (Diversion), INTENT(INOUT) :: div !arguments with intent out: REAL(KIND = float), INTENT(OUT) :: Qout !!downstream discharge at time t + dt (m3/s) !local declarations INTEGER (KIND = short) :: doy !!day of year CHARACTER (LEN = 3) :: QdivertedDOY !-------------------------------end of declaration----------------------------- !compute day of year (used to set ecological flow and weir function) doy = DayOfYear (time, 'noleap') !set doy for weir function QdivertedDOY = ToString ( div % weirDOY (doy) ) CALL TableGetValue ( valueIn = Qin, tab = div % weir, keyIn = 'Qstream', & keyOut = QdivertedDOY, match = 'linear', & valueOut = div % QinChannel, bound = 'extendconstant' ) !check environmental flow and set Qout div % QinChannel = MIN ( Qin - div % eFlow (doy), div % QinChannel) IF ( div % QinChannel < 0. ) THEN !Qin is less than eflow div % QinChannel = 0. END IF !set Qout Qout = Qin - div % QinChannel RETURN END SUBROUTINE DivertFlow !============================================================================== !| Description: ! Discharge routing through reservoir. Given Qin, find Qout and stage ! References: SUBROUTINE LevelPool & ! (time, dt, dtrk, Qin_t, Qin_tplusdt, pool) !arguments with intent in: TYPE(DateTime), INTENT(IN) :: time INTEGER, INTENT(IN) :: dt !!main time step (s) INTEGER, INTENT(IN) :: dtrk !!time step to run runge-kutta (s) REAL(KIND = float), INTENT(IN) :: Qin_t !!Input discharge at time t (m3/s) REAL(KIND = float), INTENT(IN) :: Qin_tplusdt !!Input discharge at time t + dt (m3/s) !arguments with intent inout: TYPE(Reservoir), INTENT(INOUT) :: pool !arguments with intent out: !REAL(KIND = float), INTENT(OUT) :: Qout !!downstream discharge at time t + dt (m3/s) !local declarations INTEGER :: nstep !!number of sub-steps REAL(KIND = float) :: Qin_stepbegin, Qin_stepend !(m3/s) REAL(KIND = float) :: Qinpool_t !!Actual input dicharge to the pool at time t (m3/s) REAL(KIND = float) :: Qinpool_tplusdt !!Actual input dicharge to the pool at time t+dt (m3/s) REAL(KIND = float) :: h !!current stage (m) INTEGER :: j LOGICAL :: qinIsRising INTEGER (KIND = short) :: doy !!day of year CHARACTER (LEN = 3) :: Qdiverted CHARACTER (LEN = 3) :: QoutDOY !!doy to compute Qout !-------------------------------end of declaration----------------------------- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! ! ! * Qinpool_tplusdt ! | ! ! | ! ! Qin_stepend O | ! ! | | ! ! | | ! ! Qin_stepbegin O | | ! ! . | | | ! ! . | | | ! ! . | | | ! ! . | | | ! ! . | | | ! ! O | | | ! ! | | | | ! ! | | | | ! ! Qinpool_t * ____|...........|_____|_____| ! ! ! ! <-----> ! ! dtrk ! ! t <-----------------------------> t + dt ! ! dt ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !initialize variables to run runge-kutta SELECT CASE ( TRIM(pool % typ) ) CASE ( 'off' ) !diverted discharge by weir is the input to the pool ! time t CALL TableGetValue ( valueIn = Qin_t, tab = pool % weir, keyIn = 'Qstream', & keyOut ='Qdiverted', match = 'linear', & valueOut = Qinpool_t, bound = 'extendconstant' ) ! time t+dt CALL TableGetValue ( valueIn = Qin_tplusdt, tab = pool % weir, keyIn = 'Qstream', & keyOut ='Qdiverted', match = 'linear', & valueOut = Qinpool_tplusdt, bound = 'extendconstant' ) !if pool is full discharge can not enter anymore IF (pool % stage > pool % stageMax) THEN Qinpool_t = 0. Qinpool_tplusdt = 0. END IF CASE ('on') !discharge from stream enters pool Qinpool_t = Qin_t Qinpool_tplusdt = Qin_tplusdt IF ( Qinpool_tplusdt > Qinpool_t ) THEN qinIsRising = .TRUE. ELSE qinIsRising = .FALSE. END IF CASE DEFAULT CALL Catch ('error', 'LevelPool', 'wrong reservoir type in pool: ', & argument = TRIM(pool % name)) END SELECT !initialize stage h = pool % stage !compute day of year (used to set ecological flow) doy = DayOfYear (time, 'noleap') !set doy for computing Qout QoutDOY = ToString ( pool % geometryDOY (doy) ) ! manage high level IF ( pool % highLevel .AND. h > pool % fullReservoirLevel ) THEN IF ( pool % rising) THEN IF ( qInIsRising ) THEN SELECT CASE ( pool % qoutRule ) CASE ( 1 ) !qout = qin pool % Qout = Qinpool_tplusdt END SELECT RETURN END IF ELSE SELECT CASE ( pool % qoutRule ) CASE ( 1 ) !qout = qin pool % Qout = Qinpool_tplusdt END SELECT RETURN END IF END IF !loop throughout steps nstep = dt / dtrk DO j = 1, nstep Qin_stepbegin = LinearInterp (0., REAL(dt), Qinpool_t, & Qinpool_tplusdt, 1.*(j-1)*dtrk) Qin_stepend = LinearInterp (0., REAL(dt), Qinpool_t, & Qinpool_tplusdt, 1.*j*dtrk) !update level pool h = RungeKutta ( hini = h, Qin_begin = Qin_stepbegin, & Qin_end = Qin_stepend, dt = dtrk, & geometry = pool % geometry, & rk = pool % rk, typOut = pool % typOut, & hTarget = pool % stageTarget, & eFlow = pool % eFlow (doy), & freeFlow = pool % freeFlow, & freeFlowElevation = pool % freeFlowElevation, & QoutDOY = QoutDOY ) END DO !downstream discharge SELECT CASE ( TRIM (pool % typ ) ) CASE ('on') !on-stream reservoir IF (pool % typOut == 2 .AND. h < pool % stageTarget) THEN pool % Qout = MIN( pool % eFlow (doy), Qin_tplusdt) ELSE IF ( pool % typOut == 2 .AND. h >= pool % stageTarget) THEN CALL TableGetValue ( valueIn = h, tab = pool % geometry, keyIn = 'h', & keyOut ='Qout', match = 'linear', & valueOut = pool % Qout, bound = 'extendconstant' ) IF ( pool % eFlow (doy) > 0. ) THEN pool % Qout = MAX (pool % eFlow (doy), pool % Qout) END IF ELSE IF ( pool % stage <= pool % freeFlowElevation .AND. & Qin_tplusdt < pool % freeFlow ) THEN pool % Qout = Qin_tplusdt ELSE CALL TableGetValue ( valueIn = h, tab = pool % geometry, keyIn = 'h', & keyOut ='Qout', match = 'linear', & valueOut = pool % Qout, bound = 'extendconstant' ) END IF END IF CASE ('off') !off-stream reservoir pool % Qout = Qin_tplusdt - Qinpool_tplusdt END SELECT !update stage pool % stage = h RETURN END SUBROUTINE LevelPool !============================================================================== !| Description: ! returns the wetted area of a triangular, rectangular or trapezoidal cross ! section ! ! Reference: ! ! Todini, E. (2007) A mass conservative and water storage consistent variable ! parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659. FUNCTION WettedArea & ( y, B0, alpha ) & RESULT (a) IMPLICIT NONE ! Function arguments ! Scalar arguments with intent(in): REAL (KIND = float), INTENT (IN) :: y !!water stage (m) 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) ! Scalar arguments with intent(OUT): REAL (KIND = float) :: a REAL (KIND = float) :: calpha !!cotangent of angle alpha !------------end of declaration------------------------------------------------ calpha = 1. / TAN(alpha*degToRad) a = (B0 + y * calpha) * y RETURN END FUNCTION WettedArea !============================================================================== !| Description: ! returns the wetted perimeter of a triangular, rectangular or trapezoidal cross ! section ! ! Reference: ! ! Todini, E. (2007) A mass conservative and water storage consistent variable ! parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659. FUNCTION WettedPerimeter & ( y, B0, alpha ) & RESULT (p) IMPLICIT NONE ! Function arguments ! Scalar arguments with intent(in): REAL (KIND = float), INTENT (IN) :: y !!water stage (m) 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) ! Scalar arguments with intent(OUT): REAL (KIND = float) :: p REAL (KIND = float) :: salpha !!sine of angle alpha !------------end of declaration------------------------------------------------ salpha = SIN(alpha*degToRad) p = B0 + 2 * y / salpha RETURN END FUNCTION WettedPerimeter !============================================================================== !| Description: ! returns the topwidth of a triangular, rectangular or trapezoidal cross ! section ! ! Reference: ! ! Todini, E. (2007) A mass conservative and water storage consistent variable ! parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659. FUNCTION TopWidth & ( y, B0, alpha ) & RESULT (b) IMPLICIT NONE ! Function arguments ! Scalar arguments with intent(in): REAL (KIND = float), INTENT (IN) :: y !!water stage (m) 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) ! Scalar arguments with intent(OUT): REAL (KIND = float) :: b REAL (KIND = float) :: calpha !!cotangent of angle alpha !------------end of declaration------------------------------------------------ calpha = 1. / TAN(alpha*degToRad) b = B0 + 2 * y * calpha RETURN END FUNCTION TopWidth !============================================================================== !| Description: ! returns the normal discharge (m3/s) with Chezy equation. ! ! Reference: ! ! Todini, E. (2007) A mass conservative and water storage consistent variable ! parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659. FUNCTION Discharge & ( y, B0, alpha, slope, n ) & RESULT (Q) IMPLICIT NONE ! Function arguments ! Scalar arguments with intent(in): REAL (KIND = float), INTENT (IN) :: y !!water stage (m) 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 riughness coefficient (s m^(-1/3)) ! Scalar arguments with intent(OUT): REAL (KIND = float) :: Q !------------end of declaration------------------------------------------------ Q = slope**0.5 / n * & WettedArea ( y, B0, alpha )**(5./3.) / WettedPerimeter ( y, B0, alpha )**(2./3.) RETURN END FUNCTION Discharge !============================================================================== !| Description: ! returns the normal velocity (m/s) with Chezy equation ! ! Reference: ! ! Todini, E. (2007) A mass conservative and water storage consistent variable ! parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659. FUNCTION Velocity & ( y, B0, alpha, slope, n ) & RESULT (v) IMPLICIT NONE ! Function arguments ! Scalar arguments with intent(in): REAL (KIND = float), INTENT (IN) :: y !!water stage (m) 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 riughness coefficient (s m^(-1/3)) ! Scalar arguments with intent(OUT): REAL (KIND = float) :: v !------------end of declaration------------------------------------------------ v = slope**0.5 / n * & WettedArea ( y, B0, alpha )**(2./3.) / WettedPerimeter ( y, B0, alpha )**(2./3.) RETURN END FUNCTION Velocity !============================================================================== !| Description: ! returns the celerity (m/s) ! ! Reference: ! ! Todini, E. (2007) A mass conservative and water storage consistent variable ! parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659. FUNCTION Celerity & ( y, B0, alpha, slope, n ) & RESULT (c) IMPLICIT NONE ! Function arguments ! Arguments with intent(in): REAL (KIND = float), INTENT (IN) :: y !!water stage (m) 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 riughness coefficient (s m^(-1/3)) ! Arguments with intent(OUT): REAL (KIND = float) :: c !local declaration: REAL (KIND = float) :: A !!wetted area (m2) REAL (KIND = float) :: P !!wetted perimeter (m) REAL (KIND = float) :: B !!top width (m) REAL (KIND = float) :: salpha !!sine of angle alpha !------------end of declaration------------------------------------------------ A = WettedArea ( y, B0, alpha ) P = WettedPerimeter ( y, B0, alpha ) B = TopWidth ( y, B0, alpha ) salpha = SIN(alpha*degToRad) c = 5./3. * slope**0.5 / n * A**(2./3.) / P**(2./3.) * & ( 1. - 4./5. * A / (B * P * salpha)) RETURN END FUNCTION Celerity !============================================================================== !| Description: ! returns the correction factor to be used in MCT algorithm (-) ! ! Reference: ! ! Todini, E. (2007) A mass conservative and water storage consistent variable ! parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659. FUNCTION Beta & ( y, B0, alpha, slope, n ) & RESULT (cf) IMPLICIT NONE ! Function arguments ! Arguments with intent(in): REAL (KIND = float), INTENT (IN) :: y !!water stage (m) 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 riughness coefficient (s m^(-1/3)) ! Arguments with intent(OUT): REAL (KIND = float) :: cf !local declarations REAL (KIND = float) :: A !!wetted area (m2) REAL (KIND = float) :: P !!wetted perimeter (m) REAL (KIND = float) :: B !!top width (m) REAL (KIND = float) :: salpha !!sine of angle alpha !------------end of declaration------------------------------------------------ IF (y == 0.) THEN cf = 0. RETURN END IF A = WettedArea ( y, B0, alpha ) P = WettedPerimeter ( y, B0, alpha ) B = TopWidth ( y, B0, alpha ) salpha = SIN(alpha*degToRad) cf = 5./3. * ( 1. - 4./5. * A / (B * P * salpha)) RETURN END FUNCTION Beta !============================================================================== !| Description: ! returns the normal depth (m) using Newton-Raphson. ! ! References: ! ! Todini, E. (2007) A mass conservative and water storage consistent variable ! parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659. FUNCTION NormalDepth & ( Q, B0, alpha, slope, n ) & RESULT (nd) IMPLICIT NONE ! Function arguments ! Arguments with intent(in): REAL (KIND = float), INTENT (IN) :: Q !!discharge (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 riughness coefficient (s m^(-1/3)) ! Arguments with intent(OUT): REAL (KIND = float) :: nd !local declarations INTEGER :: i INTEGER :: maxiter = 100 REAL (KIND = float) :: nd1 REAL (KIND = float) :: Q0 REAL (KIND = float) :: dQdy !first derivative REAL (KIND = float) :: tol = 0.001 !------------end of declaration------------------------------------------------ IF (Q == 0.) THEN nd = 0. END IF !first guess, nd = 1 m nd = 1. DO i = 1, maxiter Q0 = Discharge ( nd, B0, alpha, slope, n ) dQdy = TopWidth ( nd, B0, alpha ) * Celerity ( nd, B0, alpha, slope, n ) nd1 = nd - (Q0 - Q) / dQdy !test conf IF ( ABS(nd1 - nd) / nd <= tol) THEN !root found nd = nd1 RETURN END IF nd = nd1 END DO RETURN RETURN END FUNCTION NormalDepth !============================================================================== !| Description: ! returns water depth given wetted area (m) using Newton-Raphson. ! ! References: ! ! Todini, E. (2007) A mass conservative and water storage consistent variable ! parameter Muskingum-Cunge approach, Hydrol. Earth Syst. Sci.,1645-1659. FUNCTION DepthFromArea & ( A, B0, alpha ) & RESULT (y) IMPLICIT NONE ! Function arguments ! Arguments with intent(in): REAL (KIND = float), INTENT (IN) :: A !!area (m2) 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) ! Arguments with intent(OUT): REAL (KIND = float) :: y !local declaration: INTEGER :: i INTEGER :: maxiter = 100 REAL (KIND = float) :: y1 REAL (KIND = float) :: A0 REAL (KIND = float) :: dAdy !first derivative REAL (KIND = float) :: tol = 0.001 REAL (KIND = float) :: calpha !!cotangent of angle alpha !------------end of declaration------------------------------------------------ calpha = 1. / TAN(alpha*degToRad) !first guess, y = 1 m y = 1. DO i = 1, maxiter A0 = WettedArea ( y, B0, alpha ) dAdy = calpha * y + (B0 + y * calpha) y1 = y - (A0 - A) / dAdy !test conf IF ( ABS(y1 - y) / y <= tol) THEN !root found y = y1 RETURN END IF y = y1 END DO RETURN RETURN END FUNCTION DepthFromArea !============================================================================== !| Description: ! solve mass continuity of reservoir using runge kutta method ! ! References: ! ! Chow, V. T., Maidment, D. R. &Mays, L. W. (1988) Applied Hydrology. ! McGraw-Hill, New York, USA ! ! Fenton, J.D. (1992) Reservoir routing, Hydrological Sciences Journal, ! 37(3), 233-246, DOI: 10.1080/02626669209492584 FUNCTION RungeKutta & ( hini, Qin_begin, Qin_end, dt, geometry, rk, typOut, hTarget, & eFlow, freeFlow, freeFlowElevation, QoutDOY ) & RESULT (stage) IMPLICIT NONE ! Function arguments ! Arguments with intent(in): REAL (KIND = float), INTENT(IN) :: hini !!initial stage value (m) REAL (KIND = float), INTENT(IN) :: Qin_begin !!input discharge at the beginning of time step (m3/s) REAL (KIND = float), INTENT(IN) :: Qin_end !!input discharge at the end of time step (m3/s) INTEGER, INTENT(IN) :: dt !!time step (s) TYPE (Table), INTENT(IN) :: geometry !!geometry of reservoir INTEGER, INTENT(IN) :: rk !! solution order 3 or 4 INTEGER, INTENT(IN) :: typOut !!type ofoutflow: 1=free flow 2=target level REAL (KIND = float), INTENT(IN) :: hTarget !!follow a target (observed) stage (m) REAL (KIND = float), INTENT(IN) :: eFlow !!environmental flow (m3/s) REAL (KIND = float), INTENT(IN) :: freeFlow !! free flow (m3/s) REAL (KIND = float), INTENT(IN) :: freeFlowElevation !! free flow elevation (m asl) CHARACTER (LEN = 3), INTENT(IN) :: QoutDOY !!doy to compute Qout ! Arguments with intent(OUT): REAL (KIND = float) :: stage !local declaration: REAL(KIND = float) :: h1, h2 REAL(KIND = float) :: h3 !used only if rk = 4 REAL(KIND = float) :: dh1, dh2, dh3 REAL(KIND = float) :: dh4 !used only if rk = 4 REAL(KIND = float) :: dh REAL(KIND = float) :: Q1third, Q2third REAL(KIND = float) :: Qout REAL(KIND = float) :: area !(m2) !------------end of declaration------------------------------------------------ SELECT CASE (rk) CASE(3) IF ( typOut == 2 .AND. hini < hTarget) THEN Qout = MIN (eFlow, Qin_begin) !step 1 !find area in table geometry CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', & keyOut='area', match='linear', & valueOut=area, bound = 'extendconstant' ) !compute dh1 dh1 = (Qin_begin - Qout) * dt / area !update h1 h1 = hini + dh1 / 3.0 !step 2 Q1third = LinearInterp (0., REAL(dt), Qin_begin, Qin_end, 1.0/3.0*dt) !find area in table geometry CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', & keyOut='area', match='linear', & valueOut=area, bound = 'extendconstant' ) !compute dh2 dh2 = (Q1third - Qout) * dt / area !update h2 h2 = hini + 2.0 * dh2 / 3.0 !step 3 Q2third = LinearInterp (0., real(dt), Qin_begin, Qin_end, 2.0/3.0*dt) !find area in table geometry CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', & keyOut='area', match='linear', & valueOut=area, bound = 'extendconstant' ) !compute dh3 dh3 = (Q2third - Qout) * dt / area ELSE !step 1 !find Qout IF ( hini <= freeFlowElevation .AND. Qin_begin < freeFlow ) THEN Qout = Qin_begin ELSE CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', & keyOut = QoutDOY, match = 'linear', & valueOut = Qout, bound = 'extendconstant' ) IF ( eFlow > 0.) THEN Qout = MAX (eFlow, Qout) END IF END IF !find area in table geometry CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', & keyOut = 'area', match = 'linear', & valueOut = area, bound = 'extendconstant' ) !compute dh1 dh1 = (Qin_begin - Qout) * dt / area !update h1 h1 = hini + dh1 / 3.0 !step2 Q1third = LinearInterp (0., REAL(dt), Qin_begin, Qin_end, 1.0/3.0*dt) !find Qout IF ( h1 <= freeFlowElevation .AND. Q1third < freeFlow ) THEN Qout = Q1third ELSE CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', & keyOut = QoutDOY, match = 'linear', & valueOut = Qout, bound = 'extendconstant' ) END IF !find area in table geometry CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', & keyOut = 'area', match = 'linear', & valueOut = area, bound = 'extendconstant' ) !compute dh2 dh2 = (Q1third - Qout) * dt / area !update h2 h2 = hini + 2.0 * dh2 / 3.0 !step3 Q2third = LinearInterp (0., real(dt), Qin_begin, Qin_end, 2.0/3.0*dt) !find Qout IF ( h2 <= freeFlowElevation .AND. Q2third < freeFlow ) THEN Qout = Q2third ELSE CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', & keyOut = QoutDOY, match = 'linear', & valueOut = Qout, bound = 'extendconstant' ) END IF !find area in table geometry CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', & keyOut = 'area', match = 'linear', & valueOut = area, bound = 'extendconstant' ) !compute dh3 dh3 = (Q2third - Qout) * dt / area END IF dh = dh1 / 4.0 + 3.0 * dh3 / 4.0 CASE(4) !IF ( typOut == 2 .AND. hini < hTarget) THEN ! ! Qout = MIN (eFlow, Qin_begin) ! ! !step 1 ! !find area in table geometry ! CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', & ! keyOut='area', match='linear', & ! valueOut=area, bound = 'extendlinear' ) ! ! !compute dh1 ! dh1 = (Qin_begin - Qout) * area ! !write(*,*) dh1, Qout, eflow, Qin_begin, hini, htarget ! ! !update h1 ! h1 = hini + dh1 / 3.0 ! ! !step 2 ! Q1third = LinearInterp (0., REAL(dt), Qin_begin, Qin_end, 1.0/3.0*dt) ! ! !find area in table geometry ! CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', & ! keyOut='area', match='linear', & ! valueOut=area, bound = 'extendlinear' ) ! ! !compute dh2 ! dh2 = (Q1third - Qout) * dt / area ! ! !update h2 ! h2 = hini - 1.0 / 3.0 * dh1 + dh2 ! ! !step 3 ! Q2third = LinearInterp (0., real(dt), Qin_begin, Qin_end, 2.0/3.0*dt) ! ! !find area in table geometry ! CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', & ! keyOut='area', match='linear', & ! valueOut=area, bound = 'extendlinear' ) ! ! !compute dh3 ! dh3 = (Q2third - Qout) * dt / area ! ! !update h3 ! h3 = hini + dh1 - dh2 + dh3 ! ! !step4 ! !find area in table geometry ! CALL TableGetValue ( valueIn = h3, tab = geometry, keyIn = 'h', & ! keyOut='area', match='linear', & ! valueOut=area, bound = 'extendlinear' ) ! ! !compute dh4 ! dh4 = (Qin_end - Qout) * dt / area ! !ELSE !step 1 !find Qout IF ( typOut == 2 .AND. hini < hTarget) THEN Qout = MIN (eFlow, Qin_begin) ELSE IF ( typOut == 2 .AND. hini >= hTarget) THEN CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', & keyOut = QoutDOY, match = 'linear', & valueOut = Qout, bound = 'extendconstant' ) IF ( eFlow > 0. ) THEN Qout = MAX (eFlow, Qout) END IF ELSE IF ( hini <= freeFlowElevation .AND. Qin_begin < freeFlow ) THEN Qout = Qin_begin ELSE CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', & keyOut = QoutDOY, match = 'linear', & valueOut = Qout, bound = 'extendconstant' ) END IF !find area in table geometry CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', & keyOut = 'area', match = 'linear', & valueOut = area, bound = 'extendconstant' ) !compute dh1 dh1 = (Qin_begin - Qout) * dt / area !update h1 h1 = hini + dh1 / 3.0 !step2 Q1third = LinearInterp (0., REAL(dt), Qin_begin, Qin_end, 1.0/3.0*dt) IF ( typOut == 1 ) THEN !update Qout IF ( h1 <= freeFlowElevation .AND. Q1third < freeFlow ) THEN Qout = Q1third ELSE CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', & keyOut = QoutDOY, match = 'linear', & valueOut = Qout, bound = 'extendconstant' ) END IF END IF !find area in table geometry CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', & keyOut = 'area', match = 'linear', & valueOut = area, bound = 'extendconstant' ) !compute dh2 dh2 = (Q1third - Qout) * dt / area !update h2 h2 = hini - 1.0 / 3.0 * dh1 + dh2 !step3 Q2third = LinearInterp (0., real(dt), Qin_begin, Qin_end, 2.0/3.0*dt) IF ( typOut == 1 ) THEN !update Qout IF ( h2 <= freeFlowElevation .AND. Q2third < freeFlow ) THEN Qout = Q2third ELSE CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', & keyOut = QoutDOY, match = 'linear', & valueOut = Qout, bound = 'extendconstant' ) END IF END IF !find area in table geometry CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', & keyOut = 'area', match = 'linear', & valueOut = area, bound = 'extendconstant' ) !compute dh3 dh3 = (Q2third - Qout) * dt / area !update h3 h3 = hini + dh1 - dh2 + dh3 !step 4 IF ( typOut == 1 ) THEN !update Qout IF ( h3 <= freeFlowElevation .AND. Qin_end < freeFlow ) THEN Qout = Qin_end ELSE CALL TableGetValue ( valueIn = h3, tab = geometry, keyIn = 'h', & keyOut = QoutDOY, match = 'linear', & valueOut = Qout, bound = 'extendconstant' ) END IF END IF !find area in table geometry CALL TableGetValue ( valueIn = h3, tab = geometry, keyIn = 'h', & keyOut = 'area', match = 'linear', & valueOut = area, bound = 'extendconstant' ) !compute dh4 dh4 = (Qin_end - Qout) * dt / area !END IF dh = dh1 / 8.0 + 3.0 * dh2 / 8.0 + 3.0 * dh3 / 8.0 + dh4 / 8.0 CASE DEFAULT CALL Catch ('error', 'RungeKutta', 'order Runge-Kutta not valid: ', & argument = ToString(rk)) STOP END SELECT !update reservoir stage stage = hini + dh RETURN END FUNCTION RungeKutta END MODULE RoutingModels