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
  -         
  -                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
  +          IF ( OtherState%Splitf == 0 ) THEN
  +          
  +              IF ((x%q - u%PointMesh%TranslationDisp(1,1)) < OtherState%dxc) THEN
  +          
  +                 y%PointMesh%Force(1,1) = 0
  +
  +              ELSE IF ((x%q - u%PointMesh%TranslationDisp(1,1)) >= OtherState%dxc) THEN
  +
  +                   IF (x%q - u%PointMesh%TranslationDisp(1,1) < 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 
  -             
  -              y%PointMesh%Force(1,1) = 0 
  -           
  -         ENDIF
  +                   ELSE
  +
  +                        y%PointMesh%Force(1,1) = p%Cpa * ( 2 * R *  p%h )**( p%dpa + 1 ) * 1.0e6 
  +
  +                   ENDIF
  +
  +              ENDIF
  +
  +          ELSE
  +
  +                 y%PointMesh%Force(1,1) = 0 
  +
  +          ENDIF
         
         ENDIF
         
  
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
     
  !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
     
             ENDIF
             
           ELSE
               
               force = 0
               
           ENDIF
           
           xdot%q = x%dqdt
     
           xdot%dqdt = force / p%Mice
  
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
  