anchor.diff:         mmm  = 0.0D0
          mmm  = EMass0_GL(1,1,igp)
   
  +WRITE(*,*) 'igp',igp
  +WRITE(*,*) 'gw',gw(igp)
  +WRITE(*,*) 'Jaco',Jacobian
  +WRITE(*,*) 'mmm',mmm
          elem_mass = elem_mass + gw(igp) * Jacobian * mmm
          elem_CG(1:3) = elem_CG(1:3) + gw(igp) * Jacobian * mmm * NGPpos(1:3,igp)
          elem_IN(1:3,1:3) = elem_IN(1:3,1:3) + gw(igp) * Jacobian * mmm * &
  
positive.code:        end if
  
     ENDDO
  
     blade_CG(:) = blade_CG(:) / blade_mass
  
     CALL Cleanup()
  
negative.code:            ENDDO
         ENDIF
  
  WRITE(*,*) 'nelem:',nelem
         CALL BD_ComputeElementMass(Nuu0,NGPpos,EMass0_GL,&
                                    nqp,quadrature,trapezoidal_pos,trapezoidal_w,&
                                    node_elem,dof_node,&
  