subroutine gwflow_solute !rtb gwflow

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates solute advection, dispersion, chemical
!!    reactions, and sorption for groundwater solute transport. Called from
!!    gwflow_lateral once per flow time step. Contains its own sub-timestep
!!    loop (num_ts_transport steps). Uses cell_con%latl and cell_con%sat
!!    populated during the Darcy flow calculation in the same time step.

      use gwflow_module
      use time_module

      implicit none

      integer :: i = 0                    !       |cell counter
      integer :: j = 0                    !       |general counter
      integer :: k = 0                    !       |connected-cell counter
      integer :: s = 0                    !       |solute counter
      integer :: t = 0                    !       |transport time step counter
      integer :: cell_id = 0              !       |id of connected cell
      real :: gw_trans_time_step = 0.     !days   |length of transport sub-timestep
      real :: time_fraction = 0.          !       |fraction of flow step for interpolation
      real :: gw_volume_old = 0.          !m3     |gw volume at start of flow step
      real :: gw_volume_new = 0.          !m3     |gw volume at end of flow step
      real :: gw_volume_inter = 0.        !m3     |interpolated gw volume at transport step
      real :: mass_adv(100) = 0.          !g/day  |advective mass flux per solute
      real :: mass_dsp(100) = 0.          !g/day  |dispersive mass flux per solute
      real :: m_change(100) = 0.          !g      |mass change per solute
      real :: del_no_sorp = 0.            !g      |mass change without sorption
      real :: mass_sorb(100) = 0.         !g      |mass removed by sorption per solute
      real :: mass_rct_local(100) = 0.    !g/day  |mass from reactions (unused local; module mass_rct used)
      real :: Q_cell = 0.                 !m3/day |lateral flow between cells
      real :: face_sat = 0.               !m      |saturated thickness at cell interface
      real :: area1 = 0.                  !m2     |area of connected cell
      real :: area2 = 0.                  !m2     |area of current cell
      real :: area = 0.                   !m2     |smaller of the two areas
      real :: conn_length = 0.            !m      |connection length

      !compute transport sub-timestep size
      gw_trans_time_step = gw_time_step / num_ts_transport

      !transport sub-timestep loop
      do t=1,num_ts_transport
        do i=1,ncell

          if(gw_state(i)%stat == 1) then !interior cell

            !calculate old and new groundwater volume in the cell (m3)
            if(gw_state(i)%hold > gw_state(i)%botm) then
              gw_volume_old = gw_state(i)%area * (gw_state(i)%hold - gw_state(i)%botm) * gw_state(i)%spyd
            else
              gw_volume_old = 0.
            endif
            if(gw_state(i)%head > gw_state(i)%botm) then
              gw_volume_new = gw_state(i)%area * (gw_state(i)%head - gw_state(i)%botm) * gw_state(i)%spyd
            else
              gw_volume_new = 0.
            endif

            !calculate groundwater volume for the current transport time step (via interpolation)
            time_fraction = real(t)/real(num_ts_transport)
            gw_volume_inter = gw_volume_old + ((gw_volume_new-gw_volume_old)*time_fraction)

            !advection transport
            mass_adv = 0.
            do k=1,gw_state(i)%ncon !loop through the connected cells
              cell_id = cell_con(i)%cell_id(k) !id of the connected cell
              Q_cell = cell_con(i)%latl(k) !m3/day flow between cell and connected cell
              if(Q_cell > 0) then !mass entering cell
                do s=1,gw_nsolute
                  mass_adv(s) = mass_adv(s) + (Q_cell * gwsol_state(cell_id)%solute(s)%conc) !g
                enddo
              else !mass leaving cell
                do s=1,gw_nsolute
                  mass_adv(s) = mass_adv(s) + (Q_cell * gwsol_state(i)%solute(s)%conc) !g
                enddo
              endif
            enddo !go to next connected cell

            !dispersion transport
            mass_dsp = 0.
            do k=1,gw_state(i)%ncon !loop through the connected cells
              cell_id = cell_con(i)%cell_id(k) !id of the connected cell
              face_sat = cell_con(i)%sat(k) !m saturated thickness at interface between cells
              area1 = gw_state(cell_id)%area !area of connected cell
              area2 = gw_state(i)%area !area of current cell
              area = min(area1,area2) !smaller of the two
              conn_length = sqrt(area)
              do s=1,gw_nsolute !loop through the solutes
                mass_dsp(s) = mass_dsp(s) + (gw_long_disp * &
                    ((gwsol_state(cell_id)%solute(s)%conc - gwsol_state(i)%solute(s)%conc)/conn_length) * face_sat) !g
              enddo
            enddo !go to next connected cell

            !chemical reactions --> fill in mass_rct and mass_min
            cell_id = i
            mass_rct = 0.
            mass_min = 0.
            call gwflow_chem(cell_id,gw_volume_inter)

            !calculate change in mass (g)
            do s=1,gw_nsolute !loop through the solutes
              m_change(s) = (mass_adv(s) + mass_dsp(s) + mass_rct(s) + mass_min(s) + &
                  gwsol_ss(i)%solute(s)%totl) * (gw_trans_time_step/gwsol_sorb(s))
            enddo

            !calculate mass removed due to sorption (g)
            do s=1,gw_nsolute !loop through the solutes
              del_no_sorp = (mass_adv(s) + mass_dsp(s) + mass_rct(s) + mass_min(s) + &
                  gwsol_ss(i)%solute(s)%totl) * gw_trans_time_step
              mass_sorb(s) = del_no_sorp - m_change(s)
            enddo

            !calculate new mass in the cell (g)
            do s=1,gw_nsolute !loop through the solutes
              gwsol_state(i)%solute(s)%mass = gwsol_state(i)%solute(s)%mass + m_change(s)
              if(gwsol_state(i)%solute(s)%mass < 0) then
                gwsol_state(i)%solute(s)%mass = 0.
              endif
            enddo

            !calculate new concentration (g/m3)
            if(gw_volume_inter > 0) then
              do s=1,gw_nsolute !loop through the solutes
                gwsol_state(i)%solute(s)%cnew = gwsol_state(i)%solute(s)%mass / gw_volume_inter
              enddo
            else
              do s=1,gw_nsolute !loop through the solutes
                gwsol_state(i)%solute(s)%cnew = 0.
                gwsol_state(i)%solute(s)%mass = 0.
              enddo
            endif

            !store values for mass budget analysis
            do s=1,gw_nsolute !loop through the solutes
              gwsol_ss(i)%solute(s)%advn = gwsol_ss(i)%solute(s)%advn + &
                  (mass_adv(s)*(gw_trans_time_step/gwsol_sorb(s)))
              gwsol_ss(i)%solute(s)%disp = gwsol_ss(i)%solute(s)%disp + &
                  (mass_dsp(s)*(gw_trans_time_step/gwsol_sorb(s)))
              if(mass_rct(s) > 0) then
                gwsol_ss(i)%solute(s)%rcti = gwsol_ss(i)%solute(s)%rcti + &
                    (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !produced
              else
                gwsol_ss(i)%solute(s)%rcto = gwsol_ss(i)%solute(s)%rcto + &
                    (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !consumed
              endif
              gwsol_ss(i)%solute(s)%minl = gwsol_ss(i)%solute(s)%minl + &
                  (mass_min(s)*(gw_trans_time_step/gwsol_sorb(s)))
              gwsol_ss(i)%solute(s)%sorb = gwsol_ss(i)%solute(s)%sorb + mass_sorb(s)
            enddo

            !track for annual and monthly write-out
            do s=1,gw_nsolute !loop through the solutes
              if(mass_rct(s) > 0) then
                gwsol_ss_sum(i)%solute(s)%rcti = gwsol_ss_sum(i)%solute(s)%rcti + &
                    (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !produced
                gwsol_ss_sum_mo(i)%solute(s)%rcti = gwsol_ss_sum_mo(i)%solute(s)%rcti + &
                    (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !produced
              else
                gwsol_ss_sum(i)%solute(s)%rcto = gwsol_ss_sum(i)%solute(s)%rcto + &
                    (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !consumed
                gwsol_ss_sum_mo(i)%solute(s)%rcto = gwsol_ss_sum_mo(i)%solute(s)%rcto + &
                    (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !consumed
              endif
              gwsol_ss_sum(i)%solute(s)%minl = gwsol_ss_sum(i)%solute(s)%minl + &
                  (mass_min(s)*(gw_trans_time_step/gwsol_sorb(s)))
              gwsol_ss_sum_mo(i)%solute(s)%minl = gwsol_ss_sum_mo(i)%solute(s)%minl + &
                  (mass_min(s)*(gw_trans_time_step/gwsol_sorb(s)))
              gwsol_ss_sum(i)%solute(s)%sorb = gwsol_ss_sum(i)%solute(s)%sorb + mass_sorb(s)
              gwsol_ss_sum_mo(i)%solute(s)%sorb = gwsol_ss_sum_mo(i)%solute(s)%sorb + mass_sorb(s)
            enddo

          elseif(gw_state(i)%stat == 2) then !constant concentration cell
            do s=1,gw_nsolute
              gwsol_state(i)%solute(s)%cnew = 0.
            enddo
          endif

        enddo !go to next cell

        !store new concentration values into regular array
        do i=1,ncell
          do s=1,gw_nsolute !loop through the solutes
            gwsol_state(i)%solute(s)%conc = gwsol_state(i)%solute(s)%cnew
          enddo
        enddo

      enddo !go to next transport time step

      return
      end subroutine gwflow_solute