anchor.diff:              !! (tangent to curve through this GLL point).  This is simply the
            CALL BD_CrvMatrixR(p%uu0(4:6,idx_qp,nelem),R0_temp)         ! returns R0_temp (the transpose of the DCM orientation matrix)
            p%E10(:,idx_qp,nelem) = R0_temp(:,3)                        ! unit vector tangent to curve through this GLL point (derivative with respect to z in IEC coords).
  +write(51,'(6(2x,ES20.12E3))') p%uu0(:,idx_qp,nelem)
         ENDDO
      ENDDO 
   
  +
   END SUBROUTINE BD_QuadraturePointDataAt0
   
   
  
positive.code:             !!
              !! Note: \f$ \underline{\underline{R}}\underline{\underline{R}}_0 \f$ is used to go from the material basis into the inertial basis
              !!       and the transpose for the other direction.
           rot0_T0(:) = p%uu0(4:6,idx_qp,nelem)                ! initial rotation parameter set at gauss point at T=0
        
           CALL BD_CrvCompose(cc,rot_c_xi,rot0_T0,FLAG_R1R2)   ! cc = rot_c_xi composed with rot0_T0
  
negative.code:       CALL BD_GyroForce( nelem, idx_qp, m )
        CALL BD_GravityForce( nelem, idx_qp, p, m, p%gravity )
        CALL BD_ElasticForce( nelem,idx_qp,p,m,.FALSE. )     ! Calculate Fc, Fd only
  write(53,'(6(2x,ES20.12E3))') m%qp%Fc(:,idx_qp,nelem) 
  write(54,'(6(2x,ES20.12E3))') m%qp%Fd(:,idx_qp,nelem) 
        IF(p%damp_flag .NE. 0) THEN
           CALL BD_DissipativeForce( nelem, idx_qp, p,m,.FALSE. )         ! Calculate dissipative terms on Fc, Fd
  
  