anchor.diff:     
      !Create a variable that lists Y2 and adds removed constrained nodes' dofs; we will be using it to carry out other calculations with a special indexing array
      yout =0 !Initialize and populate with Y2 data  
  -   yout(1:         p%UrbarL        ) = m%UR_bar
  +   yout(1:         p%UrbarL        ) = m%UR_bar ! TODO TODO TODO Look at U_full/U_red
      yout(p%URbarL+1:p%URbarL+p%nDOFL) = m%UL
     
      !Same for a variable that deals with Udotdot
      uddout =0 !Initialize and populate with Udotdot data
  -   uddout(1          : p%URbarL          ) = m%UR_bar_dotdot
  +   uddout(1          : p%URbarL          ) = m%UR_bar_dotdot ! TODO TODO TODO look at U_full_dotdot U_red_dotdot
      uddout(p%URbarL+1 : p%URbarL+p%nDOFL  ) = m%UL_dotdot
   
      ! TODO TODO TODO, there is a lot of similarity between the three outputs sections with some code redundency
  
positive.code: 
  !------------------------------------------------------------------------------------------------------
  !>This subroutine allocates and calculated TIreact, Matrix to go from local reactions at constrained nodes to single point reactions
  SUBROUTINE ReactMatx(Init, p, WtrDpth, ErrStat, ErrMsg)
     TYPE(SD_InitType),      INTENT(  IN)  :: Init         ! Input data for initialization routine
     TYPE(SD_ParameterType), INTENT(  INOUT)  :: p         ! Parameter data
     REAL(ReKi),                   INTENT(IN)     :: WtrDpth
     INTEGER(IntKi),               INTENT(  OUT)  :: ErrStat     ! Error status of the operation
     CHARACTER(*),                 INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
     ! local variables
     INTEGER                             :: I !counter
     INTEGER                             :: rmndr !type-column index
     INTEGER                             :: nodeID  !node ID
     INTEGER(IntKi)                      :: DOFC !  DOFC = Init%NReact*6
     INTEGER(IntKi)                      :: iiDOF, iNode, nDOFPerNode ! 
     REAL(ReKi)                          :: dx, dy, dz ! distances from reaction points to subdyn origin (mudline)
     REAL(ReKi), dimension(6)            :: Line
     ErrStat=ErrID_None
     ErrMsg=""
     
     DOFC = p%nNodes_C*6 ! bjj, this is p%DOFC    !Total DOFs at the base of structure 
     
     CALL AllocAry(p%TIreact, 6, DOFC, 'p%TIReact', ErrStat, ErrMsg )
     if ( ErrStat /= ErrID_None ) return
     
     p%TIreact=0 !Initialize
     
     DO I=1,3  !Take care of first three rows
        p%TIreact(I,I:DOFC:6)=1
     ENDDO 
  
      !Other rows done per column actually  
     DO I = 1, DOFC
        nodeID = p%Nodes_C(ceiling(I/6.0),1)  !Constrained Node ID (this works in the reordered/renumbered p%Nodes_C) 
        
        dx = Init%Nodes(nodeID, 2)
        dy = Init%Nodes(nodeID, 3)
        dz = Init%Nodes(nodeID, 4) + WtrDpth
        
        rmndr = MOD(I, 6)  !It gives me the column index among the 6 different kinds
        SELECT CASE (rmndr)
           CASE (1); p%TIreact(4:6, I) = (/0.0_ReKi  , dz        , -dy/)
           CASE (2); p%TIreact(4:6, I) = (/-dz       , 0.0_ReKi , dx/)
           CASE (3); p%TIreact(4:6, I) = (/dy        , -dx       , 0.0_ReKi/)
           CASE (4); p%TIreact(4:6, I) = (/1.0_ReKi , 0.0_ReKi , 0.0_ReKi/)
           CASE (5); p%TIreact(4:6, I) = (/0.0_ReKi , 1.0_ReKi , 0.0_ReKi/)
           CASE (0); p%TIreact(4:6, I) = (/0.0_ReKi , 0.0_ReKi , 1.0_ReKi/)
           CASE DEFAULT
              ErrStat = ErrID_Fatal
              ErrMsg  = 'Error calculating transformation matrix TIreact, wrong column index '
              RETURN
           END SELECT
           !print*,'TIr',p%TIreact(:, I)
     ENDDO
  
  END SUBROUTINE ReactMatx
  
  
negative.code:          ENDDO
        ENDDO
        !Store the matrix that will let me calculate single point reaction at the base of structure
        CALL ReactMatx(Init, WtrDpth, p, ErrStat, ErrMsg)
     ENDIF
   
     ! These variables are to help follow the framework template, but the data in them is simply a copy of data
  