soil_temperature Subroutine

  subroutine soil_temperature(zz, j, je, k, bcv, ida, mo, pit, preinf, swe, tmax_tmp, tmin_tmp, tx_tmp)
    !**** PURPOSE: THIS SUBROUTINE ESTIMATES DAILY AVERAGE TEMPERATURE
    !              AT THE BOTTOM OF EACH SOIL LAYER
    !              =f(air temp, precip, residue, snow cover)
    !**** CALLED IN: HYDROTOP
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      abd(k) = sum of porosity for soil type, mm
    !      amp(j) = annual 1/2 amplitude in daily average temp, calc readwet,
    !                   degree C
    !      avt(j) = average annual air temperature, calc readwet, degree C
    !      bcv(j, je) = lag factor for residue and snow effects on surf.temp.
    !      ida = current day
    !      isolt = switch code to print from solt()
    !      isosb = number of subbasin to print from solt(), if isolt = 1
    !      iy = current year
    !      mo = current month
    !      nn = number of soil layers, calc in subbasin, cycle 100
    !block pit = 58.13
    !      precip = precipitation, mm, calc in readcli
    !      swe(j, je) = soil water, mm, calc in hydrotop (previous day)
    !      te(j, je, l) = daily ave temp at the bottom of each layer, degree C
    !      te0(j, je) = bare soil surface temp, degree C
    !      tmn(j) = daily min temp. for subbasin, readcli, degree C
    !      tmx(j) = daily max temp. for subbasin, readcli, degree C
    !      tx(j) = daily aver. temperature, readcli, degree C
    !      wft(m, j) = monthly prob. of rainy day, calc readwet
    !      z(l, k) = soil depth, read in readsol, mm
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      alx = local par
    !      b = local par
    !      dd = local par
    !      mdp = max damping depth for the soil
    !      dt = local par
    !      f = local par
    !      l = local par
    !      st0 = bare soil surfac temp
    !      ta = local par
    !      wc = local par
    !      ww = local par
    !      xi = local par
    !      xx = local par
    !      zd = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), dimension(:, :), intent(in) :: bcv
    integer, intent(in) :: ida
    integer, intent(in) :: mo
    real(dp), intent(in) :: pit
    real(dp), dimension(:, :), intent(in) :: preinf
    real(dp), dimension(:, :), intent(in) :: swe
    real(dp), intent(in) :: tmax_tmp, tmin_tmp, tx_tmp

    integer j, je, k, l
    real(dp) zz, alx, b, dd, mdp, dt, f, st0, ta, wc, ww, xi, xx, zd

    xi = ida
    alx = (xi - 200.) / pit
    f = abd(k) / (abd(k) + 686. * exp(- 5.63 * abd(k)))
    mdp = 1000. + 2500. * f
    ww = .356 - .144 * abd(k)
    b =log(500. / mdp)
    wc = swe(j, je) / (ww * z(nn, k))
    f = exp(b * ((1. - wc) / (1. + wc)) ** 2.)
    dd = f * mdp

    !**** CALC st0 - bare soil surfac temp
    if (preinf(j, je) .le. 0.) then
      st0 = wft(mo, j) * (tmax_tmp - tx_tmp) + tx_tmp
    else
      st0 = wft(mo, j) * (tx_tmp - tmin_tmp) + tmin_tmp
    end if

    ta = avt(j) + amp(j) * cos(alx)
    xx = bcv(j, je) * te0(j, je) + (1. - bcv(j, je)) * st0
    dt = xx - ta
    te0(j, je) = st0
    zz = 2. * dd
    xx = 0.

    !**** CALC te()
    do l = 1, nn
      zd = - z(l, k) / dd
      te(j, je, l) = avt(j) + (amp(j) * cos(alx + zd) + dt) * exp(zd)
    end do

    return
  end subroutine soil_temperature