anchor.diff:        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      CALL AllocAry( indx, n,  'IPIV', ErrStat2, ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  -   if (ErrStat >= AbortErrLev) return
  -   
  -   ! compute K, the coefficient matrix, based on the z-component of the entered key points:
  -   ! all of the coefficients will depend on kp_zr
  -   K(:,:) = 0.0_BDKi
  -
  -   K(1,3) = 2.0_BDKi
  -   K(1,4) = 6.0_BDKi*kp_coordinate(1,3)
  -   DO j=1,kp_member-1
  -      temp_id1 = (j-1)*4
  -
  -      K(temp_id1+2,temp_id1+1) = 1.0_BDKi
  -      K(temp_id1+2,temp_id1+2) = kp_coordinate(j,3)
  -      K(temp_id1+2,temp_id1+3) = kp_coordinate(j,3)**2
  -      K(temp_id1+2,temp_id1+4) = kp_coordinate(j,3)**3
  -
  -      K(temp_id1+3,temp_id1+1) = 1.0_BDKi
  -      K(temp_id1+3,temp_id1+2) = kp_coordinate(j+1,3)
  -      K(temp_id1+3,temp_id1+3) = kp_coordinate(j+1,3)**2
  -      K(temp_id1+3,temp_id1+4) = kp_coordinate(j+1,3)**3
  -   END DO
  -
  -    DO j=1,kp_member-2
  -       temp_id1 = (j-1)*4
  -
  -       K(temp_id1+4,temp_id1+2) = 1.0_BDKi
  -       K(temp_id1+4,temp_id1+3) = 2.0_BDKi*kp_coordinate(j+1,3)
  -       K(temp_id1+4,temp_id1+4) = 3.0_BDKi*kp_coordinate(j+1,3)**2
  -
  -       K(temp_id1+4,temp_id1+6) = -1.0_BDKi
  -       K(temp_id1+4,temp_id1+7) = -2.0_BDKi*kp_coordinate(j+1,3)
  -       K(temp_id1+4,temp_id1+8) = -3.0_BDKi*kp_coordinate(j+1,3)**2
  -
  -       K(temp_id1+5,temp_id1+3) = 2.0_BDKi
  -       K(temp_id1+5,temp_id1+4) = 6.0_BDKi*kp_coordinate(j+1,3)
  -
  -       K(temp_id1+5,temp_id1+7) = -2.0_BDKi
  -       K(temp_id1+5,temp_id1+8) = -6.0_BDKi*kp_coordinate(j+1,3)
  -    ENDDO
  -
  -    temp_id1 = (kp_member-2)*4
  -    K(n,temp_id1+3) = 2.0_BDKi
  -    K(n,temp_id1+4) = 6.0_BDKi*kp_coordinate(kp_member,3)
  -
  -       ! compute the factored K matrix so we can use it to solve for the coefficients later
  -
  -    CALL LAPACK_getrf( n, n, K,indx, ErrStat2, ErrMsg2)
  -       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  -
  -
  -    DO i=1,4 ! one for each column of kp_coordinate
  -
  -          ! compute the right hand side for the cubic spline fit
  -       RHS(:) = 0.0_BDKi
  -       DO j=1,kp_member-1
  +
  +   if (ErrStat < AbortErrLev) then
  +     ! note that if we return here instead, we could have a memory leak unless we deallocate the local arrays
  +   
  +      ! compute K, the coefficient matrix, based on the z-component of the entered key points:
  +      ! all of the coefficients will depend on kp_zr
  +      K(:,:) = 0.0_BDKi
  +
  +      K(1,3) = 2.0_BDKi
  +      K(1,4) = 6.0_BDKi*kp_coordinate(1,3)
  +      DO j=1,kp_member-1
  +         temp_id1 = (j-1)*4
  +
  +         K(temp_id1+2,temp_id1+1) = 1.0_BDKi
  +         K(temp_id1+2,temp_id1+2) = kp_coordinate(j,3)
  +         K(temp_id1+2,temp_id1+3) = kp_coordinate(j,3)**2
  +         K(temp_id1+2,temp_id1+4) = kp_coordinate(j,3)**3
  +
  +         K(temp_id1+3,temp_id1+1) = 1.0_BDKi
  +         K(temp_id1+3,temp_id1+2) = kp_coordinate(j+1,3)
  +         K(temp_id1+3,temp_id1+3) = kp_coordinate(j+1,3)**2
  +         K(temp_id1+3,temp_id1+4) = kp_coordinate(j+1,3)**3
  +      END DO
  +
  +       DO j=1,kp_member-2
             temp_id1 = (j-1)*4
   
  -          RHS(temp_id1+2) = kp_coordinate(j,i)
  -          RHS(temp_id1+3) = kp_coordinate(j+1,i)
  +          K(temp_id1+4,temp_id1+2) = 1.0_BDKi
  +          K(temp_id1+4,temp_id1+3) = 2.0_BDKi*kp_coordinate(j+1,3)
  +          K(temp_id1+4,temp_id1+4) = 3.0_BDKi*kp_coordinate(j+1,3)**2
  +
  +          K(temp_id1+4,temp_id1+6) = -1.0_BDKi
  +          K(temp_id1+4,temp_id1+7) = -2.0_BDKi*kp_coordinate(j+1,3)
  +          K(temp_id1+4,temp_id1+8) = -3.0_BDKi*kp_coordinate(j+1,3)**2
  +
  +          K(temp_id1+5,temp_id1+3) = 2.0_BDKi
  +          K(temp_id1+5,temp_id1+4) = 6.0_BDKi*kp_coordinate(j+1,3)
  +
  +          K(temp_id1+5,temp_id1+7) = -2.0_BDKi
  +          K(temp_id1+5,temp_id1+8) = -6.0_BDKi*kp_coordinate(j+1,3)
          ENDDO
   
  -          ! solve for the cubic-spline coefficients
  -       CALL LAPACK_getrs( 'N', n, K, indx, RHS, ErrStat2, ErrMsg2)
  +       temp_id1 = (kp_member-2)*4
  +       K(n,temp_id1+3) = 2.0_BDKi
  +       K(n,temp_id1+4) = 6.0_BDKi*kp_coordinate(kp_member,3)
  +
  +          ! compute the factored K matrix so we can use it to solve for the coefficients later
  +
  +       CALL LAPACK_getrf( n, n, K,indx, ErrStat2, ErrMsg2)
             CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   
  -          ! convert cubic-spline coefficients in RHS to output array, Coef
  -       DO j=1,kp_member-1
  -          DO m=1,4
  -             SP_Coef(j,m,i) = RHS( (j-1)*4 + m )
  +
  +       DO i=1,4 ! one for each column of kp_coordinate
  +
  +             ! compute the right hand side for the cubic spline fit
  +          RHS(:) = 0.0_BDKi
  +          DO j=1,kp_member-1
  +             temp_id1 = (j-1)*4
  +
  +             RHS(temp_id1+2) = kp_coordinate(j,i)
  +             RHS(temp_id1+3) = kp_coordinate(j+1,i)
  +          ENDDO
  +
  +             ! solve for the cubic-spline coefficients
  +          CALL LAPACK_getrs( 'N', n, K, indx, RHS, ErrStat2, ErrMsg2)
  +             CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +
  +             ! convert cubic-spline coefficients in RHS to output array, Coef
  +          DO j=1,kp_member-1
  +             DO m=1,4
  +                SP_Coef(j,m,i) = RHS( (j-1)*4 + m )
  +             ENDDO
             ENDDO
          ENDDO
  -    ENDDO
  -
  +   end if
  +
  +      ! this is the cleanup() routine:
      if (allocated(K   )) deallocate(K)
      if (allocated(RHS )) deallocate(RHS)
      if (allocated(indx)) deallocate(indx)
  
positive.code: 
        ! Incorporate boundary conditions (overwrite first node of continuous states and Acc array at t+dt)
     CALL BD_BoundaryGA2(x,p,u_interp,OtherState, ErrStat2, ErrMsg2)
     CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
     if (ErrStat >= AbortErrLev) return
  
  
        ! find x, acc, and xcc at t+dt
  
negative.code:       ! initialize states, given parameters and initial inputs (in BD coordinates)
     CALL BD_CalcIC_Position(u_tmp,p,x, ErrStat2, ErrMsg2)
       CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
     CALL BD_CalcIC_Velocity(u_tmp,p,x)
     CALL Cleanup()
  
  