anchor.diff:                     
                     EXIT   !We found the element for that node, exit loop on elements
                 ENDIF
  -              
  -              
             ENDDO
  -          
         ENDDO
  -
       ENDDO    
  -    
  +   ENDIF
    
  - ENDIF
  +   !_____________________________________REACTIONS_____________________________________________
  +   p%OutReact = .FALSE.
  +   DO I=1,p%NumOuts
  +      if ( ANY( p%OutParam(I)%Indx == ReactSS) ) THEN ! bjj: removed check of first 5 characters being "React" because (1) cases matter and (2) we can also ask for "-React*" or "mREACT"
  +         p%OutReact   =.TRUE.  
  +         EXIT
  +      ENDIF
  +   ENDDO
    
  - !_____________________________________REACTIONS_____________________________________________
  -p%OutReact = .FALSE.
  -DO I=1,p%NumOuts
  -   if ( ANY( p%OutParam(I)%Indx == ReactSS) ) THEN ! bjj: removed check of first 5 characters being "React" because (1) cases matter and (2) we can also ask for "-React*" or "mREACT"
  -      p%OutReact   =.TRUE.  
  -      EXIT
  +   IF (p%OutReact) THEN  !I need to store all constrained forces and moments; WE do not allow more than one member to be connected at a constrained joint for the time being
  +
  +      ALLOCATE ( p%MOutLst3(p%NReact), STAT = ErrStat2)     !this list contains different arrays for each of its elements
  +      ErrMsg2 = 'Error allocating p%MOutLst3 array in SDOut_Init'
  +      if(Failed()) return
  +
  +      DO I=1,p%NReact  !For all constrained node
  +         p%MOutLst3(I)%Noutcnt=p%Reacts(I,1) !Assign nodeID for list I, I am using Noutcnt as a temporary holder for it, since nodeID is n array
  +         NconEls=Init%NodesConnE(p%MoutLst3(I)%Noutcnt,1) !Number of elements connecting to the joint
  +         ! ElmIDs: element IDs connecting to the joint; (1,NconEls) and not (NconEls) as the same meshauxtype is used with other Moutlst
  +         ! Me: has for each selected joint, and for each element attached to that node, a 12x12 matrix (extra dimension redundant)
  +         ! Ke: has for each selected joint, and for each element attached to that node  a 12x12 matrix
  +         CALL AllocAry(p%MOutLst3(I)%ElmIDs,      1, NconEls, ' p%MOutLst3(I)%ElmIds', ErrStat2, ErrMsg2); if(Failed()) return
  +         CALL AllocAry(p%MOutLst3(I)%ElmNds,      1, NconEls, ' p%MOutLst3(I)%ElmNds', ErrStat2, ErrMsg2); if(Failed()) return
  +         CALL AllocAry(p%MOutLst3(I)%Me, 12, 12 , 1, NconEls, ' p%MOutLst3(I)%Me'    , ErrStat2, ErrMsg2); if(Failed()) return
  +         CALL AllocAry(p%MOutLst3(I)%Ke, 12, 12 , 1, NconEls, ' p%MOutLst3(I)%Ke'    , ErrStat2, ErrMsg2); if(Failed()) return
  +         CALL AllocAry(p%MOutLst3(I)%Fg,     12 , 1, NconEls, ' p%MOutLst3(I)%Fg'    , ErrStat2, ErrMsg2); if(Failed()) return
  +
  +         DO K=1, NconEls 
  +            L=Init%NodesConnE(p%MoutLst3(I)%Noutcnt,k+1)  !k-th Element Number in the set of elements attached to the selected node 
  +
  +            p%MoutLst3(I)%ElmIDs(1,K)=L     !This array has for each joint requested  the elements' ID to get results for  
  +
  +            M=p%Elems(L,2:3) !1st and 2nd node of the k-th element
  +
  +            !Select whether the joint is the 1st or second node of the element
  +            p%MoutLst3(I)%ElmNds(1,K)=1                        !store whether first or second node of element  
  +            IF (M(2) .EQ. p%MoutLst3(I)%Noutcnt ) p%MoutLst3(I)%ElmNds(1,K)=2 !store whether first or second node of element  
  +
  +            !Calculate Ke, Me to be used for output
  +            CALL ElemK_Beam( p%elemprops(L)%Area, p%elemprops(L)%Length, p%elemprops(L)%Ixx, p%elemprops(L)%Iyy, &
  +               p%elemprops(L)%Jzz, p%elemprops(L)%Shear, p%elemprops(L)%kappa, p%elemprops(L)%YoungE,  & 
  +               p%elemprops(L)%ShearG, p%elemprops(L)%DirCos, p%MoutLst3(I)%Ke(:,:,1,K) )
  +            CALL ElemM_Beam( p%elemprops(L)%Area, p%elemprops(L)%Length, p%elemprops(L)%Ixx, p%elemprops(L)%Iyy,&
  +               p%elemprops(L)%Jzz,  p%elemprops(L)%rho,  p%elemprops(L)%DirCos, p%MoutLst3(I)%Me(:,:,1,K) )   
  +            CALL ElemG( p%elemprops(L)%Area, p%elemprops(L)%Length, p%elemprops(L)%rho,  p%elemprops(L)%DirCos, p%MoutLst3(I)%Fg(:,1,K), Init%g )   
  +         ENDDO
  +      ENDDO
  +      !Store the matrix that will let me calculate single point reaction at the base of structure
  +      CALL ReactMatx(Init, p, WtrDpth, ErrStat, ErrMsg)
      ENDIF
  -ENDDO
    
  - IF (p%OutReact) THEN  !I need to store all constrained forces and moments; WE do not allow more than one member to be connected at a constrained joint for the time being
  -    
  -    ALLOCATE ( p%MOutLst3(p%NReact), STAT = ErrStat )     !this list contains different arrays for each of its elements
  -    IF ( ErrStat /= ErrID_None ) THEN
  -        ErrStat = ErrID_Fatal
  -         ErrMsg  = 'Error allocating p%MOutLst3 array in SDOut_Init'
  -         RETURN
  -    END IF
  -    
  -     
  -    DO I=1,p%NReact  !For all constrained node
  - 
  -      p%MOutLst3(I)%Noutcnt=p%Reacts(I,1) !Assign nodeID for list I, I am using Noutcnt as a temporary holder for it, since nodeID is n array
  -      
  -      !Next stuff to be eliminated
  -      !ALLOCATE( p%MOutLst3(I)%NodeIDs(Init%Ndiv+1), STAT = ErrStat )  !
  -      ! IF ( ErrStat/= 0 ) THEN
  -      !   ErrStat = ErrID_Fatal
  -      !   ErrMsg  = 'Error allocating p%MOutLst3(I)%NodeIDs arrays in SDOut_Init'
  -      !   RETURN
  -      ! END IF
  -      ! 
  -      !p%MOutLst3(I)%NodeIDs=Init%MemberNodes(I,1:Init%Ndiv+1)  !We are storing  the actual node numbers in the member
  -      !Now I need to find out which elements are attached to those nodes and still belong to the member I
  -      !ElmID2s could contain the same element twice if Ndiv=1
  -      !p%MOutLst3(I)%ElmID2s=0  !Initialize to 0
  -      
  -          !I need to get ALL the elements that belong to the same joint
  -          !make use of MemberNodes and NodesConnE
  -          
  -      NconEls=Init%NodesConnE(p%MoutLst3(I)%Noutcnt,1) !Number of elements connecting to the joint
  -           
  -       ALLOCATE( p%MOutLst3(I)%ElmIDs(1,NconEls), STAT = ErrStat )  !element IDs connecting to the joint; (1,NconEls) and not (NconEls) as the same meshauxtype is used with other Moutlst
  -       IF ( ErrStat/= 0 ) THEN
  -         ErrStat = ErrID_Fatal
  -         ErrMsg  = 'Error allocating p%MOutLst3(I)%ElmIDs arrays in SDOut_Init'
  -         RETURN
  -       END IF
  -       
  -       ALLOCATE( p%MOutLst3(I)%ElmNds(1,NconEls), STAT = ErrStat )  !This array contains for each element connected to the joint, a flag that tells me whether the joint is node 1 or 2 for the elements
  -       IF ( ErrStat/= 0 ) THEN
  -         ErrStat = ErrID_Fatal
  -         ErrMsg  = 'Error allocating p%MOutLst3(I)%ElmIDs arrays in SDOut_Init'
  -         RETURN
  -       END IF
  -       
  -    ALLOCATE( p%MOutLst3(I)%Me(12,12,1,NconEls), STAT = ErrStat ) !Me has for each selected joint, and for each element attached to that node, a 12x12 matrix (extra dimension redundant)
  -   IF ( ErrStat/= 0 ) THEN
  -      ErrStat = ErrID_Fatal
  -      ErrMsg  = 'Error allocating p%MOutLst3(I)%Me arrays in SDOut_Init'
  -      RETURN
  -   END IF 
  -   ALLOCATE( p%MOutLst3(I)%Ke(12,12,1,NconEls), STAT = ErrStat ) !Ke has for each selected joint, and for each element attached to that node  a 12x12 matrix
  -   IF ( ErrStat/= 0 ) THEN
  -      ErrStat = ErrID_Fatal
  -      ErrMsg  = 'Error allocating p%MOutLst3(I)%Ke arrays in SDOut_Init'
  -      RETURN
  -   END IF 
  -    ALLOCATE( p%MOutLst3(I)%Fg(12,1,NconEls), STAT = ErrStat ) !Fg has for each selected joint, and for each element attached to that node  a 12x1 vector
  -   IF ( ErrStat/= 0 ) THEN
  -      ErrStat = ErrID_Fatal
  -      ErrMsg  = 'Error allocating p%MOutLst3(I)%Ke arrays in SDOut_Init'
  -      RETURN
  -    END IF 
  -    
  -         DO K=1, NconEls 
  -              L=Init%NodesConnE(p%MoutLst3(I)%Noutcnt,k+1)  !k-th Element Number in the set of elements attached to the selected node 
  -              
  -              p%MoutLst3(I)%ElmIDs(1,K)=L     !This array has for each joint requested  the elements' ID to get results for  
  -              
  -              M=p%Elems(L,2:3) !1st and 2nd node of the k-th element
  -             
  -              !Select whether the joint is the 1st or second node of the element
  -              p%MoutLst3(I)%ElmNds(1,K)=1                        !store whether first or second node of element  
  -              IF (M(2) .EQ. p%MoutLst3(I)%Noutcnt ) p%MoutLst3(I)%ElmNds(1,K)=2 !store whether first or second node of element  
  -             
  -              !Calculate Ke, Me to be used for output
  -              CALL ElemK_Beam( p%elemprops(L)%Area, p%elemprops(L)%Length, p%elemprops(L)%Ixx, p%elemprops(L)%Iyy, &
  -                              p%elemprops(L)%Jzz, p%elemprops(L)%Shear, p%elemprops(L)%kappa, p%elemprops(L)%YoungE,  & 
  -                              p%elemprops(L)%ShearG, p%elemprops(L)%DirCos, p%MoutLst3(I)%Ke(:,:,1,K) )
  -              CALL ElemM_Beam( p%elemprops(L)%Area, p%elemprops(L)%Length, p%elemprops(L)%Ixx, p%elemprops(L)%Iyy,&
  -                              p%elemprops(L)%Jzz,  p%elemprops(L)%rho,  p%elemprops(L)%DirCos, p%MoutLst3(I)%Me(:,:,1,K) )   
  -              CALL ElemG( p%elemprops(L)%Area, p%elemprops(L)%Length, p%elemprops(L)%rho,  p%elemprops(L)%DirCos, p%MoutLst3(I)%Fg(:,1,K), Init%g )   
  -              
  -             
  -        ENDDO
  -          
  -    ENDDO
  -
  -    !Store the matrix that will let me calculate single point reaction at the base of structure
  -    CALL ReactMatx(Init, p, WtrDpth, ErrStat, ErrMsg)
  - ENDIF
  -
  - 
  - 
  - ! These variables are to help follow the framework template, but the data in them is simply a copy of data
  - ! already available in the OutParam data structure
  -
  -  ALLOCATE ( InitOut%WriteOutputHdr(p%NumOuts+p%OutAllint*p%OutAllDims), STAT = ErrStat )
  -  IF ( ErrStat /= ErrID_None ) THEN
  -   ErrMsg  = ' Error allocating space for WriteOutputHdr array.'
  -   ErrStat = ErrID_Fatal
  +   ! These variables are to help follow the framework template, but the data in them is simply a copy of data
  +   ! already available in the OutParam data structure
  +   CALL AllocAry(InitOut%WriteOutputHdr, p%NumOuts+p%OutAllint*p%OutAllDims, 'WriteOutputHdr', ErrStat2, ErrMsg2); if(Failed()) return
  +   CALL AllocAry(InitOut%WriteOutputUnt, p%NumOuts+p%OutAllint*p%OutAllDims, 'WriteOutputUnt', ErrStat2, ErrMsg2); if(Failed()) return
  +   DO I = 1,p%NumOuts+p%OutAllint*p%OutAllDims
  +      InitOut%WriteOutputHdr(I) = TRIM( p%OutParam(I)%Name  )
  +      InitOut%WriteOutputUnt(I) = TRIM( p%OutParam(I)%Units )      
  +   END DO  
  +   
      RETURN
  -  END IF
  -
  -  ALLOCATE ( InitOut%WriteOutputUnt(p%NumOuts+p%OutAllint*p%OutAllDims), STAT = ErrStat )
  -  IF ( ErrStat /= ErrID_None ) THEN
  -   ErrMsg  = ' Error allocating space for WriteOutputHdr array.'
  -   ErrStat = ErrID_Fatal
  -   RETURN
  -  END IF
  -
  -   
  -   DO I = 1,p%NumOuts+p%OutAllint*p%OutAllDims
  -    InitOut%WriteOutputHdr(I) = TRIM( p%OutParam(I)%Name  )
  -    InitOut%WriteOutputUnt(I) = TRIM( p%OutParam(I)%Units )      
  -   END DO  
  -   
  -    
  -RETURN
   
   CONTAINS
      LOGICAL FUNCTION Failed()
  
positive.code:       
        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, z, -y/)
              
           CASE (2)
              p%TIreact(4:6, I) = (/-z, 0.0_ReKi, x/)
              
           CASE (3)
              p%TIreact(4:6, I) = (/y, -x, 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 '
              RETURN
           END SELECT
        
     ENDDO
     
     
  END SUBROUTINE ReactMatx
  
  !====================================================================================================
  
negative.code:              IF (M(1) .EQ. p%MoutLst2(I)%NodeIDs(J)) Junk=M(2)
               
               IF (ANY(Init%MemberNodes(p%MoutLst2(I)%MemberID,:) .EQ. Junk)) THEN  !This means we are in the selected member
                    p%MoutLst2(I)%ElmID2s(K2)=L     !This array has for each node requested NODEID(J), for each member I, the element to get results for 
                    p%MoutLst2(I)%ElmNd2s(K2)=1                        !store whether first or second node of element  
                    IF (M(2) .EQ. p%MoutLst2(I)%NodeIDs(J) ) p%MoutLst2(I)%ElmNd2s(K2)=2 !store whether first or second node of element  
                          
  