groundwater.f95 Source File

! FILE gwat.f
!
! SUBROUTINES IN THIS FILE CALLED FROM
! subroutine gwmod(j) subbasin
module groundwater

  use utilities, only : dp

  implicit none

  ! subsurface flow
  real(dp), save, dimension(:), allocatable :: gw_bff ! baseflow factor for basin, used to cal
  ! groundwater height, mm
  real(dp), save, dimension(:), allocatable :: gw_gwht
  ! initial groundwater flow contribution to streamflow, mm
  real(dp), save, dimension(:), allocatable :: gw_gwq
  ! alpha factor for groundwater
  real(dp), save, dimension(:), allocatable :: gw_abf
  ! specific yield for g - w
  real(dp), save, dimension(:), allocatable :: gw_syld
  ! alpha factor for groundwater
  real(dp), save, dimension(:), allocatable :: gw_delay ! groundwater delay, days
  ! fraction of root zone perc that goes to revap
  real(dp), save, dimension(:), allocatable :: gw_revapc
  ! fraction (0 - 1) of root zone perc that percolates
  real(dp), save, dimension(:), allocatable :: gw_rchrgc
  ! revap storage, mm
  real(dp), save, dimension(:), allocatable :: gw_revapmn
  ! evapotranspiration (actual), mm
  real(dp), save :: xet
  real(dp), save :: gwseep = 0. ! groundwater seepage, mm
  real(dp), save :: revap = 0. ! revaporation from groundwater, mm
  real(dp), save :: gwchrg = 0. ! groundwater recharge, mm
  ! groundwater height, mm
  real(dp), save, dimension(:), allocatable :: gwht
  ! groundwater contribution to stream, in mm
  real(dp), save, dimension(:), allocatable :: gwq
  ! alpha factor for g - w = reaction factor * delta(T)
  real(dp), save, dimension(:), allocatable :: abf
  ! specific yield for g - w
  real(dp), save, dimension(:), allocatable :: syld
  ! groundwater delay, days
  real(dp), save, dimension(:), allocatable :: delay
  ! fraction of root zone perc that goes to revap
  real(dp), save, dimension(:), allocatable :: revapc
  ! fraction (0 - 1) of root zone perc that percolates
  real(dp), save, dimension(:), allocatable :: rchrgc
  ! revap storage, mm
  real(dp), save, dimension(:), allocatable :: revapmn
  ! exp function of alpha factor for groundwater (abf)
  real(dp), save, dimension(:), allocatable :: abf1
  ! groundwater recharge to the aquifer, mm
  real(dp), save, dimension(:), allocatable :: rchrg
  !transmission losses in the reach, calc in route, m3
  real(dp), save, dimension(:), allocatable :: revapst
  ! initial groundwater flow contribution to streamflow, mm (1 for basin)
  real(dp), save :: gwq0 = 0.2
  ! alpha factor for groundwater
  real(dp), save :: abf0 = 0.5
  !
  real(dp), save, dimension(:), allocatable :: abf2
  !
  real(dp), save, dimension(:), allocatable :: gwqs
  !
  real(dp), save, dimension(:), allocatable :: gwqd
  real(dp), save, dimension(:), allocatable :: additionalGwUptake
  real(dp), save, dimension(:), allocatable :: bff
  ! baseflow factor perennial rivers = 0, ephemeral rivers = 1.
  real(dp) :: bff0 = 0.


  ! Output variable ids
  integer :: &
    groundwater_recharge_output_id = 0, &
    baseflow_output_id = 0

  namelist / GROUNDWATER_PARAMETERS / &
    abf0, &
    gwq0, &
    bff0

