      subroutine gwflow_canal_div !rtb gwflow

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates the water exchange volume between irrigation canals and connected grid cells
!!    for canals that are connected to a point source diversion (recall object); seepage water should be
!!    removed from the diverted water volume.

      use gwflow_module
      use hydrograph_module, only : irrig
      use time_module
      use constituent_mass_module
			use hru_module, only : hru
			!div_conc_cs and div_conc_salt now in gwflow_module (temporary, pending wallo integration)
			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 :: i = 0                       !       |cell counter
	  integer :: j = 0
      character(len=18) :: canal_name = ''   !       |constructed canal name
      integer :: s = 0                       !       |counter of groundwater solutes
      integer :: cell_id = 0                 !       |cell in connection with the canal
	  integer :: irec = 0                    !       |recall ID for canal diversion
	  integer :: sol_index = 0
	  integer :: ics = 0
	  integer :: isalt = 0
	  integer :: hru_id = 0
      integer :: canal_id = 0
			integer :: wetland = 0                 !       |wetland flag
			integer :: dum = 0
      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 :: solmass(100) = 0.               !g      |solute mass transferred
      real :: heat_flux = 0.                  !J      |heat in groundawter-canal exchange water
			real :: canal_area = 0.								 !m2     |total area irrigated by canal
			real :: irrig_depth = 0.                !mm     |depth of irrigation water for the HRUs irrigated by a canal
			real :: irrig_volm = 0.
			real :: irrig_conc = 0.
			real :: irrig_mass = 0.
			real :: canal_conc = 0.
			real :: mass_div = 0.
            real :: mass_stor = 0.
            real :: mass_pond = 0.
            real :: mass_seep = 0.
            real :: mass_irrg = 0.
            real :: mass_ret = 0.



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


			  !first: calculate seepage from each canal ---------------------------------------------------------------------
			  !(seepage only occurs if the canal has water; based on point diversions)

        !loop through the canal cells
        do i=1,gw_canal_ncells_div

				  !canal associated with grid cell
				  canal_id = gw_canl_div_cell(i)%canal_id

					!only proceed if the canal has water
					if(gw_canl_div_info(canal_id)%stor > 0) then

            cell_id = gw_canl_div_cell(i)%cell_id
            if (gw_state(cell_id)%stat == 1) then !if cell is active
              !canal attributes
						  width = gw_canl_div_info(canal_id)%width
              depth = gw_canl_div_info(canal_id)%depth
              thick = gw_canl_div_info(canal_id)%thick
							bed_K = gw_canl_div_info(canal_id)%bed_K
						  !cell attributes
              length = gw_canl_div_cell(i)%leng
              canal_bed = gw_canl_div_cell(i)%elev
							stage = canal_bed + depth
              !calculate exchange rate Q (m3/day)
              flow_area = length * width !m2 = area of seepage
              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 the seepage is positive, then water is taken from the volume of water already diverted from
							!channels (point source diversion = recall object)
              if(Q < 0) then !groundwater --> canal
                if (-Q .ge.gw_state(cell_id)%stor) then !can only remove what is there
                  Q = -gw_state(cell_id)%stor
                endif
							elseif(Q > 0) then !canal --> groundwater
							  !remove seepage water from the diversion volume
								if(Q > gw_canl_div_info(canal_id)%stor) then
								  Q = gw_canl_div_info(canal_id)%stor
								endif
								gw_canl_div_info(canal_id)%stor = gw_canl_div_info(canal_id)%stor - Q
								gw_canl_div_info(canal_id)%out_seep = gw_canl_div_info(canal_id)%out_seep + Q
							endif
							gw_state(cell_id)%stor = gw_state(cell_id)%stor + Q !update available groundwater in the cell
              gw_hyd_ss(cell_id)%canl = gw_hyd_ss(cell_id)%canl + Q
              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
                if(Q < 0) then
                  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)%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
                endif
              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 diverted solute mass)
								  !calculate mass (g)
								  solmass(1) = Q * canal_out_conc(1) !no3
									solmass(2) = Q * canal_out_conc(2) !p
									sol_index = 2
                  !salts
                  if (gwsol_salt == 1) then
									  irec = gw_canl_div_info(canal_id)%divr
									  if(irec > 0) then
										  do isalt=1,cs_db%num_salts
                        sol_index = sol_index + 1
                        solmass(sol_index) = Q * div_conc_salt(isalt,irec)
                      enddo
										else
                      do isalt=1,cs_db%num_salts
                        sol_index = sol_index + 1
                        solmass(sol_index) = Q * canal_out_conc(sol_index)
                      enddo
										endif
                  endif
                  !constituents
                  if (gwsol_cons == 1) then
									  irec = gw_canl_div_info(canal_id)%divr
									  if(irec > 0) then
										  do ics=1,cs_db%num_cs
                        sol_index = sol_index + 1
                        solmass(sol_index) = Q * div_conc_cs(ics,irec)
                      enddo
										else
                      do ics=1,cs_db%num_cs
                        sol_index = sol_index + 1
                        solmass(sol_index) = Q * canal_out_conc(sol_index)
                      enddo
										endif
                  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
          endif !check if canal has water

        enddo !go to next canal cell


				!(irrigation removed -- handled by wallo water allocation system)

        !write out daily volumes and fluxes for each canal
        do i=1,gw_ncanal
          write(canal_name,'(a6,i4.4)') 'canal_',i
          if(gwflag_flux == 1) then
					write(out_canal_bal,8100) time%day,time%mo,time%day_mo,time%yrc, &
					    i,i,canal_name, &
					    gw_canl_div_info(i)%div, &
					    gw_canl_div_info(i)%stor, &
					    gw_canl_div_info(i)%out_pond, &
					    gw_canl_div_info(i)%out_seep
          endif !gwflag_flux

					!write solute mass balance (seepage mass only, irrigation handled by wallo)
					if (gw_solute_flag == 1) then
						irec = gw_canl_div_info(i)%divr
						sol_index = 2
						if (gwsol_salt == 1) then
							do isalt=1,cs_db%num_salts
                sol_index = sol_index + 1
								canal_conc = div_conc_salt(isalt,irec)
								mass_div = (gw_canl_div_info(i)%div * canal_conc) / 1000.
								mass_stor = (gw_canl_div_info(i)%stor * canal_conc) / 1000.
								mass_pond = (gw_canl_div_info(i)%out_pond * canal_conc) / 1000.
								mass_seep = (gw_canl_div_info(i)%out_seep * canal_conc) / 1000.
								if(gwflag_flux == 1) then
								write(out_canal_sol,8101) time%day,time%mo,time%day_mo,time%yrc, &
								    i,i,canal_name,gwsol_nm(sol_index), &
								    mass_div,mass_stor,mass_pond,mass_seep
								endif
              enddo
            endif
						if (gwsol_cons == 1) then
							do ics=1,cs_db%num_cs
                sol_index = sol_index + 1
								canal_conc = div_conc_cs(ics,irec)
								mass_div = (gw_canl_div_info(i)%div * canal_conc) / 1000.
								mass_stor = (gw_canl_div_info(i)%stor * canal_conc) / 1000.
								mass_pond = (gw_canl_div_info(i)%out_pond * canal_conc) / 1000.
								mass_seep = (gw_canl_div_info(i)%out_seep * canal_conc) / 1000.
								if(gwflag_flux == 1) then
								write(out_canal_sol,8101) time%day,time%mo,time%day_mo,time%yrc, &
								    i,i,canal_name,gwsol_nm(sol_index), &
								    mass_div,mass_stor,mass_pond,mass_seep
								endif
              enddo
            endif
					endif

				enddo !go to next canal

      endif !check if canal-cell exchange is active



      !standard SWAT+ data row: water balance (6 data cols)
8100  format(4i6,2i8,a18,50e13.4)
      !standard SWAT+ data row: solute balance (solute name + 6 data cols)
8101  format(4i6,2i8,a18,a13,50e13.4)

      return
      end subroutine gwflow_canal_div
