cdifs_obj_AdvanceSolutionCorrector Subroutine

    subroutine cdifs_obj_AdvanceSolutionCorrector(this,it)
      !> Perform corrector step: compute pressure at n+1 and divergence-free
      ! velocity at n+1.
      implicit none
      class(cdifs_obj), intent(inout) :: this                                  !! CDIFS solver
      integer,          intent(in)    :: it                                    !! Subiterations
      ! Work variables
      integer :: i,j,k
      real(wp):: buffr,vol,dV
      real(wp):: intRHS
      integer :: dir

      ! Compute the right-hand side of the pressure poisson equation
      associate (lo => this%block%lo, hi => this%block%hi, op => this%op, &
        rho => this%rho, dt => this%timer%dt)

          this%rhs = 0.0_wp
          do k=lo(3),hi(3)
            do j=lo(2),hi(2)
              do i=lo(1),hi(1)

                this%rhs%cell(i,j,k) = this%rhs%cell(i,j,k) + (rho/dt)*( &
                      dot_product(this%DIV(:,i,j,k,1),this%V(1)%cell(i-op%st+1:i+op%st,j,k)) &
                    + dot_product(this%DIV(:,i,j,k,2),this%V(2)%cell(i,j-op%st+1:j+op%st,k)) &
                    + dot_product(this%DIV(:,i,j,k,3),this%V(3)%cell(i,j,k-op%st+1:k+op%st)) )

              end do
            end do
          end do

      end associate

      ! Compute integral of RHS -- needs to be close to machine precision
      associate ( lo => this%block%lo, hi => this%block%hi )

        intRHS = 0.0_wp
        vol    = 0.0_wp
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              dV     = this%block%axis(1)%dxm(i) &
                     * this%block%axis(2)%dxm(j) &
                     * this%block%axis(3)%dxm(k)
              vol    = vol + dV
              intRHS = intRHS + this%rhs%cell(i,j,k)*dV
            end do
          end do
        end do

        call this%parallel%sum(vol,   buffr); vol    = buffr
        call this%parallel%sum(intRHS,buffr); intRHS = buffr/vol
      end associate

      ! Compute the pressure correction
      !--------------------------------
      associate ( dP=>this%DP, rhs=>this%rhs, hypre=>this%hypre, bcs=>this%bcs)
        call hypre%SetRHS(rhs)
        call hypre%Solve(dP)
        dP%cell = dP%cell - dP%mean()

        call bcs%UpdateBoundary(dP)
      end associate


      ! Update pressure
      !--------------------------------
      associate ( dP=>this%DP, P=>this%P, bcs=>this%bcs)
        P = P + dP
        call bcs%UpdateBoundary(P)
      end associate

      associate ( V=>this%V, dP=>this%dP, rho=>this%rho,    &
        lo=>this%block%lo, hi=>this%block%hi, op=>this%op,  &
        bcs=>this%bcs, dt=>this%timer%dt)

        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              V(1)%cell(i,j,k) = V(1)%cell(i,j,k) - (dt/rho) * dot_product(this%pGRAD(:,i,j,k,1),dP%cell(i-op%st:i+op%st-1,j,k))
              V(2)%cell(i,j,k) = V(2)%cell(i,j,k) - (dt/rho) * dot_product(this%pGRAD(:,i,j,k,2),dP%cell(i,j-op%st:j+op%st-1,k))
              V(3)%cell(i,j,k) = V(3)%cell(i,j,k) - (dt/rho) * dot_product(this%pGRAD(:,i,j,k,3),dP%cell(i,j,k-op%st:k+op%st-1))
            end do
          end do
        end do

        do dir=1,3
          call V(dir)%UpdateGhostCells
          call bcs%UpdateBoundary(V(dir))
        end do
      end associate

      ! Send diagnostics to pressure monitor
      call this%pmonitor%Set('Pressure', 1, this%timer%iter)
      call this%pmonitor%Set('Pressure', 2, this%timer%time)
      call this%pmonitor%Set('Pressure', 3, it             )
      call this%pmonitor%Set('Pressure', 4, this%hypre%rel )
      call this%pmonitor%Set('Pressure', 5, intRHS         )
      call this%pmonitor%Set('Pressure', 6, this%hypre%it  )

      ! Print pressure monitor
      call this%pmonitor%print()

      return
    end subroutine cdifs_obj_AdvanceSolutionCorrector