contains

  subroutine groundwater_initialise(mb, subbasin_input_file_id)
    use input, only: get_config_fid
    use output, only: output_register_hydrotope_var, output_register_subbasin_var
    integer, intent(in) :: mb, subbasin_input_file_id

    read(get_config_fid(), GROUNDWATER_PARAMETERS)
    groundwater_recharge_output_id = output_register_hydrotope_var("groundwater_recharge", .false.)
    baseflow_output_id = output_register_subbasin_var("baseflow", .false.)

    call groundwater_allocate(mb)
    call groundwater_read_input(subbasin_input_file_id)
  end subroutine groundwater_initialise

  subroutine groundwater_correct_params
    ! Called in swim:initialised after the subcatch parameters have been read
    abf1 = exp(- abf)
    delay = exp(-1. / (delay + 1.e-6))
    rchrg = gwq
    revapmn = - 1. * revapmn
  end subroutine groundwater_correct_params

  subroutine groundwater_allocate(mb)
    integer, intent(in) :: mb
    allocate(abf(mb))
    allocate(abf1(mb))
    allocate(abf2(mb))
    allocate(additionalGwUptake(mb))
    allocate(bff(mb))
    allocate(delay(mb))
    allocate(gwht(mb))
    allocate(gwq(mb))
    allocate(gwqd(mb))
    allocate(gwqs(mb))
    allocate(rchrg(mb))
    allocate(rchrgc(mb))
    allocate(revapc(mb))
    allocate(revapmn(mb))
    allocate(revapst(mb))
    allocate(syld(mb))
    abf = abf0
    abf1 = 0.
    abf2 = 0.
    additionalGwUptake = 0.
    bff = bff0
    delay = 200.0
    gwht = 1.
    gwq = gwq0
    gwqd = 0.
    gwqs = 0.
    rchrg = 0.
    rchrgc = 0.05
    revapc = 0.200
    revapmn = 0.
    revapst = 0.
    syld = 0.003
  end subroutine groundwater_allocate

  subroutine groundwater_allocate_subcatch(n_subcatch)
    integer, intent(in) :: n_subcatch

    allocate(gw_gwht(n_subcatch))
    gw_gwht = 0.
    allocate(gw_gwq(n_subcatch))
    gw_gwq = 0.
    allocate(gw_abf(n_subcatch))
    gw_abf = 0.
    allocate(gw_syld(n_subcatch))
    gw_syld = 0.
    allocate(gw_delay(n_subcatch))
    gw_delay = 0.
    allocate(gw_revapc(n_subcatch))
    gw_revapc = 0.
    allocate(gw_rchrgc(n_subcatch))
    gw_rchrgc = 0.
    allocate(gw_revapmn(n_subcatch))
    gw_revapmn = 0.
    allocate(gw_bff(n_subcatch))
    gw_bff = 0.
  end subroutine groundwater_allocate_subcatch


  subroutine groundwater_read_input(subbasin_input_file_id)
    use input, only : read_real_column
    integer, intent(in) :: subbasin_input_file_id

    call read_real_column(subbasin_input_file_id, "gwht", gwht, 1.0_dp)
    call read_real_column(subbasin_input_file_id, "gwq", gwq, 0.5_dp)
    call read_real_column(subbasin_input_file_id, "abf", abf, 0.048_dp)
    call read_real_column(subbasin_input_file_id, "delay", delay, 200.0_dp)
    call read_real_column(subbasin_input_file_id, "revapc", revapc, 0.2_dp)
    call read_real_column(subbasin_input_file_id, "syld", syld, 0.003_dp)
    call read_real_column(subbasin_input_file_id, "rchrgc", rchrgc, 0.05_dp)
    call read_real_column(subbasin_input_file_id, "revapmn", revapmn, 0.0_dp)

  end subroutine groundwater_read_input


  subroutine dealloc_groundwater
   deallocate(abf)
   deallocate(abf1)
   deallocate(abf2)
   deallocate(additionalGwUptake)
   deallocate(delay)
   deallocate(gwqd)
   deallocate(gwqs)
   deallocate(rchrg)
   deallocate(rchrgc)
   deallocate(revapc)
   deallocate(revapmn)
   deallocate(revapst)
  end subroutine dealloc_groundwater

  subroutine groundwater_process(j, percolation, ETact, ETpot)
    !**** PURPOSE: TO ESTIMATE GROUNDWATER CONTRIBUTION TO STREAMFLOW
    !**** CALLED IN: SUBBASIN
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    ! PARAMETERS & VARIABLES
    !
    !  >>>>> COMMON PARAMETERS & VARIABLES
    !  abf(j) = alpha factor for g-w = reaction factor * delta(T)
    !  abf1(j) = exp fun of alpha factor: abf1(j) = exp(-abf(j))
    !  delay(j) = groundwater delay, days
    !  gwchrg = groundwater recharge, mm
    !  gwht(j) = groundwater height, mm
    !  gwq(j) = groundwaterw contribution to stream, in mm
    !  gwseep = groundwater seepage, mm
    !  rchrg(j) = groundwater recharge to the aquifer, mm
    !  rchrgc(j) = fraction of root zone perc that goes into deep g-w
    !  revap = revaporation from groundwater, mm
    !  revapc(j) = fraction of root zone perc that goes to revap
    !  revapst(j) = revap storage, mm
    !  syld(j) = specific yield for g-w
    !  xet = evapotranspiration (actual), mm
    !  xsep = subbasin percolation, mm
    !  >>>>>

    !  >>>>> STATIC PARAMETERS
    !  rchrg1 = groundwater recharge to the aquifer
    !  >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    integer, intent(in) :: j
    ! [mm]
    real(dp), intent(in) :: percolation, ETact, ETpot
    ! = rchrg of previous time step
    real(dp) :: rchrg1

    rchrg1 = rchrg(j)

    !  Calc recharge
    rchrg(j) = percolation !xsep

    !  Calc g-w seepage
    gwseep = rchrg(j) * rchrgc(j)

    !! compute revap to soil profile/plant roots
    revap = revapc(j) * ETpot ! according to SWAT revap depends on ETpot not ETact
    if (revap > (ETpot - ETact) ) revap = ETpot - ETact
    revap = max(0.0_dp, revap)

    !! compute gw recharge level
    rchrg(j) = (1. - delay(j)) * rchrg(j) + delay(j) * rchrg1
    if (rchrg(j) < 1.e-6) rchrg(j) = 0.

    gwchrg = rchrg(j)

    !  Calc g-w height
    gwht(j) = gwht(j) * abf1(j) + rchrg(j) * (1. - abf1(j)) &
              / (800. * syld(j) + 1.e-6 * abf(j) + 1.e-6)
    gwht(j) = max(real(1.e-6, dp), gwht(j))

    !  Calc g-w contribution to streamflow
    gwq(j) = gwq(j) * abf1(j) + (rchrg(j)) * (1. - abf1(j))

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

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

    ! 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 revap to soil profile (capillary rise) from shallow aquifer storage
    if (revapst(j) < revapmn(j) ) then
      revap = 0.
    else
      revapst(j) = revapst(j) - revap
      if (revapst(j) < revapmn(j) ) then
        revap = revap + revapst(j) - revapmn(j)
        revapst(j) = revapmn(j)
      end if
    end if
    ! add revap to actual evapotranspiration otherwise it does not occur in water balance output
    xet = max(0.0_dp, ETact + revap)

  end subroutine groundwater_process

end module groundwater