anchor.diff:                                EStif0_GL,EMass0_GL,gravity,DistrLoad_GL,&
                                 damp_flag,beta,                          &
                                 ngp,node_elem,dof_node,elk,elf,elm,elg)
  -!WRITE(*,*) 'elf:',nelem
  -!DO j=1,18
  -!   WRITE(*,*) elf(j)
  -!ENDDO
   
          CALL BD_AssembleStiffK(nelem,node_elem,dof_elem,dof_node,elk,StifK)
          CALL BD_AssembleStiffK(nelem,node_elem,dof_elem,dof_node,elm,MassM)
          CALL BD_AssembleStiffK(nelem,node_elem,dof_elem,dof_node,elg,DampG)
          CALL BD_AssembleRHS(nelem,dof_elem,node_elem,dof_node,elf,RHS)
      ENDDO
  -!WRITE(*,*) 'RHS:'
  -!DO j=1,30
  -!   WRITE(*,*) RHS(j)
  -!ENDDO
   
      DEALLOCATE(Nuu0)
      DEALLOCATE(Nuuu)
  
positive.code:    temp_L2(1:3,1) = p%uuN0(1:3,1)
     temp_L2(1:3,p%ngp*p%elem_total+2) = p%uuN0(temp_int-5:temp_int-3,p%elem_total)
     DEALLOCATE(temp_GLL)
  
  ! Compute sectional propertities ( 6 by 6 stiffness and mass matrices)
  ! at Gauss points
     CALL AllocAry(temp_ratio,p%ngp,p%elem_total,'temp_ratio',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     DO i=1,p%ngp
         temp_GL(i) = (temp_GL(i) + 1.0D0)/2.0D0
     ENDDO
  
     DO i=1,p%elem_total
         IF(i .EQ. 1) THEN
             DO j=1,p%ngp
  
negative.code:    ! temp_L2: physical coordinates of Gauss points and two end points
     CALL AllocAry(temp_L2,3,p%ngp*p%elem_total+2,'temp_L2',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     temp_L2(:,:) = 0.0D0
     ! Temporary Gauss point intrinsic coordinates array
     CALL AllocAry(temp_GL,p%ngp,'temp_GL',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     temp_GL(:) = 0.0D0
     ! Temporary Gauss weight function array
     CALL AllocAry(temp_w,p%ngp,'GL weight array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     temp_w(:) = 0.0D0
     CALL BD_GaussPointWeight(p%ngp,temp_GL,temp_w,ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     DEALLOCATE(temp_w)
  