anchor.diff:             temp_id = (j-1)*p%dof_node
              TmpPos(1:3) = p%GlbPos(1:3) + MATMUL(p%GlbRot,p%uuN0(temp_id+1:temp_id+3,i))
              temp_CRV(:) = MATMUL(p%GlbRot,p%uuN0(temp_id+4:temp_id+6,i))
  -           CALL BD_CrvMatrixR(temp_CRV,TmpDCM)
  +           CALL BD_CrvMatrixR(temp_CRV,TmpDCM,ErrStat2,ErrMsg2)
              TmpDCM(:,:) = TRANSPOSE(TmpDCM)
              temp_id = (i-1)*p%node_elem+j
              CALL MeshPositionNode ( Mesh    = y%BldMotion  &
  
positive.code:            DistrLoad_GL(1:3,j)  = u%DistrLoad%Force(1:3,temp_id+j+1)
             DistrLoad_GL(4:6,j)  = u%DistrLoad%Moment(1:3,temp_id+j+1)
         ENDDO
         CALL BD_NodalRelRot(Nuu0,node_elem,dof_node,Nrr0)
         CALL BD_NodalRelRot(Nuuu,node_elem,dof_node,Nrrr)
         CALL BD_ElemNodalDisp(vvNf,node_elem,dof_node,nelem,Nvvv,ErrStat2,ErrMsg2)
         CALL BD_ElemNodalDisp(aaNf,node_elem,dof_node,nelem,Naaa,ErrStat2,ErrMsg2)
  
         CALL BD_ElementMatrixGA2(Nuu0,Nuuu,Nrr0,Nrrr,Nvvv,Naaa,           &
                                EStif0_GL,EMass0_GL,gravity,DistrLoad_GL,&
                                damp_flag,beta,                          &
                                ngp,node_elem,dof_node,elk,elf,elm,elg)
  
         CALL BD_AssembleStiffK(nelem,node_elem,dof_elem,dof_node,elk,StifK)
         CALL BD_AssembleStiffK(nelem,node_elem,dof_elem,dof_node,elm,MassM)
  
negative.code:    INTEGER(IntKi)               :: n
     INTEGER(IntKi)               :: temp_id1
     INTEGER(IntKi)               :: temp_id2
     INTEGER(IntKi)               :: ErrStat2                     ! Temporary Error status
     CHARACTER(ErrMsgLen)         :: ErrMsg2                      ! Temporary Error message
     CHARACTER(*), PARAMETER      :: RoutineName = 'BD_ElementMatrixGA2'
  
     ErrStat  = ErrID_None
     ErrMsg   = ""
     elk(:,:) = 0.0D0
     elf(:)   = 0.0D0
     elg(:,:) = 0.0D0
     elm(:,:) = 0.0D0
  
     CALL AllocAry(gp,ngp,'Gauss piont array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(gw,ngp,'Gauss piont weight function array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(hhx,node_elem,'Shape function array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(hpx,node_elem,'Derivative of shape function array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(GLL_temp,node_elem,'Gauss-Lobatto-Legendre (GLL) point array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(w_temp,node_elem,'GLL weight function array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     if (ErrStat >= AbortErrLev) then
        call Cleanup()
        return
     end if
     gp(:)       = 0.0D0
     gw(:)       = 0.0D0
     hhx(:)      = 0.0D0
     hpx(:)      = 0.0D0
     GLL_temp(:) = 0.0D0
     w_temp(:)   = 0.0D0
  
  
     CALL BD_GenerateGLL(node_elem-1,GLL_temp,w_temp,ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL BD_GaussPointWeight(ngp,gp,gw,ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     DO igp=1,ngp
         gpr = gp(igp)
         CALL BD_ComputeJacobian(gpr,Nuu0,node_elem,dof_node,gp,GLL_temp,ngp,igp,&
                                 hhx,hpx,Jacobian,ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         CALL BD_GaussPointDataAt0(hhx,hpx,Nuu0,Nrr0,node_elem,dof_node,&
                                   uu0,E10,ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         Stif(:,:) = 0.0D0
         Stif(1:6,1:6) = EStif0_GL(1:6,1:6,igp)
         CALL BD_GaussPointData(hhx,hpx,Nuuu,Nrrr,uu0,E10,node_elem,dof_node,&
                                uuu,uup,E1,RR0,kapa,Stif,cet,ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         CALL BD_ElasticForce(E1,RR0,kapa,Stif,cet,Fc,Fd,Oe,Pe,Qe,ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  
         mmm  = 0.0D0
         mEta = 0.0D0
  