!! 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