subroutine Std_ori(natom,rx,ry,rz,a,b,c,rmass)
  !c
  !c      rotational constants for all replicas by inline diagonalization
  !c      of inertial tensor
  !c
!!!!!!! Input
! natom    : number of atoms
! rx,ry,rz : arrays natom with cartesian coordinates of natom particles
! rmass    : vector of length natom with particle masses
!
!!!!!!! Output
! a,b,c    : rotational constants a, b, and c
!
!!!!!! Dependencies
! tensor: calculate the inertial tensor
! Jacobi, trie : diaganalization subroutines

  implicit none

  INTEGER, PARAMETER :: KINT=KIND(1)
  INTEGER, PARAMETER :: KREAL=KIND(1.0D0)

  INTEGER(KINT), INTENT(IN) :: Natom
  REAL(KREAL), INTENT(INOUT) :: rx(Natom),ry(Natom),rz(Natom)
  REAL(KREAL), INTENT(IN) :: rmass(Natom)
  REAL(KREAL), INTENT(OUT) :: a,b,c

  real(KREAL) :: Tinert(3,3),D(3),V(3,3),CoordG(3)
  real(KREAL) :: ca,sa,na,cb,sb,nb

  INTEGER(KINT) :: I
  REAL(KREAL) :: Vt,Rxt,Ryt

  call tensor(natom,rx,ry,rz,Tinert,CoordG,rmass)
  call Jacobi(Tinert,3,D,V,3)
  call Trie(3,D,V,3)

!1000 FORMAT(3F15.3)
!1010 FORMAT(4F15.3)

  a=D(1)
!     WRITE(*,1010) a,(V(i,1),i=1,3)
  b=D(2)
!     WRITE(*,1010) b,(V(i,2),i=1,3)
  c=D(3)
!     WRITE(*,1010) c,(V(i,3),i=1,3)

!C On commence par mettre le centre de masse a l'origine

  Do I=1,natom
!     WRITE(*,*) 'DBG Tensor', I, rx(I),rY(I),rz(I)
     rx(i)=rx(i)-CoordG(1)
     ry(i)=ry(i)-CoordG(2)
     rz(i)=rz(i)-CoordG(3)
!     WRITE(*,*) 'DBG Tensor', I, rx(I),rY(I),rz(I)
  END DO

  ! C Ensuite, on va faire coincider les axes d'inertie avec x,y,z
  ! C Pour etre coherent, et bien que ce ne soit pas le plus logique,
  ! C on classe les vecteurs dans le sens croissant, et on fait coincider
  ! C les axes dans l'ordre x,y,z pour l'instant.

  ! C Coincidation du premier axe avec x
  ! C Rotation autour de z pour l'amener dans le plan xOz
  ! C puis Rotation autour de y pour l'amener suivant x
  na=sqrt(V(1,1)*V(1,1)+V(2,1)*V(2,1))
  nb=1
  if (na.gt.1D-10) then
     ca=V(1,1)/na
     sa=-V(2,1)/na
  else
     ca=1
     sa=0
  END IF
  cb=na/nb
  sb=-V(3,1)/nb
  DO I=1,natom
     rxt=rx(i)
     rx(i)=rx(i)*ca-ry(i)*sa
     ry(i)=rxt*sa+ry(i)*ca
     rxt=rx(i)
     rx(i)=rx(i)*cb-rz(i)*sb
     rz(i)=rxt*sb+rz(i)*cb
  END DO

  !C on applique aussi la rotation aux deux autres axes...
  Vt=V(1,2)
  V(1,2)=V(1,2)*ca-V(2,2)*sa
  V(2,2)=Vt*sa+V(2,2)*ca
  Vt=V(1,2)
  V(1,2)=V(1,2)*cb-V(3,2)*sb
  V(3,2)=Vt*sb+V(3,2)*cb

  !c     WRITE(IOOUT,1010) b,(V(i,2),i=1,3)

  !c      Vt=V(1,3)
  !c      V(1,3)=V(1,3)*ca-V(2,3)*sa
  !c      V(2,3)=Vt*sa+V(2,3)*ca
  !c      Vt=V(1,3)
  !c      V(1,3)=V(1,3)*cb-V(3,3)*sb
  !c      V(3,3)=Vt*sb+V(3,3)*cb
  !c     WRITE(IOOUT,1010) c,(V(i,3),i=1,3)

  !C et on repart pour le 2e axe sur y
  !C commes les vecteurs sont orthonormes, il est
  !C deja dans le plan yOz d'ou une seule rotation
  !c      na=sqrt(V(2,2)*V(2,2)+V(3,2)*V(3,2))
  na=1
  nb=1
  ca=V(2,2)/na
  sa=-V(3,2)/na
  !c      write(IOOUT,*) 'ca:',ca,sa,cb,sb,na
  DO I=1,natom
     ryt=ry(i)
     ry(i)=ry(i)*ca-rz(i)*sa
     rz(i)=ryt*sa+rz(i)*ca
  END DO
  !C on applique aussi la rotation au dernier axe...
  !c      Vt=V(2,3)
  !c      V(2,3)=V(2,3)*ca-V(3,3)*sa
  !c      V(3,3)=Vt*sa+V(3,3)*ca
  !c      WRITE(IOOUT,1010) c,(V(i,3),i=1,3)

  !C et le 3e axe est deja sur z car  ils sont orthogonaux...

  return
end subroutine Std_ori
