vegetation_forcing_sveg Subroutine

  subroutine vegetation_forcing_sveg
    use modglobal, only : pref0, rlv, cp, rv, rd, rhoa
    use modfields, only : thlm, qtm, qtp, thlp, um, vm, wm
    implicit none
    integer :: i, j, k, m
    real :: ladv, q_av_leaf
    real :: e_sat, e_vap, d_vap, slope_sat, r_a, omega, qe, qh, gam
    real :: lsizev, rsv, wind2

    gam = (cp*pref0*rv)/(rlv*rd)
    do m = 1, veg%npts
      i = veg%ijk(m,1)
      j = veg%ijk(m,2)
      k = veg%ijk(m,3)

      ladv = veg%lad(m)
      q_av_leaf = sveg(m) / max(ladv, 1.0e-12)

      lsizev = max(veg%lsize(m), 1.0e-6)
      rsv = max(veg%rs(m), 1.0e-6)

      e_sat = 610.8*exp((17.27*(thlm(i,j,k)-273.15))/(thlm(i,j,k)-35.85))
      e_vap = (qtm(i,j,k) * pref0) / (0.378 * qtm(i,j,k) + 0.622)
      d_vap = max(e_sat - e_vap, 0.)
      slope_sat = (4098*e_sat)/((thlm(i,j,k)-35.85)**2)

      wind2 = max((0.5*(um(i,j,k)+um(i+1,j,k)))**2 &
                +(0.5*(vm(i,j,k)+vm(i,j+1,k)))**2 &
                +(0.5*(wm(i,j,k)+wm(i,j,k+1)))**2, 1.0e-12)
      r_a = 130*sqrt(lsizev / sqrt(wind2))

      omega = 1/(1 + 2*(gam/(slope_sat+2*gam)) * (rsv/r_a))
      qe = omega*(slope_sat/(slope_sat+2*gam))*q_av_leaf + (1-omega)*(1/(gam*rsv))*rhoa*cp*d_vap
      qh = q_av_leaf - qe

      vegp%omega(m) = omega
      vegp%qtR(m) = ladv*(omega*(slope_sat/(slope_sat+2*gam))*q_av_leaf)/(rhoa*rlv)
      vegp%qtA(m) = ladv*((1-omega)*(1/(gam*rsv))*rhoa*cp*d_vap)/(rhoa*rlv)
      vegp%qt(m) = vegp%qt(m) + ladv*qe/(rhoa*rlv)
      vegp%thl(m) = vegp%thl(m) + ladv*qh/(rhoa*cp)
      qtp(i,j,k) = qtp(i,j,k) + vegp%qt(m)
      thlp(i,j,k) = thlp(i,j,k) + vegp%thl(m)
    end do
  end subroutine vegetation_forcing_sveg