subroutine gwflow_pond !rtb gwflow

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates the volume of seepage from recharge ponds;
!!    removes water from source channel or canal diversion;
!!    writes out recharge pond water balance


      use gwflow_module
      use hydrograph_module, only : ch_stor
      use time_module
      use constituent_mass_module
      use water_allocation_module
      use climate_module
      !div_conc_cs and div_conc_salt now in gwflow_module (temporary, pending wallo integration)

      implicit none

      character(len=18) :: pond_name = ''     !       |constructed pond name
      integer :: r = 0                       !       |counter for number of recharge ponds
	    integer :: k = 0                       !       |counter for the number of cells connected to the recharge pond
	    integer :: s = 0						 !       |solute counter
	    integer :: year = 0
      integer :: day = 0
      integer :: month = 0
	    integer :: chan_id = 0                 !       |channel id
      integer :: rec_id = 0                  !       |point source if (canal diversion)
      integer :: cell_id = 0				 !       |cell id
      integer :: iwst = 0					 !       |weather station id
      integer :: dum = 0
      integer :: sol_index = 0               !       |solute counter (no3, p, salt ions, constituents)
      integer :: isalt = 0                   !       |salt ion counter
      integer :: ics = 0					 !       |constituent counter
      integer :: canal_id = 0
      real :: chan_volume = 0.               !m3     |starting volume in the source channel
      real :: cell_recharge = 0.             !m3     |recharge from the pond to the aquifer, for a single cell
      real :: pond_recharge = 0.			 !m3     |total recharge from the pond to the aquifer
      real :: div_specified = 0.             !m3     |specified diversion volume in gwflow.ponds file
      real :: div_added = 0.                 !m3     |actual volume added to the recharge pond
      real :: pond_evap = 0.                 !m3     |water evaporated from the pond during the day
      real :: pond_rain = 0.                 !m3     |rainfall added to the pond during the day
      real :: pond_volume = 0.               !m3     |pond volume before recharge occurs
      real :: sol_conc = 0.                  !g/m3   |solute concentration in the source water
      real :: sol_mass = 0.                  !kg     |solute mass removed from the source added --> added to recharge pond
      real :: div_mass(20) = 0.              !kg     |solute mass added to the recharge pond
      real :: rech_mass(20) = 0.             !kg     |solute mass leaching from the pond to the water table
      real :: rech_mass_cell(20) = 0.        !g      |solute mass leaching from the pond to an individual cell




      !only proceed if the recharge pond flag has been activated ----------------------------------------------------------------
      if (gw_pond_flag == 1) then

			  !read the diverted volumes (m3) for the current day (from pond_div.gw, if provided) ----------------------------
			  if (gw_pond_div_flag == 1) then
				  read(in_ponds,*) year,month,day,(gw_pond_info(r)%div,r=1,gw_npond)
			  else
				  do r=1,gw_npond
					  gw_pond_info(r)%div = 0.
				  enddo
			  endif

        !loop through the recharge ponds ------------------------------------------------------------------------------
			  do r=1,gw_npond

				  !zero out values
				  pond_rain = 0.
				  div_added = 0.
					pond_evap = 0.
					pond_recharge = 0.
					gw_pond_info(r)%div_uns = 0.
					div_mass = 0.
					rech_mass = 0.

				  !only proceed if the recharge pond is in operation
					if(gw_daycount .ge. gw_pond_info(r)%dy_start) then

				  !add diverted water to recharge pond storage (m3) ---------------------------------------
				  div_specified = gw_pond_info(r)%div
					if(div_specified > 0) then
					  dum = 10
					endif
					div_added = 0.
					!channel source
					if(gw_pond_info(r)%chan > 0) then
					  chan_id = gw_pond_info(r)%chan
            chan_volume = ch_stor(chan_id)%flo !current channel storage (m3)
						!determine water volume added to the recharge pond
						div_added = div_specified
						if(div_specified > ch_stor(chan_id)%flo) then
						  div_added = ch_stor(chan_id)%flo
						endif
						gw_pond_info(r)%div_uns = div_specified - div_added !unsatisfied diversion
						!remove from the channel --> add to the recharge pond
	          ch_stor(chan_id)%flo = ch_stor(chan_id)%flo - div_added
						gw_pond_info(r)%stor = gw_pond_info(r)%stor + div_added
						!remove salt mass from the channel --> add to the recharge pond
						if (gw_solute_flag == 1) then
							if(chan_volume > 10.) then
								!no3
								sol_conc = (ch_stor(chan_id)%no3 * 1000.) / chan_volume !no3 g/m3 in channel
								sol_mass = (sol_conc * div_added) / 1000. !kg removed from channel
								if(sol_mass > ch_stor(chan_id)%no3) then
								  sol_mass = ch_stor(chan_id)%no3
								endif
								ch_stor(chan_id)%no3 = ch_stor(chan_id)%no3 - sol_mass
								gw_pond_info(r)%sol_mass(1) = gw_pond_info(r)%sol_mass(1) + sol_mass !kg added to pond
								div_mass(1) = sol_mass
								!p
								sol_conc = (ch_stor(chan_id)%solp * 1000.) / chan_volume !p g/m3 in channel
								sol_mass = (sol_conc * div_added) / 1000. !kg removed from channel
								if(sol_mass > ch_stor(chan_id)%solp) then
								  sol_mass = ch_stor(chan_id)%solp
								endif
								ch_stor(chan_id)%solp = ch_stor(chan_id)%solp - sol_mass
								gw_pond_info(r)%sol_mass(2) = gw_pond_info(r)%sol_mass(2) + sol_mass !kg added to pond
								div_mass(2) = sol_mass
								!salt
								sol_index = 2
								if (gwsol_salt == 1) then
									do isalt=1,cs_db%num_salts
										sol_index = sol_index + 1
										sol_conc = (ch_water(chan_id)%salt(isalt) * 1000.) / chan_volume !g/m3 in channel water
										sol_mass = (sol_conc * div_added) / 1000. !kg removed from channel
										if(sol_mass > ch_water(chan_id)%salt(isalt)) then
										  sol_mass = ch_water(chan_id)%salt(isalt)
										endif
										ch_water(chan_id)%salt(isalt) = ch_water(chan_id)%salt(isalt) - sol_mass
										gw_pond_info(r)%sol_mass(sol_index) = gw_pond_info(r)%sol_mass(sol_index) + sol_mass !kg added to pond
										div_mass(sol_index) = sol_mass
									enddo
								endif
								!constituents
								if (gwsol_cons == 1) then
									do ics=1,cs_db%num_cs
										sol_index = sol_index + 1
										sol_conc = (ch_water(chan_id)%cs(ics) * 1000.) / chan_volume !g/m3 in channel water
										sol_mass = (sol_conc * div_added) / 1000. !kg removed from channel
										if(sol_mass > ch_water(chan_id)%cs(ics)) then
										  sol_mass = ch_water(chan_id)%cs(ics)
										endif
										ch_water(chan_id)%cs(ics) = ch_water(chan_id)%cs(ics) - sol_mass
										gw_pond_info(r)%sol_mass(sol_index) = gw_pond_info(r)%sol_mass(sol_index) + sol_mass !kg added to pond
										div_mass(sol_index) = sol_mass
									enddo
								endif
							endif
						endif
					!canal diversion source
					elseif(gw_pond_info(r)%canal > 0) then
					  !determine the canal id
					  canal_id = gw_pond_info(r)%canal
						if(gw_canl_div_info(canal_id)%stor > 0) then
						  div_added = div_specified
						  if(div_specified > gw_canl_div_info(canal_id)%stor) then
							  div_added = gw_canl_div_info(canal_id)%stor
							endif
							gw_pond_info(r)%div_uns = div_specified - div_added !unsatisfied diversion
							!remove from the canal diversion --> add to the recharge pond
							gw_canl_div_info(canal_id)%stor = gw_canl_div_info(canal_id)%stor - div_added
							gw_canl_div_info(canal_id)%out_pond = gw_canl_div_info(canal_id)%out_pond + div_added
						  gw_pond_info(r)%stor = gw_pond_info(r)%stor + div_added
							if (gw_solute_flag == 1) then
							  rec_id = gw_canl_div_info(canal_id)%divr
								!salt
								sol_index = 2
								if (gwsol_salt == 1) then
									do isalt=1,cs_db%num_salts
										sol_index = sol_index + 1
										sol_mass = (div_conc_salt(isalt,rec_id) * div_added) / 1000. !kg removed from channel
										gw_pond_info(r)%sol_mass(sol_index) = gw_pond_info(r)%sol_mass(sol_index) + sol_mass !kg added to pond
										div_mass(sol_index) = sol_mass
									enddo
								endif
								!constituents
								if (gwsol_cons == 1) then
									do ics=1,cs_db%num_cs
										sol_index = sol_index + 1
										sol_mass = (div_conc_cs(ics,rec_id) * div_added) / 1000. !kg removed from channel
										gw_pond_info(r)%sol_mass(sol_index) = gw_pond_info(r)%sol_mass(sol_index) + sol_mass !kg added to pond
										div_mass(sol_index) = sol_mass
									enddo
								endif
							endif

						endif
					!outside water source
					elseif(gw_pond_info(r)%unl == 1) then
					  gw_pond_info(r)%stor = gw_pond_info(r)%stor + div_specified
						!solute mass
						if (gw_solute_flag == 1) then
						  !no3
						  sol_conc = gw_pond_info(r)%unl_conc(1) !g/m3
							sol_mass = (sol_conc * div_specified) / 1000. !kg
							gw_pond_info(r)%sol_mass(1) = gw_pond_info(r)%sol_mass(1) + sol_mass !kg added to pond
							!p
							sol_conc = gw_pond_info(r)%unl_conc(2) !g/m3
							sol_mass = (sol_conc * div_specified) / 1000. !kg
							gw_pond_info(r)%sol_mass(2) = gw_pond_info(r)%sol_mass(2) + sol_mass !kg added to pond
							!salt
							sol_index = 2
							if (gwsol_salt == 1) then
								do isalt=1,cs_db%num_salts
									sol_index = sol_index + 1
									sol_conc = gw_pond_info(r)%unl_conc(sol_index) !g/m3
									sol_mass = (sol_conc * div_specified) / 1000. !kg
									gw_pond_info(r)%sol_mass(sol_index) = gw_pond_info(r)%sol_mass(sol_index) + sol_mass !kg added to pond
									div_mass(sol_index) = sol_mass
								enddo
							endif
							!constituents
							if (gwsol_cons == 1) then
								do ics=1,cs_db%num_cs
									sol_index = sol_index + 1
									sol_conc = gw_pond_info(r)%unl_conc(sol_index) !g/m3
									sol_mass = (sol_conc * div_specified) / 1000. !kg
									gw_pond_info(r)%sol_mass(sol_index) = gw_pond_info(r)%sol_mass(sol_index) + sol_mass !kg added to pond
									div_mass(sol_index) = sol_mass
								enddo
							endif
						endif
					else !no water source - no water added to recharge pond
					  goto 10
					endif

					!rainfall
					iwst = gw_pond_info(r)%wsta
					pond_rain = (wst(iwst)%weat%precip/1000.) * gw_pond_info(r)%area !m3
					gw_pond_info(r)%stor = gw_pond_info(r)%stor + pond_rain

					!evaporation
					iwst = gw_pond_info(r)%wsta
					pond_evap = gw_pond_info(r)%evap_co * (wst(iwst)%weat%pet/1000.) * gw_pond_info(r)%area !m3
					if(pond_evap > gw_pond_info(r)%stor) then
					  pond_evap = gw_pond_info(r)%stor
					endif
					gw_pond_info(r)%stor = gw_pond_info(r)%stor - pond_evap

				  !loop through the number of cells connected to the recharge pond
				  pond_recharge = 0.
				  do k=1,gw_pond_info(r)%ncell
					  !calculate recharge to the cell
					  cell_id = gw_pond_info(r)%cells(k)
						if(gw_state(cell_id)%stat == 1) then !only proceed if the cell is active

						  !store initial pond storage volume
						  pond_volume = gw_pond_info(r)%stor !m3

						  !calculate pond seepage (recharge)
						  cell_recharge = gw_pond_info(r)%bed_k * gw_pond_info(r)%conn_area(k) !m3/day
							!check against pond volume - can only remove what is available in storage
							if(cell_recharge > gw_pond_info(r)%stor) then
								cell_recharge = gw_pond_info(r)%stor
							endif

							!store for water balance and output
							gw_hyd_ss(cell_id)%pond = gw_hyd_ss(cell_id)%pond + cell_recharge
							gw_hyd_ss_yr(cell_id)%pond = gw_hyd_ss_yr(cell_id)%pond + cell_recharge
							gw_hyd_ss_mo(cell_id)%pond = gw_hyd_ss_mo(cell_id)%pond + cell_recharge
							!add to pond recharge --> remove from pond storage
							pond_recharge = pond_recharge + cell_recharge !add to total pond recharge
							gw_pond_info(r)%stor = gw_pond_info(r)%stor - cell_recharge !remove from storage

							!solute mass in recharge water (added to cell)
							if(pond_volume> 0) then
							!no3
							if (gw_solute_flag == 1) then
								sol_conc = (gw_pond_info(r)%sol_mass(1)*1000.) / pond_volume !g/m3
								sol_mass = (sol_conc*cell_recharge)/1000. !kg/day
								if(sol_mass > gw_pond_info(r)%sol_mass(1)) then
								  sol_mass = gw_pond_info(r)%sol_mass(1)
								endif
								gw_pond_info(r)%sol_mass(1) = gw_pond_info(r)%sol_mass(1) - sol_mass !kg leached from pond
								rech_mass_cell(1) = sol_mass * 1000. !g
								rech_mass(1) = rech_mass(1) + sol_mass !add to total from pond
								!p
								sol_conc = (gw_pond_info(r)%sol_mass(2)*1000.) / pond_volume !g/m3
								sol_mass = (sol_conc*cell_recharge)/1000. !kg/day
								if(sol_mass > gw_pond_info(r)%sol_mass(2)) then
								  sol_mass = gw_pond_info(r)%sol_mass(2)
								endif
								gw_pond_info(r)%sol_mass(2) = gw_pond_info(r)%sol_mass(2) - sol_mass !kg leached from pond
								rech_mass_cell(2) = sol_mass * 1000. !g
								rech_mass(2) = rech_mass(2) + sol_mass !add to total from pond
								!salt
								sol_index = 2
								if (gwsol_salt == 1) then
									do isalt=1,cs_db%num_salts
									  sol_index = sol_index + 1
									  sol_conc = (gw_pond_info(r)%sol_mass(sol_index)*1000.) / pond_volume !g/m3
										sol_mass = (sol_conc*cell_recharge)/1000. !kg/day
										if(sol_mass > gw_pond_info(r)%sol_mass(sol_index)) then
											sol_mass = gw_pond_info(r)%sol_mass(sol_index)
										endif
										gw_pond_info(r)%sol_mass(sol_index) = gw_pond_info(r)%sol_mass(sol_index) - sol_mass !kg leached from pond
										rech_mass_cell(sol_index) = sol_mass * 1000. !g
										rech_mass(sol_index) = rech_mass(sol_index) + sol_mass !add to total from pond
									enddo
								endif
								!constituents
								if (gwsol_cons == 1) then
									do ics=1,cs_db%num_cs
										sol_index = sol_index + 1
									  sol_conc = (gw_pond_info(r)%sol_mass(sol_index)*1000.) / pond_volume !g/m3
										sol_mass = (sol_conc*cell_recharge)/1000. !kg/day
										if(sol_mass > gw_pond_info(r)%sol_mass(sol_index)) then
											sol_mass = gw_pond_info(r)%sol_mass(sol_index)
										endif
										gw_pond_info(r)%sol_mass(sol_index) = gw_pond_info(r)%sol_mass(sol_index) - sol_mass !kg leached from pond
										rech_mass_cell(sol_index) = sol_mass * 1000. !g
										rech_mass(sol_index) = rech_mass(sol_index) + sol_mass !add to total from pond
									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)%pond = rech_mass_cell(s) !g/day
									gwsol_ss_sum(cell_id)%solute(s)%pond = gwsol_ss_sum(cell_id)%solute(s)%pond + rech_mass_cell(s) !g/day
									gwsol_ss_sum_mo(cell_id)%solute(s)%pond = gwsol_ss_sum_mo(cell_id)%solute(s)%pond + rech_mass_cell(s) !g/day
								enddo
							endif
							endif !if pond has water

						endif
					enddo !go to next cell

					endif !check if recharge pond is in operation

					if(div_added < 0) then
					  dum = 10
					endif

					!recharge pond water balance - output
          write(pond_name,'(a5,i4.4)') 'pond_',r
          if(gwflag_flux == 1) then
					write(out_pond_bal,8100) time%day,time%mo,time%day_mo,time%yrc, &
					    r,r,pond_name, &
					    gw_pond_info(r)%area, &
					    gw_pond_info(r)%stor, &
					    pond_rain, &
					    div_added, &
					    pond_evap, &
					    pond_recharge, &
					    gw_pond_info(r)%div, &
					    gw_pond_info(r)%div_uns

					!recharge pond solute mass balance - output
					do s=1,gw_nsolute
					  write(out_pond_sol,8101) time%day,time%mo,time%day_mo,time%yrc, &
					      r,r,pond_name, &
					      gw_pond_info(r)%area, &
					      gw_pond_info(r)%stor, &
					      gwsol_nm(s), &
					      gw_pond_info(r)%sol_mass(s), &
					      div_mass(s), &
					      rech_mass(s)
					enddo
          endif !gwflag_flux

