anchor.diff:  
     end function compute_temperature
   
  +  !! Advance mpcd particles
  +  !
  +  ! If the cell system has a wall in the z direction, a bounce-back collision is used.
  +  subroutine mpcd_stream(particles, cells, dt)
  +    type(particle_system_t), intent(inout) :: particles
  +    type(cell_system_t), intent(in) :: cells
  +    double precision, intent(in) :: dt
  +
  +    integer :: i
  +    double precision :: pos_min(3), pos_max(3), delta
  +
  +    pos_min(3) = 0
  +    pos_max(3) = cells% L(3) * cells% a
  +
  +    do i = 1, particles% Nmax
  +       particles% pos(:,i) = particles% pos(:,i) + particles% vel(:,i)*dt
  +       particles% pos(1,i) = modulo( particles% pos(1,i) , cells% L(1)*cells% a )
  +       particles% pos(2,i) = modulo( particles% pos(2,i) , cells% L(2)*cells% a )
  +       if (cells% has_walls) then
  +          if (particles% pos(3,i) < pos_min(3)) then
  +             ! bounce position
  +             delta = pos_min(3) - particles% pos(3,i)
  +             particles% pos(3,i) = pos_min(3) + delta
  +             ! bounce velocity
  +             particles% vel(3,i) = -particles% vel(3,i)
  +          else if (particles% pos(3,i) > pos_max(3)) then
  +             ! bounce position
  +             delta = particles% pos(3,i) - pos_max(3)
  +             particles% pos(3,i) = pos_max(3) - delta
  +             ! bounce velocity
  +             particles% vel(3,i) = -particles% vel(3,i)
  +          end if
  +       else
  +          particles% pos(2,i) = modulo( particles% pos(2,i) , cells% L(2)*cells% a )
  +       end if
  +    end do
  +
  +  end subroutine mpcd_stream
  +
   end module mpcd
  
positive.code: 
    end subroutine simple_mpcd_step
  
    function compute_temperature(particles, cells) result(te)
      class(particle_system_t), intent(in) :: particles
      class(cell_system_t), intent(in) :: cells
  
negative.code:   private
  
    public :: compute_temperature, simple_mpcd_step
    public :: wall_mpcd_step
    public :: mpcd_stream
  
  contains
  
  