anchor.diff:     !
      call CleanUp()
   
  +   print*,'DOF_I  ',p%IDI
  +   print*,'DOF_L  ',p%IDL
  +   print*,'DOF_C  ',p%IDC
  +   print*,'Nodes_C',p%Nodes_C(:,1)
  +   print*,'Nodes_L',p%Nodes_L
      print*,'Nodes_I',p%Nodes_I(:,1)
      print*,'Nodes_C',p%Nodes_C(:,1)
      print*,'Nodes_L',p%Nodes_L
  
positive.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                                :: DOFvar !placeholder used to get both PhiL or PhiM into 1 process
     INTEGER                                :: ipiv(p%nDOFL) !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  = ''
     
     CALL WrScr('   Calculating Internal Modal Eigenvectors')
          
     IF (p%SttcSolve) THEN ! STATIC TREATMENT IMPROVEMENT
        DOFvar=p%nDOFL
     ELSE
        DOFvar=nDOFM !Initialize for normal cases, dynamic only      
     ENDIF  
     
     !....................................................
     ! Set OmegaL and PhiL from Eq. 2
     !....................................................
     IF ( DOFvar > 0 ) THEN ! Only time this wouldn't happen is if no modes retained and no static improvement...
        ! NOTE: 
        !   bRemoveConstraint = False
        !   bCheckSingularity = True
        CALL EigenSolveWrap(KLL, MLL, p%nDOFL, DOFvar, .True., PhiL(:,1:DOFvar), OmegaL(1:DOFvar),  ErrStat2, ErrMsg2); if(Failed()) return
  
        ! --- Normalize PhiL
        ! bjj: break up this equation to avoid as many tenporary variables on the stack
        ! MU = MATMUL ( MATMUL( TRANSPOSE(PhiL), MLL ), PhiL )
        CALL AllocAry( Temp , p%nDOFL , p%nDOFL , 'Temp' , ErrStat2 , ErrMsg2); if(Failed()) return
        CALL AllocAry( MU   , p%nDOFL , p%nDOFL , '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, DOFvar
           PhiL(:,I) = PhiL(:,I) / SQRT( MU(I, I) )
        ENDDO    
        DO I=DOFvar+1, p%nDOFL !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)
        
        !....................................................
        ! Set p%PhiL_T and p%PhiLInvOmgL2 for static improvement
        !....................................................
        IF (p%SttcSolve) THEN   
           p%PhiL_T=TRANSPOSE(PhiL) !transpose of PhiL for static improvement
           DO I = 1, p%nDOFL
              p%PhiLInvOmgL2(:,I) = PhiL(:,I)* (1./OmegaL(I)**2)
           ENDDO 
        END IF
        
     ! ELSE .not. p%SttcSolve .and. nDOFM < p%nDOFL (in this case, PhiL, OmegaL aren't used)      
     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( p%nDOFL, p%nDOFL, 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=p%nDOFL,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,  p%nDOFR, p%nDOFL, '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 ( nDOFM .EQ. 0) THEN
        MBM = 0.0_ReKi
     ELSE
        MBM = MATMUL( PhiR_T_MLL, PhiL(:,1:nDOFM))  ! last half of operation
        MBM = MATMUL( MRL, PhiL(:,1:nDOFM) ) + MBM    !This had PhiM      
     ENDIF
     DEALLOCATE( PhiR_T_MLL )
     
  
negative.code:    !    CBparams%PhiR from Eq. 3
     !    CBparams%MBB, CBparams%MBM, and CBparams%KBB from Eq. 4.
     !................................
     CALL WrScr('   Performing Craig-Bampton decomposition')
     IF (p%SttcSolve) THEN ! STATIC TREATMENT IMPROVEMENT
        nM_Out=p%nDOFL ! Selecting all CB modes for outputs to the function below
     ELSE
        nM_Out=p%nDOFM ! Selecting only the requrested number of CB modes
     ENDIF  
     CALL CBMatrix(MRR, MLL, MRL, KRR, KLL, KRL, p%nDOFR, p%nDOFL, p%nDOFM, nM_Out, &  ! < inputs
                   CBparams%MBB, CBparams%MBM, CBparams%KBB, CBparams%PhiL, CBparams%PhiR, CBparams%OmegaL, ErrStat2, ErrMsg2)  ! <- outputs
  
     ! Set p%PhiL_T and p%PhiLInvOmgL2 for static improvement
     IF (p%SttcSolve) THEN   
        p%PhiL_T=TRANSPOSE(CBparams%PhiL) !transpose of PhiL for static improvement
        DO I = 1, nM_Out
           p%PhiLInvOmgL2(:,I) = CBparams%PhiL(:,I)* (1./CBparams%OmegaL(I)**2)
        ENDDO 
     END IF
  
     ! TODO TODO TODO DAMPING MATRIX 
     if(Failed()) return
        
  