anchor.diff:         rho(1:3,1:3) = EMass0_GL(4:6,4:6,igp)
          CALL BD_GaussPointDataMass(hhx,hpx,Nvvv,Naaa,RR0,node_elem,dof_node,vvv,aaa,vvp,mmm,mEta,rho)
          CALL BD_ElasticForce(E1,RR0,kapa,Stif,cet,Fc,Fd,Oe,Pe,Qe)
  -!WRITE(*,*) 'Fc'
  -!WRITE(*,*) Fc
  -!WRITE(*,*) 'Fd'
  -!WRITE(*,*) Fd
  -!WRITE(*,*) 'aaa'
  -!WRITE(*,*) aaa
          IF(damp_flag .EQ. 1) THEN
              CALL BD_DissipativeForce(beta,Stif,vvv,vvp,E1,Fc,Fd,Sd,Od,Pd,Qd,betaC,Gd,Xd,Yd)
          ENDIF
  
positive.code:        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)
  !WRITE(*,*) 'uuu'
  !WRITE(*,*) uuu
  !WRITE(*,*) 'E1'
  !WRITE(*,*) E1
  !WRITE(*,*) 'kapa'
  !WRITE(*,*) kapa
         mmm  = 0.0D0
         mEta = 0.0D0
         rho  = 0.0D0
  
negative.code:        CALL BD_GravityForce(mmm,mEta,gravity,Fg)
         Fg(:) = Fg(:) + DistrLoad_GL(:,igp)
  
         DO i=1,node_elem
             DO j=1,dof_node
                 temp_id1 = (i-1) * dof_node+j
  