! 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 :: disLastDay real(dp), save, dimension(:), allocatable :: bff ! baseflow factor perennial rivers = 0.1, ephemeral rivers = 1. real(dp) :: bff0 = 0.1 ! 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(disLastDay(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. disLastDay = 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(disLastDay) 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