subroutine tensor(natom,rx,ry,rz,Tinert,CoordG,rmass)
  ! c
  ! c      elements of inertial tensor for a set of atomic systems.
  ! c
  ! c      rx,ry,rz : arrays natom with cartesian coordinates of natom particles; input
  ! c      txx,tyy,tzz,txy,txz,tyz: elements of inertial tensors ; output
  ! c      a,b,c    :  returning center of mass coordinates.
  ! c      rmass    : vector of length natom with particle masses; input
  ! c      natom    : number of particles per system; input
  ! c      IOOUT     : unit number for messages, silent for IOOUT.le.0; input
  ! c      subroutines called: cmshft                   m. lewerenz jul/92
  ! c

  use Io_module

  IMPLICIT NONE

  REAL(KREAL), parameter :: zero=0.d0,one=1.d0

  INTEGER(KINT), INTENT(IN) :: Natom
  real(KREAL) :: rx(Natom),ry(Natom),rz(Natom),rmass(Natom)
  real(KREAL) :: a,b,c,Tinert(3,3),CoordG(3)

  INTEGER(KINT) :: jover, I, J
  REAL(KREAL) :: Txx, Tyy, Tzz, Txy, Txz, Tyz, Txxi
  REAL(KREAL) :: Ryb, Rzc, Rxa, Rmj2, ryb2, rzc2, rxa2, Rmj


  if(natom.le.0) then
     if(IOOUT.gt.0) write(IOOUT,'(/2a/)')   &
          '  *** error, illegal array dimension(s) in tensor,', &
          ' return without action ***'
  else if(natom.eq.1) then
     a=rx(1)
     b=ry(1)
     c=rz(1)
     txx=zero
     tyy=zero
     tzz=zero
     txy=zero
     txz=zero
     tyz=zero
  else
     !c
     !c      first center of mass coordinates, then tensor elements
     !c
     call cmshft(natom,rx,ry,rz,rmass,a,b,c,0)
     CoordG(1)=a
     CoordG(2)=b
     CoordG(3)=c
1000 FORMAT(3F15.3)
     !c       WRITE(IOOUT,1000) (CoordG(i),i=1,3)
     jover=mod(natom,2)
     if(jover.eq.0) then
        txx=zero
        tyy=zero
        tzz=zero
        txy=zero
        txz=zero
        tyz=zero
     else
        rmj=rmass(1)
        ryb=b-ry(1)
        rzc=c-rz(1)
        tyz=-rmj*ryb*rzc
        rxa=a-rx(1)
        txz=(-rmj*rxa)*rzc
        txy=(-rmj*rxa)*ryb
        txx=rmj*((ryb*ryb)+(rzc*rzc))
        tyy=rmj*((rxa*rxa)+(rzc*rzc))
        tzz=rmj*((rxa*rxa)+(ryb*ryb))
     end if
!c
!c      force vectorization of inner loop if necessary (ibm-3090)
!c
     do j=jover+1,natom,2
        rmj=rmass(j)
        rmj2=rmass(j+1)
        ryb=b-ry(j)
        rzc=c-rz(j)
        txxi=txx+rmj*((ryb*ryb)+(rzc*rzc))
        ryb2=b-ry(j+1)
        rzc2=c-rz(j+1)
        txx=txxi+rmj2*((ryb2*ryb2)+(rzc2*rzc2))
        tyz=tyz-rmj*ryb*rzc-rmj2*ryb2*rzc2
        rxa=a-rx(j)
        rxa2=a-rx(j+1)
        txy=txy-(rmj*rxa)*ryb-(rmj2*rxa2)*ryb2
        txz=txz-(rmj*rxa)*rzc-(rmj2*rxa2)*rzc2
        tyy=tyy+rmj*((rxa*rxa)+(rzc*rzc))  &
             +rmj2*((rxa2*rxa2)+(rzc2*rzc2))
        tzz=tzz+rmj*((rxa*rxa)+(ryb*ryb))   &
             +rmj2*((rxa2*rxa2)+(ryb2*ryb2)) 
     END DO
  end if
  Tinert(1,1)=txx
  Tinert(1,2)=txy
  Tinert(1,3)=txz
  Tinert(2,1)=txy
  Tinert(2,2)=tyy
  Tinert(2,3)=tyz
  Tinert(3,1)=txz
  Tinert(3,2)=tyz
  Tinert(3,3)=tzz
  return
end subroutine tensor
