reservoir_process Subroutine

  subroutine reservoir_process(ih, in1, in2, ecal, humi, ida, iyr, mo, nc, omega, ra, subp, tx, varoute)

    real(dp), dimension(:), intent(in) :: ecal
    real(dp), dimension(:), intent(in) :: humi
    integer, intent(in) :: ida
    integer, intent(in) :: iyr
    integer, intent(in) :: mo
    integer, dimension(13), intent(in) :: nc
    real(dp), dimension(12), intent(in) :: omega
    ! real(dp), intent(inout) :: precip
    real(dp), dimension(:), intent(in) :: ra
    real(dp), dimension(:), intent(in) :: subp
    real(dp), dimension(:), intent(in) :: tx
    real(dp), dimension(:, :), intent(inout) :: varoute
    ! CALLED FROM mainpro.f90 during route command
    ! actual hydrograph storage location
    integer, intent(in) :: ih
    ! inum1: subbasin number of this reservoir
    integer, intent(in) :: in1
    ! inum2: the input (upstream) subbasin or hydrograph storage location
    integer, intent(in) :: in2
    ! subbasin number
    integer :: sub
    ! reservoir number
    integer :: res
    ! potential evapotranspiration of water surface [mm]
    real(dp) :: ETpot_water_mm
    ! potential evapotranspiration of land surface [mm]
    real(dp) :: ETpot_land_mm
    ! potential evapotranspiration of land surface [m^3 / d]
    real(dp) :: ETpot_water_m3
    ! potential evapotranspiration of land surface [m^3 / d]
    real(dp) :: ETpot_land_m3
    ! area weighted pot. evapotranspiration (water / land)
    real(dp) :: ETpot_weighted_mm
    ! fraction of flooded surface area total reservoir area
    real(dp) :: a_frac_wet
    ! A = flooded surface area [km^2]
    real(dp) :: A
    ! Area of Reservoir which is dry [km^2]
    real(dp) :: ResA_dry
    ! Area of Reservoir which is wet [km^2]
    real(dp) :: ResA_wet
    ! climatic water balance [m^3 / d]
    real(dp) :: cwb
    ! precipation [mm / d] current day on reservoir subbasin
    real(dp) :: precip
    ! volume of Reservoir not filled [m3]
    real(dp) :: Res_space
    ! delta Quantity [m^3 / d]
    real(dp) :: d_Q
    ! Reservoir release [m3 / d]
    real(dp) :: release
    ! water level [m.a.s.l.]
    real(dp) :: SH
    real(dp) :: excess
    ! file unit of current reservoir
    integer :: funit
  !   TYPE (TSubbasin), POINTER :: pS ! pointer on subbasin !#### WATER MANAGEMENT MODULE ####

    !-----------------------------------------------------------------
    ! initialise
    !-----------------------------------------------------------------
    sub = in1
    res = reservoir_get(sub)
    Res_space = 0.
    d_Q = 0.
    release = 0.
    precip = subp(sub) ! Precipitation on reservoir subbasin [mm]
    rsv_Outflow(1, res) = 0.
    rsv_Outflow(2, res) = 0.

    ! set actual reservoir volume to previous day volume
    Rsv(1, res) = rsv_B(res)
    !-----------------------------------------------------------------

    !-----------------------------------------------------------------
    ! check if reservoir is active
    ! if not: set active if actual storage volume exceeds given threshold
    !-----------------------------------------------------------------
    if (.NOT. rsv_active(res) ) then
      if (Rsv(1, res) >= (rsv_Capac_Max(res) - rsv_dead_stor_capac(res)) * rsv_activate_thresh(res) ) &
        rsv_active(res) = .true.
    end if
    !-----------------------------------------------------------------

    !-----------------------------------------------------------------
    ! calculate daily data from monthly data
    ! rsv_Day_Cap_Act, rsv_Day_Fill_Min, Rsv_ann_cycle, rsv_Day_Disch_Min
    !-----------------------------------------------------------------
    call reservoir_daily_from_monthly(ida, iyr, mo, nc)
    !-----------------------------------------------------------------

    !-----------------------------------------------------------------
    ! calculate inflow volume to reservoir subbasin [m^3/d]
    !-----------------------------------------------------------------
    rsv_Inflow(res) = RSV_calc_inflow(in2, res, sub, ida, varoute)

    !    !#################################
    !    !#### WATER MANAGEMENT MODULE ####
    !    !#################################
    !    ! Add inflows from transfers
    !    if ( bWAM_Module ) then
    !       ! check if water transfers take place in current subbasin
    !       if ( wam_is_transfer_subbasin(sub) ) then
    !          pS => wam_get_pointer_subbasin(sub) ! get pointer on current subbasin
    !          ! add subbasins' %inflow(daycounter-1) of previous day
    !          if (daycounter > 1) rsv_Inflow(res) = rsv_Inflow(res) + pS%inflow(daycounter-1) * 86400.
    !       end if
    !    end if
    !    !#################################
    !-----------------------------------------------------------------

    ! calculate fraction of surface runoff to total inflow
    !rsv_frac_sr(res) = varoute(2, in2) / (varoute(2, in2) + varoute(8, in2))

    !-----------------------------------------------------------------
    ! calculate seepage from reservoir [m^3/d]
    !-----------------------------------------------------------------
    rsv_seepage(res) = RSV_calc_seepage(res)

    ! remove seepage from actual and/or dead storage
    ! re-calculate actual storage volume Rsv(2, res) after seepage in [m^3]
    ! re-calculate dead storage if necessary
    call reservoir_remove_seepage(res)
    !-----------------------------------------------------------------

    !-----------------------------------------------------------------
    ! calculate water level [m.a.s.l.] and flooded surface area [km^2]
    !-----------------------------------------------------------------
    SH = 0. ! water level [m.a.s.l.]
    SH = RSV_POL(20, rsv_pol_V(res, :), rsv_pol_L(res, :), Rsv(2, res) + rsv_dead_stor_act(res))
    A = 0. ! flooded surface area
    A = RSV_POL(20, rsv_pol_L(res, :), rsv_pol_A(res, :), SH)

    ! calculate dry and wet (flooded) reservoir surface area [km^2]
    ResA_dry = rsv_pol_A(res, 20) - A
    ResA_dry = max(real(1.e-6, dp), ResA_dry)
    ResA_wet = A
    ResA_wet = max(real(1.e-6, dp), ResA_wet)
    ! fraction of flooded reservoir subbasin area
    a_frac_wet = ResA_wet / rsv_pol_A(res, 20)
    !-----------------------------------------------------------------

    !-----------------------------------------------------------------
    ! calculate ETpot in [mm] and [m^3/d]
    !-----------------------------------------------------------------
    ETpot_water_mm = 0.
    ETpot_water_m3 = 0.
    ETpot_land_mm = 0.
    ETpot_land_m3 = 0.
    ETpot_water_mm = ET_Turc(sub, humi, mo, omega, ra, tx) * ecal(sub) * 1.3 * rsv_evapc(res)
    ETpot_land_mm = ET_Turc(sub, humi, mo, omega, ra, tx) * ecal(sub) * rsv_evapc(res)
    ETpot_land_m3 = ETpot_land_mm / 10 ** 3 * ResA_dry * 10 ** 6
    ETpot_water_m3 = ETpot_water_mm / 10 ** 3 * A * 10 ** 6
    ! area weighted evapotranspiration from reservoir subbasin
    ! ToDo: here, we assume that ET_a from land surface equals ET_pot from land surface which may lead to overestimations
    ETpot_weighted_mm = ETpot_water_mm * a_frac_wet + ETpot_land_mm * (1. - a_frac_wet)

    ! calculate climatic water balance [m^3/d]
    cwb = (precip / 10 ** 3 * rsv_pol_A(res, 20) * 10 ** 6) - (ETpot_land_m3 + ETpot_water_m3)
    !-----------------------------------------------------------------

    !-----------------------------------------------------------------
    ! remove evaporation from actual and/or dead storage
    ! re-calculate actual storage volume Rsv(3, res) after cwb, evaporation in [m^3]
    ! re-calculate dead storage if necessary
    !-----------------------------------------------------------------
    if (cwb < 0.) then
      call reservoir_remove_cwb(res)
    else
      call reservoir_add_cwb(res)
      ! determine the volume of water that must be released in order to reduce Rsv to daily maximum active capacity
      excess = cwb - (rsv_Day_Cap_Act(res) - Rsv(2, res))
      if (excess > 0.) then
        rsv_Outflow(1, res) = max(excess, 0.0_dp) ! simulated reservoir discharge
        rsv_Outflow(2, res) = max(excess, 0.0_dp) ! available reservoir discharge
      end if
    end if
    !-----------------------------------------------------------------

    !*****************************************************************
    !*****************************************************************
    ! If the active storage is empty, filling and outflow is determined
    ! by the status of the dead storage.
    !*****************************************************************
    !*****************************************************************
    if (.NOT. reservoir_is_full_dead_storage(res) ) then ! if dead storage is not full

      !-----------------------------------------------------------------
      ! as long as dead storage is not full, calc filling of dead storage
      ! and discharges downstream
      !-----------------------------------------------------------------
      excess = 0.
      excess = rsv_Inflow(res) - rsv_Day_Disch_Min(res) * 86400.

      ! if inflow is larger than required minimal discharges
      if (excess > 0.) then
        ! fill dead storage
        rsv_dead_stor_act(res) = rsv_dead_stor_act(res) + excess

        ! dead storage is full after filling
        if (reservoir_is_full_dead_storage(res) ) then
          ! set active storage
          Rsv(4, res) = max(rsv_dead_stor_act(res) - rsv_dead_stor_capac(res), 0.0_dp)
          ! set dead storage to max. dead storage capacity
          rsv_dead_stor_act(res) = rsv_dead_stor_capac(res)
          rsv_Outflow(1, res) = rsv_Day_Disch_Min(res) * 86400.

        ! dead storage is not full after filling
        else
          Rsv(4, res) = 0. ! active storage
          rsv_outflow(1, res) = rsv_Day_Disch_Min(res) * 86400.
        end if

      ! if inflow is not larger than required minimal discharges
      else
        ! route inflow through without filling the dead storage
        rsv_Outflow(1, res) = rsv_Inflow(res)
        Rsv(4, res) = 0.

      end if
      !-----------------------------------------------------------------

      ! calculate water level in reservoir [m.a.s.l.]
      call reservoir_mean_water_level(Rsv(4, res), rsv_dead_stor_act(res), res)

      ! set current reservoir volume for next day
      ! if Rsv(4, res) got part of the filling, ignore it this time step
      ! volumes are certainly negligible
      if (Rsv(4, res) == 0. ) Rsv(:, res) = 0.
      rsv_B(res) = Rsv(4, res)

      rsv_Prod_HPP(res) = 0.

    !*****************************************************************
    else
      ! if dead storage was full, continue here...
      !*****************************************************************

      !-----------------------------------------------------------------
      ! calculate reservoir outflow
      !Todo: Warum eigentlich hier?
      !-----------------------------------------------------------------
      call reservoir_outflow
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! if reservoir is operating
      !-----------------------------------------------------------------
      if (rsv_active(res) ) then

        select case(rsv_Mngmt(res)) ! Reservoir Management Option
          case(1) ! Min (rsv_Fill_Min) & max (rsv_Cap_Act) capacity and minimum discharge (Rsv_Dis_Min)
            ! continue below
            release = rsv_Day_Disch_Min(res) ! => for Option 1 no calculation of release but read in values applied, changed H. Koch, 09.07.2019

          case(2) ! safe energy yield (calculation of release depending on hydropower production)
            ! calculate release
            release = RSV_release_firm_energy_yield(res) ! release in [m3 / d]
          case(3) ! release depending on filling of reservoir (might be rising or falling release with incresed filling)
            ! calculate release
            release = RSV_release_dep_on_filling(res) ! [m3 / d]

          case default
            call log_error("reservoir_process", "No valid reservoir management option "// &
              "was selected! Valid setting of rsv_Mngmt (1-3). Check reservoir input file")
          ! do nothing
        end select

        ! calculate minimal discharges
        rsv_Day_Disch_Min(res) = release / 86400. ! [m3 / d - > m3 / s]
        rsv_Day_Fill_Min(res) = max(Rsv(4, res) - release, 0.0_dp)

      end if ! if (rsv_active(res) )
      !-----------------------------------------------------------------

      !if ( .NOT. RSV_is_full_dead_storage(res) ) then
      !   call RSV_fill_dead_storage(res), wie bisher in RSV_release_discharge_refilling
      !   set values of Rsv(4, 5, 6, res)
      ! else: call RSV_release_discharge_refilling

      call reservoir_release_refilling

      ! calculate direct withdrawal from reservoir (irrigation, drinking water etc.)
      call reservoir_withdrawal(res, mo)

      ! calculate water level in reservoir [m.a.s.l.]
      call reservoir_mean_water_level(Rsv(7, res), rsv_dead_stor_act(res), res)

      ! calculate firm energy yield, release is calculated depending on water level
      rsv_Prod_HPP(res) = RSV_produced_energy(res)

      ! set current reservoir volume for next day
      rsv_B(res) = Rsv(7, res)

      ! if (rsv_active(res) > 0) RsvE(res) = RsvE(res) + rsv_dead_stor_act(res)

      rsv_Day_Cap_Act(res) = rsv_Day_Cap_Act(res) + rsv_dead_stor_act(res) ! incl. dead storage
      rsv_Day_Fill_Min(res) = rsv_Day_Fill_Min(res) + rsv_dead_stor_act(res) ! incl. dead storage

    end if ! if dead storage is empty...
    !*****************************************************************
    !*****************************************************************

    !varoute(1, ih) = precip
    varoute(2, ih) = rsv_Outflow(1, res) * 0.5 ! rsv_frac_sr(res) ! [m^3]
    varoute(8, ih) = rsv_Outflow(1, res) * 0.5 !(1 - rsv_frac_sr(res)) ! [m^3]
    ! Add ground water contribution to gw flow component varoute(8, )
    ! Ground water contribution is computed in subroutine RSV_subbasin
    ! sl 2015-08-18
    !varoute(8, ih) = varoute(8, ih) + (rsv_gwq(res) / 1000. * ResA_wet*10**6) ! mm-->m^3
    varoute(8, ih) = varoute(8, ih) + (rsv_gwq(res) / 1000. * rsv_tot_area(res)) ! mm - - > m^3

    ! set previous day (pd) variables
    pd_outflow(res) = rsv_Outflow(1, res) ![m3]
    !pd_seepage(res) = (rsv_seepage(res)/(ResA_wet*10**6))*1000. ! [mm]
    pd_seepage(res) = (rsv_seepage(res) / (rsv_tot_area(res))) * 1000. ! [mm]
    pd_area_wet(res) = ResA_wet ! [km^2]
    pd_et(res) = ETpot_weighted_mm ! [mm]

    ! #### write output files
    funit = rsv_funit(res) ! get file unit
    write(funit, fmt="(i4, i4)", advance="NO") iyr, ida
    write(funit, fmt="(2f16.2)", advance="NO") (precip/10**3 * rsv_pol_A(res, 20)*10**6)/86400., precip ! m^3/s and mm
    write(funit, fmt="(2f16.2)", advance="NO") cwb/86400., (cwb/(rsv_pol_A(res, 20)*10**6))*10**3 ! m^3/s and mm
    write(funit, fmt="(2f16.4)", advance="NO") ResA_dry, ResA_wet ! km^2
    write(funit, fmt="(f18.2)", advance="NO") Rsv(7, res)/10**6 ! [million m^3]
    write(funit, fmt="(f18.2)", advance="NO") rsv_dead_stor_act(res)/10**6 ! [million m^3]
    write(funit, fmt="(f16.2)", advance="NO") rsv_water_level(res) ! [m.a.s.l.]
    write(funit, fmt="(2f16.2)", advance="NO") rsv_Inflow(res)/86400., rsv_Outflow(1, res)/86400. ! [m^3/s]
    write(funit, fmt="(2f16.2)", advance="NO") rsv_act_withdrawal(res)/86400., rsv_seepage(res)/86400. ! [m^3 and m^3/s]
    write(funit, fmt="(2f16.2)", advance="NO") ETpot_water_mm, ETpot_land_mm ! ETpot water and land [mm]
    write(funit, fmt="(f16.2)", advance="NO") ETpot_weighted_mm ! actual ET [mm]
    write(funit, fmt="(f16.2)", advance="NO") (ETpot_land_m3+ETpot_water_m3)/86400.
    write(funit, fmt="(f16.2)", advance="YES") rsv_Prod_HPP(res) ! produced energy in MW


    !***************************
    ! subroutine RSV_Reservoir_processes
  contains
    ! following subroutines and functions
    !***************************

    !-----------------------------------------------------------------------
    subroutine reservoir_daily_from_monthly(ida, iyr, mo, nc)
      use utilities, only: days_in_month, day_of_month

      integer, intent(in) :: ida
      integer, intent(in) :: iyr
      integer, intent(in) :: mo
      integer, dimension(13), intent(in) :: nc
      ! number of days in current month
      integer :: d_im
      ! day of month
      integer :: d_om
                      ! mo = month (global variable in SWIM)
                      ! iyr = year (global variable in SWIM)

      ! calculate date stuff
      d_im = days_in_month(mo, iyr)
      d_om = day_of_month(ida, mo, iyr, nc)

      ! maximum daily active capacity calculated from monthly values [m^3]
      ! re-calculated in RSV_release_discharge_refilling: (calculation of reservoir re-filling)
      rsv_Day_Cap_Act(res) = (rsv_Cap_Act(res, mo) - rsv_Cap_Act(res, mo - 1)) / d_im * d_om
      rsv_Day_Cap_Act(res) = rsv_Day_Cap_Act(res) + rsv_Cap_Act(res, mo - 1)

      ! daily minimum filling (minimum daily active capacity) calculated from monthly values [m^3]
      ! volume of dead storage is added at the end of the day
      ! possibly re-calculated in calculation of firm energy yield (mgt option 2)
      ! possibly re-calculated in calculation of reservoir release (mgt option 3)
      rsv_Day_Fill_Min(res) = (rsv_Fill_Min(res, mo) - rsv_Fill_Min(res, mo - 1)) / d_im * d_om
      rsv_Day_Fill_Min(res) = rsv_Day_Fill_Min(res) + rsv_Fill_Min(res, mo - 1)

      ! daily minimum filling calculated from monthly values [share of filling]
      !Todo: probably not necessary to do this every day (only once in the beginning)
      rsv_Day_ann_cycle(res) = (rsv_ann_cycle(res, mo) - rsv_ann_cycle(res, mo - 1)) / d_im * d_om
      rsv_Day_ann_cycle(res) = rsv_Day_ann_cycle(res) + rsv_ann_cycle(res, mo - 1)

      !Todo: das könnte wahrscheinlich optimiert werden (vielleicht nicht jeden Tag aufrufen)
      ! [m^3/s]

      ! set mimimal discharges according to status of reservoir (filling or operational)
      if (rsv_active(res) ) then
        rsv_Disch_Min(res, mo) = rsv_Dis_Min_Act(res, mo) ! is operational
      else
        rsv_Disch_Min(res, mo) = rsv_Dis_Min_Fill(res, mo) ! during period of filling
      end if

      ! daily minimum discharge calculated from monthly values [m^3/s]
      ! re-calculated if mgt option > 1
      rsv_Day_Disch_Min(res) = (rsv_Disch_Min(res, mo) - rsv_Disch_Min(res, mo - 1)) / d_im * d_om
      rsv_Day_Disch_Min(res) = rsv_Day_Disch_Min(res) + rsv_Disch_Min(res, mo - 1)

    end subroutine reservoir_daily_from_monthly
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    real(dp) function rsv_calc_inflow(in2, res, sub, ida, varoute) ! returns inflow in [m^3 / d]

      real(dp), dimension(:, :), intent(in) :: varoute
      integer, intent(in) :: in2, res, sub, ida

      RSV_calc_inflow = 0.
      RSV_calc_inflow = varoute(2, in2) + varoute(8, in2) ! [m^3 / d]
      if (RSV_calc_inflow < 0.) then
        RSV_calc_inflow = 0.
        call log_warn("RSV_calc_inflow", "Reservoir inflow is negative!!! (res, sub, day)", &
          ints=(/res, sub, ida/))
      end if
    end function RSV_calc_inflow
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    real(dp) function rsv_calc_seepage(res) ! returns seepage in m^3 / d
      integer, intent(in) :: res
      !Todo: rsv_loss_seepage is probably too high /10^6
      RSV_calc_seepage = (Rsv(1, res) + rsv_dead_stor_act(res)) * rsv_loss_seepage(res)
    end function RSV_calc_seepage
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    subroutine reservoir_remove_seepage(res)
    ! remove seepage from storage(s)
      ! actual reservoir
      integer, intent(in) :: res
      ! volume, m^3
      real(dp) :: vol

      vol = Rsv(1, res) - rsv_seepage(res)
      if (vol > 0.) then
        Rsv(2, res) = vol
      else ! if seepage > actual storage
        Rsv(:, res) = 0.
        ! remove negative seepage value from dead storage
        rsv_dead_stor_act(res) = max(rsv_dead_stor_act(res) + vol, 0.0_dp)
      end if

    end subroutine reservoir_remove_seepage
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    subroutine reservoir_remove_cwb(res)
      ! remove evaporation from storage(s)
      ! actual reservoir
      integer, intent(in) :: res
      ! volume, m^3
      real(dp) :: vol

      vol = Rsv(2, res) + cwb ! cwb is negative!
      if (vol > 0.) then
        Rsv(3, res) = vol
      else ! if evaporation > actual storage
        Rsv(:, res) = 0.
        ! remove negative cwb value from dead storage
        rsv_dead_stor_act(res) = max(rsv_dead_stor_act(res) + vol, 0.0_dp)
      end if

    end subroutine reservoir_remove_cwb
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    subroutine reservoir_add_cwb(res)
    ! remove evaporation from storage(s)
      ! actual reservoir
      integer, intent(in) :: res
      real(dp) :: dx

      if (.NOT. reservoir_is_full_dead_storage(res) ) then
        ! if dead storage is not full
        dx = rsv_dead_stor_capac(res) - (rsv_dead_stor_act(res) + cwb)
        if (dx < 0.) then
          ! add to dead storage
          rsv_dead_stor_act(res) = rsv_dead_stor_capac(res)
          ! add to active storage
          Rsv(3, res) = min(Rsv(2, res) + dx, rsv_Day_Cap_Act(res))
        else
          rsv_dead_stor_act(res) = rsv_dead_stor_act(res) + cwb
        end if

      else ! if dead storage was already full
        ! add cwb to reservoir storage
        Rsv(3, res) = min(Rsv(2, res) + cwb, rsv_Day_Cap_Act(res))
      end if

    end subroutine reservoir_add_cwb
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    subroutine reservoir_outflow ! [m^3 / d]
      !real(dp) :: Res_free ! [m^3]
      ! Volume of reservoir active capacity not filled [m^3]
      real(dp) :: Res_space
      !real(dp) :: dx ! The volume [m^3/d] that should be released from the reservoir if Rsv > rsv_Day_Cap_Act
      ! Volume of water surpassing Reservoir capacity [m^3 / d]
      real(dp) :: Rsv_Spill

      !       rsv_Outflow(1, res) = 0.
      !       rsv_Outflow(2, res) = 0.
      !       dx = 0.
      !       Res_free = rsv_Day_Cap_Act(res) - Rsv(2, res)
      !
      !       if (cwb > 0.) then
      !          ! if cwb positive, add to reservoir storage
      !          Rsv(3, res) = min(Rsv(2, res) + cwb, rsv_Day_Cap_Act(res))
      !          dx = cwb - Res_free
      !          if (dx > 0.) then
      !             ! determine the volume of water that must be released in order to reduce Rsv to daily maximum active capacity
      !             rsv_Outflow(1, res) = dx ! simulated reservoir discharge
      !             rsv_Outflow(2, res) = dx ! available reservoir discharge
      !          end if
      !       else
      !          ! if cwb < 0. remove evaporated volume from actual storage volume
      !          Rsv(3, res) = max(Rsv(2, res) + cwb, 0.)
      !       end if

      ! calculate reservoir volume (not filled)
      Res_space = rsv_Day_Cap_Act(res) - Rsv(3, res)
      if (Res_space < 0.) then
        rsv_Outflow(1, res) = rsv_Outflow(1, res) + Res_space * (- 1.) + rsv_Inflow(res)
        rsv_Outflow(2, res) = rsv_Outflow(2, res) + Res_space * (- 1.) + rsv_Inflow(res)
        Rsv(4, res) = rsv_Day_Cap_Act(res)
      else
        ! Calculation of reservoir spill
        Rsv_Spill = rsv_Inflow(res) - Res_space
        if (Rsv_Spill > 0.) then
          rsv_Outflow(1, res) = rsv_Outflow(1, res) + Rsv_Spill
          rsv_Outflow(2, res) = rsv_Outflow(2, res) + Rsv_Spill
          Rsv(4, res) = rsv_Day_Cap_Act(res)
        else
          Rsv(4, res) = Rsv(3, res) + rsv_Inflow(res)
        end if
      end if

    end subroutine reservoir_outflow
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    subroutine reservoir_release_refilling
      ! Volume of water available in the Reservoir [m^3 / d]
      real(dp) :: release
      ! [m^3]
      real(dp) :: Res_free
      ! Volume of water surpassing Reservoir capacity [m^3]
      real(dp) :: Rsv_Spill
      ! additional active capacity [m3]
      real(dp) :: Rsv_Add_Cap

      ! calculate reservoir release
      release = 0.
      release = release + Rsv(4, res) - rsv_Day_Fill_Min(res)
      if (release > 0.) then
        rsv_Outflow(1, res) = rsv_Outflow(1, res) + release
        rsv_Outflow(2, res) = rsv_Outflow(2, res) + release
        Rsv(5, res) = max(Rsv(4, res) - release, 0.0_dp)
      else
        Rsv(5, res) = Rsv(4, res)
      end if

      ! calculate discharge downstream
      d_Q = rsv_Outflow(1, res) - max(rsv_Day_Disch_Min(res) * 86400.0_dp, 0.0_dp)
      ! calculate reservoir re-filling
      ! Rsv(5, res) = actual storage volume where release (for hydropower or minimal requirements downstream) is already abstracted
      rsv_Outflow(2, res) = 0.
      if (d_Q > 0.) then
        rsv_Outflow(1, res) = max(rsv_Day_Disch_Min(res) * 86400.0_dp, 0.0_dp)
        Res_free = Rsv(5, res) - rsv_Day_Cap_Act(res)
        Rsv_Spill = d_Q + Res_free
        if (Rsv_Spill > 0.) then
          release = rsv_Outflow(1, res) + Rsv_Spill
          Rsv_Spill = MIN(release, rsv_cap_hpp(res) * 86400.0_dp) ! rsv_cap_hpp is in [m^3 / s]
          Res_space = release - Rsv_Spill
          Rsv_Spill = release - rsv_Outflow(1, res)
          if (rsv_Day_Cap_Act(res) < (rsv_Capac_Max(res) - rsv_dead_stor_act(res)) ) then
            if (Res_space > 0.) then
              release = MAX(((rsv_cap_hpp(res) * 86400.0_dp) - rsv_Outflow(1, res)), 0.0_dp)
              Res_space = MAX((Rsv_Spill - release), 0.0_dp)
              Rsv_Add_Cap = rsv_Capac_Max(res) - rsv_dead_stor_act(res) - rsv_Day_Cap_Act(res)
              Rsv_Add_Cap = MIN(Res_space, Rsv_Add_Cap)
              rsv_Day_Cap_Act(res) = rsv_Day_Cap_Act(res) + Rsv_Add_Cap
              Rsv_Spill = MAX((Rsv_Spill - Rsv_Add_Cap), 0.0_dp)
            end if
          end if
          rsv_Outflow(1, res) = rsv_Outflow(1, res) + Rsv_Spill
          rsv_Outflow(2, res) = rsv_Outflow(2, res) + Rsv_Spill
          Rsv(6, res) = Rsv_Day_Cap_Act(res)
        else ! if (Rsv_Spill > 0.)
          Rsv(6, res) = Rsv(5, res) + d_Q
        end if ! if (Rsv_Spill > 0.)
      else ! if (d_Q > 0.)
        Rsv(6, res) = Rsv(5, res)
      end if ! if (d_Q > 0.)

    end subroutine reservoir_release_refilling
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    subroutine reservoir_withdrawal(res, mo)
      integer, intent(in) :: mo
      ! calculate direct water withdrawal from reservoir [m^3/d]
      ! (drinking water abstraction, irrigation etc.)
      ! reservoir number
      integer, intent(in) :: res
      ! TYPE (TSubbasin), POINTER :: pS ! water management module, pointer on subbasin

      rsv_act_withdrawal(res) = 0.

      !rsv_Day_Fill_Min = the actual filling (current release already removed)

      ! if (rsv_Mngmt(res) > 1) then ! => also valid for Option 1, changed H. Koch, 09.07.2019
      rsv_Day_Fill_Min(res) = rsv_Day_Fill_Min(res) - (Rsv(6, res) * rsv_shr_withdr(res))
      rsv_Day_Fill_Min(res) = max(rsv_Day_Fill_Min(res), 0.0_dp)
      ! end if ! => also valid for Option 1, changed H. Koch, 09.07.2019

      Res_space = Rsv(6, res) - max(rsv_Day_Fill_Min(res), 0.0_dp)
      Res_space = max(Res_space, 0.0_dp)

      !       !#################################
      !       !#### WATER MANAGEMENT MODULE ####
      !       !#################################
      !       if ( bWAM_Module ) then
      !          ! check if water transfers take place in current subbasin
      !          if ( wam_is_transfer_subbasin(sub) ) then
      !             pS => wam_get_pointer_subbasin(sub) ! get pointer on current subbasin
      !             !!!! overwrite withdrawals from reservoir input file!!!
      !             rsv_Withdr_Mon(res, mo) = pS%totDemand(daycounter) ! [m3/s]
      !          end if
      !       end if
      !       !#################################

      if (Res_space > (rsv_Withdr_Mon(res, mo) * 86400.) ) then
        Rsv(7, res) = Rsv(6, res) - (rsv_Withdr_Mon(res, mo) * 86400.)
        rsv_act_withdrawal(res) = rsv_Withdr_Mon(res, mo) * 86400.
      else
        d_Q = rsv_Outflow(2, res)
        if ((Res_space + d_Q) < (rsv_Withdr_Mon(res, mo) * 86400.) ) then
          rsv_Outflow(1, res) = rsv_Outflow(1, res) - d_Q
          rsv_Outflow(2, res) = rsv_Outflow(2, res) - d_Q
          Rsv(7, res) = Rsv(6, res) - Res_space
          rsv_act_withdrawal(res) = Res_space + d_Q
        else
          d_Q = (rsv_Withdr_Mon(res, mo) * 86400.) - Res_space
          rsv_Outflow(1, res) = rsv_Outflow(1, res) - d_Q
          rsv_Outflow(2, res) = rsv_Outflow(2, res) - d_Q
          Rsv(7, res) = Rsv(6, res) - Res_space
          rsv_act_withdrawal(res) = Res_space + d_Q
        end if
      end if

      !       !#################################
      !       !#### WATER MANAGEMENT MODULE ####
      !       !#################################
      !       if ( bWAM_Module ) then
      !          ! check if water transfers take place in current subbasin
      !          if ( wam_is_transfer_subbasin(sub) ) then
      !             pS => wam_get_pointer_subbasin(sub) ! get pointer on current subbasin
      !             ! assign actual withdrawal to subbasins' water user(s)
      !             if ( rsv_act_withdrawal(res) > 0. ) &
      !                call wam_distrib_water_frm_rsrvoir (pS, rsv_act_withdrawal(res)) ! [m3]
      !          end if
      !       end if
      !       !#################################

    end subroutine reservoir_withdrawal
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    subroutine reservoir_mean_water_level(vol_act, vol_dead, res)
      ! calculate water level in [m.a.s.l.]
      integer, intent(in) :: res
      ! volume of active storage
      real(dp), intent(in) :: vol_act
      ! volume of dead storage
      real(dp), intent(in) :: vol_dead
      real(dp) :: filling

      filling = 0.
      filling = (vol_act + rsv_B(res)) / 2. ! mean of actual storage t = 0 and t = - 1
      SH = RSV_POL(20, rsv_pol_V(res, :), rsv_pol_L(res, :), filling + vol_dead)
      SH = max(SH, rsv_level_min(res))
      rsv_water_level(res) = min(rsv_level_max(res), SH)
      rsv_height_hpp(res) = min(rsv_water_level(res) - rsv_level_min(res), rsv_level_hpp(res))

    end subroutine reservoir_mean_water_level
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    real(dp) function rsv_produced_energy(res)
      integer, intent(in) :: res
      real(dp) :: release_hpp

      RSV_produced_energy = 0.
      release_hpp = 0.

      release_hpp = min(rsv_Outflow(1, res) / 86400., rsv_cap_hpp(res))
      RSV_produced_energy = rsv_height_hpp(res) * release_hpp * 9.81 * 1.0 * rsv_eff_hpp(res) / 1000. ! Efficiency of HPP; included 2015 - 07 - 07 (H. Koch)
      !     MW = H [m] * Q [m³/s] * g [m/s]* ρ [kg/m³] * η [-] / 1000.
    end function RSV_produced_energy
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    real(dp) function rsv_release_firm_energy_yield(res)
      ! Reservoir management option 2
      ! Release for firm energy yield depends on water level.
      ! Returns reservoir release in [m3/d]
      ! reservoir number
      integer, intent(in) :: res
      ! [m.a.s.l.]
      real(dp) :: water_level
      ! [m3 / s]
      real(dp) :: release_hpp

      water_level = 0.
      RSV_release_firm_energy_yield = 0.
      release_hpp = 0.
      SH = 0.
      ! water level [m.a.s.l.]
      SH = RSV_POL(20, rsv_pol_V(res, :), rsv_pol_L(res, :), Rsv(4, res) + rsv_dead_stor_act(res))
      water_level = min(rsv_level_max(res), SH)

      ! [m3/s]
      release_hpp = RSV_POL(20, rsv_pol_L2(res, :), rsv_pol_HP(res, :), water_level)

      ! Maximum Release depending on Reservoir filling
      release_hpp = MIN(release_hpp, rsv_cap_hpp(res)) ! [m3 / s]

      ! return [m3/d]
      RSV_release_firm_energy_yield = min(release_hpp * 86400., Rsv(4, res) * rsv_Day_ann_cycle(res))

    end function RSV_release_firm_energy_yield
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    real(dp) function rsv_release_dep_on_filling(res)
      ! Reservoir management option 3
      ! Release depends on filling.
      ! reservoir number
      integer, intent(in) :: res
      ! [m.a.s.l.]
      real(dp) :: water_level
      ! [m3 / d]
      real(dp) :: release

      water_level = 0.
      release = 0.
      RSV_release_dep_on_filling = 0.

      ! water level [m.a.s.l.]
      SH = RSV_POL(20, rsv_pol_V(res, :), rsv_pol_L(res, :), Rsv(4, res) + rsv_dead_stor_act(res))
      water_level = min(rsv_level_max(res), SH)

      release = RSV_POL(20, rsv_pol_L2(res, :), rsv_pol_HP(res, :), water_level) ![m3 / s]

      ! release [m3/d]
      RSV_release_dep_on_filling = release * 86400.

    end function RSV_release_dep_on_filling
    !-----------------------------------------------------------------------

  end subroutine reservoir_process