cdifs_obj_AdvanceSolutionIB Subroutine

    subroutine cdifs_obj_AdvanceSolutionIB(this)
      !> Performs intermediate IB step>
      implicit none
      class(cdifs_obj), intent(inout) :: this                                  !! CDIFS solver
      ! Work variable
      type(eulerian_obj_r):: tmp(3)
      integer             :: i,j,k

      ! Initialize temporary fields
      call tmp(1)%Initialize('V1',this%block,this%parallel,0)
      call tmp(2)%Initialize('V2',this%block,this%parallel,0)
      call tmp(3)%Initialize('V3',this%block,this%parallel,0)

      ! Build forcing
      associate ( V => this%V, rho => this%rho, ibS => this%ibS, &
          ibF => this%ibF, op => this%op, dt => this%timer%dt )

        ! Update SDF
        ibS = 0.0_wp
        if (this%use_RP) call this%RP%UpdateSDF(ibS)
        if (this%use_IB) call this%IB%UpdateSDF(ibS)
        call this%bcs%UpdateBoundary(this%ibS)

        ! Interpolate velocities
        tmp(1) = op%intrp1(V(1))
        tmp(2) = op%intrp2(V(2))
        tmp(3) = op%intrp3(V(3))
        call this%bcs%UpdateBoundary(tmp)

        ! Recompute the forcing
        ibF(1) = 0.0_wp; ibF(2) = 0.0_wp; ibF(3) = 0.0_wp
        if (this%use_RP) call this%RP%GetIBForcing(tmp(1),tmp(2),tmp(3),rho,ibS,ibF,dt)
        if (this%use_IB) call this%IB%GetIBForcing(tmp(1),tmp(2),tmp(3),rho,ibS,ibF,dt)
        call this%bcs%UpdateBoundary(this%ibF)

      end associate

      call tmp(1)%Finalize()
      call tmp(2)%Finalize()
      call tmp(3)%Finalize()

      ! Apply forcing
      associate ( V => this%V, ibF => this%ibF, 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)
              if (this%maskV(1)%cell(i,j,k)>0) &
                V(1)%cell(i,j,k) = V(1)%cell(i,j,k) + dt/rho * dot_product(op%c_intrp1m(:,i),ibF(1)%cell(i-op%st:i+op%st-1,j,k))
              if (this%maskV(2)%cell(i,j,k)>0) &
                V(2)%cell(i,j,k) = V(2)%cell(i,j,k) + dt/rho * dot_product(op%c_intrp2m(:,j),ibF(2)%cell(i,j-op%st:j+op%st-1,k))
              if (this%maskV(3)%cell(i,j,k)>0) &
                V(3)%cell(i,j,k) = V(3)%cell(i,j,k) + dt/rho * dot_product(op%c_intrp3m(:,k),ibF(3)%cell(i,j,k-op%st:k+op%st-1))
            end do
          end do
        end do

        call bcs%UpdateBoundary(V)
      end associate

      return
    end subroutine cdifs_obj_AdvanceSolutionIB