      subroutine gwflow_soil(hru_id) !rtb gwflow

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates the water exchange volume between the aquifer and the soil profile
!!    (exchange volumes are used in gwflow_simulate, in groundwater balance equations)
      
      use gwflow_module
      use soil_module, only : soil
      use hydrograph_module, only : ob
      use hru_module, only : gwsoilq
      
      implicit none

      integer, intent (in) :: hru_id		 !       |hru number
      integer :: k                       !       |counter
      integer :: s                       !       |solute counter
      integer :: jj                      !       |soil layer counter
      integer :: cell_id                 !       |cell in connection with the channel
      real :: Q                          !m3     |volume transferred from cell
      real :: hru_Q                      !m3     |total volume transferred to soil profile
      real :: hru_soilz                  !m      |thickness of HRU soil profile
      real :: vadose_z    							 !m			 |thickness of cell vadose zone
      real :: poly_area                  !m2     |area of cell within the HRU
      real :: solmass(100)               !g      |solute mass transferred from cell
      real :: hru_mass(100)              !g      |total solute mass transferred to soil profile
      real :: water_depth(100)           !m      |depth of groundwater in each soil layer
      real :: water_depth_tot            !m      |total depth of groundwater in the soil profile
      real :: sol_thick                  !m      |thickness of soil layer
      real :: layer_fraction             !       |fraction of saturated soil profile within the layer
      real :: layer_transfer             !       |amount of water and solute transferred to the soil layer
      real :: hru_area_m2                !m2     |surface area of the hru
      

      
      !only proceed if gw-->soil exchange is active
      if (gw_soil_flag == 1) then
      
        !HRU soil thickness
        hru_soilz = soil(hru_id)%phys(soil(hru_id)%nly)%d / 1000. !m
        
        !loop through the grid cells connected to the HRU -------------------------------------------------------------
        hru_Q = 0.
        hru_mass = 0.
        do k=1,hru_num_cells(hru_id)
          
          !cell in connection with the HRU
          cell_id = hru_cells(hru_id,k)
        
          !only proceed if the cell is active
          if(gw_state(cell_id)%stat == 1) then 
          
            !calculate the current thickness of the vadose zone (m)
            vadose_z = gw_state(cell_id)%elev - gw_state(cell_id)%head !thickness of vadose zone (m)
            
            !if water table is within the soil profile --> calculate groundwater volume (Q) to transfer
            Q = 0.
            if(vadose_z < hru_soilz) then !water table is within the soil profile
					    poly_area = gw_state(cell_id)%area * cells_fract(hru_id,k) !area of cell within HRU
              Q = (hru_soilz - vadose_z) * poly_area * gw_state(cell_id)%spyd !m3 of groundwater to transfer to the soil profile 
              !store for water balance calculations (in gwflow_simulate)
              gw_ss(cell_id)%soil = gw_ss(cell_id)%soil + (Q*(-1)) !negative = leaving the aquifer
              gw_ss_sum(cell_id)%soil = gw_ss_sum(cell_id)%soil + (Q*(-1))
              !sum for HRU
              hru_Q = hru_Q + Q
              !solutes
              if (gw_solute_flag == 1) then
                !solute mass transferred from aquifer to soil
                do s=1,gw_nsolute
                  solmass(s) = Q * gwsol_state(cell_id)%solute(s)%conc !g
                  if(solmass(s) > gwsol_state(cell_id)%solute(s)%mass) then !can only remove what is there
                    solmass(s) = gwsol_state(cell_id)%solute(s)%mass
                  endif
                  gwsol_ss(cell_id)%solute(s)%soil = solmass(s) * (-1) !negative = leaving the aquifer
                  gwsol_ss_sum(cell_id)%solute(s)%soil = gwsol_ss_sum(cell_id)%solute(s)%soil + (solmass(s)*(-1))
                  hru_mass(s) = hru_mass(s) + solmass(s) !g
                enddo  
              endif !end solutes
            endif !check for water table in soil profile
                      
          endif !check if cell is active
          
        enddo !go to next connected cell
        
        !add water and solute mass to HRU soil profile ----------------------------------------------------------------
        if(hru_Q > 0) then
          hru_area_m2 = ob(hru_id)%area_ha * 10000. !area of the HRU in m2     
          !determine which soil layers are to receive groundwater
          water_depth = 0.
          water_depth_tot = 0.
          do jj=1,soil(hru_id)%nly
            sol_thick = soil(hru_id)%phys(jj)%thick / 1000. !mm --> m
            if((soil(hru_id)%phys(jj)%d/1000.) > vadose_z) then
              if(jj == 1) then !top layer
                water_depth(jj) = (soil(hru_id)%phys(jj)%d/1000.) - vadose_z
              else
                if(vadose_z > (soil(hru_id)%phys(jj-1)%d/1000.)) then !water table is within this layer
                  water_depth(jj) = (soil(hru_id)%phys(jj)%d/1000.) - vadose_z
                else !water table is above this layer
                  water_depth(jj) = sol_thick
                endif
              endif
              water_depth_tot = water_depth_tot + water_depth(jj)
            endif  
          enddo !next soil layer
          !transfer the groundwater and the solute mass
          do jj=1,soil(hru_id)%nly
            if(water_depth_tot > 0) then
              layer_fraction = water_depth(jj) / water_depth_tot
            else
              layer_fraction = 0.
            endif
            layer_transfer = (hru_Q*layer_fraction) / hru_area_m2 * 1000. !m3 --> mm
            soil(hru_id)%phys(jj)%st = soil(hru_id)%phys(jj)%st + layer_transfer !mm
            gwsoilq(hru_id) = gwsoilq(hru_id) + layer_transfer !store for hru output
            if (gw_solute_flag == 1) then
              do s=1,gw_nsolute !loop through the solutes
                layer_transfer = (hru_mass(s)*layer_fraction) / 1000. / ob(hru_id)%area_ha !g --> kg/ha
                hru_soil(hru_id,jj,s) = hru_soil(hru_id,jj,s) + layer_transfer !kg/ha (mass added to soil profile in nut_nlch, nut_solp, salt_lch, cs_lch)
              enddo
            endif
          enddo  
        endif
        
      endif !check if gw-->soil transfer is active
      
      return
      end subroutine gwflow_soil     