anchor.diff:    end if
   
     call enforce_BCs
  +
  +  fx = 0.0_WP
  +  fy = 0.0_WP
  +  fz = 0.0_WP
   
     return
   end subroutine posvel_update
  
positive.code: ! One step of SPH                       !
  ! ===================================== !
  subroutine step
    implicit none
  
    call neighbor_find
    call compute_density
    call compute_pressure
    call compute_visc
    call ext_forces
    call posvel_update
  
  end subroutine step
  
  
negative.code:       dz = pz(i) - pz(j)
        r2 = dx*dx + dy*dy + dz*dz
        h2 = h*h
  
        if(r2 .gt. h2) then
          visckernel = 0.0_WP
        else
          c = 45.0_WP*mass/PI/h**5/rho(j)*(1.0_WP-sqrt(r2/h2))
          visckernel = -c * VISCOSITY
        end if
  
        return
  