soil_curve_number_runoff Subroutine

  subroutine soil_curve_number_runoff(j, je, alai, blai, canstor, igro, nucr, preinf, cnmx)
    !**** PURPOSE: THIS SUBROUTINE COMPUTES DAILY RUNOFF GIVEN DAILY PRECIPITATION
    !     AND SNOW MELT USING A MODIFIED SCS CURVE NUMBER APPROACH
    !**** CALLED IN: HYDROTOP
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      alai(j, je) = Leaf Area Index (LAI)
    !      blai(icr) = max LAI for crop
    !      canmax(n) = canopy maximum storage for interception, mm, calc in init
    !      canstor(j, je) = canopy water storage, mm
    !      cn = Curve Number, current
    !      icurn = switch code to print from curn()
    !      icursb = number of subbasin to print from curn(), if icurn = 1
    !      ida = current day
    !      igro(j, je) = vegetation index, =1 if vegetation is growing
    !      nn = number of soil layers, calc in subbasin, cycle 100
    !      nucr(j, je) = crop number (database)
    !      precip = precipitation, mm, read in readcli
    !      preinf(j, je) = precipitation adjusted for canopy storage, mm
    !      qd = daily surface runoff, mm
    !      smx(j, je) = retention coef, calc in curno
    !      ste(j, je, l) = water storage in a layer, mm, calc in hydrotop & purk
    !      te(j, je, l) = soil temperature, degree C, calc in solt
    !      wf(2, j, je) = shape parameters eq.6, calc in curno
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      bb = local par
    !      canmxl = local par
    !      l = local par
    !      pb = local par
    !      r2 = local par
    !      sum = soil water content in all layers
    !      xx = local par
    !      xx1 = local par
    !      xx3 = local par
    !      xx4 = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), dimension(:, :), intent(in) :: alai
    real(dp), dimension(:), intent(in) :: blai
    real(dp), dimension(:, :), intent(inout) :: canstor
    integer, dimension(:, :), intent(in) :: igro
    integer, dimension(:, :), intent(in) :: nucr
    real(dp), dimension(:, :), intent(inout) :: preinf
    real(dp), intent(in) :: cnmx

    integer j, je, l
    real(dp) bb, canmxl, pb, r2, sum, xx, xx1, xx3, xx4

    sum = 0.
    do l = 1, nn
      sum = sum + ste(j, je, l)
    end do

      !**** Canopy interception
    xx1 = 0.
    canmxl = 0.
    xx1 = preinf(j, je)

      !**** CALC canopy storage
    if (igro(j, je) .ge. 1) then
      canmxl = cnmx * alai(j, je) / blai(nucr(j, je))
    else
      canmxl = 0.
    endif

    xx3 = preinf(j, je) - canmxl
    if (xx3 < 0.) then
      canstor(j, je) = xx1
    else
      canstor(j, je) = canmxl
    endif

    preinf(j, je) = preinf(j, je) - canstor(j, je)
    xx4 = preinf(j, je) - canstor(j, je)

    xx = wf(1, j, je) - wf(2, j, je) * sum
    if (xx .lt. - 20.) xx = - 20.
    if (xx .gt. 20.) xx = 20.

    r2 = smx(j, je) * (1. - sum / (sum + exp(xx)))
    if (te(j, je, 2) .le. 0.) r2 = smx(j, je) * (1. - exp(- .000862 * r2))
    cn = 25400. / (r2 + 254.)
    r2 = 25400. / cn - 254.
    bb = .2 * r2
    pb = xx4 - bb

    !**** CALC daily surface runoff qd
    if (pb .gt. 0.) then
      qd = pb * pb / (xx4 + .8 * r2)
    else
      qd = 0.
    end if

    return
  end subroutine soil_curve_number_runoff