RungeKutta – FEST

private function RungeKutta(hini, Qin_begin, Qin_end, dt, geometry, rk, typOut, hTarget, eFlow, freeFlow, freeFlowElevation, QoutDOY, hmin) result(stage)

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

Arguments

Type	Intent		Name
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
real(kind=float),	intent(in)	::	hmin	minimum stage value (m)

Return Value real(kind=float)

Variables

Type	Visibility		Name
real(kind=float),	public	::	Q1third
real(kind=float),	public	::	Q2third
real(kind=float),	public	::	Qout
real(kind=float),	public	::	S1	used to compute deltaS (m3)
real(kind=float),	public	::	S2	used to compute deltaS (m3)
real(kind=float),	public	::	area
real(kind=float),	public	::	deltaS	delta storage volume (m3)
real(kind=double),	public	::	dh
real(kind=double),	public	::	dh1
real(kind=double),	public	::	dh2
real(kind=double),	public	::	dh3
real(kind=double),	public	::	dh4
real(kind=float),	public	::	h1
real(kind=float),	public	::	h2
real(kind=float),	public	::	h3

Source Code

FUNCTION RungeKutta &
( hini, Qin_begin, Qin_end, dt, geometry, rk, typOut, hTarget, &
  eFlow, freeFlow, freeFlowElevation, QoutDOY, hmin ) &
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 
REAL (KIND = float), INTENT(IN)  :: hmin !!minimum stage value (m)


! 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 = double) :: dh1, dh2, dh3
REAL(KIND = double) :: dh4	!used only if rk = 4
REAL(KIND = double) :: dh
REAL(KIND = float) :: Q1third, Q2third
REAL(KIND = float) :: Qout
REAL(KIND = float) :: area !(m2)
REAL(KIND = float) :: S1, S2  !!used to compute deltaS (m3)
REAL(KIND = float) :: deltaS  !!delta storage volume (m3)


!------------end of declaration------------------------------------------------


SELECT CASE (rk)
        
CASE(2)
    
    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 = DBLE ( (Qin_begin - Qout) * dt / area )
            
        !update h1
        h1 = DBLE ( hini + dh1 )
            
        !step 2
                  
        !find area in table geometry
        CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', &
                                keyOut='area', match='linear', &
                                valueOut=area, bound = 'extendconstant' )
            
        !compute dh2
        dh2 = DBLE ( (Qin_end - 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 = DBLE ( (Qin_begin - Qout) * dt / area )
            
        !update h1
        h1 = DBLE ( hini + dh1 )
  
        !step2
        !find Qout
        IF ( h1 <= freeFlowElevation .AND. Qin_begin < freeFlow ) THEN
            Qout = Qin_begin
        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 = DBLE ( (Qin_end - Qout) * dt / area )
            
	END IF	  
  
    dh = DBLE ( 0.5 * (dh1 + dh2) )
    
    
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 = DBLE ( (Qin_begin - Qout) * dt / area )
            
        !update h1
        h1 = DBLE ( 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 = DBLE ( (Q1third - Qout) * dt / area )

        !update h2
        h2 = DBLE ( 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 = DBLE ( (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 = DBLE ( (Qin_begin - Qout) * dt / area )
            
        !update h1
        h1 = DBLE ( 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 = DBLE ( (Q1third - Qout) * dt / area )
            
        !update h2
        h2 = DBLE ( 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 = DBLE ( (Q2third - Qout) * dt / area )
	END IF	  
  
    dh = DBLE ( dh1 / 4.0 + 3.0 * dh3 / 4.0 )

CASE(4)
        
            
    !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 = DBLE( (Qin_begin - Qout) * dt / area)
            
    !update h1
    h1 = DBLE(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 = DBLE( (Q1third - Qout) * dt / area)
            
    !update h2
    h2 = DBLE(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 = DBLE ((Q2third - Qout) * dt / area)
            
    !update h3
    h3 =DBLE( 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 = DBLE((Qin_end - Qout) * dt / area)
            
!END IF
        
dh = DBLE(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 = DBLE (hini + dh)

!detect and filter instability
IF (QfilterApply) THEN
    deltaS = Qfilter * dt !volume for a Qfilter (m3/s) net input discharge

    !compute S1 volume at previous step stage
    CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', &
                                    keyOut = 'volume', match = 'linear', &
                                    valueOut = S1, bound = 'extendconstant' )

    !compute S2 volume at current step stage
    CALL TableGetValue ( valueIn = stage, tab = geometry, keyIn = 'h', &
                                    keyOut = 'volume', match = 'linear', &
                                    valueOut = S2, bound = 'extendconstant' )

    IF ( (S2 - S1) > deltaS) THEN
        !too rapid volume change, possible instability detected
        stage = hini
    END IF
END IF


IF ( stage < hmin ) THEN
    stage = hmin
END IF

RETURN
END FUNCTION RungeKutta

RungeKutta Function

Contents

Variables

Source Code

private function RungeKutta(hini, Qin_begin, Qin_end, dt, geometry, rk, typOut, hTarget, eFlow, freeFlow, freeFlowElevation, QoutDOY, hmin) result(stage)

Arguments

Return Value real(kind=float)

Variables

Source Code