anchor.diff:  c... current x
               x_tmp(:) = x(i,:) + disp_snap(i,:)
   c... correct x
  -            call sim_closest_pos (disp_snap(i,1), disp_snap(i,2),
  -     &                            disp_snap(i,3), i)
  +c            call sim_closest_pos (disp_snap(i,1), disp_snap(i,2),
  +c     &                            disp_snap(i,3), i)
   c... hardcoding; we should use simmetrix V_movedParamPoint to do this
  -c            mag = sqrt(x_tmp(2)*x_tmp(2) + x_tmp(3)*x_tmp(3)) ! cylinder axis is x
  -c            rad = 0.05
  -c            x_crt(2) = x_tmp(2) / mag * rad
  -c            x_crt(3) = x_tmp(3) / mag * rad
  +            mag = sqrt(x_tmp(2)*x_tmp(2) + x_tmp(3)*x_tmp(3)) ! cylinder axis is x
  +            rad = 0.05
  +            x_crt(2) = x_tmp(2) / mag * rad
  +            x_crt(3) = x_tmp(3) / mag * rad
   c... end hardcoding
   c... update iBC and BC
               iBC(i) = ibset(iBC(i), 14)
  
positive.code:         subroutine resetSnapBC ( x,       disp_snap,
       &                           iBC,     BC,     surfID_snap)
  c
          use iso_c_binding
          use simmetrix_snap
  c
          include "common.h"
  c
  
negative.code: c      module simmetrix_snap
  c      use iso_c_binding
  c
  c      interface
  c        subroutine sim_closest_pos ( dx, dy, dz, id )
  c     &    bind(C, NAME='sim_closest_pos')
  c        use iso_c_binding
  c          real(c_double), value :: dx, dy, dz
  c          integer(c_int), value :: id
  c        end subroutine
  c
  c        subroutine get_model_velocity (v)
  c     &    bind(C, NAME='get_model_velocity')
  c        use iso_c_binding
  c          real(c_double) :: v
  c        end subroutine
  c      end interface
  c
  c      end module
  c
  c----------------------------------------------------------------------
  c
  