LevelPool – FEST

public subroutine LevelPool(time, dt, dtrk, Qin_t, Qin_tplusdt, pool)

Discharge routing through reservoir. Given Qin, find Qout and stage References:

Arguments

Type	Intent		Name
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)
type(Reservoir),	intent(inout)	::	pool

Variables

Type	Visibility		Name
real(kind=float),	public	::	Q1	used to compute deltaQ (m3/s)
real(kind=float),	public	::	Q2	used to compute deltaQ (m3/s)
character(len=3),	public	::	Qdiverted
real(kind=float),	public	::	Qin_stepbegin
real(kind=float),	public	::	Qin_stepend
real(kind=float),	public	::	Qinpool_t	Actual input dicharge to the pool at time t (m3/s)
real(kind=float),	public	::	Qinpool_tplusdt	Actual input dicharge to the pool at time t+dt (m3/s)
character(len=3),	public	::	QoutDOY	doy to compute 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	::	deltaQ	delta Qout (m3/s)
real(kind=float),	public	::	deltaS	delta storage volume (m3)
real(kind=float),	public	::	deltah	delta h to compute dtvar (m)
integer(kind=short),	public	::	doy	day of year
integer,	public	::	dtcalc	actual time step for runge kutta (s)
real(kind=float),	public	::	dtvar	variable time step (s)
real(kind=float),	public	::	h	current stage (m)
integer,	public	::	j
integer,	public	::	nstep	number of sub-steps
logical,	public	::	qinIsRising

Source Code

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
INTEGER                         :: dtcalc !!actual time step for runge kutta (s) 
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) 
REAL(KIND = float)              :: dtvar !!variable time step (s) 
REAL(KIND = float)              :: deltah !!delta h to compute dtvar (m)
REAL(KIND = float)              :: Q1, Q2  !!used to compute deltaQ (m3/s)
REAL(KIND = float)              :: S1, S2  !!used to compute deltaS (m3)
REAL(KIND = float)              :: deltaQ  !!delta Qout (m3/s)
REAL(KIND = float)              :: deltaS  !!delta storage volume (m3)
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


    !check variable time step
IF ( dtAdjust) THEN
    deltah = 0.01 !1 cm drop to compute dtvar

    !compute S1 volume at stage = h
    CALL TableGetValue ( valueIn = h, tab = pool % geometry, keyIn = 'h', &
                                    keyOut = 'volume', match = 'linear', &
                                    valueOut = S1, bound = 'extendlinear' )


    !compute S2 volume at stage = h - deltah
    CALL TableGetValue ( valueIn = h - deltah, tab = pool % geometry, keyIn = 'h', &
                                    keyOut = 'volume', match = 'linear', &
                                    valueOut = S2, bound = 'extendlinear' )
    
    deltaS = ABS (S2 - S1)

    IF ( deltaS == 0. ) THEN
        CALL Catch ('error', 'LevelPool', 'detected constant volume in reservoir: ', &
	                          argument = TRIM(pool % name))
    END IF


    !compute Q1 at stage = h
    CALL TableGetValue ( valueIn = h, tab = pool % geometry, keyIn = 'h', &
                                    keyOut = QoutDOY, match = 'linear', &
                                    valueOut = Q1, bound = 'extendlinear' )

    !compute Q2 at stage = h - deltah
    CALL TableGetValue ( valueIn = h - deltah, tab = pool % geometry, keyIn = 'h', &
                                    keyOut = QoutDOY, match = 'linear', &
                                    valueOut = Q2, bound = 'extendlinear' )

    deltaQ = ABS (Q1 - Q2)

    IF (deltaQ == 0.) THEN
        deltaQ = 1.
    END IF

    !dtvar = 2. * deltaS / deltaQ
    dtvar = 1. * deltaS / deltaQ
    IF ( dtMin > 0 .AND. dtvar < dtMin) THEN
        dtvar = dtMin
        nstep = dt / dtvar
    ELSE
        IF ( dtvar < dtrk) THEN
          nstep = dt / dtvar + 1
        ELSE
          nstep = dt / dtrk
        END IF
    END IF

    dtcalc = dt / nstep
    
ELSE
    nstep = dt / dtrk
    dtcalc = dtrk
END IF

!loop throughout steps
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 = dtcalc, &
                    geometry = pool % geometry, &
                    rk = pool % rk, typOut = pool % typOut, &
                    hTarget = pool % stageTarget, &
                    eFlow = pool % eFlow (doy), &
                    freeFlow = pool % freeFlow, &
                    freeFlowElevation = pool % freeFlowElevation, &
                    QoutDOY = QoutDOY, hmin = pool % stageMin  )
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 = QoutDOY, 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 = QoutDOY, 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

LevelPool Subroutine

Contents

Variables

Source Code

public subroutine LevelPool(time, dt, dtrk, Qin_t, Qin_tplusdt, pool)

Arguments

Variables

Source Code