soil_percolation Subroutine

  subroutine soil_percolation(j, je, k, j1, j2, swe)
    !**** PURPOSE: THIS SUBROUTINE COMPUTES PERCOLATION AND LATERAL SUBSURFACE
    !              FLOW FROM A SOIL LAYER WHEN FIELD CAPACITY IS EXCEEDED.
    !**** METHOD: SWRRB
    !**** CALLED IN: PURK
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !inp j1 - layer, from purk
    !inp j2 - layer, from purk
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      hk(l, k) = beta coef. to estimate hydr. cond., used in perc
    !      hrtc(j, je, l) = return flow travel time, calc. hrflowtt, h
    !      hwss(2, j, je) = fun(field cap), calc in subbasin from hsumfc(j, jea)
    !      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
    !      prk = lateral subsurface flow, calc in perc for 4mm slugs, mm
    !      sc(l, k) = saturated conductivity, mm/h, calc, if isc = 1
    !      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
    !      swe(j, je) = soil water content, mm
    !      te(j, je, l) = daily average temp at the bottom of each layer, degree C
    !      ul(l, k) = upper limit water content in layer, mm
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      adjf = local par
    !      cr = local par
    !      fx = local par
    !      iconfig = local par
    !      rtf = local par
    !      rtw = local par
    !      stu = local par
    !      stu1 = local par
    !      stu2 = local par
    !      stz = local par
    !      sup = local par
    !      sup1 = local par
    !      xx = local par
    !      zz = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !**** Include common parameters

    real(dp), dimension(:, :), intent(in) :: swe

    integer j, je, k, j1, j2
    real(dp) cr, fx, rtf, rtw, stu, stu1, stu2, stz, sup, sup1, xx, zz

    !**** IF TEMP OF LAYER IS BELOW 0 C - NO WATER FLOW
    if (te(j, je, j1) .le. 0.) then
      sep = 0.
      prk = 0.
      rtf = 0.
      sup1 = 0.
      rtw = 0.
    else
      sup = su - sep
      sup1 = sup

      !**** CALC SEEPAGE TO NEXT LAYER
      xx = hwss(1, j, je) - hwss(2, j, je) * swe(j, je)
      if (xx .lt. - 20.) xx = - 20.
      rtw = 10. * (1. - swe(j, je) / (swe(j, je) + exp(xx)))
      rtf = hrtc(j, je, j1) * rtw
      if (rtf .gt. 0.) rtf = 1. - exp(- 1. / rtf)
      stz = ste(j, je, j1) / ul(j1, k)
      if (stz .ge. 1.) then
        fx = 1.
      else
        fx = stz ** hk(j1, k)
      end if
      stu = ste(j, je, j2) / ul(j2, k)
      stu1 = ul(j1, k) / ste(j, je, j1)
      stu2 = ste(j, je, j1) / ul(j1, k)

      if (stu .ge. 1.) then
        sep = 0.
        prk = rtf * sup
      else
        cr = sqrt(1. - stu)
        ! SL ***
        zz = 24. * fx * cr * sc(j1, k) * sccor(j) / sup
        if (zz .gt. 10.) then
          sep = sup
          prk = 0.
        else
          sep = sup * (1. - exp(- zz))
          sup = sup - sep
          prk = rtf * sup
        endif
      end if
    end if

    !**** Write
    if (j .eq. 1 .and. je .eq. 1) then
      if (sep <= 0.) sep = .000001
    end if

    return
  end subroutine soil_percolation