subroutine gwflow_pump_allo(ob_id_dmd,demand_vol,extracted,dmd_unmet) !rtb gwflow !! ~ ~ ~ PURPOSE ~ ~ ~ !! this subroutine determines the volume of groundwater that is extracted !! from gwflow grid cells !! (pumping volume are used in gwflow_simulate, in groundwater balance equations) !! adapted for swatplusGW water allocation (4-arg call from wallo_withdraw) use gwflow_module use organic_mineral_mass_module, only : soil1 use hru_module, only : hru,irrn,irrp use hydrograph_module, only : wet,ob use soil_module, only : soil use constituent_mass_module use res_salt_module, only : wetsalt_d use res_cs_module, only : wetcs_d use salt_module, only : hsaltb_d use cs_module, only : hcsb_d implicit none integer, intent (in) :: ob_id_dmd ! |object that has a water demand (HRU ID) real, intent (in) :: demand_vol !m3 |volume of water demand real, intent (inout) :: extracted !m3 |volume of groundwater extracted from aquifer for irrigation real, intent (inout) :: dmd_unmet !m3 |volume of irrigation demand not met by aquifer integer :: ob_id_src = 0 ! |source object (0 = HRU self-supply) character*4 :: demand_type = "hru" ! |demand type (default: HRU irrigation) integer :: i = 0 ! |counter integer :: s = 0 ! |solute counter integer :: cell_id = 0 ! |gwflow cell integer :: wetland = 0 ! |wetland flag integer :: isalt = 0 ! |salt ion counter integer :: ics = 0 ! |constituent counter integer :: sol_index = 0 integer :: hru_id = 0 real :: cell_demand = 0. !m3 |volume of irrigation water demand per gwflow cell real :: gwvol_avail = 0. !m3 |groundwater storage in the gwflow cell real :: gwvol_removed = 0. !m3 |groundwater removed from the gwflow cell for irrigation demand real :: gwvol_unmet = 0. !m3 |groundwater not available to meet irrigation demand real :: gw_mass = 0. !kg |mass of solute in groundwater of the gwflow cell real :: irr_mass(100) = 0. !kg |mass of solute removed from aquifer for irrigation real :: mass_diff = 0. !kg |difference between irrigation mass and actual groundwater mass real :: sum_pump = 0. !m3 |total pumping for the HRU real :: hru_area_m2 = 0. !m2 |HRU area real :: heat_flux = 0. !J |heat flux in pumped groundwater real :: soil_volm = 0. !m3 |volume of soil water in soil layer real :: soil_heat = 0. !J |heat in soil water of soil layer !zero out the total met and unmet demand extracted = 0. dmd_unmet = 0. !first option: hru irrigation from all cells geographically connected to the hru !hru demand --> irrigation ------------------------------------------------------------------------------------------------ if(demand_type == "hru" .and. ob_id_src == 0) then !only proceed if the HRU is connected (intersected with) gwflow grid cells sum_pump = 0. hru_id = ob_id_dmd if(hru_num_cells(hru_id) > 0) then !groundwater volume to remove from each cell connected to the HRU cell_demand = demand_vol / hru_num_cells(hru_id) !groundwater to remove from each cell connected to the HRU !loop through the cells that are connected to the HRU do i=1,hru_num_cells(hru_id) !id of the current cell cell_id = hru_cells(hru_id,i) !check for available groundwater in the cell if(gw_state(cell_id)%head > gw_state(cell_id)%botm) then !if water table is above bedrock gwvol_avail = ((gw_state(cell_id)%head-gw_state(cell_id)%botm) * gw_state(cell_id)%area) * gw_state(cell_id)%spyd !m3 else gwvol_avail = 0. endif !remove irrigation demand from storage; if storage is less than demand, remove all groundwater gwvol_removed = 0. gwvol_unmet = 0. if(gwvol_avail < cell_demand) then gwvol_removed = gwvol_avail gwvol_unmet = cell_demand - gwvol_avail !amount that is not available for irrigation else gwvol_removed = cell_demand endif extracted = extracted + gwvol_removed dmd_unmet = dmd_unmet + gwvol_unmet !save the pumping volume (m3), for use in gwflow_simulate gw_hyd_ss(cell_id)%ppag = gwvol_removed * (-1) !m3 negative = leaving the aquifer gw_hyd_ss_yr(cell_id)%ppag = gw_hyd_ss_yr(cell_id)%ppag + (gwvol_removed * (-1)) !store for annual water gw_hyd_ss_mo(cell_id)%ppag = gw_hyd_ss_mo(cell_id)%ppag + (gwvol_removed * (-1)) !store for monthly water !sum the pumping for the current HRU sum_pump = sum_pump + gwvol_removed !save the unsatisfied pumping volume (m3), for output gw_hyd_ss(cell_id)%ppdf = gwvol_unmet gw_hyd_ss_yr(cell_id)%ppdf = gw_hyd_ss_yr(cell_id)%ppdf + gwvol_unmet !add heat in irrigation water to soil profile if(gw_heat_flag == 1) then !remove heat from groundwater heat_flux = gwheat_state(cell_id)%temp * gw_rho * gw_cp * gwvol_removed !J if(heat_flux >= gwheat_state(cell_id)%stor) then heat_flux = gwheat_state(cell_id)%stor endif gw_heat_ss(cell_id)%ppag = heat_flux * (-1) gw_heat_ss_yr(cell_id)%ppag = gw_heat_ss_yr(cell_id)%ppag + heat_flux * (-1) endif !add solute mass in irrigation water to HRU soil profile !(mass is removed from the aquifer via mass balance equation in gwflow_simulate.f) if (gw_solute_flag == 1) then !determine the mass (kg) removed for irrigation water do s=1,gw_nsolute !loop through the solutes gw_mass = gwsol_state(cell_id)%solute(s)%mass / 1000. !kg in cell irr_mass(s) = (gwsol_state(cell_id)%solute(s)%conc * gwvol_removed) / 1000. !kg removed in irrigation water mass_diff = 0. if(irr_mass(s).gt.gw_mass) then mass_diff = irr_mass(s) - gw_mass endif irr_mass(s) = irr_mass(s) - mass_diff if(irr_mass(s).lt.0) irr_mass(s) = 0. enddo !add solute mass to soil profile of demand object (hru) wetland = hru(hru_id)%dbs%surf_stor !check if HRU is a wetland if(wetland > 0) then !add to wetland !nutrients wet(hru_id)%no3 = wet(hru_id)%no3 + irr_mass(1) !kg wet(hru_id)%solp = wet(hru_id)%solp + irr_mass(2) !kg sol_index = 2 !salts if (gwsol_salt == 1) then do isalt=1,cs_db%num_salts sol_index = sol_index + 1 wet_water(hru_id)%salt(isalt) = wet_water(hru_id)%salt(isalt) + irr_mass(sol_index) !kg wetsalt_d(hru_id)%salt(isalt)%irrig = wetsalt_d(hru_id)%salt(isalt)%irrig + irr_mass(sol_index) !kg enddo endif !constituents if (gwsol_cons == 1) then do ics=1,cs_db%num_cs sol_index = sol_index + 1 wet_water(hru_id)%cs(ics) = wet_water(hru_id)%cs(ics) + irr_mass(sol_index) !kg wetcs_d(hru_id)%cs(ics)%irrig = wetcs_d(hru_id)%cs(ics)%irrig + irr_mass(sol_index) !kg enddo endif else !add to soil profile !nutrients soil1(hru_id)%mn(1)%no3 = soil1(hru_id)%mn(1)%no3 + (irr_mass(1)/hru(hru_id)%area_ha) !kg/ha - add to soil layer soil1(hru_id)%mp(1)%lab = soil1(hru_id)%mp(1)%lab + (irr_mass(2)/hru(hru_id)%area_ha) !kg/ha - add to soil layer sol_index = 2 !salts if (gwsol_salt == 1) then do isalt=1,cs_db%num_salts sol_index = sol_index + 1 cs_soil(hru_id)%ly(1)%salt(isalt) = cs_soil(hru_id)%ly(1)%salt(isalt) + (irr_mass(sol_index)/hru(hru_id)%area_ha) !kg/ha - add to soil layer hsaltb_d(hru_id)%salt(isalt)%irgw = hsaltb_d(hru_id)%salt(isalt)%irgw + (irr_mass(sol_index)/hru(hru_id)%area_ha) !kg/ha - include in soil salt balance enddo endif !constituents if (gwsol_cons == 1) then do ics=1,cs_db%num_cs sol_index = sol_index + 1 cs_soil(hru_id)%ly(1)%cs(ics) = cs_soil(hru_id)%ly(1)%cs(ics) + (irr_mass(sol_index)/hru(hru_id)%area_ha) !kg/ha - add to soil layer hcsb_d(hru_id)%cs(ics)%irgw = hcsb_d(hru_id)%cs(ics)%irgw + (irr_mass(sol_index)/hru(hru_id)%area_ha) !kg/ha - include in soil constituent balance enddo endif endif !add to mass balance arrays (to be used in gwflow_simulate) do s=1,gw_nsolute !loop through the solutes gwsol_ss(cell_id)%solute(s)%ppag = (irr_mass(s)*1000.) * (-1) !g; negative = leaving the aquifer gwsol_ss_sum(cell_id)%solute(s)%ppag = gwsol_ss_sum(cell_id)%solute(s)%ppag - (irr_mass(s)*1000.) gwsol_ss_sum_mo(cell_id)%solute(s)%ppag = gwsol_ss_sum_mo(cell_id)%solute(s)%ppag - (irr_mass(s)*1000.) enddo endif enddo !go to next gwflow cell endif hru_pump(hru_id) = sum_pump !store pumping (m3) for the HRU !irrigation pumping from a single cell ------------------------------------------------------------------------------------ elseif(demand_type == "hru" .and. ob_id_src > 0) then hru_id = ob_id_dmd sum_pump = 0. cell_demand = demand_vol !m3 groundwater volume to remove from the cell if(grid_type == "structured") then cell_id = cell_id_list(ob_id_src) endif !check for available groundwater in the cell if(gw_state(cell_id)%head > gw_state(cell_id)%botm) then !if water table is above bedrock gwvol_avail = ((gw_state(cell_id)%head-gw_state(cell_id)%botm) * gw_state(cell_id)%area) * gw_state(cell_id)%spyd !m3 else gwvol_avail = 0. endif !remove irrigation demand from storage; if storage is less than demand, remove all groundwater gwvol_removed = 0. gwvol_unmet = 0. if(gwvol_avail < cell_demand) then gwvol_removed = gwvol_avail gwvol_unmet = cell_demand - gwvol_avail !amount that is not available for irrigation else gwvol_removed = cell_demand endif extracted = extracted + gwvol_removed dmd_unmet = dmd_unmet + gwvol_unmet !save the pumping volume (m3), for use in gwflow_simulate gw_hyd_ss(cell_id)%ppag = gwvol_removed * (-1) !m3 negative = leaving the aquifer gw_hyd_ss_yr(cell_id)%ppag = gw_hyd_ss_yr(cell_id)%ppag + (gwvol_removed * (-1)) !store for annual water gw_hyd_ss_mo(cell_id)%ppag = gw_hyd_ss_mo(cell_id)%ppag + (gwvol_removed * (-1)) !store for monthly water !sum the pumping for the current HRU sum_pump = sum_pump + gwvol_removed !save the unsatisfied pumping volume (m3), for output gw_hyd_ss(cell_id)%ppdf = gwvol_unmet gw_hyd_ss_yr(cell_id)%ppdf = gw_hyd_ss_yr(cell_id)%ppdf + gwvol_unmet !add heat in irrigation water to soil profile if(gw_heat_flag == 1) then !remove heat from groundwater heat_flux = gwheat_state(cell_id)%temp * gw_rho * gw_cp * gwvol_removed !J if(heat_flux >= gwheat_state(cell_id)%stor) then heat_flux = gwheat_state(cell_id)%stor endif gw_heat_ss(cell_id)%ppag = heat_flux * (-1) gw_heat_ss_yr(cell_id)%ppag = gw_heat_ss_yr(cell_id)%ppag + heat_flux * (-1) endif !add solute mass in irrigation water to HRU soil profile !(mass is removed from the aquifer via mass balance equation in gwflow_simulate.f) if (gw_solute_flag == 1) then !determine the mass (kg) removed for irrigation water do s=1,gw_nsolute !loop through the solutes gw_mass = gwsol_state(cell_id)%solute(s)%mass / 1000. !kg in cell irr_mass(s) = (gwsol_state(cell_id)%solute(s)%conc * gwvol_removed) / 1000. !kg removed in irrigation water mass_diff = 0. if(irr_mass(s).gt.gw_mass) then mass_diff = irr_mass(s) - gw_mass endif irr_mass(s) = irr_mass(s) - mass_diff if(irr_mass(s).lt.0) irr_mass(s) = 0. enddo !add solute mass to soil profile of demand object (hru) wetland = hru(hru_id)%dbs%surf_stor !check if HRU is a wetland if(wetland > 0) then !add to wetland !nutrients wet(hru_id)%no3 = wet(hru_id)%no3 + irr_mass(1) !kg wet(hru_id)%solp = wet(hru_id)%solp + irr_mass(2) !kg sol_index = 2 !salts if (gwsol_salt == 1) then do isalt=1,cs_db%num_salts sol_index = sol_index + 1 wet_water(hru_id)%salt(isalt) = wet_water(hru_id)%salt(isalt) + irr_mass(sol_index) !kg wetsalt_d(hru_id)%salt(isalt)%irrig = wetsalt_d(hru_id)%salt(isalt)%irrig + irr_mass(sol_index) !kg enddo endif !constituents if (gwsol_cons == 1) then do ics=1,cs_db%num_cs sol_index = sol_index + 1 wet_water(hru_id)%cs(ics) = wet_water(hru_id)%cs(ics) + irr_mass(sol_index) !kg wetcs_d(hru_id)%cs(ics)%irrig = wetcs_d(hru_id)%cs(ics)%irrig + irr_mass(sol_index) !kg enddo endif else !add to soil profile !nutrients soil1(hru_id)%mn(1)%no3 = soil1(hru_id)%mn(1)%no3 + (irr_mass(1)/hru(hru_id)%area_ha) !kg/ha - add to soil layer soil1(hru_id)%mp(1)%lab = soil1(hru_id)%mp(1)%lab + (irr_mass(2)/hru(hru_id)%area_ha) !kg/ha - add to soil layer sol_index = 2 !salts if (gwsol_salt == 1) then do isalt=1,cs_db%num_salts sol_index = sol_index + 1 cs_soil(hru_id)%ly(1)%salt(isalt) = cs_soil(hru_id)%ly(1)%salt(isalt) + (irr_mass(sol_index)/hru(hru_id)%area_ha) !kg/ha - add to soil layer hsaltb_d(hru_id)%salt(isalt)%irgw = hsaltb_d(hru_id)%salt(isalt)%irgw + (irr_mass(sol_index)/hru(hru_id)%area_ha) !kg/ha - include in soil salt balance enddo endif !constituents if (gwsol_cons == 1) then do ics=1,cs_db%num_cs sol_index = sol_index + 1 cs_soil(hru_id)%ly(1)%cs(ics) = cs_soil(hru_id)%ly(1)%cs(ics) + (irr_mass(sol_index)/hru(hru_id)%area_ha) !kg/ha - add to soil layer hcsb_d(hru_id)%cs(ics)%irgw = hcsb_d(hru_id)%cs(ics)%irgw + (irr_mass(sol_index)/hru(hru_id)%area_ha) !kg/ha - include in soil constituent balance enddo endif endif !add to mass balance arrays (to be used in gwflow_simulate) do s=1,gw_nsolute !loop through the solutes gwsol_ss(cell_id)%solute(s)%ppag = (irr_mass(s)*1000.) * (-1) !g; negative = leaving the aquifer gwsol_ss_sum(cell_id)%solute(s)%ppag = gwsol_ss_sum(cell_id)%solute(s)%ppag - (irr_mass(s)*1000.) gwsol_ss_sum_mo(cell_id)%solute(s)%ppag = gwsol_ss_sum_mo(cell_id)%solute(s)%ppag - (irr_mass(s)*1000.) enddo endif hru_pump(hru_id) = sum_pump !store pumping (m3) for the HRU !municipal demand from all cells connected geographically to the hru ------------------------------------------------------ elseif(demand_type == "muni" .and. ob_id_src == 0) then hru_id = ob_id_dmd if(hru_num_cells(hru_id) > 0) then !groundwater volume to remove from each cell connected to the HRU cell_demand = demand_vol / hru_num_cells(hru_id) !groundwater to remove from each cell connected to the HRU !loop through the cells that are connected to the HRU do i=1,hru_num_cells(hru_id) !id of the current cell cell_id = hru_cells(hru_id,i) !check for available groundwater in the cell if(gw_state(cell_id)%head > gw_state(cell_id)%botm) then !if water table is above bedrock gwvol_avail = ((gw_state(cell_id)%head-gw_state(cell_id)%botm) * gw_state(cell_id)%area) * gw_state(cell_id)%spyd !m3 else gwvol_avail = 0. endif !remove irrigation demand from storage; if storage is less than demand, remove all groundwater gwvol_removed = 0. gwvol_unmet = 0. if(gwvol_avail < cell_demand) then gwvol_removed = gwvol_avail gwvol_unmet = cell_demand - gwvol_avail !amount that is not available for irrigation else gwvol_removed = cell_demand endif extracted = extracted + gwvol_removed dmd_unmet = dmd_unmet + gwvol_unmet !save the pumping volume (m3), for use in gwflow_simulate gw_hyd_ss(cell_id)%ppex = gwvol_removed * (-1) !m3 negative = leaving the aquifer gw_hyd_ss_yr(cell_id)%ppex = gw_hyd_ss_yr(cell_id)%ppex + (gwvol_removed * (-1)) !store for annual water gw_hyd_ss_mo(cell_id)%ppex = gw_hyd_ss_mo(cell_id)%ppex + (gwvol_removed * (-1)) !store for monthly water !save the unsatisfied pumping volume (m3), for output gw_hyd_ss(cell_id)%ppdf = gwvol_unmet gw_hyd_ss_yr(cell_id)%ppdf = gw_hyd_ss_yr(cell_id)%ppdf + gwvol_unmet !remove heat from aquifer if(gw_heat_flag == 1) then heat_flux = gwheat_state(cell_id)%temp * gw_rho * gw_cp * gwvol_removed !J if(heat_flux >= gwheat_state(cell_id)%stor) then heat_flux = gwheat_state(cell_id)%stor endif gw_heat_ss(cell_id)%ppex = heat_flux * (-1) gw_heat_ss_yr(cell_id)%ppex = gw_heat_ss_yr(cell_id)%ppex + heat_flux * (-1) endif !determine solute mass to be from the aquifer via mass balance equation (in gwflow_simulate.f) if (gw_solute_flag == 1) then !determine the mass (kg) removed for irrigation water do s=1,gw_nsolute !loop through the solutes gw_mass = gwsol_state(cell_id)%solute(s)%mass / 1000. !kg in cell irr_mass(s) = (gwsol_state(cell_id)%solute(s)%conc * gwvol_removed) / 1000. !kg removed in irrigation water mass_diff = 0. if(irr_mass(s).gt.gw_mass) then mass_diff = irr_mass(s) - gw_mass endif irr_mass(s) = irr_mass(s) - mass_diff if(irr_mass(s).lt.0) irr_mass(s) = 0. enddo !add to mass balance arrays (to be used in gwflow_simulate) do s=1,gw_nsolute !loop through the solutes gwsol_ss(cell_id)%solute(s)%ppex = (irr_mass(s)*1000.) * (-1) !g; negative = leaving the aquifer gwsol_ss_sum(cell_id)%solute(s)%ppex = gwsol_ss_sum(cell_id)%solute(s)%ppex - (irr_mass(s)*1000.) gwsol_ss_sum_mo(cell_id)%solute(s)%ppex = gwsol_ss_sum_mo(cell_id)%solute(s)%ppex - (irr_mass(s)*1000.) enddo endif enddo !go to next gwflow cell endif !municipal demand from a single cell -------------------------------------------------------------------------------------- elseif(demand_type == "muni" .and. ob_id_src > 0) then cell_demand = demand_vol !m3 groundwater volume to remove from the cell if(grid_type == "structured") then cell_id = cell_id_list(ob_id_src) endif !check for available groundwater in the cell if(gw_state(cell_id)%head > gw_state(cell_id)%botm) then !if water table is above bedrock gwvol_avail = ((gw_state(cell_id)%head-gw_state(cell_id)%botm) * gw_state(cell_id)%area) * gw_state(cell_id)%spyd !m3 else gwvol_avail = 0. endif !remove irrigation demand from storage; if storage is less than demand, remove all groundwater gwvol_removed = 0. gwvol_unmet = 0. if(gwvol_avail < cell_demand) then gwvol_removed = gwvol_avail gwvol_unmet = cell_demand - gwvol_avail !amount that is not available for irrigation else gwvol_removed = cell_demand endif extracted = extracted + gwvol_removed dmd_unmet = dmd_unmet + gwvol_unmet !save the pumping volume (m3), for use in gwflow_simulate gw_hyd_ss(cell_id)%ppex = gwvol_removed * (-1) !m3 negative = leaving the aquifer gw_hyd_ss_yr(cell_id)%ppex = gw_hyd_ss_yr(cell_id)%ppex + (gwvol_removed * (-1)) !store for annual water gw_hyd_ss_mo(cell_id)%ppex = gw_hyd_ss_mo(cell_id)%ppex + (gwvol_removed * (-1)) !store for monthly water !save the unsatisfied pumping volume (m3), for output gw_hyd_ss(cell_id)%ppdf = gwvol_unmet gw_hyd_ss_yr(cell_id)%ppdf = gw_hyd_ss_yr(cell_id)%ppdf + gwvol_unmet !remove heat from aquifer if(gw_heat_flag == 1) then !remove heat from groundwater heat_flux = gwheat_state(cell_id)%temp * gw_rho * gw_cp * gwvol_removed !J if(heat_flux >= gwheat_state(cell_id)%stor) then heat_flux = gwheat_state(cell_id)%stor endif gw_heat_ss(cell_id)%ppex = heat_flux * (-1) gw_heat_ss_yr(cell_id)%ppex = gw_heat_ss_yr(cell_id)%ppex + heat_flux * (-1) endif !determine solute mass to be from the aquifer via mass balance equation (in gwflow_simulate.f) if(gw_solute_flag == 1) then !determine the mass (kg) removed by pumping do s=1,gw_nsolute !loop through the solutes gw_mass = gwsol_state(cell_id)%solute(s)%mass / 1000. !kg in cell irr_mass(s) = (gwsol_state(cell_id)%solute(s)%conc * gwvol_removed) / 1000. !kg removed in irrigation water mass_diff = 0. if(irr_mass(s).gt.gw_mass) then mass_diff = irr_mass(s) - gw_mass endif irr_mass(s) = irr_mass(s) - mass_diff if(irr_mass(s).lt.0) irr_mass(s) = 0. enddo !add to mass balance arrays do s=1,gw_nsolute !loop through the solutes gwsol_ss(cell_id)%solute(s)%ppex = (irr_mass(s)*1000.) * (-1) !g; negative = leaving the aquifer gwsol_ss_sum(cell_id)%solute(s)%ppex = gwsol_ss_sum(cell_id)%solute(s)%ppex - (irr_mass(s)*1000.) gwsol_ss_sum_mo(cell_id)%solute(s)%ppex = gwsol_ss_sum_mo(cell_id)%solute(s)%ppex - (irr_mass(s)*1000.) enddo endif endif !check: hru or muni end subroutine gwflow_pump_allo