      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