anchor.diff:      gwmag = sxGW*sxGW + syGW*syGW + szGW+szGW
       gwmag = sqrt(gwmag)
   
  -    fx(i) = fx(i) - STC*mass/rho(i)*LW*sxGW/gwmag
  -    fy(i) = fy(i) - STC*mass/rho(i)*LW*syGW/gwmag
  -    fz(i) = fz(i) - STC*mass/rho(i)*LW*szGW/gwmag
  +    fx(i) = fx(i) - STC*mass*LW*sxGW/gwmag
  +    fy(i) = fy(i) - STC*mass*LW*syGW/gwmag
  +    fz(i) = fz(i) - STC*mass*LW*szGW/gwmag
   
     end do
   
  
positive.code:   use precision
    implicit none
    real(WP),dimension(:),allocatable :: fx, fy, fz  ! Total force
  end module forces
  
  ! ===================================== !
  
negative.code:   real(WP),dimension(:),allocatable ::  vx,  vy,  vz  ! Velocity
    real(WP),dimension(:),allocatable ::  vlx,  vly,  vlz  ! Velocity half time step behind
    real(WP),dimension(:),allocatable ::  rho           ! Density
    !dir$ attributes align:64 :: px, py, pz
    !dir$ attributes align:64 :: vx, vy, vz
    !dir$ attributes align:64 :: vlxx, vly, vlz
    !dir$ attributes align:64 :: rho
    !DIR$ ASSUME_ALIGNED px: 64
    !DIR$ ASSUME_ALIGNED py: 64
    !DIR$ ASSUME_ALIGNED pz: 64
    !DIR$ ASSUME_ALIGNED vx: 64
    !DIR$ ASSUME_ALIGNED vy: 64
    !DIR$ ASSUME_ALIGNED vz: 64
    !DIR$ ASSUME_ALIGNED vlx: 64
    !DIR$ ASSUME_ALIGNED vly: 64
    !DIR$ ASSUME_ALIGNED vlz: 64
    !DIR$ ASSUME_ALIGNED rho: 64
  
  end module
  
  ! ===================================== !
  