cdifs_AdvanceSolutionMomentumRHS3 Function

    function cdifs_AdvanceSolutionMomentumRHS3(this) result(res)
      !> Computes Momentum operator RHS for velocity component in x3-direction.
      implicit none
      class(cdifs_obj),     intent(inout):: this                               !! CDIFS solver
      type(eulerian_obj_r)               :: res                                !! result
      ! Work variables
      integer :: i,j,k
      integer :: n,kk
      integer :: is, js, ks
      real(wp):: Flux_c
      real(wp):: CONV
      real(wp):: Vinterp(3)
      integer :: stm

      call res%initialize('rhs',this%block,this%parallel,stag=3)

      ! Build the right-hand side of the momentum equation with delta formulation
      ! and store it in RES
      associate ( &
        lo => this%block%lo, hi => this%block%hi, x => this%block%x, &
        xm => this%block%xm, y => this%block%y, ym => this%block%ym, &
        z  => this%block%z,  zm => this%block%zm, rho => this%rho,   &
        visc => this%visc, V => this%V, P => this%P, op => this%op )

        stm = 2*op%st-1
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)

              res%cell(i,j,k) = 0.0_wp

              ! Convective part
              ! -- Flux in 1-dir
              do n = -op%st+1,op%st
                ! Shift
                is = i+n; js = j; ks = k

                ! Convective: conv13
                Vinterp(1) = dot_product(op%c_intrp3m(:,ks),    V(1)%cell(is,js,ks-op%st:ks+op%st-1))
                CONV = 0.0_wp
                do kk = 1, op%st
                  Flux_c  = Vinterp(1) * dot_product(op%morinishi_int(:,kk),V(3)%cell(is-op%st:is+op%st-1,js,ks))
                  CONV = CONV + op%alpha(kk)*op%morinishi_ddx(n+op%st,kk)/(x(i+1) - x(i)) * Flux_c
                end do

                ! Add convective, viscous, and pressure terms
                res%cell(i,j,k) = res%cell(i,j,k) - CONV
              end do
              ! -- Flux in 2-dir
              do n= -op%st+1,op%st
                ! Shift
                is = i; js = j+n; ks = k

                ! Convective: conv23
                Vinterp(2) = dot_product(op%c_intrp3m(:,ks),    V(2)%cell(is,js,ks-op%st:ks+op%st-1))
                CONV = 0.0_wp
                do kk = 1, op%st
                  Flux_c  = Vinterp(2) * dot_product(op%morinishi_int(:,kk),V(3)%cell(is,js-op%st:js+op%st-1,ks))
                  CONV = CONV + op%alpha(kk)*op%morinishi_ddx(n+op%st,kk)/(y (j+1)-y (j  )) * Flux_c
                end do

                ! Add convective, viscous, and pressure terms
                res%cell(i,j,k) = res%cell(i,j,k) - CONV
              end do
              ! -- Flux in 3-dir
              do n= -op%st,op%st-1
                ! Shift
                is = i; js = j; ks = k+n

                ! Convective: conv33
                Vinterp(3) = dot_product(op%c_intrp3 (:,ks),V(3)%cell(is,js,ks-op%st+1:ks+op%st))
                CONV = 0.0_wp
                do kk = 1, op%st
                  Flux_c  = Vinterp(3) * dot_product(op%morinishi_int(:,kk),V(3)%cell(is,js,ks-op%st+1:ks+op%st))
                  CONV = CONV + op%alpha(kk)*op%morinishi_ddx(n+op%st+1,kk)/(zm(k)-zm(k-1)) * Flux_c
                end do

                ! Add convective, viscous, and pressure terms
                res%cell(i,j,k) = res%cell(i,j,k) - CONV
              end do

              ! Include pressure part
              res%cell(i,j,k) = res%cell(i,j,k) -(1.0_wp/rho)*                               &
                                  dot_product(this%pGRAD(:,i,j,k,3),P%cell(i,j,k-op%st:k+op%st-1))
              ! Include viscous part
              res%cell(i,j,k) = res%cell(i,j,k) + visc/rho*(                                 &
                                  dot_product(this%LAP(:,i,j,k,1,3),V(3)%cell(i-stm:i+stm,j,k)) &
                                + dot_product(this%LAP(:,i,j,k,2,3),V(3)%cell(i,j-stm:j+stm,k)) &
                                + dot_product(this%LAP(:,i,j,k,3,3),V(3)%cell(i,j,k-stm:k+stm)) )
            end do
          end do
        end do
      end associate

      return
    end function cdifs_AdvanceSolutionMomentumRHS3