vegetation_forcing Subroutine

  subroutine vegetation_forcing
    use modglobal,  only : ib,ie,jb,je,kb,ke,lmoist,ltempeq,nsv,itree_mode,TREE_MODE_DRAG_ONLY,TREE_MODE_SVEG,TREE_MODE_LEGACY_SEB
    use modfields,  only : um,vm,wm,up,vp,wp,svp,svm
    implicit none
    integer :: i,j,k,m,n,npts
    integer :: mf
    real :: dcoefv, ladv, udv

    if (.not. vegetation_ready) return

    npts = veg%npts

    if (npts <= 0) return

    call reset_vegetation_sources()

    ! ========================================================================
    ! Loop 1: Momentum drag forces (precomputed staggered faces, no branching)
    ! ========================================================================
    do mf = 1, npts_u
      i = ijk_u(mf,1)
      j = ijk_u(mf,2)
      k = ijk_u(mf,3)
      dcoefv = dcoef_u(mf)
      veg_up(mf) = - dcoefv * um(i,j,k) * &
                         sqrt( um(i,j,k)**2 &
                         + (0.25*(vm(i,j,k) + vm(i,j+1,k) + vm(i-1,j,k) + vm(i-1,j+1,k)))**2 &
                         + (0.25*(wm(i,j,k) + wm(i,j,k+1) + wm(i-1,j,k) + wm(i-1,j,k+1)))**2 )
      up(i,j,k) = up(i,j,k) + veg_up(mf)
    end do

    do mf = 1, npts_v
      i = ijk_v(mf,1)
      j = ijk_v(mf,2)
      k = ijk_v(mf,3)
      dcoefv = dcoef_v(mf)
      veg_vp(mf) = - dcoefv * vm(i,j,k) * &
                         sqrt( vm(i,j,k)**2 &
                         + (0.25*(um(i,j,k) + um(i+1,j,k) + um(i,j-1,k) + um(i+1,j-1,k)))**2 &
                         + (0.25*(wm(i,j,k) + wm(i,j,k+1) + wm(i,j-1,k) + wm(i,j-1,k+1)))**2 )
      vp(i,j,k) = vp(i,j,k) + veg_vp(mf)
    end do

    do mf = 1, npts_w
      i = ijk_w(mf,1)
      j = ijk_w(mf,2)
      k = ijk_w(mf,3)
      dcoefv = dcoef_w(mf)
      veg_wp(mf) = - dcoefv * wm(i,j,k) * &
                         sqrt( wm(i,j,k)**2 &
                         + (0.25*(um(i,j,k) + um(i+1,j,k) + um(i,j,k-1) + um(i+1,j,k-1)))**2 &
                         + (0.25*(vm(i,j,k) + vm(i,j+1,k) + vm(i,j,k-1) + vm(i,j+1,k-1)))**2 )
      wp(i,j,k) = wp(i,j,k) + veg_wp(mf)
    end do

    ! ========================================================================
    ! Loop 2: Canopy energy balance (heat and moisture fluxes)
    ! ========================================================================
    if (lmoist .and. ltempeq) then
      if (itree_mode == TREE_MODE_LEGACY_SEB) then
        call vegetation_forcing_legacy
      else if (itree_mode == TREE_MODE_SVEG) then
        call vegetation_forcing_sveg
      end if
    end if

    ! ========================================================================
    ! Loop 3: Scalar deposition
    ! ========================================================================
    if (nsv > 0) then
      do m = 1, npts
        i = veg%ijk(m,1)
        j = veg%ijk(m,2)
        k = veg%ijk(m,3)

        ladv = veg%lad(m)
        udv  = veg%ud(m)

        do n = 1, nsv
          vegp%sv(m,n) = vegp%sv(m,n) - svm(i,j,k,n) * ladv * udv
          svp(i,j,k,n) = svp(i,j,k,n) + vegp%sv(m,n)
        end do
      end do
    end if

  end subroutine vegetation_forcing