anchor.diff:     CALL AllocAry(TempMembers,     p%NMembers,    MembersCol, 'TempMembers',     ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'SD_Discrt') 
      CALL AllocAry(TempProps,       MaxNProp,      PropSetsCol,'TempProps',       ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'SD_Discrt') 
      CALL AllocAry(TempReacts,      p%NReact,      ReactCol,   'TempReacts',      ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'SD_Discrt')
  -   
   
      IF ( ErrStat >= AbortErrLev ) THEN
         CALL CleanUp_Discrt()
  
positive.code:       y2  = Init%Nodes(N2, 3)
        z2  = Init%Nodes(N2, 4)
  
        CALL GetDirCos(X1, Y1, Z1, X2, Y2, Z2, DirCos, L, ErrStat2, ErrMsg2)
           CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'AssembleKM' )
           IF (ErrStat >= AbortErrLev) THEN
              CALL CleanUp_AssembleKM()
              RETURN
           END IF
        
           
  !BJJ: TODO: for efficiency, this if check should be OUTSIDE the DO loop. 
        ! 1: uniform Euler-Bernouli
        ! 3: uniform Timoshenko
        IF ( (Init%FEMMod == 1).OR.(Init%FEMMod == 3)) THEN ! uniform element 
           r1 = 0.25*(D1 + D2)
           t  = 0.5*(t1+t2)
           
           IF ( EqualRealNos(t, 0.0_ReKi) ) THEN
              r2 = 0
           ELSE
              r2 = r1 - t
           ENDIF
           
           A = Pi_D*(r1*r1-r2*r2)
           Ixx = 0.25*Pi_D*(r1**4-r2**4)
           Iyy = Ixx
           Jzz = 2.0*Ixx
           
           IF( Init%FEMMod == 1 ) THEN ! uniform Euler-Bernoulli
              Shear = .false.
              kappa = 0
           ELSEIF( Init%FEMMod == 3 ) THEN ! uniform Timoshenko
              Shear = .true.
            ! kappa = 0.53            
              
                 ! equation 13 (Steinboeck et al) in SubDyn Theory Manual 
              nu = E / (2.0_ReKi*G) - 1.0_ReKi
              D_outer = 2.0_ReKi * r1  ! average (outer) diameter
              D_inner = D_outer - 2*t  ! remove 2x thickness to get inner diameter
              ratioSq = ( D_inner / D_outer)**2
              kappa =   ( 6.0 * (1.0 + nu) **2 * (1.0 + ratioSq)**2 ) &
                      / ( ( 1.0 + ratioSq )**2 * ( 7.0 + 14.0*nu + 8.0*nu**2 ) + 4.0 * ratioSq * ( 5.0 + 10.0*nu + 4.0 *nu**2 ) )
                          
           ENDIF
           
           p%ElemProps(i)%Area = A
           p%ElemProps(i)%Length = L
           p%ElemProps(i)%Ixx = Ixx
           p%ElemProps(i)%Iyy = Iyy
           p%ElemProps(i)%Jzz = Jzz
           p%ElemProps(i)%Shear = Shear
           p%ElemProps(i)%kappa = kappa
           p%ElemProps(i)%YoungE = E
           p%ElemProps(i)%ShearG = G
           p%ElemProps(i)%Rho = rho
           p%ElemProps(i)%DirCos = DirCos
           
           
           CALL ElemK(A, L, Ixx, Iyy, Jzz, Shear, kappa, E, G, DirCos, Ke)
           CALL ElemM(A, L, Ixx, Iyy, Jzz, rho, DirCos, Me)
           CALL ElemG(A, L, rho, DirCos, FGe, Init%g)
           
        ELSEIF  (Init%FEMMod == 2) THEN ! tapered Euler-Bernoulli
           CALL SetErrStat ( ErrID_Fatal, 'FEMMod = 2 is not implemented.', ErrStat, ErrMsg, 'AssembleKM' )
           CALL CleanUp_AssembleKM()
           RETURN
           
        ELSEIF  (Init%FEMMod == 4) THEN ! tapered Timoshenko
           CALL SetErrStat ( ErrID_Fatal, 'FEMMod = 4 is not implemented.', ErrStat, ErrMsg, 'AssembleKM' )
           CALL CleanUp_AssembleKM()
           RETURN
           
        ELSE
           CALL SetErrStat ( ErrID_Fatal, 'FEMMod is not valid. Please choose from 1, 2, 3, and 4. ', ErrStat, ErrMsg, 'AssembleKM' )
           CALL CleanUp_AssembleKM()
           RETURN
           
        ENDIF  
                                                                                                                                                                 
  
        
        ! assemble element matrices to global matrices
           
        DO J = 1, NNE
           jn = nn(j)
           
           Init%FG( (jn*6-5):(jn*6) ) = Init%FG( (jn*6-5):(jn*6) ) &
                                      + FGe( (J*6-5):(J*6) )
           
           DO K = 1, NNE
              kn = nn(k)
              
              Init%K( (jn*6-5):(jn*6), (kn*6-5):(kn*6) ) = Init%K( (jn*6-5):(jn*6), (kn*6-5):(kn*6) ) &
                                                    + Ke( (J*6-5):(J*6), (K*6-5):(K*6) )
                    
              Init%M( (jn*6-5):(jn*6), (kn*6-5):(kn*6) ) = Init%M( (jn*6-5):(jn*6), (kn*6-5):(kn*6) ) &
                                                    + Me( (J*6-5):(J*6), (K*6-5):(K*6) )
                       
                       
           ENDDO !K
                       
        ENDDO !J
                       
                       
     ENDDO ! I end loop over elements
                 
        
        ! add concentrated mass 
     DO I = 1, Init%NCMass
        DO J = 1, 3
            r = ( NINT(Init%CMass(I, 1)) - 1 )*6 + J
            Init%M(r, r) = Init%M(r, r) + Init%CMass(I, 2)
            
        ENDDO
        DO J = 4, 6
            r = ( NINT(Init%CMass(I, 1)) - 1 )*6 + J
            Init%M(r, r) = Init%M(r, r) + Init%CMass(I, J-1)
        ENDDO
  
     ENDDO ! I concentrated mass
   
        ! add concentrated mass induced gravity force
     DO I = 1, Init%NCMass
        
        r = ( NINT(Init%CMass(I, 1)) - 1 )*6 + 3
        Init%FG(r) = Init%FG(r) - Init%CMass(I, 2)*Init%g 
  
     ENDDO ! I concentrated mass induced gravity
     
     CALL CleanUp_AssembleKM()
     RETURN
     
  CONTAINS
  !..............
     SUBROUTINE CleanUp_AssembleKM()
  ! deallocate temp matrices
        IF (ALLOCATED(Ke)) DEALLOCATE(Ke)
        IF (ALLOCATED(Me)) DEALLOCATE(Me)
        IF (ALLOCATED(FGe)) DEALLOCATE(FGe)
        IF (ALLOCATED(nn)) DEALLOCATE(nn)
     END SUBROUTINE CleanUp_AssembleKM
     
  END SUBROUTINE AssembleKM
  !------------------------------------------------------------------------------------------------------
  !------------------------------------------------------------------------------------------------------
  
  SUBROUTINE GetDirCos(X1, Y1, Z1, X2, Y2, Z2, DirCos, L, ErrStat, ErrMsg)
     !This should be from local to global -RRD
     ! bjj: note that this is the transpose of what is normally considered the Direction Cosine Matrix  
     !      in the FAST framework. It seems to be used consistantly in the code (i.e., the transpose 
     !      of this matrix is used later).
     
     
     REAL(ReKi) ,      INTENT(IN   )  :: x1, y1, z1, x2, y2, z2  ! (x,y,z) positions of two nodes making up an element
     REAL(ReKi) ,      INTENT(  OUT)  :: DirCos(3, 3)            ! calculated direction cosine matrix
     REAL(ReKi) ,      INTENT(  OUT)  :: L                       ! length of element
     
     INTEGER(IntKi),   INTENT(  OUT)  :: ErrStat                 ! Error status of the operation
     CHARACTER(*),     INTENT(  OUT)  :: ErrMsg                  ! Error message if ErrStat /= ErrID_None
     
     REAL(ReKi)                       ::  Dx,Dy,Dz, Dxy          ! distances between nodes
  !real(reki) :: dxyz         
     ErrMsg  = ""
     ErrStat = ErrID_None
     
  
negative.code:             EXIT  !Exit J loop if node found -RRD
           ENDIF
        ENDDO
        IF (.not. found) THEN
           CALL Abort(' React has node not in the node list !')
           RETURN
        ENDIF
        DO J = 1, 6
           Init%BCs( (I-1)*6+J, 1) = (Node2-1)*6+J;
           Init%BCs( (I-1)*6+J, 2) = p%Reacts(I, J+1);
        ENDDO
     ENDDO
     p%Reacts=TempReacts   !UPDATED REACTS
        
     ! Initialize interface constraint vector
     ! Change the node number
     Init%IntFc = 0
     DO I = 1, Init%NInterf
        Node1 = Init%Interf(I, 1);
  