subroutine gwflow_heat !rtb gwflow

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates heat advection and dispersion for groundwater
!!    heat transport. Called from gwflow_lateral once per flow time step, after
!!    the Darcy head update. Uses cell_con%latl and cell_con%sat populated
!!    during the flow calculation in the same time step.

      use gwflow_module
      use time_module

      implicit none

      integer :: i = 0                    !       |cell counter
      integer :: k = 0                    !       |connected-cell counter
      integer :: cell_id = 0              !       |id of connected cell
      real :: heat_adv = 0.               !J/day  |total advective heat flux for current cell
      real :: cell_adv = 0.               !J/day  |advective heat flux from one connection
      real :: heat_dsp = 0.               !J/day  |total dispersive heat flux for current cell
      real :: Q_heat = 0.                 !J/day  |dispersive heat flux from one connection
      real :: heat_change = 0.            !J      |net change in heat storage for current cell
      real :: face_thmc = 0.              !J/d-m-T|harmonic mean thermal conductivity at interface
      real :: Q_cell = 0.                 !m3/day |lateral flow between cells
      real :: dist_x = 0.                 !m      |x distance between cell centroids
      real :: dist_y = 0.                 !m      |y distance between cell centroids
      real :: grad_distance = 0.          !m      |distance between cell centroids
      real :: gradient = 0.               !T/m    |temperature gradient
      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
      real :: face_sat = 0.               !m      |saturated thickness at interface
      real :: flow_area = 0.              !m2     |cross-section area at interface

      !store current heat (J) - used to track heat storage during cell loop
      do i=1,ncell
        heat_cell(i) = gwheat_state(i)%stor !J
      enddo

      !loop through cells - update heat and temperature
      do i=1,ncell
        if(gw_state(i)%stat == 1) then !interior cell

          !starting heat in the cell
          heat_cell(i) = gwheat_state(i)%stor !J

          !advection heat transport (transported with groundwater flow)
          heat_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 !heat entering cell
              cell_adv = Q_cell * gw_rho * gw_cp * gwheat_state(cell_id)%temp !J
              !check against available heat in contributing cell
              if(cell_adv > heat_cell(cell_id)) then
                cell_adv = heat_cell(cell_id)
                heat_cell(cell_id) = heat_cell(cell_id) - cell_adv
              endif
            else !heat leaving cell
              cell_adv = Q_cell * gw_rho * gw_cp * gwheat_state(i)%temp !J
              if((cell_adv*-1) > heat_cell(i)) then
                cell_adv = heat_cell(i) * -1 !can only take what is there
              endif
            endif
            heat_adv = heat_adv + cell_adv !J
            heat_cell(i) = heat_cell(i) + cell_adv
            if(gw_state(cell_id)%stat == 2) then !boundary flow
              gw_heat_ss(i)%bndr = gw_heat_ss(i)%bndr + (cell_adv*gw_time_step)
            else
              gw_heat_ss(i)%latl = gw_heat_ss(i)%latl + (cell_adv*gw_time_step)
              gw_heat_ss_yr(i)%latl = gw_heat_ss_yr(i)%latl + (cell_adv*gw_time_step)
            endif
          enddo !go to next connected cell

          !dispersive heat transport (Fourier's Law)
          heat_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
            !temperature gradient (del_T/del_m)
            dist_x = gw_state(i)%xcrd - gw_state(cell_id)%xcrd !m
            dist_y = gw_state(i)%ycrd - gw_state(cell_id)%ycrd !m
            grad_distance = sqrt((dist_x)**2 + (dist_y)**2)
            gradient = (gwheat_state(cell_id)%temp - gwheat_state(i)%temp) / grad_distance
            !interface cross-section area (m2)
            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) !connection length
            face_sat = cell_con(i)%sat(k) !m saturated thickness at interface between cells
            flow_area = face_sat * conn_length
            !heat flux (J/day) (harmonic mean of connected cells)
            face_thmc = conn_length / (((conn_length/2.)/gwheat_state(cell_id)%thmc) + &
                        ((conn_length/2.)/gwheat_state(i)%thmc))
            Q_heat = flow_area * face_thmc * gradient !m2 * J/d-m-T * T/m = J/d
            if(Q_heat > 0) then
              !check against available heat in contributing cell
              if(Q_heat > heat_cell(cell_id)) then
                Q_heat = heat_cell(cell_id)
                heat_cell(cell_id) = heat_cell(cell_id) - Q_heat
              endif
            endif
            if(Q_heat < 0) then
              !check against available heat in current cell
              if((Q_heat*-1) > heat_cell(i)) then
                Q_heat = heat_cell(i) * -1 !can only take what is there
              endif
            endif
            heat_dsp = heat_dsp + Q_heat
            heat_cell(i) = heat_cell(i) + Q_heat
          enddo !go to next connected cell
          gw_heat_ss(i)%disp = heat_dsp * gw_time_step
          gw_heat_ss_yr(i)%disp = gw_heat_ss_yr(i)%disp + (heat_dsp*gw_time_step)

          !change in heat (J); new heat and temperature
          heat_change = (heat_adv + heat_dsp + gw_heat_ss(i)%totl) * gw_time_step
          gwheat_state(i)%stor = gwheat_state(i)%stor + heat_change
          if(gw_state(i)%stor > 0) then
            gwheat_state(i)%tnew = gwheat_state(i)%stor / (gw_rho * gw_cp * gw_state(i)%stor) !deg C
          else
            gwheat_state(i)%tnew = 0.
          endif

        endif !check for interior cell

      enddo !go to the next cell

      !store new temperature values into regular array
      do i=1,ncell
        gwheat_state(i)%told = gwheat_state(i)%temp
        gwheat_state(i)%temp = gwheat_state(i)%tnew
      enddo

      return
      end subroutine gwflow_heat