reservoir_subbasin Subroutine

  subroutine reservoir_subbasin(in1, da, aff, dart, flu, neap, precip, qtl, sbp, sda, smq, smsq, snowval, sub, subp, sumcn, susb, tx, varoute, wysb, xeo, xpercn, xqi, xsep, xssf, xssfn, xswind, xyno3, xysp, yd, yon, yph)
    ! CALLED FROM: mainpro.f90 during subbasin cycle in routing scheme.
    ! This is a modified version of the original 'subbasin' subroutine in SWIM.
    ! It was basically implemented to maintain output statistics.
    ! NOTE: reservoir subbasin estimates should not be added in subroutine add (route.f90)
    use evapotranspiration, only : eta_output_id, etp_output_id, ra, tmn
    use groundwater, only : &
      abf1, &
      baseflow_output_id, &
      delay, &
      groundwater_recharge_output_id, &
      gwchrg, &
      gwq, &
      gwseep, &
      rchrg, &
      revap, &
      revapst, &
      xet
    use hydrotope, only : &
      precipitation_output_id, &
      tmax_output_id, &
      tmin_output_id, &
      water_yield_output_id
    use output, only : &
      output_store_hydrotope_value, &
      output_store_subbasin_value, &
      output_store_subbasin_values
    use snow, only : bSnowModule, sml, snowfall_weq_output_id, tmx, tsnfall
    use soil, only : sep, surface_runoff_output_id, xqd

    real(dp), intent(in) :: da
    real(dp), intent(out) :: aff
    real(dp), dimension(:), intent(in) :: dart
    real(dp), dimension(:), intent(in) :: flu
    integer, dimension(:), intent(in) :: neap
    real(dp), intent(inout) :: precip
    real(dp), intent(inout) :: qtl
    real(dp), dimension(:), intent(inout) :: sbp
    real(dp), dimension(:, :), intent(inout) :: sda
    real(dp), dimension(:), intent(inout) :: smq
    real(dp), dimension(:), intent(inout) :: smsq
    real(dp), intent(out) :: snowVal
    real(dp), dimension(30), intent(inout) :: sub
    real(dp), dimension(:), intent(in) :: subp
    real(dp), intent(in) :: sumcn
    real(dp), dimension(:, :), intent(inout) :: susb
    real(dp), dimension(:), intent(in) :: tx
    real(dp), dimension(:, :), intent(out) :: varoute
    real(dp), intent(inout) :: wysb
    real(dp), intent(inout) :: xeo
    real(dp), intent(inout) :: xpercn
    real(dp), intent(out) :: xqi
    real(dp), intent(inout) :: xsep
    real(dp), intent(inout) :: xssf
    real(dp), intent(inout) :: xssfn
    real(dp), intent(inout) :: xswind
    real(dp), intent(inout) :: xyno3
    real(dp), intent(inout) :: xysp
    real(dp), intent(inout) :: yd
    real(dp), intent(inout) :: yon
    real(dp), intent(inout) :: yph
    integer, intent(in) :: in1
    integer :: i, j, res, hru
    ! = rchrg of previous time step
    real(dp) :: rchrg1
    real(dp) :: SH, A, ResA_dry, ResA_wet
    real(dp) :: out = 0.

    j = in1
    res = reservoir_get(j)

    aff = da * 1000. * flu(j)

    !**** COMPUTE TOTALS for precip:
    precip = subp(j)
    susb(1, j) = susb(1, j) + precip
    sub(1) = sub(1) + precip * flu(j)
    sbp(j) = sbp(j) + precip

    !#### COMPUTE snow, recalc precip
    !     ATTN: Snow accumulation, melt, and snow evaporation
    !        are calculated at hydrotope level!
    snowVal = 0.
    if (tx(j) .le. tsnfall(j)) then
      snowVal = precip
      sub(2) = sub(2) + precip * flu(j)
      call output_store_subbasin_value(snowfall_weq_output_id, j, precip)
      susb(2, j) = susb(2, j) + precip
      precip = 0.
    end if
    susb(3, j) = susb(3, j) + sml
    sub(3) = sub(3) + sml * flu(j)

    !#### Hydrotope loop is not called for reservoir subbasins
    !#### This means that hydrotope output files are not written if a
    !#### reservoir hydrotope was chosen for this subbasin.

    !**** COMPUTE TOTALS for ground water: sub(15...19)
    !### sl 2013-11-08
    !call gwmod(j, pd_seepage(res), pd_et(res))
    pd_seepage(res) = max(0.0_dp, pd_seepage(res))
    rchrg1 = rchrg(j)
    rchrg(j) = pd_seepage(res)
    gwseep = pd_seepage(res) * (1. - rsv_gwc(res)) !rchrgc(j)
    rchrg(j) = (1. - delay(j)) * rchrg(j) + delay(j) * rchrg1
    if (rchrg(j) < 1.e-6) rchrg(j) = 0.
    gwchrg = rchrg(j)
    gwq(j) = gwq(j) * abf1(j) + (rchrg(j)) * (1. - abf1(j))

    !  Calc revap storage
    revapst(j) = revapst(j) + rchrg(j)

    ! remove seepage to deep aquifer from revap storage (shallow aquifer)
    revapst(j) = revapst(j) - gwseep
    if (revapst(j) < 0. ) then
      gwseep = gwseep + revapst(j)
      revapst(j) = 0.
    end if

    !! remove ground water flow from shallow aquifer storage
    revapst(j) = revapst(j) - gwq(j)
    if (revapst(j) < 0. ) then
      gwq(j) = gwq(j) + revapst(j)
      revapst(j) = 0.
    end if

    !! Ground water contribution is added to varoute(8, ) in RSV_Reservoir_processes
    rsv_gwq(res) = max(0.0_dp, gwq(j))


    !#### values of previous day are being used here !!!
    ! 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)
    ! sl 2015-08-18
    !rsv_tot_area(res) = (ResA_wet+ResA_dry)*10**6 ! km2 --> m2 total resevoir area

    !! pd_outflow is the total outflow from the reservoir in [mm]
    !! For surface and subsurface (subbasin) output it must be corrected on ground water contribution
    !! Otherwise, some water balance output files are accounting double
    !! The half is distributed to surface runoff and the other half to subsurface
    ! sl 2015-08-18
    !upstream_area = rsv_upstream_area(res)*10**6 - tot_area(res)
    !if ( upstream_area <= 0. ) upstream_area = tot_area(res) ! if reservoir is a headwater
    !   out = pd_outflow(res) - rsv_gwq(res) - (rsv_Inflow(res) / upstream_area * 1000.) ! mm
    out = (pd_outflow(res) - rsv_Inflow(res)) / rsv_tot_area(res) * 1000. ! mm
    xqd = out * 0.5 ! rsv_frac_sr(res) ! surface runoff in [mm]
    xssf = out * 0.5 !(1 - rsv_frac_sr(res)) ! subsurface runoff in [mm]
    !    end if

    qtl = 0. ! transmission losses in [mm]
    xeo = pd_et(res) ! pot. evapotranspiration in [mm]
    xet = xeo ! act. evapotranspiration in [mm]

    do hru = 1, neap(j) !nbr_hru_sub(j)
      !**** CALC variables for GIS output
      call output_store_hydrotope_value(precipitation_output_id, j, hru, precip)
      call output_store_hydrotope_value(surface_runoff_output_id, j, hru, xqd + xssf)
      call output_store_hydrotope_value(groundwater_recharge_output_id, j, hru, sep)
      call output_store_hydrotope_value(eta_output_id, j, hru, pd_et(res))
      call output_store_hydrotope_value(etp_output_id, j, hru, pd_et(res))

    end do
    !### sl 2013-11-08

    !####

    !**** COMPUTE TOTALS for water: sub(4...13)
    !### sl
    !! flu(j) replaced by tot_area / da

    sub(4) = sub(4) + tmx(j) * flu(j) ! daily max temp. for subbasin, readcli, degree C
    call output_store_subbasin_value(tmax_output_id, j, tmx(j))
    sub(5) = sub(5) + tmn(j) * flu(j) ! daily min temp. for subbasin, readcli, degree C
    call output_store_subbasin_value(tmin_output_id, j, tmn(j))
    sub(6) = sub(6) + ra(j) * flu(j) ! solar radiation
    sub(7) = sub(7) + sumcn * flu(j) ! weighted aver. Cur.Num. in subbasin (SUM(cn))
    susb(8, j) = susb(8, j) + xqd
    sub(8) = sub(8) + xqd * flu(j)
    smq(j) = smq(j) + xqd
    susb(9, j) = susb(9, j) + xssf
    sub(9) = sub(9) + xssf * flu(j)
    smsq(j) = smsq(j) + xssf
    susb(10, j) = susb(10, j) + qtl
    sub(10) = sub(10) + qtl * flu(j)
    revapst(j) = revapst(j) + qtl

    xsep = pd_seepage(res) ! percolation [mm]
    susb(11, j) = susb(11, j) + xsep
    sub(11) = sub(11) + xsep * flu(j)
    susb(12, j) = susb(12, j) + xeo !pot. evapotranspiration
    sub(12) = sub(12) + xeo * flu(j)
    !  actual is equal to potential evapotranspiration in reservoir subbasin
    susb(13, j) = susb(13, j) + xet
    sub(13) = sub(13) + xet * flu(j)

    susb(15, j) = susb(15, j) + gwq(j)
    sub(15) = sub(15) + gwq(j) * flu(j)
    call output_store_subbasin_value(baseflow_output_id, j, gwq(j))
    susb(16, j) = susb(16, j) + 0. !revap
    sub(16) = sub(16) + revap * flu(j) ! revap * tot_area / da
    susb(17, j) = susb(17, j) + gwseep
    sub(17) = sub(17) + gwseep * flu(j)
    susb(18, j) = susb(18, j) + gwchrg
    sub(18) = sub(18) + gwchrg * flu(j)
    xswind = 1.
    susb(19, j) = susb(19, j) + xswind
    sub(19) = sub(19) + xswind * flu(j)
    xxswind = xxswind + xswind * flu(j)

    !**** COMPUTE TOTALS for water yield (wysb) & sediments: sub(20...25)
    !wysb = xqi + xssf + gwq(j) - qtl
    !wysb = (pd_outflow(res) + gwq(j)) ! - (rsv_Inflow(res)*1000./(rsv_pol_A(res, 20)*10**6)) ! [mm]
    ! sl 2015-08-18
    !wysb = pd_outflow(res)
    ! wysb can be negative if inflow is larger than outflow + gwq
    wysb = pd_outflow(res) / rsv_tot_area(res) * 1000. + rsv_gwq(res) - qtl
    !wysb = max(wysb, 0.)

    susb(20, j) = susb(20, j) + wysb
    sub(20) = sub(20) + wysb * flu(j)
    call output_store_subbasin_value(water_yield_output_id, j, wysb)
    !Todo: sediment
    !Todo: water quality parameters
    !    sub(21) = sub(21) + yd/(100.*da*tot_area / da)
    !    sym(j) = sym(j) + yd
    !    susb(21, inum1) = susb(21, inum1) + yd/(100.*da*tot_area / da)
    !    sub(22) = sub(22) + yon * tot_area / da
    !    susb(22, inum1) = susb(22, inum1) + yon
    !    sub(23) = sub(23) + yph * tot_area / da
    !    susb(23, inum1) = susb(23, inum1) + yph
    !    sub(24) = sub(24) + ysp * tot_area / da
    !    susb(24, inum1) = susb(24, inum1) + ysp


    !**** COMPUTE SUBBASIN OUTPUTS FOR ROUTING: sda(), varoute()
    !     ATTN: sda(6, j) = sum of 3 fluxes after retention, new version
    !     ATTN: coef (dart()*1000) to transform from mm to m3 (42, 48)
    !     ATTN: coef (dart()*100) to transform kg/ha to kg (44-47)
    !     ATTN: wysb(mm)*dart(km2)*1000 = wysb*dart*1000 (m3)
    !     ATTN: sda(2, j), sda(8, j) in m3
    !     ATTN: xyno3 in kg/ha, sda(6) & varoute(6, .) in kg

    xqi = xqd !pd_outflow(res) * 0.5 ! rsv_frac_sr(res) ! [mm]
    yd = 0.
    yon = 0.
    yph = 0.
    xyno3 = 0.
    xssfn = 0.
    xpercn = 0.
    xysp = 0.
    !  * dart * 1000. replaced by tot_area/1000.
    !sda(2, j) = xqi * dart(j) * 1000.
    ! sl 2015-08-18
    !sda(2, j) = xqi * rsv_tot_area(res) / 1000.
    sda(2, j) = pd_outflow(res) * 0.5 / rsv_tot_area(res) * 1000.

    sda(3, j) = yd
    sda(4, j) = yon * dart(j) * 100
    sda(5, j) = yph * dart(j) * 100
    sda(6, j) = (xyno3 + xssfn + xpercn) * dart(j) * 100
    sda(7, j) = xysp * dart(j) * 100
    ! sl 2015-08-18
    !sda(8, j) = (wysb - xqi) * rsv_tot_area(res) / 1000.
    sda(8, j) = (wysb - pd_outflow(res) * 0.5) / rsv_tot_area(res) * 1000.
    xwysb = xwysb + wysb / rsv_tot_area(res) * 1000.

    do i = 1, 8
      varoute(i, j) = sda(i, j)
    end do

  end subroutine reservoir_subbasin