10  	  enddo !go to next recharge pond

        !write out mass (kg) and concentration (g/m3) for each pond
        if(gwflag_flux == 1) then
				!no3
				write(out_pond_mass,102) time%yrc,time%mo,time%day,gwsol_nm(1),(gw_pond_info(r)%sol_mass(1),r=1,gw_npond)
				do r=1,gw_npond
					if(gw_pond_info(r)%stor > 0) then
						gw_pond_info(r)%sol_conc(1) = (gw_pond_info(r)%sol_mass(1)*1000.) / gw_pond_info(r)%stor
					else
						gw_pond_info(r)%sol_conc(1) = 0.
					endif
				enddo
				write(out_pond_conc,102) time%yrc,time%mo,time%day,gwsol_nm(1),(gw_pond_info(r)%sol_conc(1),r=1,gw_npond)
				!p
				write(out_pond_mass,102) time%yrc,time%mo,time%day,gwsol_nm(2),(gw_pond_info(r)%sol_mass(2),r=1,gw_npond)
				do r=1,gw_npond
					if(gw_pond_info(r)%stor > 0) then
						gw_pond_info(r)%sol_conc(2) = (gw_pond_info(r)%sol_mass(2)*1000.) / gw_pond_info(r)%stor
					else
						gw_pond_info(r)%sol_conc(2) = 0.
					endif
				enddo
				write(out_pond_conc,102) time%yrc,time%mo,time%day,gwsol_nm(2),(gw_pond_info(r)%sol_conc(2),r=1,gw_npond)
				!salt ions
				sol_index = 2
				if (gwsol_salt == 1) then
					do isalt=1,cs_db%num_salts
						sol_index = sol_index + 1
						write(out_pond_mass,102) time%yrc,time%mo,time%day,gwsol_nm(sol_index),(gw_pond_info(r)%sol_mass(sol_index),r=1,gw_npond)
						do r=1,gw_npond
							if(gw_pond_info(r)%stor > 0) then
								gw_pond_info(r)%sol_conc(sol_index) = (gw_pond_info(r)%sol_mass(sol_index)*1000.) / gw_pond_info(r)%stor
							else
								gw_pond_info(r)%sol_conc(sol_index) = 0.
							endif
						enddo
						write(out_pond_conc,102) time%yrc,time%mo,time%day,gwsol_nm(sol_index),(gw_pond_info(r)%sol_conc(sol_index),r=1,gw_npond)
					enddo
				endif
				!constituents
				if (gwsol_cons == 1) then
					do ics=1,cs_db%num_cs
						sol_index = sol_index + 1
						write(out_pond_mass,102) time%yrc,time%mo,time%day,gwsol_nm(sol_index),(gw_pond_info(r)%sol_mass(sol_index),r=1,gw_npond)
						do r=1,gw_npond
							if(gw_pond_info(r)%stor > 0) then
								gw_pond_info(r)%sol_conc(sol_index) = (gw_pond_info(r)%sol_mass(sol_index)*1000.) / gw_pond_info(r)%stor
							else
								gw_pond_info(r)%sol_conc(sol_index) = 0.
							endif
						enddo
						write(out_pond_conc,102) time%yrc,time%mo,time%day,gwsol_nm(sol_index),(gw_pond_info(r)%sol_conc(sol_index),r=1,gw_npond)
					enddo
				endif
        endif !gwflag_flux


      endif !check if recharge pond seepage is active


      !standard SWAT+ data row: water balance
8100  format(4i6,2i8,a18,50e13.4)
      !standard SWAT+ data row: solute balance (area, storage real; solute name string; mass values)
8101  format(4i6,2i8,a18,2e13.4,a13,50e13.4)
      !wide format: pond mass/conc per pond
102   format (3i8,5x,a10,5000(e15.4))


      return
      end subroutine gwflow_pond