anchor.diff:  
      ! --- Additional Damping and stiffness at pin/ball/universal joints
      CALL InsertJointStiffDamp(p, Init, ErrStat2, ErrMsg2); if(Failed()) return
  +
  +
  +   ! Allocate the output matrices and the index mapping array
  +!    DOF_reduced = Init%TDOF-p%DOFCr   
  +!    CALL AllocAry(p%IDCr,  p%DOFCr, 'IDCR',  ErrStat2, ErrMsg2 )
  +!    CALL AllocAry(RetDOFs,  DOF_reduced, 'RetDOFs',  ErrStat2, ErrMsg2 )
  +!    !........
  +!    ! set the indices we want to keep (i.e., a mapping from reduced to full matrix)
  +!    J  = 1   
  +!    J2 = 1
  +!    DO I=1,Init%TDOF
  +!       idx(J) = idx(I)      
  +!       IF ( idx(J) /= 0 ) THEN 
  +!           J = J + 1   
  +!       ELSE 
  +!           p%IDCR(J2) = I
  +!           j2=j2+1
  +!       ENDIF
  +!                                 
  +!    END DO
  +!    !idx retains the rows (column) indices that are still active and has 0-padding (from elementTDOF-DOFCR to TDOF) to account for the removed ones, and everywhere else      
  +!    !idx2 retains the rows (column) indices that are removed (fixed DOFs, there are DOFCR of them)   
  +!    RetDOFs=idx(1:DOF_reduced)  
  +
   
      ! --------------------------------------------------------------------------------
      ! --- CB, Misc  
  
positive.code:    !    CBparams%PhiR from Eq. 3
     !    CBparams%MBB, CBparams%MBM, and CBparams%KBB from Eq. 4.
     !................................
     CALL CBMatrix(MRR, MLL, MRL, KRR, KLL, KRL, p%nDOFM, Init, &  ! < inputs
                   CBparams%MBB, CBparams%MBM, CBparams%KBB, CBparams%PhiL, CBparams%PhiR, CBparams%OmegaL, ErrStat2, ErrMsg2, p)  ! <- outputs (p is also input )
     ! TODO TODO TODO DAMPING MATRIX 
     if(Failed()) return
        
  
negative.code:    REAL(ReKi) , allocatable               :: Temp(:, :)        ! temp matrix for intermediate steps [bjj: made allocatable to try to avoid stack issues]
     REAL(ReKi) , allocatable               :: PhiR_T_MLL(:,:)   ! PhiR_T_MLL(p%nDOFR,p%nDOFL) = transpose of PhiR * MLL (temporary storage)
     INTEGER                                :: I !, lwork !counter, and varibales for inversion routines
     INTEGER                                :: ipiv(nL) !the integer vector ipvt of length min(m,n), containing the pivot indices. 
                                                         !Returned as: a one-dimensional array of (at least) length min(m,n), containing integers,
                                                         !where 1 <= less than or equal to ipvt(i) <= less than or equal to m.
     INTEGER(IntKi)                         :: ErrStat2                                                                    
     CHARACTER(ErrMsgLen)                   :: ErrMsg2
     CHARACTER(*), PARAMETER                :: RoutineName = 'CBMatrix'
     ErrStat = ErrID_None 
     ErrMsg  = ''
     
     if (nM_out>nL) then
        ErrMsg2='Cannot request more modes than internal degrees of Freedom'; ErrStat2=ErrID_Fatal; 
        if(Failed()) return;
     endif
     
     !....................................................
     ! Set OmegaL and PhiL from Eq. 2
     !....................................................
     IF ( nM_out > 0 ) THEN ! Only time this wouldn't happen is if no modes retained and no static improvement...
        ! bCheckSingularity = True
        CALL EigenSolveWrap(KLL, MLL, nL, nM_out, .True., PhiL(:,1:nM_out), OmegaL(1:nM_out),  ErrStat2, ErrMsg2); if(Failed()) return
  
        ! --- Normalize PhiL
        ! bjj: break up this equation to avoid as many temporary variables on the stack
        ! MU = MATMUL ( MATMUL( TRANSPOSE(PhiL), MLL ), PhiL )
        CALL AllocAry( Temp , nL , nL , 'Temp' , ErrStat2 , ErrMsg2); if(Failed()) return
        CALL AllocAry( MU   , nL , nL , 'Mu'   , ErrStat2 , ErrMsg2); if(Failed()) return
        MU   = TRANSPOSE(PhiL)
        Temp = MATMUL( MU, MLL )
        MU   = MATMUL( Temp, PhiL )
        DEALLOCATE(Temp)
        ! PhiL = MATMUL( PhiL, MU2 )  !this is the nondimensionalization (MU2 is diagonal)   
        DO I = 1, nM_out
           PhiL(:,I) = PhiL(:,I) / SQRT( MU(I, I) )
        ENDDO    
        DO I=nM_out+1, nL !loop done only if .not. p%SttcSolve .and. nDOFM < p%nDOFL (and actually, in that case, these values aren't used anywhere anyway)
           PhiL(:,I) = 0.0_ReKi
           OmegaL(I) = 0.0_ReKi
        END DO     
        DEALLOCATE(MU)
     END IF
        
     !....................................................
     ! Set PhiR from Eq. 3:
     !....................................................   
     ! now factor KLL to compute PhiR: KLL*PhiR=-TRANSPOSE(KRL)
     ! ** note this must be done after EigenSolveWrap() because it modifies KLL **
     CALL LAPACK_getrf( nL, nL, KLL, ipiv, ErrStat2, ErrMsg2); if(Failed()) return
     
     PhiR = -1.0_ReKi * TRANSPOSE(KRL) !set "b" in Ax=b  (solve KLL * PhiR = - TRANSPOSE( KRL ) for PhiR)
     CALL LAPACK_getrs( TRANS='N', N=nL, A=KLL, IPIV=ipiv, B=PhiR, ErrStat=ErrStat2, ErrMsg=ErrMsg2); if(Failed()) return
     
     !....................................................
     ! Set MBB, MBM, and KBB from Eq. 4:
     !....................................................
     CALL AllocAry( PhiR_T_MLL,  nR, nL, 'PhiR_T_MLL', ErrStat2, ErrMsg2); if(Failed()) return
        
     PhiR_T_MLL = TRANSPOSE(PhiR)
     PhiR_T_MLL = MATMUL(PhiR_T_MLL, MLL)
     MBB = MATMUL(MRL, PhiR)
     MBB = MRR + MBB + TRANSPOSE( MBB ) + MATMUL( PhiR_T_MLL, PhiR )
        
     IF ( nM == 0) THEN
        MBM = 0.0_ReKi
     ELSE
        MBM = MATMUL( PhiR_T_MLL, PhiL(:,1:nM))  ! last half of operation
        MBM = MATMUL( MRL, PhiL(:,1:nM) ) + MBM    !This had PhiM      
     ENDIF
     DEALLOCATE( PhiR_T_MLL )
     
  