anchor.diff:  
           R = p%StrWd/2 
           
  -        IF ( OtherState%Splitf == 0 ) THEN
  -   
  -           IF ((x%q - u%PointMesh%TranslationDisp(1,1)) >= OtherState%dxc .AND. (x%q - u%PointMesh%TranslationDisp(1,1)) < OtherState%dxc+R ) THEN
  -   
  -!bjj: force is NaN sometimes because  R**2 - (R - x%q + u%PointMesh%TranslationDisp(1,1))**2 < 0             
  -                force = -p%Cpa * ( 2 * p%h * ( R**2 - (R - x%q + u%PointMesh%TranslationDisp(1,1))**2 )**0.5 )**( p%dpa + 1 ) * 1.0e6 + p%FdrN
  -         
  -           ELSE IF (  (x%q - u%PointMesh%TranslationDisp(1,1)) >= OtherState%dxc+R ) THEN
  -             
  -                force = -p%Cpa * ( 2 * R *  p%h )**( p%dpa + 1 ) * 1.0e6  + p%FdrN
  -            
  -           ELSE
  -   
  -               force = 0. + p%FdrN
  -   
  +       IF ( OtherState%Splitf == 0 ) THEN
  +          
  +           IF ((x%q - u%PointMesh%TranslationDisp(1,1)) < OtherState%dxc) THEN
  +          
  +               force = 0 + p%FdrN
  +
  +           ELSE IF ((x%q - u%PointMesh%TranslationDisp(1,1)) >= OtherState%dxc) THEN
  +
  +                IF (x%q - u%PointMesh%TranslationDisp(1,1) < R) THEN
  +      
  +                    force = -p%Cpa * ( 2 * p%h * ( R**2 - (R - x%q + u%PointMesh%TranslationDisp(1,1))**2 )**0.5 )**( p%dpa + 1 ) * 1.0e6 + p%FdrN
  +           
  +                ELSE
  +
  +                    force = -p%Cpa * ( 2 * R *  p%h )**( p%dpa + 1 ) * 1.0e6 + p%FdrN
  +
  +                ENDIF
  + 
              ENDIF
  -           
  -         ELSE
  -             
  -             force = 0
  -             
  -         ENDIF
  -         
  +
  +        ELSE
  +
  +              force = 0 
  + 
  +        ENDIF
  +        
            xdot%q = x%dqdt
      
            xdot%dqdt = force / p%Mice
  
positive.code:       
             R = p%StrWd/2
        
           IF (OtherState%Splitf == 0 ) THEN
               
             IF ((x%q - u%PointMesh%TranslationDisp(1,1)) >= OtherState%dxc .AND. (x%q - u%PointMesh%TranslationDisp(1,1)) < OtherState%dxc+R ) THEN
           
                  y%PointMesh%Force(1,1) =  p%Cpa * ( 2 * p%h * ( R**2 - (R - x%q + u%PointMesh%TranslationDisp(1,1))**2 )**0.5 )**( p%dpa + 1 ) * 1.0e6 
              
             ELSE IF (  (x%q - u%PointMesh%TranslationDisp(1,1)) >= OtherState%dxc+R ) THEN
               
                 y%PointMesh%Force(1,1) = p%Cpa * ( 2 * R *  p%h )**( p%dpa + 1 ) * 1.0e6 
               
             ELSE
               
                 y%PointMesh%Force(1,1) = 0 
                
             ENDIF
             
           ELSE 
               
                y%PointMesh%Force(1,1) = 0 
             
           ENDIF
        
        ENDIF
        
  
negative.code:       
        IF ( OtherState%Splitf == 0 ) THEN
            
         IF ((x%q - u%PointMesh%TranslationDisp(1,1)) < OtherState%dxc) THEN
            
             IceForce = 0
  
         ELSE IF ((x%q - u%PointMesh%TranslationDisp(1,1)) >= OtherState%dxc) THEN
  
             IF (x%q - u%PointMesh%TranslationDisp(1,1) < R) THEN
        
                IceForce =  p%Cpa * ( 2 * p%h * ( R**2 - (R - x%q + u%PointMesh%TranslationDisp(1,1))**2 )**0.5 )**( p%dpa + 1 ) * 1.0e6
             
             ELSE
  
                IceForce = p%Cpa * ( 2 * R *  p%h )**( p%dpa + 1 ) * 1.0e6 
  
             ENDIF
  
         ENDIF
  
      ELSE
  
        IceForce = 0 
  
      ENDIF
        
        IF ( IceForce >= p%Fsp ) THEN
            OtherState%Splitf = 1
  