anchor.diff:                    ,ErrStat           = ErrStat                             &
                     ,ErrMess           = ErrMsg                              &
                     ,Force             = .TRUE.                              &
  -                  ,Moment            = .TRUE.                              &
  -                  )
  -      
  -   
  +                  ,Moment            = .TRUE.                              )
      !---------------------------------------------------------------------
      !    Interface nodes
      !---------------------------------------------------------------------
  -   
      DO I = 1,NNodes_I 
  -        
  -         ! Create the node on the mesh
  +      ! Create the node on the mesh
         nodeIndx = IDI(I*6) / 6     !integer division gives me the actual node index, is it true? Yes it is not the nodeID
  -      
         CALL MeshPositionNode (   inputMesh           &
                                 , I                   &
                                 , Nodes(nodeIndx,2:4) &  ! position
                                 , ErrStat             &
  -                              , ErrMsg                )
  -      
  +                              , ErrMsg              )
         IF ( ErrStat /= ErrID_None ) RETURN
  -       
  -      
  -         ! Create the mesh element
  -         
  +
  +      ! Create the mesh element
         CALL MeshConstructElement (   inputMesh          &
                                     , ELEMENT_POINT      &                         
                                     , ErrStat            &
                                     , ErrMsg             &
  -                                  , I                  &
  -                                              )
  -         
  +                                  , I                  )
      END DO
  -     
      
      !---------------------------------------------------------------------
      !    Interior nodes
      !---------------------------------------------------------------------
  -   
      DO I = 1,NNodes_L 
  -        
  -         ! Create the node on the mesh
  +      ! Create the node on the mesh
         nodeIndx = IDL(I*6) / 6     !integer division gives me the actual node index, is it true? Yes it is not the nodeID of the input file that may not be sequential, but the renumbered list of nodes
  -      
         CALL MeshPositionNode (   inputMesh           &
                                 , I + NNodes_I        &
  -                              , Nodes(nodeIndx,2:4) &  
  +                              , Nodes(nodeIndx,2:4) &
                                 , ErrStat             &
  -                              , ErrMsg                )
  -      
  +                              , ErrMsg              )
         IF ( ErrStat /= ErrID_None ) RETURN
  -       
  -      
  -         ! Create the mesh element
  -         
  +
  +      ! Create the mesh element
         CALL MeshConstructElement (   inputMesh          &
                                     , ELEMENT_POINT      &                         
                                     , ErrStat            &
                                     , ErrMsg             &
  -                                  , I + NNodes_I       &
  -                                              )
  -         
  +                                  , I + NNodes_I       )
      END DO
  -   
      
      !---------------------------------------------------------------------
      !    Base Reaction nodes
      !---------------------------------------------------------------------
  -   
      DO I = 1,NNodes_C 
  -        
  -         ! Create the node on the mesh
  +      ! Create the node on the mesh
         nodeIndx = IDC(I*6) / 6     !integer division gives me the actual node index, is it true? Yes it is not the nodeID
  -      
         CALL MeshPositionNode (   inputMesh                 &
                                 , I + NNodes_I + NNodes_L   &
                                 , Nodes(nodeIndx,2:4)       &  
                                 , ErrStat                   &
  -                              , ErrMsg                )
  -      
  +                              , ErrMsg                    )
         IF ( ErrStat /= ErrID_None ) RETURN
  -       
  -      
  -         ! Create the mesh element
  -         
  +      
  +      ! Create the mesh element
         CALL MeshConstructElement (   inputMesh                 &
                                     , ELEMENT_POINT             &                         
                                     , ErrStat                   &
                                     , ErrMsg                    &
  -                                  , I + NNodes_I + NNodes_L   &
  -                                              )
  -         
  +                                  , I + NNodes_I + NNodes_L   )
      END DO
  -     
  -   CALL MeshCommit ( inputMesh   &
  -                   , ErrStat     &
  -                   , ErrMsg       )
  -   
  +   CALL MeshCommit ( inputMesh, ErrStat, ErrMsg )
      IF ( ErrStat /= ErrID_None ) RETURN
  -      
  -         
  -         ! Create the Interior Points output mesh as a sibling copy of the input mesh
  -         
  +         
  +   ! Create the Interior Points output mesh as a sibling copy of the input mesh
      CALL MeshCopy (    SrcMesh      = inputMesh              &
                        ,DestMesh     = outputMesh             &
                        ,CtrlCode     = MESH_SIBLING           &
  
positive.code:       END DO
        
        !Repeat for the acceleration, there should be a way to combine into 1 loop
        
        L1 = p%NNodes_I+1
        L2 = p%NNodes_I+p%NNodes_L
        junk=   RESHAPE(m%UL_dotdot,(/6  ,p%NNodes_L/)) 
        y%Y2mesh%TranslationAcc (  :,L1:L2)   = junk(1:3,:) 
        y%Y2mesh%RotationAcc    (  :,L1:L2)   = junk(4:6,:) 
   
   
        
        ! ---------------------------------------------------------------------------------
        ! Base reaction nodes
        ! ---------------------------------------------------------------------------------
        
        L1 = p%NNodes_I+p%NNodes_L+1   
        L2 = p%NNodes_I+p%NNodes_L+p%NReact
  
        y%Y2mesh%TranslationDisp(  :,L1:L2)   = 0.0
        CALL Eye( y%Y2mesh%Orientation(:,:,L1:L2), ErrStat2, ErrMsg2 )  ! set this orientation to the identity matrix
           CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_CalcOutput')
  
        y%Y2mesh%TranslationVel (  :,L1:L2)   = 0.0
        y%Y2mesh%RotationVel    (  :,L1:L2)   = 0.0
        y%Y2mesh%TranslationAcc (  :,L1:L2)   = 0.0
        y%Y2mesh%RotationAcc    (  :,L1:L2)   = 0.0
           
  
     !________________________________________
     ! Set loads outputs on y%Y1Mesh
     !________________________________________
        
        
        ! ---------------------------------------------------------------------------------
        !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 
       
           startDOF = (I-1)*6 + 1
                    
           !Take care of Hydrodynamic Forces that will go into INterface Forces later
           HydroForces(startDOF:startDOF+5 ) =  (/u%LMesh%Force(:,I),u%LMesh%Moment(:,I)/)  !(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%qml > 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_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%MBB, 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
        
  
        
        ! values on the interface mesh are Y1 (SubDyn forces) + Hydrodynamic forces
        y%Y1Mesh%Force (:,1) = Y1(1:3) 
        y%Y1Mesh%Moment(:,1) = Y1(4:6) 
  
              
     !________________________________________
     ! CALCULATE OUTPUT TO BE WRITTEN TO FILE 
     !________________________________________
  
        
           ! OutSwtch determines whether or not to actually output results via the WriteOutput array
           !    0 = No one needs the SubDyn outputs provided via the WriteOutput array.
           !    1 = SubDyn will generate an output file of its own.  
           !    2 = the caller will handle the outputs, but SubDyn needs to provide them.
           !    3 = Both 1 and 2
  
        IF ( p%OutSwtch > 0 ) THEN
                      
           !_________________________________________________!
           ! call CalcContStateDeriv one more time to store these qmdotdot for debugging purposes in the output file
           !find xdot at t
           
           IF ( p%NModes > 0 ) THEN
              ! note that this re-sets m%udotdot_TP and m%UFL, but they are the same values as earlier in this routine so it doesn't change results in SDOut_MapOutputs()
              CALL SD_CalcContStateDeriv( t, u, p, x, xd, z, OtherState, m, dxdt, ErrStat2, ErrMsg2 )
                 CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_CalcOutput')
  
                 IF ( ErrStat >= AbortErrLev ) THEN
                    CALL Cleanup
                    RETURN
                 END IF
  
              !Assign the acceleration to the x variable since it will be used for output file purposes for SSqmdd01-99, and dxdt will disappear
              m%qmdotdot=dxdt%qmdot
              ! Destroy dxdt because it is not necessary for the rest of the subroutine
              CALL SD_DestroyContState( dxdt, ErrStat2, ErrMsg2)
                 CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_CalcOutput')
                 
           END IF
           !_________________________________________________!
            
                               
              ! Write the previous output data into the output file           
           IF ( ( p%OutSwtch == 1 .OR. p%OutSwtch == 3 ) .AND. ( t > m%LastOutTime ) ) THEN
              IF ((m%Decimat .EQ. p%OutDec) .OR. (m%Decimat .EQ. 0))  THEN
                 m%Decimat=1  !reset counter
                 CALL SDOut_WriteOutputs( p%UnJckF, m%LastOutTime, m%SDWrOutput, p, ErrStat2, ErrMsg2 )   
                    CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_CalcOutput')
                    IF ( ErrStat >= AbortErrLev ) THEN
                       CALL Cleanup
                       RETURN
                    END IF                        
              ELSE      
                 m%Decimat=m%Decimat+1
              ENDIF
           END IF        
           
              ! Map calculated results into the AllOuts Array + perform averaging and all necessary extra calculations
           CALL SDOut_MapOutputs(t, u,p,x,y, m, AllOuts, ErrStat2, ErrMsg2)
              CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_CalcOutput')
              IF ( ErrStat >= AbortErrLev ) THEN
                 CALL Cleanup
                 RETURN
              END IF                        
              
           ! Put the output data in the WriteOutput array
           DO I = 1,p%NumOuts+p%OutAllInt*p%OutAllDims
  
negative.code:             ENDIF    
          ENDDO
      ENDDO
  
     !--------------------------------------
     ! write discretized data to a txt file
     WRITE(*, '(A,e15.6)')  'Matrix SYmmetry Check: Largest (abs) relative error is:', MaxErr
     WRITE(*, '(A,I4,I4)')  'Matrix SYmmetry Check: (I,J)=', imax,jmax
  
  END SUBROUTINE SymMatDebug
  
  End Module SubDyn
  