soil_process Subroutine

  subroutine soil_process(j, je, k, rain, swe)
    !**** PURPOSE: THIS SUBROUTINE DIVIDES EACH LAYER'S FLOW INTO 4 MM SLUGS
    !              and CALLS PERC to CALC PERCOLATION and LATERAL SUBSURFACE FLOW
    !**** CALLED IN: HYDROTOP
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      fc(l, k) = field capacity, mm, calc readsol
    !      flate(j, je, l) = lateral subsurface flow, sum of prk for layer, mm
    !      ida = current day
    !      ipehd = number of hydrotope to print from perc(), if iperc = 1
    !      iperc = switch code to print from perc()
    !      ipesb = number of subbasin to print from perc(), if iperc = 1
    !      nn = number of soil layers, calc in subbasin
    !      poe(j, je, l) = percolation, sum of sep for layer, mm
    !      prk = lateral subsurface flow, calc in perc for 4mm slugs, mm
    !      rain = preinf(j, je)-qd, mm
    !      sep = percolation, calc in perc for 4mm slugs, recalc here, mm
    !      ste(j, je, l) = water storage in a layer, recalc here, mm
    !      su = water excess above FC, mm
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      add = local par
    !      amt = local par
    !      j1 = local par
    !      j2 = local par
    !      l = local par
    !      n1 = local par
    !      n2 = local par
    !      n3 = local par
    !      sum = local par
    !      tot = local par
    !      vvv = local par
    !      vvvs = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !**** Include common parameters

    real(dp), intent(in) :: rain
    real(dp), dimension(:, :), intent(in) :: swe
    integer j, je, k, j1, j2, l, n1, n2, n3
    real(dp) river_route_add, sum, tot, vvv
    dimension vvv(11)

    sum = 0
    n1 = 1
    n3 = 1

    !**** Calc vvv()
    n2 = 0
    prk = 0.
    do l = 1, nn
      flate(j, je, l) = 0.
      poe(j, je, l) = 0.
    end do
    vvv(1) = ste(j, je, 1) - fc(1, k)
    if (vvv(1) .ge. 0.) then
      ste(j, je, 1) = fc(1, k)
    else
      vvv(1) = 0.
    end if
    vvv(1) = vvv(1) + rain
    if (vvv(1) .gt. 0.) then
      n1 = 1
      n2 = 1
    end if
    do l = 2, nn
      vvv(l) = ste(j, je, l) - fc(l, k)
      if (vvv(l) .gt. 0.) then
        ste(j, je, l) = fc(l, k)
        if (n2 .gt. 0) CYCLE
        n2 = 1
        n1 = l
      else
        vvv(l) = 0.
      end if
    end do
    vvv(nn + 1) = 0.
    tot = 0.
    river_route_add = 0.

    !#### CALL PERC to calculate percolation and subsurface flow
    do while (n2 .ne. 0)
      n2 = 0
      do j1 = n1, nn
        sum = 0.
        j2 = j1 + 1
        if (j1 .eq. nn) j2 = nn
        if (vvv(j1) .gt. 0.) then
          if (n2 .le. 0) then
              n2 = 1
              n3 = j1
          end if
          su = vvv(j1)
          if (su .gt. amt) su = amt
          ste(j, je, j1) = ste(j, je, j1) + su
          sep = 0.
          ste(j, je, j1) = ste(j, je, j1) - sep
          sum = sep
          if (su .ne. sep) then
            if (ste(j, je, j1) .gt. fc(j1, k)) then
              call soil_percolation(j, je, k, j1, j2, swe)
              ste(j, je, j1) = ste(j, je, j1) - sep - prk
              sum = sum + sep
              river_route_add = river_route_add + prk
              flate(j, je, j1) = flate(j, je, j1) + prk
              poe(j, je, j1) = poe(j, je, j1) + sep
            endif
          end if
          vvv(j1) = vvv(j1) - su
          j2 = j1 + 1
          vvv(j2) = vvv(j2) + sum
        end if
      end do
      tot = tot + sum
      n1 = n3
    end do

    sep = tot
    prk = river_route_add

    return
  end subroutine soil_process