anchor.diff:     ENDDO
   
   
  +write(51,*) 'p%uu0'
      DO nelem = 1,p%elem_total
          DO idx_qp = 1,p%nqp
               !> ### Calculate the the initial displacement fields in an element
  
positive.code:    m%elf(:,:,:)      = 0.0_BDKi
     m%elm(:,:,:,:)    = 0.0_BDKi
  
     DO idx_qp=1,p%nqp
  
  
  
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
  
  