      subroutine gwflow_canal(chan_id) !rtb gwflow

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates the water exchange volume between irrigation canals and connected grid cells
!!    (exchange volumes are used in gwflow_simulate, in groundwater balance equations), for canals
!!    that remove water from a specified channel.

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

      implicit none

      integer, intent (in) :: chan_id		 !       |channel number
      integer :: c = 0											 !       |counter for canals connected to the channel
      integer :: k = 0                       !       |counter for cells connected to a canal
      integer :: s = 0                       !       |counter of groundwater solutes
      integer :: canal_id = 0                !       |canal in connection with the channel
      integer :: cell_id = 0                 !       |cell in connection with the canal
      integer :: day_beg = 0                 !       |beginning day (julian) of active canal
      integer :: day_end = 0                 !       |ending day (julian) of active canal
      integer :: isalt = 0                   !       |salt ion counter
      integer :: ics = 0                     !       |constituent counter
      integer :: sol_index = 0
      integer :: dum = 0
      real :: chan_volume = 0.                !m3     |water volume in channel before groundwater exchange occurs
      real :: width = 0.                      !m      |canal width
      real :: depth = 0.                      !m      |canal depth
      real :: thick = 0.                      !m      |canal bed thickness
      real :: length = 0.                     !m      |length of canal in the cell
      real :: stage = 0.                      !m      |stage of canal in the cell
      real :: bed_K = 0.											 !m/day  |hydraulic conductivity of canal bed in the cell
	  real :: reduc = 0.
	  real :: daycount_real = 0.
      real :: flow_area = 0.                  !m2     |groundwater flow area of water exchange, in cell
      real :: canal_bed = 0.                  !m      |canal bed elevation in the cell
      real :: head_diff = 0.                  !m      |head difference between canal stage and groundwater head
      real :: Q = 0.                          !m3/day |water exchange flow rate, calculated by Darcy's Law
      real :: chan_csol(100) = 0.             !g/m3   |solute concentration in channel water
      real :: solmass(100) = 0.               !g      |solute mass transferred
      real :: conc_nh3 = 0.
      real :: conc_no2 = 0.
      real :: conc_dox = 0.
      real :: conc_orgn = 0.
      real :: mass_nh3 = 0.
      real :: mass_no2 = 0.
      real :: mass_dox = 0.
      real :: mass_orgn = 0.
      real :: heat_flux = 0.                  !J      |heat flux in groundwater-canal exchange water
      real :: chan_heat = 0.                  !J      |heat in channel water



      !only proceed if canal-cell exchange is active
      if (gw_canal_flag == 1) then

        !record starting channel volume (m3) and solute concentrations (g/m3)
        chan_csol = 0.
        conc_nh3 = 0.
        conc_no2 = 0.
        conc_dox = 0.
        conc_orgn = 0.
        chan_volume = ch_stor(chan_id)%flo
        if (chan_volume > 10.) then
          chan_csol(1) = (ch_stor(chan_id)%no3 * 1000.) / chan_volume
          chan_csol(2) = (ch_stor(chan_id)%solp * 1000.) / chan_volume
          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
            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
            enddo
          endif
          conc_nh3 = (ch_stor(chan_id)%nh3 * 1000.) / chan_volume
          conc_no2 = (ch_stor(chan_id)%no2 * 1000.) / chan_volume
          conc_dox = (ch_stor(chan_id)%dox * 1000.) / chan_volume
          conc_orgn = (ch_stor(chan_id)%orgn * 1000.) / chan_volume
        endif

        !loop through the canals that are connected to the channel
        do c=1,gw_chan_canl_info(chan_id)%ncanal

          !attributes of the canal
          canal_id = gw_chan_canl_info(chan_id)%canals(c)
          width = gw_chan_canl_info(chan_id)%wdth(c)
          depth = gw_chan_canl_info(chan_id)%dpth(c)
          thick = gw_chan_canl_info(chan_id)%thck(c)
          day_beg = gw_chan_canl_info(chan_id)%dayb(c)
          day_end = gw_chan_canl_info(chan_id)%daye(c)

          !only proceed if canal is "on"
          if (time%day.ge.day_beg .and. time%day.le.day_end) then

            !loop through the cells that are connected to the current canal
            do k=1,gw_canl_info(canal_id)%ncon

              !id of the gwflow cell
              cell_id = gw_canl_info(canal_id)%cells(k)

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

                !attributes of the cell
                length = gw_canl_info(canal_id)%leng(k)
                stage = gw_canl_info(canal_id)%elev(k)
								bed_K = gw_canl_info(canal_id)%hydc(k)

                !calculate exchange rate Q (m3/day)
                flow_area = length * width !m2 = area of seepage
                canal_bed = stage - depth !m
                Q = 0.
                if(gw_state(cell_id)%head < canal_bed) then
                  head_diff = depth
                  Q =  bed_K * (head_diff / thick) * flow_area !canal seepage (positive Q: entering aquifer)
                elseif (gw_state(cell_id)%head > stage) then
                  head_diff = gw_state(cell_id)%head - stage
                  Q = bed_K * (head_diff / thick) * flow_area * (-1) !gw discharge (negative Q: leaving aquifer)
                elseif (gw_state(cell_id)%head > canal_bed .and. gw_state(cell_id)%head < stage) then
                  head_diff = stage - gw_state(cell_id)%head
                  Q = bed_K * (head_diff / thick) * flow_area !canal seepage (positive Q: entering aquifer)
                endif

                !store values in gwflow source/sink arrays
                if(Q < 0) then !groundwater --> canal
                  !if ((Q*-1 == 1).ge.gw_state(cell_id)%stor) then !can only remove what is there
                  if (-Q .ge. gw_state(cell_id)%stor) then !can only remove what is there
                    Q = -gw_state(cell_id)%stor
                  endif
                  gw_hyd_ss(cell_id)%canl = gw_hyd_ss(cell_id)%canl + Q
                  gw_state(cell_id)%stor = gw_state(cell_id)%stor + Q !update available groundwater in the cell
                else !canal --> groundwater (seepage)
                  if(Q > 0) then !canal seepage; remove water from channel
                    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)%canl = gw_hyd_ss(cell_id)%canl + Q !store for water balance calculations
                    ch_stor(chan_id)%flo = ch_stor(chan_id)%flo - Q !remove water from channel
                  endif
                endif
                gw_hyd_ss_yr(cell_id)%canl = gw_hyd_ss_yr(cell_id)%canl + Q !store for annual water
                gw_hyd_ss_mo(cell_id)%canl = gw_hyd_ss_mo(cell_id)%canl + Q !store for monthly water

                !heat
                if(gw_heat_flag == 1) then
                  chan_heat = ch_stor(chan_id)%temp * gw_rho * gw_cp * chan_volume !J in channel
                  if(Q < 0) then !mass is leaving the cell --> canal
                    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
                    gwheat_state(cell_id)%stor = gwheat_state(cell_id)%stor + heat_flux !update heat in the cell
                  else
                    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
                      gwheat_state(cell_id)%stor = gwheat_state(cell_id)%stor + heat_flux !update heat in the cell
                    endif
                  endif
                  !update channel heat and temperature
                  chan_heat = chan_heat + (heat_flux*(-1))
                  if(ch_stor(chan_id)%flo > 0) then
                    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 heat balance
                  gw_heat_ss(cell_id)%canl = gw_heat_ss(cell_id)%canl + heat_flux
                  gw_heat_ss_yr(cell_id)%canl = gw_heat_ss_yr(cell_id)%canl + heat_flux !J
                endif

                !calculate solute mass (g/day) transported between cell and channel
                if (gw_solute_flag == 1) then
                  if(Q < 0) then !mass is leaving the cell --> canal
                    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
                    enddo
                  else !mass entering the cell from the canal (i.e., from the channel that provides the canal water)
                    !calculate mass (g)
                    do s=1,gw_nsolute
                      solmass(s) = Q * chan_csol(s)
                    enddo
                    mass_nh3 = Q * conc_nh3
                    mass_no2 = Q * conc_no2
                    mass_dox = Q * conc_dox
                    mass_orgn = Q * conc_orgn
                    !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
                    !nh3
                    if((mass_nh3/1000.) > ch_stor(chan_id)%nh3) then
                      mass_nh3 = ch_stor(chan_id)%nh3 * 1000. !g
                    endif
                    ch_stor(chan_id)%nh3 = ch_stor(chan_id)%nh3 - (mass_nh3/1000.) !kg
                    !no2
                    if((mass_no2/1000.) > ch_stor(chan_id)%no2) then
                      mass_no2 = ch_stor(chan_id)%no2 * 1000. !g
                    endif
                    ch_stor(chan_id)%no2 = ch_stor(chan_id)%no2 - (mass_no2/1000.) !kg
                    !dox
                    if((mass_dox/1000.) > ch_stor(chan_id)%dox) then
                      mass_dox = ch_stor(chan_id)%dox * 1000. !g
                    endif
                    ch_stor(chan_id)%dox = ch_stor(chan_id)%dox - (mass_dox/1000.) !kg
                    !orgn
                    if((mass_orgn/1000.) > ch_stor(chan_id)%orgn) then
                      mass_orgn = ch_stor(chan_id)%orgn * 1000. !g
                    endif
                    ch_stor(chan_id)%orgn = ch_stor(chan_id)%orgn - (mass_orgn/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
                  endif
                  !store in mass balance arrays
                  do s=1,gw_nsolute !loop through the solutes
                    gwsol_ss(cell_id)%solute(s)%canl = gwsol_ss(cell_id)%solute(s)%canl + solmass(s)
                    gwsol_ss_sum(cell_id)%solute(s)%canl = gwsol_ss_sum(cell_id)%solute(s)%canl + solmass(s)
										gwsol_ss_sum_mo(cell_id)%solute(s)%canl = gwsol_ss_sum_mo(cell_id)%solute(s)%canl + solmass(s)
                  enddo
                endif !end solutes

              endif !check if cell is active

            enddo !go to next cell connected to the canal

          endif !check if canal is "on"

        enddo !go to next canal connected to the channel

      endif !check if canal-cell exchange is active

      return
      end subroutine gwflow_canal
