anchor.diff:     ! NOTE: for debugging, add preprocessor definition SD_SUMMARY_DEBUG
      !       this will add additional matrices to the SubDyn summary file.
      !............................
  -
      TYPE(ProgDesc), PARAMETER  :: SD_ProgDesc = ProgDesc( 'SubDyn', '', '' )
         
      ! ..... Public Subroutines ...................................................................................................
  -
      PUBLIC :: SD_Init                           ! Initialization routine
  -
      PUBLIC :: SD_End                            ! Ending routine (includes clean up)
  -
      PUBLIC :: SD_UpdateStates                   ! Loose coupling routine for solving for constraint states, integrating
  -                                                 
      PUBLIC :: SD_CalcOutput                     ! Routine for computing outputs
  -
      PUBLIC :: SD_CalcContStateDeriv             ! Tight coupling routine for computing derivatives of continuous states
      
  -   
   CONTAINS
  -   
   
   SUBROUTINE CreateTPMeshes( TP_RefPoint, inputMesh, outputMesh, ErrStat, ErrMsg )
  -
      REAL(ReKi),                INTENT( IN    ) :: TP_RefPoint(3)
      TYPE(MeshType),            INTENT( INOUT ) :: inputMesh
      TYPE(MeshType),            INTENT( INOUT ) :: outputMesh
      INTEGER(IntKi),            INTENT(   OUT)  :: ErrStat     ! Error status of the operation
      CHARACTER(*),              INTENT(   OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
      
  -   
      ! NOTE: The initialization of the fields for these meshes is to be handled by FAST/Driver
  -   
  -      CALL MeshCreate( BlankMesh        = inputMesh         &
  -                     ,IOS               = COMPONENT_INPUT   &
  -                     ,Nnodes            = 1                 &
  -                     ,ErrStat           = ErrStat           &
  -                     ,ErrMess           = ErrMsg            &
  -                     ,TranslationDisp   = .TRUE.            &
  -                     ,Orientation       = .TRUE.            &
  -                     ,TranslationVel    = .TRUE.            &
  -                     ,RotationVel       = .TRUE.            &
  -                     ,TranslationAcc    = .TRUE.            &
  -                     ,RotationAcc       = .TRUE.            &
  -                     )
  -      
  -      
  -         ! Create the node on the mesh
  -            
  -      CALL MeshPositionNode (   inputMesh           &
  -                              , 1                   &
  -                              , TP_RefPoint         &  
  -                              , ErrStat             &
  -                              , ErrMsg                ) !note: assumes identiy matrix as reference orientation
  -      
  -      IF ( ErrStat >= AbortErrLev ) RETURN
  -       
  -         
  -         ! Create the mesh element
  -         
  -      CALL MeshConstructElement (   inputMesh          &
  -                                  , ELEMENT_POINT      &                         
  -                                  , ErrStat            &
  -                                  , ErrMsg             &
  -                                  , 1                  &
  -                                              )
  -      CALL MeshCommit ( inputMesh   &
  -                      , ErrStat            &
  -                      , ErrMsg             )
  -   
  -      IF ( ErrStat >= AbortErrLev ) RETURN
  -      
  -         
  -         ! Create the Transition Piece reference point output mesh as a sibling copy of the input mesh
  -         
  -      CALL MeshCopy ( SrcMesh      = inputMesh              &
  -                     ,DestMesh     = outputMesh             &
  -                     ,CtrlCode     = MESH_SIBLING           &
  -                     ,IOS          = COMPONENT_OUTPUT       &
  -                     ,ErrStat      = ErrStat                &
  -                     ,ErrMess      = ErrMsg                 &
  -                     ,Force        = .TRUE.                 &
  -                     ,Moment       = .TRUE.                 &
  -                     ) 
  -     
  +   CALL MeshCreate( BlankMesh        = inputMesh         &
  +                  ,IOS               = COMPONENT_INPUT   &
  +                  ,Nnodes            = 1                 &
  +                  ,ErrStat           = ErrStat           &
  +                  ,ErrMess           = ErrMsg            &
  +                  ,TranslationDisp   = .TRUE.            &
  +                  ,Orientation       = .TRUE.            &
  +                  ,TranslationVel    = .TRUE.            &
  +                  ,RotationVel       = .TRUE.            &
  +                  ,TranslationAcc    = .TRUE.            &
  +                  ,RotationAcc       = .TRUE.            )
  +   
  +   ! Create the node on the mesh
  +   CALL MeshPositionNode (   inputMesh           &
  +                           , 1                   &
  +                           , TP_RefPoint         &
  +                           , ErrStat             &
  +                           , ErrMsg              ) !note: assumes identiy matrix as reference orientation
  +   IF ( ErrStat >= AbortErrLev ) RETURN
  +      
  +   ! Create the mesh element
  +   CALL MeshConstructElement (   inputMesh          &
  +                               , ELEMENT_POINT      &                         
  +                               , ErrStat            &
  +                               , ErrMsg             &
  +                               , 1                  )
  +   CALL MeshCommit ( inputMesh, ErrStat, ErrMsg )
  +   IF ( ErrStat >= AbortErrLev ) RETURN
  +   
  +   ! Create the Transition Piece reference point output mesh as a sibling copy of the input mesh
  +   CALL MeshCopy ( SrcMesh      = inputMesh              &
  +                  ,DestMesh     = outputMesh             &
  +                  ,CtrlCode     = MESH_SIBLING           &
  +                  ,IOS          = COMPONENT_OUTPUT       &
  +                  ,ErrStat      = ErrStat                &
  +                  ,ErrMess      = ErrMsg                 &
  +                  ,Force        = .TRUE.                 &
  +                  ,Moment       = .TRUE.                 ) 
   END SUBROUTINE CreateTPMeshes
   
   SUBROUTINE CreateY2Meshes( NNode, Nodes, NNodes_I, IDI, NNodes_L, IDL, NNodes_C, IDC, inputMesh, outputMesh, ErrStat, ErrMsg )
  -
      INTEGER(IntKi),            INTENT( IN    ) :: NNode                     !total number of nodes in the structure, used to size the array Nodes, i.e. its rows
      REAL(ReKi),                INTENT( IN    ) :: Nodes(NNode, JointsCol)
      INTEGER(IntKi),            INTENT( IN    ) :: NNodes_I                  ! number interface nodes   i.e. Y2 stuff at the beginning
  
positive.code:       Init%RootName = TRIM(InitInput%RootName)//'.SD'
     END IF
     
     
     !.................................
     ! Parse the SubDyn inputs 
     !.................................
              
     CALL SD_Input(InitInput%SDInputFile, Init, p, ErrStat2, ErrMsg2)
        CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )
        IF ( ErrStat >= AbortErrLev ) THEN
           CALL CleanUp()
           RETURN
        END IF
  
        
     
     !.................................
     !-------------  Discretize the structure according to the division size -----------------
     !.................................
           
     CALL SD_Discrt(Init,p, ErrStat2, ErrMsg2) ! sets Init%NNode, Init%NElm
        CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )
        IF ( ErrStat >= AbortErrLev ) THEN
           CALL CleanUp()
           RETURN
        END IF
     
        
        
     CALL AssembleKM(Init,p, ErrStat2, ErrMsg2)
        CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )
        IF ( ErrStat >= AbortErrLev ) THEN
           CALL CleanUp()
           RETURN
        END IF
      
  
     !.................................
     ! Calculate values for FEMparams (for summary file output only
     !.................................
     
           ! Solve dynamics problem
           
     FEMparams%NOmega = Init%TDOF - p%Nreact*6 !removed an extra "-6"  !Note if fixity changes at the reaction points, this will need to change
       
     CALL AllocAry(FEMparams%Omega,            FEMparams%NOmega, 'FEMparams%Omega', ErrStat2, ErrMsg2 ); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )   
     CALL AllocAry(FEMparams%Modes, Init%TDOF, FEMparams%NOmega, 'FEMparams%Modes', ErrStat2, ErrMsg2 ); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )
        IF ( ErrStat >= AbortErrLev ) THEN
           CALL CleanUp()
           RETURN
        END IF
     
     
        ! We call the EigenSolver here only so that we get a print-out the eigenvalues from the full system (minus Reaction DOF)
        ! The results, Phi is not used in the remainder of this Init subroutine, Omega goes to outsummary.
        
     CALL EigenSolve( Init%K, Init%M, Init%TDOF, FEMparams%NOmega, .True., Init, p, FEMparams%Modes, FEMparams%Omega, ErrStat2, ErrMsg2 )
        CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )
        IF ( ErrStat >= AbortErrLev ) THEN
           CALL CleanUp()
           RETURN
        END IF 
     
  
     !.................................
     ! Craig-Bampton reduction
     !.................................
        
           
  
     CALL Craig_Bampton(Init, p, CBparams, ErrStat2, ErrMsg2) !sets many parameters...
        CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )
        IF ( ErrStat >= AbortErrLev ) THEN
           CALL CleanUp()
           RETURN
        END IF
  
         
     !.................................
     ! Define initial system states here:
     !.................................
     IF ( p%qmL > 0 ) THEN
     
        CALL AllocAry(x%qm,       p%qmL, 'x%qm',       ErrStat2, ErrMsg2 ); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )
        CALL AllocAry(x%qmdot,    p%qmL, 'x%qmdot',    ErrStat2, ErrMsg2 ); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )
        CALL AllocAry(m%qmdotdot, p%qmL, 'm%qmdotdot', ErrStat2, ErrMsg2 ); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_Init' )
           IF ( ErrStat >= AbortErrLev ) THEN
              CALL CleanUp()
              RETURN
           END IF   
        x%qm      = 0.0_ReKi   
        x%qmdot   = 0.0_ReKi
        m%qmdotdot= 0.0_ReKi
        
     END IF
     
     
     xd%DummyDiscState  = 0.0_ReKi
     z%DummyConstrState = 0.0_ReKi
     
     
           ! Allocate OtherState%xdot if using multi-step method; initialize n
     IF ( ( p%IntMethod .eq. 2) .OR. ( p%IntMethod .eq. 3)) THEN
        !bjj: note that the way SD_UpdateStates is implemented, "n" doesn't need to be initialized here
        Allocate( OtherState%xdot(4), STAT=ErrStat2 )
  
