anchor.diff:     REAL(ReKi)               :: r1, r2, t, Iyy, Jzz, Ixx, A, kappa, nu, ratioSq, D_inner, D_outer
      LOGICAL                  :: shear
      REAL(ReKi), ALLOCATABLE  :: Ke(:,:), Me(:, :), FGe(:) ! element stiffness and mass matrices gravity force vector
  -   INTEGER, ALLOCATABLE     :: nn(:)                     ! node number in element 
  +   INTEGER, DIMENSION(NNE)  :: nn                        ! node number in element 
      INTEGER                  :: r
  -   
  -   
      INTEGER(IntKi)           :: ErrStat2
      CHARACTER(1024)          :: ErrMsg2
  -
  -   
  -      ! for current application
  -   IF ( (Init%FEMMod .LE. 3) .and. (Init%FEMMod .GE. 0)) THEN
  -      NNE = 2
  -   ELSE
  -      ErrStat = ErrID_Fatal
  -      ErrMsg = 'Invalid FEMMod in AssembleKM'
  -      RETURN
  -   ENDIF                              
  +   
  +   ! for current application
  +   if    (Init%FEMMod == 2) THEN ! tapered Euler-Bernoulli
  +       CALL Abort ('FEMMod = 2 is not implemented.')
  +       return
  +   elseif (Init%FEMMod == 4) THEN ! tapered Timoshenko
  +       CALL Abort ('FEMMod = 2 is not implemented.')
  +       return
  +   elseif ((Init%FEMMod == 1) .or. (Init%FEMMod == 3)) THEN !
  +      ! 1: uniform Euler-Bernouli,  3: uniform Timoshenko
  +   else
  +       CALL Abort('FEMMod is not valid. Please choose from 1, 2, 3, and 4. ')
  +       return
  +   endif
      
      ! total degrees of freedom of the system 
      Init%TDOF = 6*Init%NNode
      
  -         ! Assemble system stiffness and mass matrices with gravity force vector
  -   
      ALLOCATE( p%ElemProps(Init%NElem), STAT=ErrStat2)
  -      IF (ErrStat2 /= 0) THEN
  -         CALL SetErrStat ( ErrID_Fatal, 'Error allocating p%ElemProps', ErrStat, ErrMsg, 'AssembleKM' )
  -         CALL CleanUp_AssembleKM()
  -      RETURN
  +   IF (ErrStat2 /= 0) THEN
  +       CALL Abort('Error allocating p%ElemProps')
  +       return
      ENDIF
   
  -   CALL AllocAry( Ke,     NNE*6,         NNE*6 , 'Ke',      ErrStat2, ErrMsg2); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'AssembleKM' )    ! element stiffness matrix
  -   CALL AllocAry( Me,     NNE*6,         NNE*6 , 'Me',      ErrStat2, ErrMsg2); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'AssembleKM' )    ! element mass matrix 
  -   CALL AllocAry( FGe,    NNE*6,                 'FGe',     ErrStat2, ErrMsg2); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'AssembleKM' )    ! element gravity force vector 
  -   CALL AllocAry( nn,     NNE,                   'nn',      ErrStat2, ErrMsg2); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'AssembleKM' )    ! node number in element array 
  -   
  -   CALL AllocAry( Init%K, Init%TDOF, Init%TDOF , 'Init%K',  ErrStat2, ErrMsg2); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'AssembleKM' )    ! system stiffness matrix 
  -   CALL AllocAry( Init%m, Init%TDOF, Init%TDOF , 'Init%M',  ErrStat2, ErrMsg2); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'AssembleKM' )    ! system mass matrix 
  -   CALL AllocAry( Init%FG,Init%TDOF,             'Init%FG', ErrStat2, ErrMsg2); CALL SetErrStat ( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'AssembleKM' )    ! system gravity force vector 
  -    
  -   IF (ErrStat >= AbortErrLev) THEN
  -      CALL CleanUp_AssembleKM()
  -      RETURN
  -   ENDIF
  -   
  -   Init%K = 0.0_ReKi
  -   Init%M = 0.0_ReKi
  +   CALL AllocAry( Ke,     NNE*6,         NNE*6 , 'Ke',      ErrStat2, ErrMsg2); if(Failed()) return; ! element stiffness matrix
  +   CALL AllocAry( Me,     NNE*6,         NNE*6 , 'Me',      ErrStat2, ErrMsg2); if(Failed()) return; ! element mass matrix 
  +   CALL AllocAry( FGe,    NNE*6,                 'FGe',     ErrStat2, ErrMsg2); if(Failed()) return; ! element gravity force vector 
  +   CALL AllocAry( Init%K, Init%TDOF, Init%TDOF , 'Init%K',  ErrStat2, ErrMsg2); if(Failed()) return; ! system stiffness matrix 
  +   CALL AllocAry( Init%m, Init%TDOF, Init%TDOF , 'Init%M',  ErrStat2, ErrMsg2); if(Failed()) return; ! system mass matrix 
  +   CALL AllocAry( Init%FG,Init%TDOF,             'Init%FG', ErrStat2, ErrMsg2); if(Failed()) return; ! system gravity force vector 
  +   Init%K  = 0.0_ReKi
  +   Init%M  = 0.0_ReKi
      Init%FG = 0.0_ReKi
   
      
  -      ! loop over all elements
  +   ! loop over all elements
      DO I = 1, Init%NElem
      
         DO J = 1, NNE
  
