cdifs_AdvanceSolutionMomentumRHS1 Function

    function cdifs_AdvanceSolutionMomentumRHS1(this) result(res)
      !> Computes Momentum operator RHS for velocity component in x1-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=1)

      ! 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,op%st-1
                ! Shift
                is = i+n; js = j; ks = k

                ! Convective term: conv11
                Vinterp(1) = dot_product(op%c_intrp1 (:,is),    V(1)%cell(is-op%st+1:is+op%st,js,ks))
                CONV = 0.0_wp
                do kk = 1, op%st
                  Flux_c  = Vinterp(1) * dot_product(op%morinishi_int(:,kk),V(1)%cell(is-op%st+1:is+op%st,js,ks))
                  CONV = CONV + op%alpha(kk)*op%morinishi_ddx(n+op%st+1,kk)/(xm(i) - xm(i-1)) * 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 term: conv21
                Vinterp(2) = dot_product(op%c_intrp1m(:,is),    V(2)%cell(is-op%st:is+op%st-1,js,ks))
                CONV = 0.0_wp
                do kk = 1, op%st
                  Flux_c  = Vinterp(2) * dot_product(op%morinishi_int(:,kk),V(1)%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+1,op%st
                ! Shift
                is = i; js = j; ks = k+n

                ! Convective: conv31
                Vinterp(3) = dot_product(op%c_intrp1m(:,is),    V(3)%cell(is-op%st:is+op%st-1,js,ks))
                CONV = 0.0_wp
                do kk = 1, op%st
                  Flux_c  = Vinterp(3) * dot_product(op%morinishi_int(:,kk),V(1)%cell(is,js,ks-op%st:ks+op%st-1))
                  CONV = CONV + op%alpha(kk)*op%morinishi_ddx(n+op%st,kk)/(z (k+1)-z (k  )) * 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,1),P%cell(i-op%st:i+op%st-1,j,k))

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

      return
    end function cdifs_AdvanceSolutionMomentumRHS1