      subroutine gwflow_channel_exch(chan_id) !rtb gwflow

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates the water exchange volume between the channel and the connected grid cells
!!    (exchange volumes are used in gwflow_simulate, in groundwater balance equations)

      use gwflow_module
      use hydrograph_module, only : ch_stor,ch_out_d
      use sd_channel_module, only : sd_ch
      use constituent_mass_module
      use time_module

      implicit none

      integer, intent (in) :: chan_id    !       |channel number
      integer :: k = 0                   !       |counter
      integer :: s = 0                   !       |solute counter
      integer :: cell_id = 0             !       |cell in connection with the channel
      integer :: isalt = 0               !       |salt ion counter
      integer :: ics = 0                 !       |constituent counter
      integer :: sol_index = 0           !
      integer :: dum = 0
      real :: chan_depth = 0.            !m      |channel depth
      real :: chan_width = 0.            !m      |channel width
      real :: chan_length = 0.           !m      |length of channel, in cell
      real :: bed_elev = 0.              !m      |elevation of channel bed, in cell
      real :: bed_K = 0.                 !m/day  |hydraulic conductivity of channel bed, in cell
      real :: bed_thick = 0.             !m      |thickness of channel bed, in cell
      real :: chan_stage = 0.            !m      |elevation of water in channel, in cell
      real :: flow_area = 0.             !m2     |groundwater flow area of water exchange, in cell
      real :: gw_head = 0.               !m      |current groundwater head in cell
      real :: Q = 0.                     !m3/day |water exchange flow rate, calculated by Darcy's Law
      real :: head_diff = 0.             !m      |head difference between channel stage and groundwater head
      real :: chan_volume = 0.           !m3     |water volume in channel before groundwater exchange occurs
      real :: stor_change = 0.
      real :: sat_change = 0.
      real :: chan_csol(100) = 0.        !g/m3   |solute concentration in channel water
      real :: solmass(100) = 0.          !g      |solute mass transferred
      real :: chan_heat = 0.             !J      |current heat in channel
      real :: heat_flux = 0.             !J/day  |heat transferred between groundwater and channel
      real :: chan_flow = 0.
			real :: chan_temp = 0.


      !current channel storage (m3) and temperature (deg C)
      chan_volume = ch_stor(chan_id)%flo

      !characteristics of channel
      chan_depth = sd_ch(chan_id)%chd !depth (m) of water in channel
      chan_width = sd_ch(chan_id)%chw !width (m) of channel

      !loop through the cells connected to the channel
      do k=1,gw_chan_info(chan_id)%ncon

        !cell in connection with the channel
        cell_id = gw_chan_info(chan_id)%cells(k)

        !only proceed if the cell is active
        if(gw_state(cell_id)%stat == 1) then

          !characteristics of cell
          chan_length = gw_chan_info(chan_id)%leng(k)
          bed_elev = gw_chan_info(chan_id)%elev(k) - gw_bed_change
          bed_K = gw_chan_info(chan_id)%hydc(k)
          bed_thick = gw_chan_info(chan_id)%thck(k)

          !derived values
					if(gw_chan_dep_flag == 1) then !depth zone for daily channel depth: specified in gwflow.chancells_depth
					  chan_depth = gw_chan_dep(gw_chan_info(chan_id)%dpzn(k))
					endif
					chan_stage = bed_elev + chan_depth !stage of water in channel (m)
          flow_area = chan_width * chan_length !water exchange flow area (m2)

          !calculate flow exchange rate (m3/day) using Darcy's Law
          !head difference --> head gradient --> flow rate
          gw_head = gw_state(cell_id)%head
          Q = 0.
          if(gw_head < bed_elev) then
            head_diff = chan_depth
            Q =  bed_K * (head_diff / bed_thick) * flow_area !stream leakage (positive Q: entering aquifer)
          elseif (gw_head > chan_stage) then
            head_diff = gw_head - chan_stage
            Q = bed_K * (head_diff / bed_thick) * flow_area * (-1) !gw discharge (negative Q: leaving aquifer)
          elseif (gw_head > bed_elev .and. gw_head < chan_stage) then
            head_diff = chan_stage - gw_head
            Q = bed_K * (head_diff / bed_thick) * flow_area !stream leakage (positive Q: entering aquifer)
          endif

          !store values in gwflow source/sink arrays
          if(Q < 0) then !aquifer --> channel
            if (-Q >= gw_state(cell_id)%stor) then !can only remove what is there
              Q = -gw_state(cell_id)%stor
            endif
						gw_hyd_ss(cell_id)%gwsw = gw_hyd_ss(cell_id)%gwsw + Q
          else !channel --> aquifer
            if(Q > ch_stor(chan_id)%flo) then !can only remove what is there
              Q = ch_stor(chan_id)%flo
            endif
            gw_hyd_ss(cell_id)%swgw = gw_hyd_ss(cell_id)%swgw + Q
          endif
		  !remove groundwater from storage; calculate change in saturated thickness and new gw head
          stor_change = Q !m3
		  gw_state(cell_id)%stor = gw_state(cell_id)%stor + stor_change !m3
          sat_change = stor_change / (gw_state(cell_id)%spyd * gw_state(cell_id)%area) !m
          gw_state(cell_id)%head = gw_state(cell_id)%head + sat_change !m
		  !annual and monthly values
          gw_hyd_ss_yr(cell_id)%gwsw = gw_hyd_ss_yr(cell_id)%gwsw + Q
          gw_hyd_ss_mo(cell_id)%gwsw = gw_hyd_ss_mo(cell_id)%gwsw + Q
          !store for channel object (positive value = water added to channel)
					chan_flow = ch_stor(chan_id)%flo
          ch_stor(chan_id)%flo = ch_stor(chan_id)%flo + (Q*(-1))

          !heat
          if(gw_heat_flag == 1) then
            !current heat in channel
						chan_temp = ch_stor(chan_id)%temp
            chan_heat = ch_stor(chan_id)%temp * gw_rho * gw_cp * chan_flow !J in channel
            if(Q < 0) then !leaving the cell (aquifer --> channel)
              heat_flux = gwheat_state(cell_id)%temp * gw_rho * gw_cp * Q !J
              if(-heat_flux >= gwheat_state(cell_id)%stor) then
                heat_flux = -gwheat_state(cell_id)%stor
              endif
              gw_heat_ss(cell_id)%gwsw = gw_heat_ss(cell_id)%gwsw + heat_flux
              gwheat_state(cell_id)%stor = gwheat_state(cell_id)%stor + heat_flux !update heat in the cell
            else !entering cell (channel --> aquifer)
              heat_flux = 0.
              if(ch_stor(chan_id)%temp > 0) then
                heat_flux = ch_stor(chan_id)%temp * gw_rho * gw_cp * Q !J
                if(heat_flux > chan_heat) then
                  heat_flux = chan_heat
                endif
                gw_heat_ss(cell_id)%swgw = gw_heat_ss(cell_id)%swgw + heat_flux
              endif
            endif
            !update channel heat and temperature
            chan_heat = chan_heat + (heat_flux*(-1))
            if(ch_stor(chan_id)%flo > 0) then
							chan_temp = chan_heat / (gw_rho * gw_cp * ch_stor(chan_id)%flo)
              ch_stor(chan_id)%temp = chan_heat / (gw_rho * gw_cp * ch_stor(chan_id)%flo)
            else
              ch_stor(chan_id)%temp = 0.
						endif
						ch_out_d(chan_id)%temp = ch_stor(chan_id)%temp
            !store for annual values
            gw_heat_ss_yr(cell_id)%gwsw = gw_heat_ss_yr(cell_id)%gwsw + heat_flux !J
          endif

          !calculate solute mass (g/day) transported between cell and channel
          if (gw_solute_flag == 1) then
            chan_csol = 0.
            solmass = 0.
            if(Q < 0) then !mass leaving the cell (aquifer --> channel)
              do s=1,gw_nsolute !loop through the solutes
                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)%gwsw = solmass(s)
                gwsol_ss_sum(cell_id)%solute(s)%gwsw = gwsol_ss_sum(cell_id)%solute(s)%gwsw + solmass(s)
								gwsol_ss_sum_mo(cell_id)%solute(s)%gwsw = gwsol_ss_sum_mo(cell_id)%solute(s)%gwsw + solmass(s)
              enddo
              !add solute mass to channel
              ch_stor(chan_id)%no3 = ch_stor(chan_id)%no3 + (solmass(1)*(-1)/1000.) !kg
              ch_stor(chan_id)%solp = ch_stor(chan_id)%solp + (solmass(2)*(-1)/1000.) !kg
              sol_index = 2
              !salts
              if (gwsol_salt == 1) then
                do isalt=1,cs_db%num_salts
                  sol_index = sol_index + 1
                  ch_water(chan_id)%salt(isalt) = ch_water(chan_id)%salt(isalt) + (solmass(sol_index)*(-1)/1000.) !kg
                enddo
              endif
              !constituents
              if (gwsol_cons == 1) then
                do ics=1,cs_db%num_cs
                  sol_index = sol_index + 1
                  ch_water(chan_id)%cs(ics) = ch_water(chan_id)%cs(ics) + (solmass(sol_index)*(-1)/1000.) !kg
                enddo
              endif
            else !mass entering cell (channel --> aquifer)
              !calculate solute concentrations in channel water
              if(chan_volume > 10.) then
                chan_csol(1) = (ch_stor(chan_id)%no3 * 1000.) / chan_volume !no3 g/m3 in channel
                chan_csol(2) = (ch_stor(chan_id)%solp * 1000.) / chan_volume !p g/m3 in channel
                sol_index = 2
                !salts
                if (gwsol_salt == 1) then
                  do isalt=1,cs_db%num_salts
                    sol_index = sol_index + 1
                    chan_csol(sol_index) = (ch_water(chan_id)%salt(isalt)*1000.) / chan_volume !g/m3 in channel water
                  enddo
                endif
                !constituents
                if (gwsol_cons == 1) then
                  do ics=1,cs_db%num_cs
                    sol_index = sol_index + 1
                    chan_csol(sol_index) = (ch_water(chan_id)%cs(ics)*1000.) / chan_volume !g/m3 in channel water
                    if(chan_csol(sol_index) < 0) then
                      dum = 10
                    endif
                  enddo
                endif
              endif
              !calculate mass (g)
              do s=1,gw_nsolute
                solmass(s) = Q * chan_csol(s)
              enddo
              !check against available solute mass in the channel
              !no3
              if((solmass(1)/1000.) > ch_stor(chan_id)%no3) then
                solmass(1) = ch_stor(chan_id)%no3 * 1000. !g
              endif
              ch_stor(chan_id)%no3 = ch_stor(chan_id)%no3 - (solmass(1)/1000.) !kg
              !p
              if((solmass(2)/1000.) > ch_stor(chan_id)%solp) then
                solmass(2) = ch_stor(chan_id)%solp * 1000. !g
              endif
              ch_stor(chan_id)%solp = ch_stor(chan_id)%solp - (solmass(2)/1000.) !kg
              sol_index = 2
              !salts
              if (gwsol_salt == 1) then
                do isalt=1,cs_db%num_salts
                  sol_index = sol_index + 1
                  if((solmass(sol_index)/1000.) > ch_water(chan_id)%salt(isalt)) then
                    solmass(sol_index) = ch_water(chan_id)%salt(isalt) * 1000. !g
                  endif
                  ch_water(chan_id)%salt(isalt) = ch_water(chan_id)%salt(isalt) - (solmass(sol_index)/1000.) !kg
                enddo
              endif
              !constituents
              if (gwsol_cons == 1) then
                do ics=1,cs_db%num_cs
                  sol_index = sol_index + 1
                  if((solmass(sol_index)/1000.) > ch_water(chan_id)%cs(ics)) then
                    solmass(sol_index) = ch_water(chan_id)%cs(ics) * 1000. !g
                  endif
                  ch_water(chan_id)%cs(ics) = ch_water(chan_id)%cs(ics) - (solmass(sol_index)/1000.) !kg
                enddo
              endif
              !store for mass balance calculations (in gwflow_simulate)
              do s=1,gw_nsolute !loop through the solutes
                gwsol_ss(cell_id)%solute(s)%swgw = solmass(s)
                gwsol_ss_sum(cell_id)%solute(s)%swgw = gwsol_ss_sum(cell_id)%solute(s)%swgw + solmass(s)
								gwsol_ss_sum_mo(cell_id)%solute(s)%swgw = gwsol_ss_sum_mo(cell_id)%solute(s)%swgw + solmass(s)
              enddo
            endif

          endif !end solutes

        endif !check if cell is active

      enddo !go to next cell

      return
      end subroutine gwflow_channel_exch
