anchor.diff:        INTEGER(IntKi), pointer      :: DOFList(:)
         INTEGER(IntKi)               :: startDOF
         REAL(ReKi)                   :: DCM(3,3)
  -      REAL(ReKi)                   :: HydroForces(6*p%nNodes_I) !  !Forces from all interface nodes listed in one big array  ( those translated to TP ref point HydroTP(6) are implicitly calculated in the equations)
  +      REAL(ReKi)                   :: F_I(6*p%nNodes_I) !  !Forces from all interface nodes listed in one big array  ( those translated to TP ref point HydroTP(6) are implicitly calculated in the equations)
         TYPE(SD_ContinuousStateType) :: dxdt        ! Continuous state derivatives at t- for output file qmdotdot purposes only
         INTEGER(IntKi)               :: ErrStat2    ! Error status of the operation (occurs after initial error)
         CHARACTER(ErrMsgLen)         :: ErrMsg2     ! Error message if ErrStat2 /= ErrID_None
  
positive.code:       !Y1= TP reaction Forces, i.e. force that the jacket exerts onto the TP and above  
        ! ---------------------------------------------------------------------------------
        ! Eq. 15: Y1 = -(C1*x + D1*u + FY)  [note the negative sign!!!!]
        !NEED TO ADD HYDRODYNAMIC FORCES AT THE Interface NODES
          !Aggregate the forces and moments at the interface nodes to the reference point
          !TODO: where are these HydroTP, HydroForces documented?
        DO I = 1, p%nNodes_I 
           iSDNode = p%Nodes_I(I,1)
           iY2Node = iSDNode
           startDOF = (I-1)*6 + 1 ! NOTE: this works since interface is assumed to be sorted like LMesh and have 6 DOF per nodes
           !Take care of Hydrodynamic Forces that will go into INterface Forces later
           HydroForces(startDOF:startDOF+5) =  (/u%LMesh%Force(1:3,iY2Node),u%LMesh%Moment(1:3,iY2Node)/)  !(6,NNODES_I)
        ENDDO
                  
        !HydroTP =  matmul(transpose(p%TI),HydroForces) ! (6,1) calculated below
        ! note: matmul( HydroForces, p%TI ) = matmul( transpose(p%TI), HydroForces) because HydroForces is 1-D            
        IF ( p%nDOFM > 0) THEN
           Y1 = -(   matmul(p%C1_11, x%qm) + matmul(p%C1_12,x%qmdot)                                    &  ! -(   C1(1,1)*x(1) + C1(1,2)*x(2)
                   + matmul(p%KBB,   m%u_TP) + matmul(p%D1_12, m%udot_TP) + matmul(p%D1_13, m%udotdot_TP) + matmul(p%D1_14, m%UFL)   &  !    + D1(1,1)*u(1) + 0*u(2) + D1(1,3)*u(3) + D1(1,4)*u(4)
                   - matmul( HydroForces, p%TI )  + p%FY )                                                                            !    + D1(1,5)*u(5) + Fy(1) )
        ELSE ! No retained modes, so there are no states
           Y1 = -( matmul(p%KBB,   m%u_TP)   + matmul(p%D1_12, m%udot_TP) + matmul(p%D1_13, m%udotdot_TP) + matmul(p%D1_14, m%UFL)   &  ! -(  0*x + D1(1,1)*u(1) + 0*u(2) + D1(1,3)*u(3) + D1(1,4)*u(4)
                   - matmul( HydroForces, p%TI )  + p%FY )                                             !    + D1(1,5)*u(5) + Fy(1) )
        END IF
        ! Computing extra moments due to lever arm introduced by interface displacement
        !               Y1(:3) = -f_TP
  
negative.code:       m%udot_TP    = (/u%TPMesh%TranslationVel( :,1), u%TPMesh%RotationVel(:,1)/)
        m%udotdot_TP = (/u%TPMesh%TranslationAcc( :,1), u%TPMesh%RotationAcc(:,1)/)
  
        ! External force on internal and interface nodes
        call GetExtForceOnInternalDOF( u, p, m, m%UFL, ErrStat2, ErrMsg2 ); if(Failed()) return
        call GetExtForceOnInterfaceDOF(p, m%Fext, F_I)
  
        !________________________________________
        ! Set motion outputs on y%Y2mesh
  