negative.code: LOGICAL                      :: Echo  
  INTEGER(IntKi)               :: UnIn
  INTEGER(IntKi)               :: IOS
  INTEGER(IntKi)               :: UnEc   !Echo file ID
  
  REAL(ReKi),PARAMETER        :: WrongNo=-9999.   ! Placeholder value for bad(old) values in JDampings
  
  INTEGER(IntKi)               :: I, J, flg, K
  REAL(ReKi)                   :: Dummy_ReAry(SDMaxInpCols) 
  INTEGER(IntKi)               :: Dummy_IntAry(SDMaxInpCols)
  INTEGER(IntKi)       :: ErrStat2
  CHARACTER(ErrMsgLen) :: ErrMsg2
  ! Initialize ErrStat
  ErrStat = ErrID_None
  ErrMsg  = ""
  
  UnEc = -1 
  Echo = .FALSE.
  
  CALL GetNewUnit( UnIn )   
    
  CALL OpenFInpfile(UnIn, TRIM(SDInputFile), ErrStat2, ErrMsg2)
  
  IF ( ErrStat2 /= ErrID_None ) THEN
     Call Abort('Could not open SubDyn input file')
     return
  END IF
  
  CALL GetPath( SDInputFile, PriPath )    ! Input files will be relative to the path where the primary input file is located.
  
  
  !-------------------------- HEADER ---------------------------------------------
  CALL ReadCom( UnIn, SDInputFile, 'SubDyn input file header line 1', ErrStat2, ErrMsg2 ); if(Failed()) return
  CALL ReadCom( UnIn, SDInputFile, 'SubDyn input file header line 2', ErrStat2, ErrMsg2 ); if(Failed()) return
  
  !-------------------------- SIMULATION CONTROL PARAMETERS ----------------------
  CALL ReadCom( UnIn, SDInputFile, ' SIMULATION CONTROL PARAMETERS ', ErrStat2, ErrMsg2 ); if(Failed()) return
  CALL ReadVar(UnIn, SDInputFile, Echo, 'Echo', 'Echo Input File Logic Variable',ErrStat2, ErrMsg2); if(Failed()) return
  
  IF ( Echo )  THEN 
     CALL OpenEcho ( UnEc, TRIM(Init%RootName)//'.ech' ,ErrStat2, ErrMsg2)
     IF ( ErrStat2 /= 0 ) THEN
        CALL Abort("Could not open SubDyn echo file")
        return
     END IF
     REWIND(UnIn)
     !bjj: note we don't need to do error checking here; it was already checked (this is just a repeat of above)
     CALL ReadCom( UnIn, SDInputFile, 'SubDyn input file header line 1', ErrStat2, ErrMsg2 )
     CALL ReadCom( UnIn, SDInputFile, 'SubDyn input file header line 2', ErrStat2, ErrMsg2 )
     CALL ReadCom( UnIn, SDInputFile, 'SIMULATION CONTROL PARAMETERS'  , ErrStat2, ErrMsg2, UnEc )
     CALL ReadVar( UnIn, SDInputFile, Echo, 'Echo', 'Echo Input File Logic Variable',ErrStat2, ErrMsg2, UnEc )
  ENDIF 
  
  ! Read time step   ("default" means use the glue-code default)
  CALL ReadVar( UnIn, SDInputFile, Line, 'SDdeltaT', 'Subdyn Time Step',ErrStat2, ErrMsg2, UnEc ); if(Failed()) return
  
  CALL Conv2UC( Line )    ! Convert Line to upper case.
  IF ( TRIM(Line) == 'DEFAULT' )  THEN   ! .TRUE. when one wants to use the default value timestep provided by the glue code.
      p%SDdeltaT=Init%DT
  ELSE                                   ! The input must have been specified numerically.
     READ (Line,*,IOSTAT=IOS)  p%SDdeltaT
     CALL CheckIOS ( IOS, SDInputFile, 'SDdeltaT', NumType, ErrStat2,ErrMsg2 ); if(Failed()) return
  
     IF ( ( p%SDdeltaT <=  0 ) )  THEN 
        call Abort('SDdeltaT must be greater than or equal to 0.')
        return         
     END IF  
  END IF
        
  CALL ReadVar ( UnIn, SDInputFile, p%IntMethod, 'IntMethod', 'Integration Method',ErrStat2, ErrMsg2, UnEc ); if(Failed()) return
  CALL ReadLVar(UnIn, SDInputFile, p%SttcSolve, 'SttcSolve', 'Solve dynamics about static equilibrium point', ErrStat2, ErrMsg2, UnEc); if(Failed()) return
  !-------------------- FEA and CRAIG-BAMPTON PARAMETERS---------------------------
  CALL ReadCom  ( UnIn, SDInputFile, ' FEA and CRAIG-BAMPTON PARAMETERS ', ErrStat2, ErrMsg2, UnEc ); if(Failed()) return
  CALL ReadIVar ( UnIn, SDInputFile, Init%FEMMod, 'FEMMod', 'FEM analysis mode'             ,ErrStat2, ErrMsg2, UnEc ); if(Failed()) return ! 0= Euler-Bernoulli(E-B); 1=Tapered E-B; 2= Timoshenko; 3= tapered Timoshenko
  CALL ReadIVar ( UnIn, SDInputFile, Init%NDiv  , 'NDiv'  , 'Number of divisions per member',ErrStat2, ErrMsg2, UnEc ); if(Failed()) return
  CALL ReadLVar ( UnIn, SDInputFile, Init%CBMod , 'CBMod' , 'C-B mod flag'                  ,ErrStat2, ErrMsg2, UnEc ); if(Failed()) return
  
  IF (Check( (p%IntMethod < 1) .OR.(p%IntMethod > 4)     , 'IntMethod must be 1 through 4.')) return
  IF (Check( (Init%FEMMod < 0 ) .OR. ( Init%FEMMod > 4 ) , 'FEMMod must be 0, 1, 2, or 3.')) return
  IF (Check( Init%NDiv < 1                               , 'NDiv must be a positive integer')) return
  
  IF (Init%CBMod) THEN
     ! Nmodes - Number of interal modes to retain.
     CALL ReadIVar ( UnIn, SDInputFile, p%Nmodes, 'Nmodes', 'Number of internal modes',ErrStat2, ErrMsg2, UnEc ); if(Failed()) return
  
     IF (Check( p%Nmodes < 0 , 'Nmodes must be a non-negative integer.')) return
     
     if ( p%Nmodes > 0 ) THEN
        ! Damping ratios for retained modes
        CALL AllocAry(Init%JDampings, p%Nmodes, 'JDamping', ErrStat2, ErrMsg2) ; if(Failed()) return
        Init%JDampings=WrongNo !Initialize
     
        CALL ReadAry( UnIn, SDInputFile, Init%JDampings, p%Nmodes, 'JDamping', 'Damping ratio of the internal modes', ErrStat2, ErrMsg2, UnEc );
        ! note that we don't check the ErrStat2 here; if the user entered fewer than Nmodes values, we will use the
        ! last entry to fill in remaining values.
        !Check 1st value, we need at least one good value from user or throw error
        IF ((Init%JDampings(1) .LT. 0 ) .OR. (Init%JDampings(1) .GE. 100.0)) THEN
              CALL Abort('Damping ratio should be larger than 0 and less than 100')
              return
        ELSE
           DO I = 2, p%Nmodes
              IF ( Init%JDampings(I) .EQ. WrongNo ) THEN
  