anchor.diff:        ! 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)
  -     if (ErrStat >= AbortErrLev) return
      CALL BD_CalcIC_Velocity(u_tmp,p,x)
      CALL Cleanup()
   
  
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:       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry( indx, n,  'IPIV', ErrStat2, ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  
     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
  
            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
               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
     end if
  
        ! this is the cleanup() routine:
     if (allocated(K   )) deallocate(K)
     if (allocated(RHS )) deallocate(RHS)
     if (allocated(indx)) deallocate(indx)
  