ResPart_set_GetHydroForces Subroutine

    impure subroutine ResPart_set_GetHydroForces(this,P,U,V,W,ibVF,visc)
      !> Computes hydrodynamic force on particle.
      implicit none
      class(ResPart_set),   intent(inout) :: this                              !! Collection of Resolved Particles
      type(Eulerian_obj_r), intent(in)    :: P                                 !! Fluid pressure field
      type(Eulerian_obj_r), intent(in)    :: U                                 !! Fluid velocity field in 1-dir
      type(Eulerian_obj_r), intent(in)    :: V                                 !! Fluid velocity field in 2-dir
      type(Eulerian_obj_r), intent(in)    :: W                                 !! Fluid velocity field in 3-dir
      type(Eulerian_obj_r), intent(in)    :: ibVF                              !! Solid volume fraction
      real(wp),             intent(in)    :: visc                              !! Fluid viscosity
      ! work variable
      integer(leapI8):: id
      integer  :: n,m
      integer  :: dir
      real(wp) :: disp(3)

      call this%GetSurfaceStresses(P,U,V,W,ibVF,visc)

      ! Bring markers that may be in other blocks
      ! to the one containing the centroid
      call this%Regroup()
      call this%UpdateLookup()

      ! Zero foces and torques
      select type (particles =>this%p)
      type is (ResPart_obj)
        do n=1,this%count_
          particles(n)%Fh=0.0_wp
          particles(n)%Th=0.0_wp
        end do
      end select

      ! Compute hydrodynamic forces exerted on particles
      select type (particles =>this%p)
      type is (ResPart_obj)
        select type(markers=>this%ib%p)
        type is (marker_obj)
          do m=1,this%ib%count_
            ! Get ID of parent centroid
            id=markers(m)%s
            ! Get local index of the centroid
            n = this%lookup(id)
            ! Update force applied on the resolved particle
            particles(n)%Fh = particles(n)%Fh + markers(m)%SA*markers(m)%f

            ! Compute displacement for torque
            disp = markers(m)%p-particles(n)%p

            ! Treatment for periodicity
            do dir=1,3
              if (this%block%periods(dir)) then
                if ( disp(dir).gt.0 .and. abs(disp(dir)).gt.0.5_wp*particles(n)%d) then
                   disp(dir) = disp(dir) - (this%block%pmax(dir)-this%block%pmin(dir))
                end if

                if ( disp(dir).lt.0 .and. abs(disp(dir)).gt.0.5_wp*particles(n)%d) then
                   disp(dir) = disp(dir) + (this%block%pmax(dir)-this%block%pmin(dir))
               end if
              end if
            end do

            ! Update torque applied on the resolved particle
            particles(n)%Th = particles(n)%Th + markers(m)%SA*cross_product(disp,markers(m)%f)
          end do
        end select
      end select

      return
    end subroutine ResPart_set_GetHydroForces