anchor.diff:  MODULE SD_FEM
     USE NWTC_Library
     USE SubDyn_Types
  -  
     IMPLICIT NONE
  -  
  -   INTEGER,         PARAMETER  :: LAKi            = R8Ki                  ! Define the kind to be used for LAPACK routines for getting eigenvalues/vectors. Apparently there is a problem with SGGEV's eigenvectors
  -  
  +
  +  INTEGER,          PARAMETER  :: LAKi            = R8Ki                  ! Define the kind to be used for LAPACK routines for getting eigenvalues/vectors. Apparently there is a problem with SGGEV's eigenvectors
  + 
     INTEGER(IntKi),   PARAMETER  :: MaxMemjnt       = 10                    ! Maximum number of members at one joint
     INTEGER(IntKi),   PARAMETER  :: MaxOutChs       = 2000                  ! Max number of Output Channels to be read in
     INTEGER(IntKi),   PARAMETER  :: TPdofL          = 6                     ! 6 degrees of freedom (length of u subarray [UTP])
  
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 SUBROUTINE NodeCon
  
  !----------------------------------------------------------------------------
  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, J, n, Node, Node1, Node2, Prop, Prop1, Prop2   
     INTEGER                       :: OldJointIndex(Init%NJoints)
     INTEGER                       :: NNE      ! number of nodes per element
  