anchor.diff:  ! Allocate continuous states
      CALL AllocAry(x%q,p%dof_total,'x%q',ErrStat2,ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +   x%q(:) = 0.0D0
      CALL AllocAry(x%dqdt,p%dof_total,'x%dqdt',ErrStat2,ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +   x%dqdt(:) = 0.0D0
   ! Allocate other states: Acceleration and algorithm accelerations 
   ! for generalized-alpha time integator
      CALL AllocAry(OtherState%acc,p%dof_total,'OtherState%acc',ErrStat2,ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +   OtherState%acc(:) = 0.0D0
      CALL AllocAry(OtherState%xcc,p%dof_total,'OtherState%xcc',ErrStat2,ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +   OtherState%xcc(:) = 0.0D0
   
      ! Define system output initializations (set up mesh) here:
      CALL MeshCreate( BlankMesh        = u%RootMotion            &
  
positive.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 )
     ! Temporary Gauss point intrinsic coordinates array
     CALL AllocAry(temp_GL,p%ngp,'temp_GL',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     ! Temporary Gauss weight function array
     CALL AllocAry(temp_w,p%ngp,'GL weight array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL BD_GaussPointWeight(p%ngp,temp_GL,temp_w,ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     DEALLOCATE(temp_w)
  
negative.code: 
     CALL AllocAry(p%Stif0_GL,6,6,p%ngp*p%elem_total,'Stif0_GL',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     p%Stif0_GL(:,:,:) = 0.0D0
     CALL AllocAry(p%Mass0_GL,6,6,p%ngp*p%elem_total,'Mass0_GL',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     p%Mass0_GL(:,:,:) = 0.0D0
     DO i=1,p%elem_total
         DO j=1,p%ngp
             temp_id = (i-1)*p%ngp+j
  