anchor.diff:        ! Column 6: member type
         p%Elems(iMem, iMType) = Init%Members(iMem, iMType) ! 
      END DO !iMem, loop through members
  +    
  +   ! TODO in theory, we shouldn't need these anymore
  +   ! deallocate(Init%Members)
  +   ! deallocate(Init%Joints)
  +CONTAINS
  +   LOGICAL FUNCTION Failed()
  +      call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SD_ReIndex_CreateNodesAndElems') 
  +      Failed =  ErrStat >= AbortErrLev
  +   END FUNCTION Failed
  +   SUBROUTINE Fatal(ErrMsg_in)
  +      CHARACTER(len=*), intent(in) :: ErrMsg_in
  +      CALL SetErrStat(ErrID_Fatal, ErrMsg_in, ErrStat, ErrMsg, 'SD_ReIndex_CreateNodesAndElems');
  +   END SUBROUTINE Fatal
  +END SUBROUTINE SD_ReIndex_CreateNodesAndElems
  +
  +!----------------------------------------------------------------------------
  +SUBROUTINE SD_Discrt(Init,p, ErrStat, ErrMsg)
  +   TYPE(SD_InitType),            INTENT(INOUT)  ::Init
  +   TYPE(SD_ParameterType),       INTENT(INOUT)  ::p
  +   INTEGER(IntKi),               INTENT(  OUT)  :: ErrStat     ! Error status of the operation
  +   CHARACTER(*),                 INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
  +   ! local variable
  +   INTEGER                       :: I, iMem, J, n, iNode, Node, Node1, Node2, Prop, Prop1, Prop2   
  +   INTEGER                       :: NNE      ! number of nodes per element
  +   INTEGER                       :: MaxNProp
  +   REAL(ReKi), ALLOCATABLE       :: TempPropsBeams(:, :)
  +   REAL(ReKi), ALLOCATABLE       :: TempPropsCable(:, :)
  +   REAL(ReKi), ALLOCATABLE       :: TempPropsRigid(:, :)
  +   INTEGER, ALLOCATABLE          :: TempMembers(:, :) ,TempReacts(:,:)         
  +   INTEGER                       :: knode, kelem, kprop, nprop
  +   REAL(ReKi)                    :: x1, y1, z1, x2, y2, z2, dx, dy, dz, dd, dt, d1, d2, t1, t2
  +   LOGICAL                       :: found, CreateNewProp
  +   INTEGER(IntKi)                :: eType !< Element Type
  +   INTEGER(IntKi)                :: mType !< Member Type
  +   CHARACTER(1024)               :: sType !< String for element type
  +   INTEGER(IntKi)                :: ErrStat2
  +   CHARACTER(1024)               :: ErrMsg2
  +   ErrStat = ErrID_None
  +   ErrMsg  = ""
  +   
  +   ! number of nodes per element
  +   IF( ( Init%FEMMod >= 0 ) .and. (Init%FEMMod <= 3) ) THEN
  +      NNE = 2 
  +   ELSE
  +      CALL Fatal('FEMMod '//TRIM(Num2LStr(Init%FEMMod))//' not implemented.')
  +      RETURN
  +   ENDIF
  +   
  +   Init%NNode = Init%NJoints + ( Init%NDiv - 1 )*p%NMembers    ! Calculate total number of nodes according to divisions 
  +   Init%NElem = p%NMembers*Init%NDiv                           ! Total number of element   
  +   MaxNProp   = Init%NPropSets + Init%NElem*NNE                ! Maximum possible number of property sets (temp): This is property set per element node, for all elements (bjj, added Init%NPropSets to account for possibility of entering many unused prop sets)
  +   
  +   ! Calculate total number of nodes and elements according to element types
  +   ! for 3-node or 4-node beam elements
  +   Init%NNode = Init%NNode + (NNE - 2)*Init%NElem
  +   !bjj: replaced with max value instead of NNE: Init%MembersCol = Init%MembersCol + (NNE - 2) 
  +   
  +   ! check the number of interior modes
  +   IF ( p%Nmodes > 6*(Init%NNode - Init%NInterf - p%NReact) ) THEN
  +      CALL Fatal(' NModes must be less than or equal to '//TRIM(Num2LStr( 6*(Init%NNode - Init%NInterf - p%NReact) )))
  +      RETURN
  +   ENDIF
  +   
  +   CALL AllocAry(Init%MemberNodes,p%NMembers,    Init%NDiv+1,'Init%MemberNodes',ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'SD_Discrt')  ! for two-node element only, otherwise the number of nodes in one element is different
  +   CALL AllocAry(Init%BCs,        6*p%NReact,    2,          'Init%BCs',        ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'SD_Discrt') !!!! RRD: THIS MAY NEED TO CHANGE IF NOT ALL NODES ARE RESTRAINED
  +   CALL AllocAry(Init%IntFc,      6*Init%NInterf,2,          'Init%IntFc',      ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'SD_Discrt')
  +   
  +   CALL AllocAry(TempMembers,     p%NMembers,    MembersCol, 'TempMembers',     ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'SD_Discrt') 
  +   CALL AllocAry(TempPropsBeams,   MaxNProp,      PropSetsCol     ,'TempPropsBeams',       ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'SD_Discrt') 
  +   CALL AllocAry(TempPropsCable,  MaxNProp,      CablePropSetsCol,'TempPropsCable',       ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'SD_Discrt') 
  +   CALL AllocAry(TempPropsRigid,  MaxNProp,      RigidPropSetsCol,'TempPropsRigid',       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()
  +      RETURN
  +   ENDIF
  +
  +   ! ---
  +   CALL SD_ReIndex_CreateNodesAndElems(Init,p, ErrStat2, ErrMsg2);  if(Failed()) return
      
      ! Initialize TempMembers
      TempMembers = p%Elems(1:p%NMembers,:)
      
      ! Initialize Temp property set, first user defined sets
  -   TempPropsBeam = 0
  -   TempPropsBeam(1:Init%NPropSets, :) = Init%PropSets   
  -   print*,'TempPropsBeam(:,1:5)'
  +   TempPropsBeams = 0
  +   TempPropsCable = 0
  +   TempPropsRigid = 0
  +   TempPropsBeams(1:Init%NPropSets, :)      = Init%PropSets   
  +   TempPropsCable(1:Init%NCablePropSets, :) = Init%CablePropSets
  +   TempPropsRigid(1:Init%NRigidPropSets, :) = Init%RigidPropSets
  +   print*,'TempPropsBeams(:,1:5)'
      DO I=1,Init%NPropSets
  -      print*,TempPropsBeam(I,1:5)
  +      print*,TempPropsBeams(I,1:5)
      enddo
  -
  -
  -   TempPropsCable = 0
  -   TempPropsCable(1:Init%NCablePropSets, :) = Init%CablePropSets
      print*,'TempPropsCable(:,:)'
      DO I=1,Init%NCablePropSets
         print*,TempPropsCable(I,1:4)
      enddo
   
  -   TempPropsRigid = 0
  -   TempPropsRigid(1:Init%NRigidPropSets, :) = Init%RigidPropSets
      print*,'TempPropsRigid(:,:)'
      DO I=1,Init%NRigidPropSets
         print*,TempPropsRigid(I,:)
  
positive.code: 
     SUBROUTINE CleanUp_Discrt()
        ! deallocate temp matrices
        IF (ALLOCATED(TempPropsBeam))    DEALLOCATE(TempPropsBeam)
        IF (ALLOCATED(TempPropsCable))   DEALLOCATE(TempPropsCable)
        IF (ALLOCATED(TempPropsRigid))   DEALLOCATE(TempPropsRigid)
        IF (ALLOCATED(TempMembers)) DEALLOCATE(TempMembers)
  
negative.code:           Init%MemberNodes(I,           1) = Node1
            Init%MemberNodes(I, Init%NDiv+1) = Node2
            
            IF  ( ( .not. EqualRealNos(TempPropsBeams(Prop1, 2),TempPropsBeams(Prop2, 2) ) ) &
             .OR. ( .not. EqualRealNos(TempPropsBeams(Prop1, 3),TempPropsBeams(Prop2, 3) ) ) &
             .OR. ( .not. EqualRealNos(TempPropsBeams(Prop1, 4),TempPropsBeams(Prop2, 4) ) ) )  THEN
            
               CALL Fatal(' Material E,G and rho in a member must be the same')
               RETURN
  