compute_sgs_fluxes Subroutine

  subroutine compute_sgs_fluxes(usgs,vsgs,wsgs,thlsgs,qtsgs,sv1sgs,sv2sgs,sv3sgs,sv4sgs)
    use modglobal,  only : ib,ie,ih,jb,je,jh,kb,ke,kh,dxf,dxhi,dyi,dzf,dzfi,dzhi,dzhiq,dzh2i,&
                           ltempeq,lmoist,nsv
    use modfields,  only : um,vm,wm,thlm,qtm,svm
    use modsubgrid, only : ekh,ekm
    implicit none

    real, intent(out) :: usgs  (ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh)
    real, intent(out) :: vsgs  (ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh)
    real, intent(out) :: wsgs  (ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh)
    real, intent(out) :: thlsgs(ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh)
    real, intent(out) :: qtsgs (ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh)
    real, intent(out) :: sv1sgs(ib:ie,jb:je,kb:ke+kh)
    real, intent(out) :: sv2sgs(ib:ie,jb:je,kb:ke+kh)
    real, intent(out) :: sv3sgs(ib:ie,jb:je,kb:ke+kh)
    real, intent(out) :: sv4sgs(ib:ie,jb:je,kb:ke+kh)

    real    :: emom
    integer :: i,j,k

    usgs = 0.; vsgs = 0.; wsgs = 0.; thlsgs = 0.; qtsgs = 0.
    sv1sgs = 0.; sv2sgs = 0.; sv3sgs = 0.; sv4sgs = 0.

    !> Momentum fluxes at the cell corners (uw and vw)
    do k=kb,ke+kh
      do j=jb,je
        do i=ib,ie

          ! interps ekm to cell corner (uw)
          emom = ( dzf(k-1) * ( ekm(i,j,k)*dxf(i-1)  + ekm(i-1,j,k)*dxf(i) )  + &
                   dzf(k)   * ( ekm(i,j,k-1)*dxf(i-1) + ekm(i-1,j,k-1)*dxf(i) ) )*dxhi(i) * dzhiq(k)
          usgs(i,j,k)  = -emom * ( (um(i,j,k)-um(i,j,k-1)) *dzhi(k) &
                      +(wm(i,j,k)-wm(i-1,j,k))  *dxhi(i))

          ! interps ekm to cell corner (vw)
          emom = ( dzf(k-1) * ( ekm(i,j,k)  + ekm(i,j-1,k) )  + &
                   dzf(k)   * ( ekm(i,j,k-1) + ekm(i,j-1,k-1) ) ) * dzhiq(k)

          vsgs(i,j,k)  = -emom * ( (vm(i,j,k)-vm(i,j,k-1)) *dzhi(k) &
                      +(wm(i,j,k)-wm(i,j-1,k))  *dyi)

        end do
      end do
    end do

    !> Normal stress tau_33 = -2*K_m*dw/dz, which lives at the cell centre
    do k=kb,ke
      do j=jb,je
        do i=ib,ie
          wsgs(i,j,k) = -2. * ekm(i,j,k) * (wm(i,j,k+1)-wm(i,j,k)) * dzfi(k)
        end do
      end do
    end do

    !> Scalar fluxes at the w-levels
    if (ltempeq) then
      do k=kb,ke+kh
        thlsgs(:,:,k) = -0.5 * (dzf(k-1)*ekh(:,:,k) + dzf(k)*ekh(:,:,k-1)) &
                        * (thlm(:,:,k)-thlm(:,:,k-1)) * dzh2i(k)
      end do
    end if

    if (lmoist) then
      do k=kb,ke+kh
        qtsgs(:,:,k) = -0.5 * (dzf(k-1)*ekh(:,:,k) + dzf(k)*ekh(:,:,k-1)) &
                        * (qtm(:,:,k)-qtm(:,:,k-1)) * dzh2i(k)
      end do
    end if

    if (nsv>0) then
      do k=kb,ke+kh
        sv1sgs(ib:ie,jb:je,k) = -0.5 * (dzf(k-1)*ekh(ib:ie,jb:je,k) + dzf(k)*ekh(ib:ie,jb:je,k-1)) &
                        * (svm(ib:ie,jb:je,k,1)-svm(ib:ie,jb:je,k-1,1)) * dzh2i(k)
      end do
    end if

    if (nsv>1) then
      do k=kb,ke+kh
        sv2sgs(ib:ie,jb:je,k) = -0.5 * (dzf(k-1)*ekh(ib:ie,jb:je,k) + dzf(k)*ekh(ib:ie,jb:je,k-1)) &
                        * (svm(ib:ie,jb:je,k,2)-svm(ib:ie,jb:je,k-1,2)) * dzh2i(k)
      end do
    end if

    if (nsv>2) then
      do k=kb,ke+kh
        sv3sgs(ib:ie,jb:je,k) = -0.5 * (dzf(k-1)*ekh(ib:ie,jb:je,k) + dzf(k)*ekh(ib:ie,jb:je,k-1)) &
                        * (svm(ib:ie,jb:je,k,3)-svm(ib:ie,jb:je,k-1,3)) * dzh2i(k)
      end do
    end if

    if (nsv>3) then
      do k=kb,ke+kh
        sv4sgs(ib:ie,jb:je,k) = -0.5 * (dzf(k-1)*ekh(ib:ie,jb:je,k) + dzf(k)*ekh(ib:ie,jb:je,k-1)) &
                        * (svm(ib:ie,jb:je,k,4)-svm(ib:ie,jb:je,k-1,4)) * dzh2i(k)
      end do
    end if

  end subroutine compute_sgs_fluxes