positive.code:    F(11) = -F(5)
     !F(12) is 0 for g along z alone
     
     
  END SUBROUTINE ElemG
  !------------------------------------------------------------------------------------------------------
  SUBROUTINE LumpForces(Area1,Area2,crat,L,rho, g, DirCos, F)
  !bjj: note this routine is a work in progress, intended for future version of SubDyn. Compare with ElemG.
  
           !This rountine calculates the lumped gravity forces at the nodes given the element geometry
           !It assumes a linear variation of the dimensions from node 1 to node 2, thus the area may be quadratically varying if crat<>1
     REAL(ReKi), INTENT( OUT)           :: F(12)
     REAL(ReKi), INTENT( IN )           :: Area1,Area2,crat !X-sectional areas at node 1 and node 2, t2/t1 thickness ratio
     REAL(ReKi), INTENT( IN )           :: g !gravity
     REAL(ReKi), INTENT( IN )           :: L !Length of element
     REAL(ReKi), INTENT( IN )           :: rho !density
     REAL(ReKi), INTENT( IN )           :: DirCos(3, 3)
  
      !LOCALS
     REAL(ReKi)                         :: TempCoeff,a0,a1,a2  !coefficients of the gravity quadratically distributed force
  
     A1 = -99999  !bjj initialized this to avoid getting warning by Intel Analyzers; needs to be set differently
     A2 = -99999  !bjj initialized this to avoid getting warning by Intel Analyzers; needs to be set differently
     
     !Calculate quadratic polynomial coefficients
     a0=A1
     a2=( (Area1+A2) - (Area1*crat+Area2/crat) )/L**2.  !*x**2
     a1= (Area2-Area1)/L -a2*L                       !*x                   
     
     !Now calculate the Lumped Forces
     F = 0
  
negative.code:             END IF
              J = J + 1
           END DO 
           IF ( .NOT. found) THEN
              CALL Abort(' Member '//TRIM(Num2LStr(I))//' has JointID'//TRIM(Num2LStr(n-1))//' = '// TRIM(Num2LStr(Node))//' which is not in the node list !')
              RETURN
           END IF
        END DO ! loop through nodes/joints
        
        ! loop through the PropSetIDs for this member and find the corresponding indices into the Joints array
        ! we're setting these two values:   
        ! p%Elems(I, 4) = property set for node 1 (note this sets the YoungE, ShearG, and MatDens columns for the ENTIRE element)   
        ! p%Elems(I, 5) = property set for node 2 (note this should be used only for the XsecD and XsecT properties in the element [for a linear distribution from node 1 to node 2 of D and T])
        DO n=4,5 ! Member column for MPropSetID1 and MPropSetID2
           Prop = Init%Members(I, n)  ! n=4 or 5
           ! ...... search for index of property set whose PropSetID matches Prop ......
           J = 1
           found = .false.      
           DO WHILE ( .NOT. found .AND. J <= Init%NPropSets )
  