C  7 Apr 06 - OQ  - changes to allow spherical harmonics
c 17 Jan 06 - MWS - GetAlphaAnisotropy,InitHyper/RamanData: math ops
C  5 Jul 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  5 Feb 05 - OQ  - I/O support for extended TDHF properties
C
C     The mu,alpha,beta and gamma tensors are described as vectors of
C     size 3, 9, 27 and 21, respectively.
C     They should be type in Fortran 90 but since it's Fortran 77...
C
C     The descriptions of each tensor element in terms of tensor indexes
C     are done throught 'IELEMM', 'IELEMA', 'IELEMB' and 'IELEMG' while
C     the reverse operations are done throught 'INDEXM', 'INDEXA',
C     'INDEXB' and 'INDEXG', respectively.  Those tables should not be
C     used directly!
C     The tensor elements should be accessed indirectly via
C         'GetMuElem', 'GetAlphaElem', 'GetBetaElem' and 'GetGammaElem',
C     and
C         'SetMuElem', 'SetAlphaElem', 'SetBetaElem' and 'SetGammaElem'.
C
C     The tensors can be printed via
C         'PrintMu', 'PrintAlpha', 'PrintBeta' and 'PrintGamma'.
C     Extended informations can be printed via
C         'PrintMuX', 'PrintAlphaX', 'PrintBetaX' and 'PrintGammaX'.
C
C*MODULE TDXPRP  *DECK NLOINI
      SUBROUTINE NLOINI
      IMPLICIT NONE
      integer i,j,k,l,jelemm(1,3),jndexm(3),jelema(2,9),jndexa(3,3),
     &        jelemb(3,27),jndexb(3,3,3),jelemg(4,21),jndexg(3,3,3,3),
     &        jelemt(2,9),jndext(3,3),jelemUA(3,27),jndexUA(3,3,3)
C     Used common blocks
      integer NELEMM,IELEMM,INDEXM
      common/NLO_MU/NELEMM,IELEMM(1,3),INDEXM(3)
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
      integer NELEMB,IELEMB,INDEXB
      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
      integer NELEMG,IELEMG,INDEXG
      common/NLO_GAMMA/NELEMG,IELEMG(4,21),INDEXG(3,3,3,3)
      integer NELEMT,IELEMT,INDEXT
      common/NLO_THETA/NELEMT,IELEMT(2,9),INDEXT(3,3)
      integer NELEMUA,IELEMUA,INDEXUA
      common/NLO_UPPERA/NELEMUA,IELEMUA(3,27),INDEXUA(3,3,3)
C     number of tensor elements for mu, alpha, beta, gamma
C     description of each tensor elements for mu
      DATA JELEMM/1,2,3/
C     index of the tensor element for a given combination of xyz indexes
      DATA JNDEXM/1,2,3/
C     description of each tensor elements for alpha
      DATA JELEMA/1,1, 2,1, 3,1,
     &            1,2, 2,2, 3,2,
     &            1,3, 2,3, 3,3/
C     index of the tensor element for a given combination of xyz indexes
      DATA JNDEXA/1,2,3,
     &            4,5,6,
     &            7,8,9/
C     description of each tensor elements for beta
      DATA JELEMB/1,1,1, 2,1,1, 3,1,1,
     &            1,2,1, 2,2,1, 3,2,1,
     &            1,3,1, 2,3,1, 3,3,1,
     &            1,1,2, 2,1,2, 3,1,2,
     &            1,2,2, 2,2,2, 3,2,2,
     &            1,3,2, 2,3,2, 3,3,2,
     &            1,1,3, 2,1,3, 3,1,3,
     &            1,2,3, 2,2,3, 3,2,3,
     &            1,3,3, 2,3,3, 3,3,3/
C     index of the tensor element for a given combination of xyz indexes
      DATA JNDEXB/ 1, 2, 3,  4, 5, 6,  7, 8, 9,
     &            10,11,12, 13,14,15, 16,17,18,
     &            19,20,21, 22,23,24, 25,26,27/
C     description of each tensor elements for gamma
      DATA JELEMG/1,1,1,1, 2,2,2,2, 3,3,3,3,
     &            1,1,2,2, 2,2,1,1, 3,3,1,1,
     &            1,1,3,3, 2,2,3,3, 3,3,2,2,
     &            1,2,1,2, 2,1,2,1, 3,1,3,1,
     &            1,3,1,3, 2,3,2,3, 3,2,3,2,
     &            1,2,2,1, 2,1,1,2, 3,1,1,3,
     &            1,3,3,1, 2,3,3,2, 3,2,2,3/
C     index of the tensor element for a given combination of xyz indexes
      DATA JNDEXG/1, 0, 0,  0, 5, 0,  0, 0, 6,
     &            0,11, 0, 16, 0, 0,  0, 0, 0,
     &            0, 0,12,  0, 0, 0, 19, 0, 0,
     &            0,17, 0, 10, 0, 0,  0, 0, 0,
     &            4, 0, 0,  0, 2, 0,  0, 0, 9,
     &            0, 0, 0,  0, 0,15,  0,20, 0,
     &            0, 0,18,  0, 0, 0, 13, 0, 0,
     &            0, 0, 0,  0, 0,21,  0,14, 0,
     &            7, 0, 0,  0, 8, 0,  0, 0, 3/
C     description of each tensor elements for theta
      DATA JELEMT/1,1, 2,1, 3,1,
     &            1,2, 2,2, 3,2,
     &            1,3, 2,3, 3,3/
C     index of the tensor element for a given combination of xyz indexes
      DATA JNDEXT/1,2,3,
     &            4,5,6,
     &            7,8,9/
C     description of each tensor elements for beta
      DATA JELEMUA/1,1,1, 2,1,1, 3,1,1,
     &             1,2,1, 2,2,1, 3,2,1,
     &             1,3,1, 2,3,1, 3,3,1,
     &             1,1,2, 2,1,2, 3,1,2,
     &             1,2,2, 2,2,2, 3,2,2,
     &             1,3,2, 2,3,2, 3,3,2,
     &             1,1,3, 2,1,3, 3,1,3,
     &             1,2,3, 2,2,3, 3,2,3,
     &             1,3,3, 2,3,3, 3,3,3/
C     index of the tensor element for a given combination of xyz indexes
      DATA JNDEXUA/ 1, 2, 3,  4, 5, 6,  7, 8, 9,
     &             10,11,12, 13,14,15, 16,17,18,
     &             19,20,21, 22,23,24, 25,26,27/
C
      nelemm  = 3
      nelema  = 9
      nelemb  = 27
      nelemg  = 21
      nelemt  = 9
      nelemua = 27
C
      do j=1,nelemm
        ielemm(1,j)=jelemm(1,j)
      enddo
      do j=1,nelema
         do i=1,2
            ielema(i,j) = jelema(i,j)
         enddo
      enddo
      do j=1,nelemb
         do i=1,3
            ielemb(i,j) = jelemb(i,j)
         enddo
      enddo
      do j=1,nelemg
         do i=1,4
            ielemg(i,j) = jelemg(i,j)
         enddo
      enddo
      do j=1,nelemt
         do i=1,2
            ielemt(i,j) = jelemt(i,j)
         enddo
      enddo
      do j=1,nelemua
         do i=1,3
            ielemua(i,j) = jelemua(i,j)
         enddo
      enddo
      do i=1,3
         indexm(i)  =jndexm(i)
         do j=1,3
            indexa(j,i)=jndexa(j,i)
            indext(j,i)=jndext(j,i)
            do k=1,3
               indexb(k,j,i)=jndexb(k,j,i)
               indexua(k,j,i)=jndexua(k,j,i)
               do l=1,3
                  indexg(l,k,j,i)=jndexg(l,k,j,i)
               enddo
            enddo
         enddo
      enddo
      END
C
C*MODULE TDXPRP  *DECK GetMuElem
      function GetMuElem(mu,i)
      implicit none
C     Arguments
      double precision GetMUELEM
      double precision mu(3)
      integer i
C     Used common block
      integer NELEMM,IELEMM,INDEXM
      common/NLO_MU/NELEMM,IELEMM(1,3),INDEXM(3)
C     Body of the subprogram
      GetMUELEM = mu(INDEXM(i))
      end
C
C*MODULE TDXPRP  *DECK GetAlphaElem
      function GetAlphaElem(alpha,i,j)
      implicit none
C     Arguments
      double precision GetALPHAELEM
      double precision alpha(9)
      integer i,j
C     Used common block
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C     Body of the subprogram
      GetALPHAELEM = alpha(INDEXA(i,j))
      end
C
C*MODULE TDXPRP  *DECK GetBetaElem
      function GetBetaElem(beta,i,j,k)
      implicit none
C     Arguments
      double precision GetBETAELEM
      double precision beta(27)
      integer i,j,k
C     Used common block
      integer NELEMB,IELEMB,INDEXB
      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
C     Body of the subprogram
      GetBETAELEM = beta(INDEXB(i,j,k))
      end
C
C*MODULE TDXPRP  *DECK GetGammaElem
      function GetGammaElem(gamma,i,j,k,l)
      implicit none
C     Arguments
      double precision GetGAMMAELEM
      double precision gamma(21)
      integer i,j,k,l
C     Used common block
      integer NELEMG,IELEMG,INDEXG
      common/NLO_GAMMA/NELEMG,IELEMG(4,21),INDEXG(3,3,3,3)
C     Body of the subprogram
      GetGAMMAELEM = gamma(INDEXG(i,j,k,l))
      end
C
C--C*MODULE TDXPRP  *DECK GetThetaElem
C--      function GetThetaElem(theta,i,j)
C--      implicit none
C--c     Arguments
C--      double precision GetThetaELEM
C--      double precision theta(9)
C--      integer i,j
C--c     Used common block
C--      integer NELEMT,IELEMT,INDEXT
C--      common/NLO_THETA/NELEMT,IELEMT(2,9),INDEXT(3,3)
C--c     Body of the subprogram
C--      GetTHETAELEM = theta(INDEXT(i,j))
C--      end
C--c
C--C*MODULE TDXPRP  *DECK GetUpperAElem
C--      function GetUpperAElem(uppera,i,j,k)
C--      implicit none
C--c     Arguments
C--      double precision GetUpperAELEM
C--      double precision uppera(27)
C--      integer i,j,k
C--c     Used common block
C--      integer NELEMUA,IELEMUA,INDEXUA
C--      common/NLO_UPPERA/NELEMUA,IELEMUA(3,27),INDEXUA(3,3,3)
C--c     Body of the subprogram
C--      GetUPPERAELEM = uppera(INDEXUA(i,j,k))
C--      end
C
C*MODULE TDXPRP  *DECK SetMuElem
      subroutine SetMuElem(mu,i,v)
      implicit none
C     Arguments
      double precision mu(3),v
      integer i
C     Used common block
      integer NELEMM,IELEMM,INDEXM
      common/NLO_MU/NELEMM,IELEMM(1,3),INDEXM(3)
C     Body of the subprogram
      mu(INDEXM(i)) = v
      end
C
C--C*MODULE TDXPRP  *DECK SetAlphaElem
C--      subroutine SetAlphaElem(alpha,i,j,v)
C--      implicit none
C--c     Arguments
C--      double precision alpha(9),v
C--      integer i,j
C--c     Used common block
C--      integer NELEMA,IELEMA,INDEXA
C--      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C--c     Body of the subprogram
C--      alpha(INDEXA(i,j)) = v
C--      end
C--c
C--C*MODULE TDXPRP  *DECK SetBetaElem
C--      subroutine SetBetaElem(beta,i,j,k,v)
C--      implicit none
C--c     Arguments
C--      double precision beta(27),v
C--      integer i,j,k
C--c     Used common block
C--      integer NELEMB,IELEMB,INDEXB
C--      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
C--c     Body of the subprogram
C--      beta(INDEXB(i,j,k)) = v
C--      end
C--c
C--C*MODULE TDXPRP  *DECK SetGammaElem
C--      subroutine SetGammaElem(gamma,i,j,k,l,v)
C--      implicit none
C--c     Arguments
C--      double precision gamma(21),v
C--      integer i,j,k,l
C--c     Used common block
C--      integer NELEMG,IELEMG,INDEXG
C--      common/NLO_GAMMA/NELEMG,IELEMG(4,21),INDEXG(3,3,3,3)
C--c     Body of the subprogram
C--      gamma(INDEXG(i,j,k,l)) = v
C--      end
C--c
C--C*MODULE TDXPRP  *DECK SetThetaElem
C--      subroutine SetThetaElem(theta,i,j,v)
C--      implicit none
C--c     Arguments
C--      double precision theta(9),v
C--      integer i,j
C--c     Used common block
C--      integer NELEMT,IELEMT,INDEXT
C--      common/NLO_THETA/NELEMT,IELEMT(2,9),INDEXT(3,3)
C--c     Body of the subprogram
C--      theta(INDEXT(i,j)) = v
C--      end
C--c
C--C*MODULE TDXPRP  *DECK SetUpperAElem
C--      subroutine SetUpperAElem(uppera,i,j,k,v)
C--      implicit none
C--c     Arguments
C--      double precision uppera(27),v
C--      integer i,j,k
C--c     Used common block
C--      integer NELEMUA,IELEMUA,INDEXUA
C--      common/NLO_UPPERA/NELEMUA,IELEMUA(3,27),INDEXUA(3,3,3)
C--c     Body of the subprogram
C--      uppera(INDEXUA(i,j,k)) = v
C--      end
C
C*MODULE TDXPRP  *DECK GetMuNorm
      function GetMuNorm(mu)
      implicit none
C     Arguments
      double precision GetMUNORM
      double precision mu(3)
C     Called subprograms
      double precision GetMUELEM
C     Body of the subprogram
      GetMUNORM = sqrt(GetMUELEM(mu,1)**2+
     &                 GetMUELEM(mu,2)**2+
     &                 GetMUELEM(mu,3)**2)
      end
C
C*MODULE TDXPRP  *DECK GetAlphaMean
      function GetAlphaMean(alpha)
      implicit none
C     Arguments
      double precision GetALPHAMEAN
      double precision alpha(9)
C     variables and constants
      double precision THREE
      parameter (THREE = 3.0d+00)
C     Called subprograms
      double precision GetALPHAELEM
C     Body of the subprogram
      GetALPHAMEAN = (GetALPHAELEM(alpha,1,1)+
     &                GetALPHAELEM(alpha,2,2)+
     &                GetALPHAELEM(alpha,3,3))/THREE
      end
C
C*MODULE TDXPRP  *DECK GetAlphaAnisotropy
      function GetAlphaAnisotropy(alpha)
      implicit none
C     Arguments
      double precision GetAlphaAnisotropy
      double precision alpha(9)
C     Variables and constants
      double precision dumXY, dumXZ,dumYZ,temp
      double precision HALF,THREE
      parameter (HALF=0.5d+00,THREE=3.0d+00)
C     Called subprograms
      double precision GetAlphaElem
C     Body of the subprogram
      dumXY = GetAlphaElem(alpha,1,1)-GetAlphaElem(alpha,2,2)
      dumXZ = GetAlphaElem(alpha,1,1)-GetAlphaElem(alpha,3,3)
      dumYZ = GetAlphaElem(alpha,2,2)-GetAlphaElem(alpha,3,3)
C
      temp= HALF*(dumXY*dumXY+dumXZ*dumXZ+dumYZ*dumYZ)+
     &     THREE*( GetAlphaElem(alpha,1,2)*GetAlphaElem(alpha,2,1)
     &            +GetAlphaElem(alpha,1,3)*GetAlphaElem(alpha,3,1)
     &            +GetAlphaElem(alpha,2,3)*GetAlphaElem(alpha,3,2))
      if(abs(temp).lt.1.0D-06) temp=0.0D+00
      GetAlphaAnisotropy = sqrt(temp)
      end
C
C*MODULE TDXPRP  *DECK GetBetaXYZ
      function GetBetaXYZ(beta,i)
      implicit none
C     Arguments
      double precision GetBetaXYZ
      double precision beta(27)
      integer i
C     Variables and constants
      integer j
      double precision dummy, ZERO
      parameter (ZERO=0.0d+00)
C     Called subprograms
      double precision GetBetaElem
C     Body of the subprogram
      dummy = ZERO
      do j = 1,3
         dummy = dummy + GetBetaElem(beta,i,j,j)
     &                 + GetBetaElem(beta,j,i,j)
     &                 + GetBetaElem(beta,j,j,i)
      end do
      GetBetaXYZ = dummy
      end
C
C*MODULE TDXPRP  *DECK GetBetaKerrPerpXYZ
      function GetBetaKerrPerpXYZ(beta,i)
      implicit none
C     Arguments
      double precision GetBetaKerrPerpXYZ
      double precision beta(27)
      integer i
C     Variables and constants
      integer j
      double precision dummy, ZERO,TWO,THREE
      parameter (ZERO=0.0d+00,TWO=2.0d+00,THREE=3.0d+00)
C     Called subprograms
      double precision GetBetaElem
C     Body of the subprogram
      dummy = ZERO
      do j = 1,3
         dummy = dummy + TWO  *GetBetaElem(beta,i,j,j)
     &                 - THREE*GetBetaElem(beta,j,i,j)
     &                 + TWO  *GetBetaElem(beta,j,j,i)
      end do
      GetBetaKerrPerpXYZ = dummy
      end
C
C*MODULE TDXPRP  *DECK GetBetaBarComp
      function GetBetaBarComp(beta,i)
      implicit none
C     Arguments
      double precision GetBetaBarComp
      double precision beta(27)
      integer i
C     Variables and constants
      double precision FIVE
      parameter (FIVE=5.0d+00)
C     Called subprograms
      double precision GetBetaXYZ
C     Body of the subprogram
      GetBetaBarComp = GetBetaXYZ(beta,i)/FIVE
      end
C
C*MODULE TDXPRP  *DECK GetBetaBar
      function GetBetaBar(beta,mu)
      implicit none
C     Arguments
      double precision GetBetaBar
      double precision beta(27),mu(3)
C     Variables and constants
      double precision norm
C     Called subprograms
      double precision GetMuNorm,GetBetaBarComp,GetMuElem
C     Body of the subprogram
      norm = GetMuNorm(mu)
      GetBetaBar =
     &     (
     &       GetBetaBarComp(beta,1)*GetMuElem(mu,1) +
     &       GetBetaBarComp(beta,2)*GetMuElem(mu,2) +
     &       GetBetaBarComp(beta,3)*GetMuElem(mu,3)
     &     ) /norm
      end
C
C*MODULE TDXPRP  *DECK GetBetaKerrPerp
      function GetBetaKerrPerp(beta,mu)
      implicit none
C     Arguments
      double precision GetBetaKerrPerp
      double precision beta(27),mu(3)
C     Variables and constants
      double precision norm,FIVE
      parameter (FIVE=5.0d+00)
C     Called subprograms
      double precision GetMuNorm,GetBetaKerrPerpXYZ,GetMuElem
C     Body of the subprogram
      norm = GetMuNorm(mu)
      GetBetaKerrPerp =
     &     (
     &       GetBetaKerrPerpXYZ(beta,1)*GetMuElem(mu,1) +
     &       GetBetaKerrPerpXYZ(beta,2)*GetMuElem(mu,2) +
     &       GetBetaKerrPerpXYZ(beta,3)*GetMuElem(mu,3)
     &     ) /norm/FIVE
      end
C
C*MODULE TDXPRP  *DECK GetBetaKerr
      function GetBetaKerr(beta,mu)
      implicit none
C     Arguments
      double precision GetBetaKerr
      double precision beta(27), mu(3)
C     Variables and constants
      double precision mult
      parameter (mult=1.5d+00)
C     Called subprograms
      double precision GetBetaBar,GetBetaKerrPerp
C     Body of the subprogram
      GetBetaKerr = mult*(GetBetaBar(beta,mu)-GetBetaKerrPerp(beta,mu))
      end
C
C*MODULE TDXPRP  *DECK GetBetaAnisotropy
      function GetBetaAnisotropy(beta)
      implicit none
C     Arguments
      double precision GetBetaAnisotropy
      double precision beta(27)
C     Variables and constants
      double precision SIX
      parameter (SIX=6.0d+00)
C     Called subprograms
      double precision GetBetaElem
C     Body of the subprogram
      GetBetaAnisotropy = (
     &     GetBetaElem(beta,1,2,3)-GetBetaElem(beta,1,3,2) +
     &     GetBetaElem(beta,2,3,1)-GetBetaElem(beta,2,1,3) +
     &     GetBetaElem(beta,3,1,2)-GetBetaElem(beta,3,2,1)
     &                    )/SIX
      end
C
C*MODULE TDXPRP  *DECK GetBetaVecComp
      function GetBetaVecComp(beta,i)
      implicit none
C     Arguments
      double precision GetBetaVecComp
      double precision beta(27)
      integer i
C     Variables and constants
      double precision THREE
      parameter (THREE=3.0d+00)
C     Called subprograms
      double precision GetBetaXYZ
C     Body of the subprogram
      GetBetaVecComp = GetBetaXYZ(beta,i)/THREE
      end
C
C*MODULE TDXPRP  *DECK GetBetaVec
      function GetBetaVec(beta)
      implicit none
C     Arguments
      double precision GetBetaVec
      double precision beta
C     Variables and constants
      double precision dumX,dumY,dumZ
C     Called subprograms
      double precision GetBetaVecComp
C     Body of the subprogram
      dumX = GetBetaVecComp(beta,1)
      dumY = GetBetaVecComp(beta,2)
      dumZ = GetBetaVecComp(beta,3)
      GetBetaVec = sqrt( dumX*dumX+dumY*dumY+dumZ*dumZ )
      end
C
C*MODULE TDXPRP  *DECK GetGammaMean
      function GetGammaMean(gamma)
      implicit none
C     Arguments
      double precision GetGammaMean
      double precision gamma(21)
C     Variables and constants
      integer i,j
      double precision dummy,ZERO,FACT
      parameter (ZERO=0.0d+00,FACT=15.0d+00)
C     Called subprograms
      double precision GetGammaElem
C     Body of the subprogram
      dummy = ZERO
      do i = 1,3
         do j = 1,3
            dummy = dummy + GetGammaElem(gamma,i,i,j,j)
     &                    + GetGammaElem(gamma,i,j,i,j)
     &                    + GetGammaElem(gamma,i,j,j,i)
         end do
      end do
      GetGammaMean = dummy/fact
      end
C
C*MODULE TDXPRP  *DECK GetGammaKerrPerp
      function GetGammaKerrPerp(gamma)
      implicit none
C     Arguments
      double precision GetGammaKerrPerp
      double precision gamma(21)
C     Variables and constants
      integer i,j
      double precision dummy,ZERO,TWO,FACT
      parameter (ZERO=0.0d+00,TWO=2.0d+00,FACT=15.0d+00)
C     Called subprograms
      double precision GetGammaElem
C     Body of the subprogram
      dummy = ZERO
      do i = 1,3
         do j = 1,3
            dummy = dummy + TWO*GetGammaElem(gamma,i,i,j,j)
     &                    - GetGammaElem(gamma,i,j,i,j)
         end do
      end do
      GetGammaKerrPerp = dummy/fact
      end
C
C*MODULE TDXPRP  *DECK GetGammaKerr
      function GetGammaKerr(gamma)
      implicit none
C     Arguments
      double precision GetGammaKerr
      double precision gamma(21)
C     Variables and constants
      double precision FACT
      parameter (FACT=1.5d+00)
C     Called subprograms
      double precision GetGammaMean,GetGammaKerrPerp
C     Body of the subprogram
      GetGammaKerr = fact*(GetGammaMean(gamma)-GetGammaKerrPerp(gamma))
      end
C
C*MODULE TDXPRP  *DECK GetGammaChi1
      function GetGammaChi1(gamma)
      implicit none
C     Arguments
      double precision GetGammaChi1
      double precision gamma(21)
C     Variables and constants
      integer i,j
      double precision dummy,ZERO,FOUR,NINETY
      parameter (ZERO=0.0d+00,FOUR=4.0d+00,NINETY=90.0d+00)
C     Called subprograms
      double precision GetGammaElem
C     Body of the subprogram
      dummy = ZERO
      do i = 1,3
         do j = 1,3
            dummy = dummy + FOUR*GetGammaElem(gamma,i,i,j,j)
     &                    - GetGammaElem(gamma,i,j,i,j)
     &                    - GetGammaElem(gamma,i,j,j,i)
         end do
      end do
      GetGammaChi1 = dummy/NINETY
      end
C
C*MODULE TDXPRP  *DECK GetGammaChi2
      function GetGammaChi2(gamma)
      implicit none
C     Arguments
      double precision GetGammaChi2
      double precision gamma(21)
C     Variables and constants
      integer i,j
      double precision dummy,ZERO,FOUR,NINETY
      parameter (ZERO=0.0d+00,FOUR=4.0d+00,NINETY=90.0d+00)
C     Called subprograms
      double precision GetGammaElem
C     Body of the subprogram
      dummy = ZERO
      do i = 1,3
         do j = 1,3
            dummy = dummy - GetGammaElem(gamma,i,i,j,j)
     &                    + FOUR*GetGammaElem(gamma,i,j,i,j)
     &                    - GetGammaElem(gamma,i,j,j,i)
         end do
      end do
      GetGammaChi2 = dummy/NINETY
      end
C
C*MODULE TDXPRP  *DECK GetGammaChi3
      function GetGammaChi3(gamma)
      implicit none
C     Arguments
      double precision GetGammaChi3
      double precision gamma(21)
C     Variables and constants
      integer i,j
      double precision dummy,ZERO,FOUR,NINETY
      parameter (ZERO=0.0d+00,FOUR=4.0d+00,NINETY=90.0d+00)
C     Called subprograms
      double precision GetGammaElem
C     Body of the subprogram
      dummy = ZERO
      do i = 1,3
         do j = 1,3
            dummy = dummy - GetGammaElem(gamma,i,i,j,j)
     &                    - GetGammaElem(gamma,i,j,i,j)
     &                    + FOUR*GetGammaElem(gamma,i,j,j,i)
         end do
      end do
      GetGammaChi3 = dummy/NINETY
      end
C
C--C*MODULE TDXPRP  *DECK TransformToTraceLessTheta
C--      subroutine TransformToTraceLessTheta(theta)
C--      implicit none
C--c     Arguments
C--      double precision theta(9)
C--c     variables and constants
C--      integer i
C--      double precision mean,dummy
C--      double precision TWO,THREE,fact
C--      parameter (TWO=2.0d+00,THREE=3.0d+00,fact=THREE/TWO)
C--c     Called subprograms
C--      double precision GetAlphaMean,GetTHETAELEM
C--c     Body of the subprogram
C--      call dscal(9,fact,theta,1)
C--      mean=GetAlphaMean(theta)
C--      do i=1,3
C--         dummy=GetTHETAELEM(theta,i,i)-mean
C--         call SetTHETAELEM(theta,i,i,dummy)
C--      end do
C--      end
C
C*MODULE TDXPRP  *DECK PrintMu
      subroutine PrintMu(iw,mu)
      implicit none
C     Arguments
      integer iw
      double precision mu(3)
C     Variables and constants
      integer i
      character*1 FIELDDIR(3)
C     Called subprograms
      double precision GetMuElem
C     Body of the subprogram
      DATA FIELDDIR /'x','y','z'/
      write(iw,9000) (FIELDDIR(i), i=1,3)
      write(iw,9001) (GetMuElem(mu,i), i=1,3)
 9000 format(6X,3(2X,7X,A1,7X))
 9001 format(6X,3(2X,F15.6))
      end
C
C*MODULE TDXPRP  *DECK PrintAlpha
      subroutine PrintAlpha(iw,alpha)
      implicit none
C     Arguments
      integer iw
      double precision alpha(9),alp1,alp2,alp3
C     Variables and constants
      integer i
      character*1 FIELDDIR(3)
C     Called subprograms
      double precision GetAlphaElem
C     Body of the subprogram
      DATA FIELDDIR /'x','y','z'/
      write(iw,9000) (FIELDDIR(i), i=1,3)
      do i = 1,3
         alp1= GetAlphaElem(alpha,i,1)
         alp2= GetAlphaElem(alpha,i,2)
         alp3= GetAlphaElem(alpha,i,3)
         write(iw,9001) FIELDDIR(i),alp1,alp2,alp3
C--      write(iw,9001) FIELDDIR(i),(GetAlphaElem(alpha,i,j), j=1,3)
      end do
 9000 format(6X,       3(2X,7X,A1,7X))
 9001 format(4X,A1,".",3(2X,F15.6))
      end
C
C*MODULE TDXPRP  *DECK PrintBeta
      subroutine PrintBeta(iw,beta)
      implicit none
C     Arguments
      integer iw
      double precision beta(27)
C     Variables and constants
      integer i,j
      character*1 FIELDDIR(3)
      double precision beta1,beta2,beta3
C     Called subprograms
      double precision GetBetaElem
C     Body of the subprogram
      DATA FIELDDIR /'x','y','z'/
      write(iw,9000) (FIELDDIR(i), i=1,3)
      do i = 1,3
         do j = 1,3
            beta1 = GetBetaElem(beta,i,j,1)
            beta2 = GetBetaElem(beta,i,j,2)
            beta3 = GetBetaElem(beta,i,j,3)
            write(iw,9001) FIELDDIR(i),FIELDDIR(j),beta1,beta2,beta3
C--  &                     (GetBetaElem(beta,i,j,k), k=1,3)
         end do
      end do
 9000 format(6X,          3(2X,7X,A1,7X))
 9001 format(3X,A1,A1,".",3(2X,F15.6))
      end
C
C*MODULE TDXPRP  *DECK printGamma
      subroutine PrintGamma(iw,gamma)
      implicit none
C     Arguments
      integer iw
      double precision gamma(21)
C     Variables and constants
      integer i,j,k
      character*1 FIELDDIR(3)
C     Used common block
      integer NELEMG,IELEMG,INDEXG
      common/NLO_GAMMA/NELEMG,IELEMG(4,21),INDEXG(3,3,3,3)
C     Body of the subprogram
      DATA FIELDDIR /'x','y','z'/
      do i = 1,NELEMG,3
         write(iw,9001) ((FIELDDIR(IELEMG(k,j)),k=1,4),gamma(j),j=i,i+2)
      end do
 9001 format(3(2X,4A1,2X,F15.6))
      end
C--c
C--C*MODULE TDXPRP  *DECK PrintTheta
C--      subroutine PrintTheta(iw,theta)
C--      implicit none
C--c     Arguments
C--      integer iw
C--      double precision theta(9),th1,th2,th3
C--c     Variables and constants
C--      integer i
C--      character*1 FIELDDIR(3)
C--c     Called subprograms
C--      double precision GetThetaElem
C--c     Body of the subprogram
C--      DATA FIELDDIR /'x','y','z'/
C--      write(iw,9000) (FIELDDIR(i), i=1,3)
C--      do i = 1,3
C--         th1 = GetThetaElem(theta,i,1)
C--         th2 = GetThetaElem(theta,i,2)
C--         th3 = GetThetaElem(theta,i,3)
C--         write(iw,9001) FIELDDIR(i),th1,th2,th3
C--ccc      write(iw,9001) FIELDDIR(i),(GetThetaElem(theta,i,j), j=1,3)
C--      end do
C-- 9000 format(6X,       3(2X,7X,A1,7X))
C-- 9001 format(4X,A1,".",3(2X,F15.6))
C--      end
C--c
C--C*MODULE TDXPRP  *DECK PrintUpperA
C--      subroutine PrintUpperA(iw,uppera)
C--      implicit none
C--c     Arguments
C--      integer iw
C--      double precision uppera(27)
C--c     Variables and constants
C--      integer i,j
C--      character*1 FIELDDIR(3)
C--      double precision gua1,gua2,gua3
C--c     Called subprograms
C--      double precision GetUpperAElem
C--c     Body of the subprogram
C--      DATA FIELDDIR /'x','y','z'/
C--      write(iw,9000) (FIELDDIR(i), i=1,3)
C--      do i = 1,3
C--         do j = 1,3
C--            gua1 = GetUpperAElem(uppera,i,j,1)
C--            gua2 = GetUpperAElem(uppera,i,j,2)
C--            gua3 = GetUpperAElem(uppera,i,j,3)
C--            write(iw,9001) FIELDDIR(i),FIELDDIR(j),gua1,gua2,gua3
C--ccc  &                     (GetUpperAElem(uppera,i,j,k), k=1,3)
C--         end do
C--         write(iw,9002)
C--      end do
C-- 9000 format(7X,              3(2X,7X,A1,7X))
C-- 9001 format(3X,A1," , ",A1,".",3(2X,F15.6))
C-- 9002 format()
C--      end
C
C*MODULE TDXPRP  *DECK PrintMuX
      subroutine PrintMuX(iw,mu,title)
      implicit none
C     Arguments
      integer iw
      double precision mu(3)
      character*(*) title
C     Called subprograms
      double precision GetMuNorm
C     Body of the subprogram
      write(iw,9001) title
      call PrintMu(iw,mu)
      write(iw,9002) GetMuNorm(mu)
 9001 format(/,2X,A)
 9002 format(/,2X,'norm : ',F15.6)
      end
C
C*MODULE TDXPRP  *DECK PrintAlphaX
      subroutine PrintAlphaX(iw,alpha,w,title)
      implicit none
C     Arguments
      integer iw
      double precision alpha(9),w(2)
      character*(*) title
C     Called subprograms
      double precision GetAlphaMean,GetAlphaAnisotropy
C     Body of the subprogram
      write(iw,9001) title,w(1),w(2)
      call PrintAlpha(iw,alpha)
      write(iw,9002) GetAlphaMean(alpha),GetAlphaAnisotropy(alpha)
 9001 format(/,2X,A,'(',F10.6,';',F10.6,')')
 9002 format(/,2X,'Mean       : ',F15.6,
     &       /,2X,'Anisotropy : ',F15.6)
      end
C
C*MODULE TDXPRP  *DECK PrintBetaX
      subroutine PrintBetaX(iw,beta,mu,w,title)
      implicit none
C     Arguments
      integer iw
      double precision beta(27),mu(3),w(3)
      character*(*) title
C     Variables and constants
      integer i
      character*1 FIELDDIR(3)
      double precision ZERO,bbar(3)
      parameter (ZERO=0.0d+00)
C     Called subprograms
      double precision GetBetaBarComp,GetBetaBar,
     &                 GetBetaKerrPerp,GetBetaKerr,
     &                 GetBetaAnisotropy,
     &                 GetBetaVecComp,GetBetaVec
C     Body of the subprogram
      DATA FIELDDIR /'x','y','z'/
      write(iw,9001) title,w(1),w(2),w(3)
      call PrintBeta(iw,beta)
      do i=1,3
         bbar(i) = GetBetaBarComp(beta,i)
      enddo
      write(iw,9002) (FIELDDIR(i),bbar(i),i=1,3),GetBetaBar(beta,mu)
C--   write(iw,9002) (FIELDDIR(i),GetBetaBarComp(beta,i),i=1,3),
C--  &               GetBetaBar(beta,mu)
      if(w(3).eq.ZERO) then
         write(iw,9003) GetBetaKerrPerp(beta,mu),GetBetaKerr(beta,mu)
      end if
      if(w(2).ne.w(3).and.-w(1).ne.w(2).and.-w(1).ne.w(3)) then
         write(iw,9004) GetBetaAnisotropy(beta)
      end if
      write(iw,9005) (FIELDDIR(i),GetBetaVecComp(beta,i),i=1,3),
     &               GetBetaVec(beta)
 9001 format(/,2X,A,'(',F10.6,';',F10.6,',',F10.6,')')
 9002 format(/,2X,A1,'   : ',F15.6,
     &         2X,'B(i)=(B(i,j,j)+B(j,i,j)+B(j,j,i))/5',
     &       /,2X,A1,'   : ',F15.6,
     &       /,2X,A1,'   : ',F15.6,
     &       /,2X,'BAR : ',F15.6,
     &         2X,'BAR=B(i)*MU(i)/|MU|')
 9003 format(  2X,'Perp : ',F15.6,
     &       /,2X,'Kerr : ',F15.6,2X,'3/2*(BAR-Perp)')
 9004 format(/,2X,'Anisotropy : ',F15.6,
     &         2X,'=(Bxyz-Bxzy+Byzx-Byxz+Bzxy-Bzyx)/6')
 9005 format(/,2X,A1,'   : ', F15.6,
     &         2X,'B(i)=(B(i,j,j)+B(j,i,j)+B(j,j,i))/3',
     &       /,2X,A1,'   : ', F15.6,
     &       /,2X,A1,'   : ', F15.6,
     &       /,2X,'VEC : ', F15.6,2X,'norm of Beta VEC')
      end
C
C*MODULE TDXPRP  *DECK PrintGammaX
      subroutine PrintGammaX(iw,gamma,w,title)
      implicit none
C     Arguments
      integer iw
      double precision gamma(21),w(4)
      character*(*) title
C     Variables and constants
      double precision ZERO
      parameter (ZERO=0.0d+00)
C     Called subprograms
      double precision GetGammaMean,GetGammaKerrPerp,GetGammaKerr
      double precision GetGammaChi1,GetGammaChi2,GetGammaChi3
C     Body of the subprogram
      write(iw,9001) title,w(1),w(2),w(3),w(4)
      call PrintGamma(iw,gamma)
      write(iw,9002) GetGammaMean(gamma)
      if(w(3).eq.ZERO.and.w(4).eq.ZERO) then
         write(iw,9003) GetGammaKerrPerp(gamma),GetGammaKerr(gamma)
      end if
      if(w(2).ne.w(3).and.
     &   w(2).ne.ZERO.and.w(3).ne.ZERO.and.w(4).eq.ZERO) then
         write(iw,9004) GetGammaChi1(gamma),GetGammaChi2(gamma),
     &                  GetGammaChi3(gamma)
      end if
 9001 format(/,2X,A,'(',F10.6,';',F10.6,',',F10.6,',',F10.6,')')
 9002 format(/,2X,'Mean : ',F15.6)
 9003 format(  2X,'Perp : ',F15.6,
     &       /,2X,'Kerr : ',F15.6,2X,'3/2*(Mean-Perp)')
 9004 format(  2X,'Chi1 : ',F15.6,2X,
     &            '(4G(i,i,j,j)-G(i,j,i,j)-G(i,j,j,i))/90',
     &       /,2X,'Chi2 : ',F15.6,2X,
     &            '(-G(i,i,j,j)+4G(i,j,i,j)-G(i,j,j,i))/90',
     &       /,2X,'Chi3 : ',F15.6,2X,
     &            '(-G(i,i,j,j)-G(i,j,i,j)+4G(i,j,j,i))/90')
      end
C
C--C*MODULE TDXPRP  *DECK PrintThetaX
C--      subroutine PrintThetaX(iw,theta,title)
C--      implicit none
C--c     Arguments
C--      integer iw
C--      double precision theta(9),thetatl(9)
C--      character*(*) title
C--c     Body of the subprogram
C--      write(iw,9001) title
C--      call PrintTheta(iw,theta)
C--      call dcopy(9,theta,1,thetatl,1)
C--      call TransformToTraceLessTheta(thetatl)
C--      write(iw,9002) title
C--      call PrintTheta(iw,thetatl)
C-- 9001 format(/,2X,A)
C-- 9002 format(/,2X,A,'(traceless, ',
C--     &              '3/2[q_i (r_i r_j - r_i^2/3 delta_ij)])')
C--      end
C--c
C--C*MODULE TDXPRP  *DECK PrintUpperAX
C--      subroutine PrintUpperAX(iw,upperA,w,title)
C--      implicit none
C--c     Arguments
C--      integer iw
C--      double precision upperA(3,9),w(2)
C--      character*(*) title
C--c     Variables
C--      integer i
C--      double precision dummy(9),upperatl(3,9)
C--c     Body of the subprogram
C--      write(iw,9001) title,w(1),w(2)
C--      call PrintUpperA(iw,upperA)
C--      do i=1,3
C--         call dcopy(9,upperA(i,1),3,dummy,1)
C--         call TransformToTraceLessTheta(dummy)
C--         call dcopy(9,dummy,1,upperatl(i,1),3)
C--      end do
C--      write(iw,9002) title,w(1),w(2)
C--      call PrintUpperA(iw,upperatl)
C--c
C-- 9001 format(/,2X,A,'(',F10.6,';',F10.6,')')
C-- 9002 format(/,2X,A,'(',F10.6,';',F10.6,')',
C--     &       /,4x,'(from traceless quadrupole: ',
C--     &            '3/2[q_i (r_i r_j - r_i^2/3 delta_ij)])')
C--      end
C
C*MODULE TDXPRP  *DECK VIBNLOINI
      subroutine VIBNLOINI
      implicit none
C
C            this sets /NLO_PERM/ if that is ever needed below
C--c     Used common block
C--      integer i,j
C--      integer IPERMA,IPERMB,IPERMG
C--      integer JPERMA(2,2),JPERMB(3,6),JPERMG(4,24)
C--      common/NLO_PERM/IPERMA(2,2),IPERMB(3,6),IPERMG(4,24)
C--c     description of the 2 permutations over the alpha indexes
C--      DATA JPERMA/1,2,
C--     &            2,1/
C--c     description of the 6 permutations over the beta indexes
C--      DATA jPERMB/1,2,3, 1,3,2,
C--     &            2,1,3, 2,3,1,
C--     &            3,1,2, 3,2,1/
C--c     description of the 24 permutations over the gamma indexes
C--      DATA jPERMG/
C--     &    1,2,3,4, 1,2,4,3, 1,3,2,4, 1,3,4,2, 1,4,2,3, 1,4,3,2,
C--     &    2,1,3,4, 2,1,4,3, 2,3,1,4, 2,3,4,1, 2,4,1,3, 2,4,3,1,
C--     &    3,1,2,4, 3,1,4,2, 3,2,1,4, 3,2,4,1, 3,4,1,2, 3,4,2,1,
C--     &    4,1,2,3, 4,1,3,2, 4,2,1,3, 4,2,3,1, 4,3,1,2, 4,3,2,1/
C--      do i=1,2
C--        do j=1,2
C--          iperma(j,i) = jperma(j,i)
C--        enddo
C--      enddo
C--      do i=1,6
C--        do j=1,3
C--          ipermb(j,i) = jpermb(j,i)
C--        enddo
C--      enddo
C--      do i=1,24
C--        do j=1,4
C--          ipermg(j,i) = jpermg(j,i)
C--        enddo
C--      enddo
C
      return
      end
C
C--c*MODULE TDXPRP  *DECK MakeMu2
C--      subroutine MakeMu2(w1)
C--      implicit none
C--c     Arguments
C--      double precision w1
C--c     Used common blocks
C--      integer soli,ngeom,nocc,nvir,nmos,naos,nxyz,nmodes,nfld,ibfld,
C--     &           iefld
C--      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,
C--     &           nfld,ibfld,iefld
C--      DOUBLE PRECISION XX
C--      COMMON /FMCOM / XX(1)
C--c     Variables
C--      double precision w(2)
C--      integer lw,ldmdq
C--      double precision dummy(1)
C--c     Called subprogram
C--      logical GetFNM,GetdMudQ
C--c     Body of the subprogram
C--      dummy(1) = 0.0d+00
C--      w(1) = -w1
C--      w(2) =  w1
C--c
C--      if(.not.GetdMudQ(dummy,.false.)) call MakedMudXNI()
C--c
C--      call NEWMEM(lW,nmodes)
C--      call NEWMEM(ldmdq,3*nmodes)
C--c
C--      if(.not.GetFNM(XX(lW),.true.,dummy,.false.,nxyz,nmodes)) then
C--         call abrtx("MAKEMU2:1")
C--      end if
C--      if(.not.GetdMudQ(XX(ldMdQ),.true.)) then
C--         call abrtx("MAKEMU2:2")
C--      end if
C--c
C--      call CompMu2(w,XX(lW),XX(ldMdQ),nmodes)
C--c
C--      call DELETEMEM(ldmdq,3*nmodes)
C--      call DELETEMEM(lW,nmodes)
C--c
C--      CALL TIMIT(1)
C--c
C--      end
C--c
C--c*MODULE TDXPRP  *DECK MakeMualpha
C--      subroutine MakeMualpha(w1,w2)
C--      implicit none
C--c     Arguments
C--      double precision w1,w2
C--c     Used common blocks
C--      integer soli,ngeom,nocc,nvir,nmos,naos,nxyz,nmodes,nfld,ibfld,
C--     &           iefld
C--      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,
C--     &           nfld,ibfld,iefld
C--      DOUBLE PRECISION XX
C--      COMMON /FMCOM / XX(1)
C--c     Variables
C--      double precision w(3)
C--      integer lw,ldmdq,ldadq
C--      double precision dummy(1)
C--      double precision ZERO
C--      parameter (ZERO=0.0d+00)
C--c     Called subprogram
C--      logical GotDip,GetFNM,GetdMudQ,GetdAlphadQ
C--c     Body of the subprogram
C--      dummy(1) = ZERO
C--      w(1) = -w1-w2
C--      w(2) =  w1
C--      w(3) =  w2
C--c
C--      call GetDipole(dummy,.false.,GotDip)
C--      if(.not.GotDip) call MakeDipole()
C--      if(.not.GetdMudQ(dummy,.false.)) call MakedMudXNI()
C--      if(.not.GetdAlphadQ(ZERO,dummy,.false.))
C--     &     call MakedalphadXNI(ZERO)
C--c
C--      call NEWMEM(lW,nmodes)
C--      call NEWMEM(ldmdq,3*nmodes)
C--      call NEWMEM(ldadq,9*nmodes)
C--c
C--      if(.not.GetFNM(XX(lW),.true.,dummy,.false.,nxyz,nmodes)) then
C--         call abrtx("MAKEMUALPHA:1")
C--      end if
C--      if(.not.GetdMudQ(XX(ldMdQ),.true.)) then
C--         call abrtx("MAKEMUALPHA:2")
C--      end if
C--      if(.not.GetdAlphadQ(ZERO,XX(ldAdQ),.true.)) then
C--         call abrtx("MAKEMUALPHA:3")
C--      end if
C--c
C--      call CompMuAlpha(w,XX(lW),XX(ldMdQ),XX(ldAdQ),nmodes)
C--c
C--      call DELETEMEM(ldadq,9*nmodes)
C--      call DELETEMEM(ldmdq,3*nmodes)
C--      call DELETEMEM(lW,nmodes)
C--c
C--      CALL TIMIT(1)
C--c
C--      end
C--c
C--C*MODULE TDXPRP  *DECK CompMu2
C--      subroutine CompMu2(Wopt,Wvib,dMdQ,nmode)
C--      implicit none
C--c     Arguments
C--      integer nmode
C--      double precision Wopt(2)
C--      double precision Wvib(nmode)
C--      double precision dMdQ(3,nmode)
C--c     Used common block
C--      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
C--      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C--      integer IPERMA,IPERMB,IPERMG
C--c         this common's initialization is commented out in vibnloini
C--      common/NLO_PERM/IPERMA(2,2),IPERMB(3,6),IPERMG(4,24)
C--      integer NELEMA,IELEMA,INDEXA
C--      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C--c     Variables and constants
C--      integer nelem
C--      parameter (nelem=9)
C--      double precision total(nelem),dummy(nelem)
C--      double precision w_a,w_i
C--      integer indexfld(2),i,j,index
C--      integer ifld,jfld,iperm,a
C--      character*50 buffer
C--      double precision HALF
C--      parameter (HALF=0.5d+00)
C--c     Called subprogram
C--      double precision GetMuElem
C--c     Body of the subprogram
C--c
C--      write(iw,3000) Wopt(1),Wopt(2)
C--c
C--      call VCLR(total,1,nelem)
C--      do a = 1,nmode
C--         w_a = Wvib(a)
C--c
C--         call VCLR(dummy,1,nelem)
C--c
C--         do index = 1,nelem
C--c
C--            ifld = IELEMA(1,index)
C--            jfld = IELEMA(2,index)
C--c
C--            indexfld(1) = ifld
C--            indexfld(2) = jfld
C--c
C--            do iperm=1,2
C--               i = indexfld(IPERMA(1,iperm))
C--               j = indexfld(IPERMA(2,iperm))
C--c
C--               w_i = Wopt(2)
C--c
C--               dummy(index) = dummy(index) +
C--     &             HALF*GetMuElem(dMdQ(1,a),i)
C--     &                 *GetMuElem(dMdQ(1,a),j)/(w_a*w_a-w_i*w_i)
C--            end do
C--            total(index)=total(index)+dummy(index)
C--         end do
C--c
C--         write(buffer,3001) a
C--         call PrintAlphaX(iw,dummy,Wopt,buffer)
C--c
C--      end do
C--c
C--      write(buffer,3002)
C--      call PrintAlphaX(iw,total,Wopt,buffer)
C--c
C-- 3000 FORMAT(/,1x,'********************',
C--     &     21('*'),'*',
C--     &       /,1x,'[mu2](0,0) to Alpha(',
C--     &     F10.6,';',F10.6,')',
C--     &       /,1x,'********************',
C--     &     21('*'),'*')
C-- 3001 FORMAT(/,1x,'[mu2](0,0) : mode ',I2)
C-- 3002 FORMAT(/,1x,'[mu2](0,0) : total')
C--      end
C--c
C--C*MODULE TDXPRP  *DECK CompMuAlpha
C--      subroutine CompMuAlpha(Wopt,Wvib,dMdQ,dAdQ,nmode)
C--      implicit none
C--c     Arguments
C--      integer nmode
C--      double precision Wopt(3)
C--      double precision Wvib(nmode)
C--      double precision dMdQ(3,nmode)
C--      double precision dAdQ(9,nmode)
C--c     Used common block
C--      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
C--      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C--      integer IPERMA,IPERMB,IPERMG
C--c         this common's initialization is commented out in vibnloini
C--      common/NLO_PERM/IPERMA(2,2),IPERMB(3,6),IPERMG(4,24)
C--      integer NELEMB,IELEMB,INDEXB
C--      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
C--c     Variables and constants
C--      integer nelem
C--      parameter (nelem=27)
C--      double precision total(nelem),dummy(nelem)
C--      double precision w_a,w_i,denom
C--      integer indexfld(3),j,k,index
C--      integer ifld,jfld,kfld,iperm,a
C--      character*50 buffer
C--      double precision HALF
C--      parameter (HALF=0.5d+00)
C--      double precision mu(3)
C--c     Called subprogram
C--      logical GotDip
C--      double precision GetMuElem,GetAlphaElem
C--c
C--c     Body of the subprogram
C--c
C--      mu(1) = 0.0d+00
C--      call GetDipole(mu,.true.,GotDip)
C--      if(.not.GotDip) then
C--         call abrtx("COMPMUALPHA")
C--      end if
C--c
C--      write(iw,3000) Wopt(1),Wopt(2)
C--c
C--      call VCLR(total,1,nelem)
C--      do a = 1,nmode
C--         w_a = Wvib(a)
C--c
C--         call VCLR(dummy,1,nelem)
C--c
C--         do index = 1,nelem
C--c
C--            ifld = IELEMB(1,index)
C--            jfld = IELEMB(2,index)
C--            kfld = IELEMB(3,index)
C--c
C--            indexfld(1) = ifld
C--            indexfld(2) = jfld
C--            indexfld(3) = kfld
C--c
C--            do iperm=1,6
C--ccc            i = indexfld(IPERMB(1,iperm))
C--               j = indexfld(IPERMB(2,iperm))
C--               k = indexfld(IPERMB(3,iperm))
C--c
C--               w_i = Wopt(IPERMA(1,iperm))
C--               denom = w_a*w_a-w_i*w_i
C--c
C--               dummy(index) = dummy(index) +
C--     &              HALF*GetMuElem(dMdQ(1,a),a)
C--     &                  *GetAlphaElem(dAdQ(1,a),j,k)/denom
C--            end do
C--            total(index)=total(index)+dummy(index)
C--         end do
C--c
C--         write(buffer,3001) a
C--         call PrintBetaX(iw,dummy,mu,Wopt,buffer)
C--c
C--      end do
C--c
C--      write(buffer,3002)
C--      call PrintBetaX(iw,total,mu,Wopt,buffer)
C--c
C-- 3000 FORMAT(/,1x,'********************',
C--     &     21('*'),'*',
C--     &       /,1x,'[mualpha](0,0) to Beta(',
C--     &     F10.6,';',F10.6,',',F10.6,')',
C--     &       /,1x,'********************',
C--     &     21('*'),'*')
C-- 3001 FORMAT(/,1x,'[mualpha](0,0) : mode ',I2)
C-- 3002 FORMAT(/,1x,'[mualpha](0,0) : total')
C--      end
C
C     Rdata and hRdata are vectors of size 7.
C     They should be type in Fortran 90 but since it's Fortran 77...
C
C     Rdata(1) : harmonic frequency in [cm^-1]
C     Rdata(2) : multiplicity
C     Rdata(3) : Ak
C     Rdata(4) : Gk
C     Rdata(5) : intensity in [Ang.^4/AMU]
C     Rdata(6) : l-depolarization ratio
C     Rdata(7) : n-depolarization ratio
C
C     hRdata(1) : harmonic frequency in [cm^-1]
C     hRdata(2) : multiplicity
C     hRdata(3) : Biii
C     hRdata(4) : Bijj
C     hRdata(5) : Bijk
C     hRdata(6) : intensity in [Ang.^6 AMU^-1 StatVolt^-2]
C     hRdata(7) : n-depolarization ratio
C     hRdata(8) : p-depolarization ratio
C
C     They are initiliazed by calling 'InitRamanData'
C     and 'InitHyperRamanData', respectively.
C
C     A summary of the Raman and hyper-Raman informations can be
C     printed via 'PrintRaman' and 'PrintHyperRaman'.
C
C*MODULE TDXPRP  *DECK InitRamanData
      subroutine InitRamanData(Rdata,dadq,freq)
      implicit none
C     Arguments
      double precision Rdata(7),dadq(9),freq
C     Variables and constants
      double precision ak,gk,intensity,rho_l,rho_n,den_l,den_n
      double precision FACTOR,FACT,ONE,THREE,FOUR,SIX,SEVEN,FORTYFIVE
      parameter (FACTOR=0.52917165d+00**4*1822.889d+00)
      parameter (FACT=219474.6313710d+00)
      parameter (ONE=1.0d+00,THREE=3.0d+00,FOUR=4.0d+00,SIX=6.0d+00,
     *           SEVEN=7.0d+00,FORTYFIVE=45.0d+00)
C     Called subprograms
      double precision GetAlphaMean,GetAlphaAnisotropy
C     Body of the subprogram
      ak = GetAlphaMean(dadq)
      gk = GetAlphaAnisotropy(dadq)
      intensity = FORTYFIVE*ak*ak+SEVEN*gk*gk
c
      den_l =  FORTYFIVE*ak*ak + FOUR*gk*gk
      if(abs(den_l).gt.1.0D-06) then
         rho_l = (THREE*gk*gk)/den_l
      else
         rho_l = 99999.999999D+00
      end if
c
      den_n = FORTYFIVE*ak*ak + SEVEN*gk*gk
      if(abs(den_n).gt.1.0D-06) then
         rho_n = (SIX*gk*gk)/den_n
      else
         rho_n = 99999.999999D+00
      end if
c
      Rdata(1) = freq*FACT
      Rdata(2) = ONE
      Rdata(3) = sqrt(factor)*ak
      Rdata(4) = sqrt(factor)*gk
      Rdata(5) = factor*intensity
      Rdata(6) = rho_l
      Rdata(7) = rho_n
      end
C
C*MODULE TDXPRP  *DECK InitHyperRamanData
      subroutine InitHyperRamanData(hRdata,dbdq,freq)
      implicit none
C     Arguments
      double precision hRdata(8),dbdq(27),freq
C     Variables and constants
      integer i,j,k
      double precision biii,bijj,bijk,intensity,rho_n,rho_p
      double precision FACTOR,FACT,ZERO,ONE,TWO
      parameter (FACTOR=(0.86392d+00/0.52917165d+00)**2*1822.889d+00)
      parameter (FACT=219474.6313710d+00)
      parameter (ZERO=0.0d+00,ONE=1.0d+00,TWO=2.0d+00)
C     Called subprograms
      double precision GetBetaElem
C     Body of the subprogram
      biii = ZERO
      bijj = ZERO
      bijk = ZERO
      do i = 1,3
         bijj = bijj
     &             +(1.0d+00/35.0d+00)*GetBetaElem(DBDQ,i,i,i)
     &                                *GetBetaElem(DBDQ,i,i,i)
         biii = biii
     &             +(1.0d+00/7.0d+00) *GetBetaElem(DBDQ,i,i,i)
     &                                *GetBetaElem(DBDQ,i,i,i)
         do j = 1,3
            if(i.ne.j) then
               bijj = bijj
     &            +(4.0d+00/105.0d+00)*GetBetaElem(DBDQ,i,i,i)
     &                                *GetBetaElem(DBDQ,i,j,j)
     &            -(4.0d+00/70.0d+00) *GetBetaElem(DBDQ,i,i,i)
     &                                *GetBetaElem(DBDQ,j,j,i)
     &            +(8.0d+00/105.0d+00)*GetBetaElem(DBDQ,i,i,j)
     &                                *GetBetaElem(DBDQ,i,i,j)
     &            +(3.0d+00/35.0d+00) *GetBetaElem(DBDQ,i,j,j)
     &                                *GetBetaElem(DBDQ,i,j,j)
     &            -(4.0d+00/70.0d+00) *GetBetaElem(DBDQ,i,i,j)
     &                                *GetBetaElem(DBDQ,j,i,i)
               biii = biii
     &            +(4.0d+00/35.0d+00) *GetBetaElem(DBDQ,i,i,j)
     &                                *GetBetaElem(DBDQ,i,i,j)
     &            +(2.0d+00/35.0d+00) *GetBetaElem(DBDQ,i,i,i)
     &                                *GetBetaElem(DBDQ,i,j,j)
     &            +(4.0d+00/35.0d+00) *GetBetaElem(DBDQ,j,i,i)
     &                                *GetBetaElem(DBDQ,i,i,j)
     &            +(4.0d+00/35.0d+00) *GetBetaElem(DBDQ,i,i,i)
     &                                *GetBetaElem(DBDQ,j,j,i)
     &            +(1.0d+00/35.0d+00) *GetBetaElem(DBDQ,j,i,i)
     &                                *GetBetaElem(DBDQ,j,i,i)
C              trick : k must different from i and j:
C              i j k i+j+k
C              1 2 3   6
C              1 3 2   6
C              2 1 3   6
C              2 3 1   6
C              3 1 2   6
C              3 2 1   6
C              therefore, i+j+k is always equal to 6
               k=6-i-j
               bijj = bijj
     &            +(1.0d+00/35.0d+00) *GetBetaElem(DBDQ,i,j,j)
     &                                *GetBetaElem(DBDQ,i,k,k)
     &            -(4.0d+00/210.0d+00)*GetBetaElem(DBDQ,i,i,k)
     &                                *GetBetaElem(DBDQ,j,j,k)
     &            -(4.0d+00/210.0d+00)*GetBetaElem(DBDQ,i,i,j)
     &                                *GetBetaElem(DBDQ,j,k,k)
     &            +(2.0d+00/35.0d+00) *GetBetaElem(DBDQ,i,j,k)
     &                                *GetBetaElem(DBDQ,i,j,k)
     &            -(4.0d+00/210.0d+00)*GetBetaElem(DBDQ,i,j,k)
     &                                *GetBetaElem(DBDQ,j,i,k)
               biii = biii
     &            +(4.0d+00/105.0d+00)*GetBetaElem(DBDQ,i,i,j)
     &                                *GetBetaElem(DBDQ,j,k,k)
     &            +(1.0d+00/105.0d+00)*GetBetaElem(DBDQ,j,i,i)
     &                                *GetBetaElem(DBDQ,j,k,k)
     &            +(4.0d+00/105.0d+00)*GetBetaElem(DBDQ,i,i,j)
     &                                *GetBetaElem(DBDQ,k,k,j)
     &            +(2.0d+00/105.0d+00)*GetBetaElem(DBDQ,i,j,k)
     &                                *GetBetaElem(DBDQ,i,j,k)
     &            +(4.0d+00/105.0d+00)*GetBetaElem(DBDQ,i,j,k)
     &                                *GetBetaElem(DBDQ,j,i,k)
            end if
         end do
      end do
      intensity = (biii+bijj) * FACTOR
      if(abs(biii+bijj).gt.1.0D-06) then
         rho_n = TWO * bijj / (biii+bijj)
      else
         rho_n = 99999.999999D+00
      endif
      if(abs(biii).gt.1.0D-06) then
         rho_p = bijj / biii
      else
         rho_n = 99999.999999D+00
      endif
      hRdata(1) = freq*FACT
      hRdata(2) = ONE
      hRdata(3)=biii
      hRdata(4)=bijj
      hRdata(5)=bijk
      hRdata(6)=intensity
      hRdata(7)=rho_n
      hRdata(8)=rho_p
      end
C
C*MODULE TDXPRP  *DECK PrintRaman
      subroutine PrintRaman(iw,listdadq,nmodes,listfreq,omega)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision listdadq(9,nmodes),listfreq(nmodes),omega
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lTable
C     Body of the subprogram
      call NEWMEM(lTable,6*nmodes)
      call PrintRamanX(iw,listdadq,nmodes,listfreq,omega,XX(lTable))
      call DELETEMEM(lTable,6*nmodes)
      end
C
C*MODULE TDXPRP  *DECK PrintRamanX
      subroutine PrintRamanX(iw,listdadq,nmodes,listfreq,omega,table)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision listdadq(9,nmodes),listfreq(nmodes),omega
      double precision table(6,nmodes)
C     Variables and constants
      logical samemode
      integer imode,i,itab,ntab
      double precision dadq(9),Rdata(7)
      double precision last_f,last_i,diff_f,diff_i,intensity,maxint
      double precision ZERO,ONE,THRS_I,THRS_F
      parameter (ZERO=0.0d+00,ONE=1.0d+00)
      parameter (THRS_I=0.1d+00,THRS_F=0.1d+00)
C     Called subprogram
      double precision abs
C     Body of the subprogram
      ntab = 0
      last_f = ZERO
      last_i = ZERO
      maxint = ZERO
      intensity = ZERO
      do imode = 1,nmodes
         do i = 1,9
            dadq(i) = listdadq(i,imode)
         end do
         call InitRamanData(Rdata,dadq,listfreq(imode))
         diff_f = abs(Rdata(1)-last_f)
         diff_i = abs(Rdata(5)-last_i)
         last_f = Rdata(1)
         last_i = Rdata(5)
         samemode = (diff_i.le.THRS_I).and.(diff_f.lt.THRS_F)
         if(samemode) then
            intensity = intensity + Rdata(5)
         else
            ntab = ntab+1
            table(1,ntab) = ZERO
            table(2,ntab) = ZERO
            table(3,ntab) = ZERO
            table(4,ntab) = ZERO
            table(5,ntab) = ZERO
            table(6,ntab) = ZERO
            intensity = Rdata(5)
         end if
         table(1,ntab) = table(1,ntab) + Rdata(1)
         table(2,ntab) = table(2,ntab) + ONE
         table(3,ntab) = table(3,ntab) + Rdata(5)
         table(4,ntab) = ZERO
         table(5,ntab) = table(5,ntab) + Rdata(6)
         table(6,ntab) = table(6,ntab) + Rdata(7)
C
         if(intensity.gt.maxint) maxint = intensity
      end do
C
      do itab = 1,ntab
         table(1,itab) = table(1,itab)/table(2,itab)
         table(4,itab) = 100.0d+00*table(3,itab)/maxint
         table(5,itab) = table(5,itab)/table(2,itab)
         table(6,itab) = table(6,itab)/table(2,itab)
      end do
C
      write(IW,2999)
      write(IW,3000) omega
      write(IW,3001)
      write(IW,3002)
      write(IW,3001)
      do itab = 1,ntab
         write(iw,3003) (table(i,itab),i=1,6)
      end do
      write(IW,3001)
      write(IW,2999)
C
 2999 FORMAT(/,1x,67('*'))
 3000 FORMAT(/,
     & 1x,'Raman Intensity at omega = ',F10.6,/,
     & 5x,'Intensity expressed in [Ang.^4/AMU]')
 3001 FORMAT(2x,64('-'))
 3002 FORMAT(
     & 2x,'  Freq  |Mult|  Intensity     (%)  |',
     &    'l-depol ratio|n-depol ratio',/,
     & 2x,' [cm^-1]|    |                     |',
     &    '             |')
 3003 FORMAT(2x,F8.1,'|',F4.0,'|',F13.2,1x,'(',F5.1,')',2('|',F13.6))
      end
C
C*MODULE TDXPRP  *DECK PrintHyperRaman
      subroutine PrintHyperRaman(iw,listdbdq,nmodes,listfreq,omega)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision listdbdq(27,nmodes),listfreq(nmodes),omega
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lTable
C     Body of the subprogram
      call NEWMEM(lTable,6*nmodes)
      call PrintHyperRamanX(iw,listdbdq,nmodes,listfreq,omega,
     &                      XX(lTable))
      call DELETEMEM(lTable,6*nmodes)
C
      end
C
C*MODULE TDXPRP  *DECK PrintHyperRamanX
      subroutine PrintHyperRamanX(iw,listdbdq,nmodes,listfreq,omega,
     &                            table)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision listdbdq(27,nmodes),listfreq(nmodes),omega
      double precision table(6,nmodes)
C     Variables and constants
      logical samemode
      integer imode,i,itab,ntab
      double precision dbdq(27),hRdata(8)
      double precision last_f,last_i,diff_f,diff_i,intensity,maxint
      double precision ZERO,ONE,THRS_I,THRS_F
      parameter (ZERO=0.0d+00,ONE=1.0d+00)
      parameter (THRS_I=0.1d+00,THRS_F=0.1d+00)
C     Called subprogram
      double precision abs
C     Body of the subprogram
      ntab = 0
      last_f = ZERO
      last_i = ZERO
      maxint = ZERO
      intensity = ZERO
      do imode = 1,nmodes
         do i = 1,27
            dbdq(i) = listdbdq(i,imode)
         end do
         call InitHyperRamanData(hRdata,dbdq,listfreq(imode))
         diff_f = abs(hRdata(1)-last_f)
         diff_i = abs(hRdata(6)-last_i)
         last_f = hRdata(1)
         last_i = hRdata(6)
         samemode = (diff_i.le.THRS_I).and.(diff_f.lt.THRS_F)
         if(samemode) then
            intensity = intensity + hRdata(6)
         else
            ntab = ntab+1
            table(1,ntab) = ZERO
            table(2,ntab) = ZERO
            table(3,ntab) = ZERO
            table(4,ntab) = ZERO
            table(5,ntab) = ZERO
            table(6,ntab) = ZERO
            intensity = hRdata(6)
         end if
         table(1,ntab) = table(1,ntab) + hRdata(1)
         table(2,ntab) = table(2,ntab) + ONE
         table(3,ntab) = table(3,ntab) + hRdata(6)
         table(4,ntab) = ZERO
         table(5,ntab) = table(5,ntab) + hRdata(7)
         table(6,ntab) = table(6,ntab) + hRdata(8)
C
         if(intensity.gt.maxint) maxint = intensity
      end do
C
      do itab = 1,ntab
         table(1,itab) = table(1,itab)/table(2,itab)
         table(4,itab) = 100.0d+00*table(3,itab)/maxint
         table(5,itab) = table(5,itab)/table(2,itab)
         table(6,itab) = table(6,itab)/table(2,itab)
      end do
C
      write(IW,2999)
      write(IW,3000) omega
      write(IW,3001)
      write(IW,3002)
      write(IW,3001)
      do itab = 1,ntab
         write(iw,3003) (table(i,itab),i=1,6)
      end do
      write(IW,3001)
      write(IW,2999)
C
 2999 FORMAT(/,1x,67('*'))
 3000 FORMAT(/,
     & 1x,'Hyper Raman Intensity at omega = ',F10.6,/,
     & 5x,'Intensity expressed in [Ang.^6 AMU^-1 StatVolt^-2]')
 3001 FORMAT(2x,64('-'))
 3002 FORMAT(
     & 2x,'  Freq  |Mult|  Intensity     (%)  |',
     &    'n-depol ratio|p-depol ratio',/,
     & 2x,' [cm^-1]|    |                     |',
     &    '             |')
 3003 FORMAT(2x,F8.1,'|',F4.0,'|',F13.2,1x,'(',F5.1,')',2('|',F13.6))
C
      end
C
C*MODULE TDXPRP  *DECK GetAtomAndDir
      subroutine GetAtomAndDir(ixyz,iat,idir,dir)
      implicit none
C     Arguments
      integer ixyz
      integer iat,idir
      character*1 dir
C     Variables and Constants
      character*1 dirs(3)
      data dirs /'x','y','z'/
C     Body of the subprogram
      iat=(ixyz-1)/3+1
      idir=(ixyz-1)-(iat-1)*3+1
      dir=dirs(idir)
      end
C
C*MODULE TDXPRP  *DECK PrintdMudX
      subroutine PrintdMudX(iw,dmdx,nxyz,title)
      implicit none
C     Arguments
      integer iw,nxyz
      double precision dmdx(3,nxyz)
      character*(*) title
C     Variables
      integer ixyz,iat,idir
      character*1 dir
C     Body of the subprogram
      write(iw,9001) title
      do ixyz=1,nxyz
         call GetAtomAndDir(ixyz,iat,idir,dir)
         write(iw,9002) ixyz,dir,iat
         call PrintMu(iw,dmdx(1,ixyz))
      end do
 9001 format(/,2X,A)
 9002 format(2X,'param.',I4,':',1('d',A1,'(',I3,')'))
      end
C
C*MODULE TDXPRP  *DECK PrintdMudQ
      subroutine PrintdMudQ(iw,wa,dmdq,nmodes,title)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision wa(nmodes),dmdq(3,nmodes)
      character*(*) title
C     Variables
      integer imode
      double precision freq
      double precision FACT
      parameter (FACT=219474.6313710d+00)
C     Body of the subprogram
      write(iw,9001) title
      do imode=1,nmodes
         freq = wa(imode)*FACT
         write(iw,9002) imode,freq
         call PrintMu(iw,dmdq(1,imode))
      end do
 9001 format(/,2X,A)
 9002 format(2X,'mode ',I4, '(w=',F6.1,' cm^-1)')
      end
C
C*MODULE TDXPRP  *DECK PrintdAlphadX
      subroutine PrintdAlphadX(iw,dadx,w,nxyz,title)
      implicit none
C     Arguments
      integer iw,nxyz
      double precision dadx(9,nxyz),w(2)
      character*(*) title
C     Variables
      integer ixyz,iat,idir
      character*1 dir
C     Body of the subprogram
      write(iw,9001) title,w(1),w(2)
      do ixyz=1,nxyz
         call GetAtomAndDir(ixyz,iat,idir,dir)
         write(iw,9002) ixyz,dir,iat
         call PrintAlpha(iw,dadx(1,ixyz))
      end do
 9001 format(/,2X,A,'(',F10.6,';',F10.6,')')
 9002 format(2X,'param.',I4,':',1('d',A1,'(',I3,')'))
      end
C
C*MODULE TDXPRP  *DECK PrintdAlphadQ
      subroutine PrintdAlphadQ(iw,wa,dadq,w,nmodes,title)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision wa(nmodes),dadq(9,nmodes),w(2)
      character*(*) title
C     Variables
      integer imode
      double precision Rdata(7)
      character*25 buffer
C     Body of the subprogram
      write(iw,9001) title,w(1),w(2)
      do imode=1,nmodes
         call InitRamanData(Rdata,dadq(1,imode),wa(imode))
         write(buffer,9002) imode,Rdata(1)
         call PrintAlphaX(iw,dadq(1,imode),w,buffer)
         write(iw,9003) Rdata(3),Rdata(4),Rdata(5),Rdata(6),Rdata(7)
C
      end do
 9001 format(/,2X,A,'(',F10.6,';',F10.6,')')
 9002 format(2X,'mode ',I4, '(',F6.1,' cm^-1)')
 9003 format(/,2X,'Raman',
     &       /,3X,'Ak=',F12.6,' Gk=',F12.6,
     &       /,3X,'Intensity=',F10.4,
     &       /,3X,'l-depolarization ratio=',F15.4,
     &       /,3X,'n-depolarization ratio=',F15.4)
      end
C
C*MODULE TDXPRP  *DECK PrintdBetadX
      subroutine PrintdBetadX(iw,dbdx,w,nxyz,title)
      implicit none
C     Arguments
      integer iw,nxyz
      double precision dbdx(27,nxyz),w(3)
      character*(*) title
C     Variables
      integer ixyz,iat,idir
      character*1 dir
C     Body of the subprogram
      write(iw,9001) title,w(1),w(2),w(3)
      do ixyz=1,nxyz
         call GetAtomAndDir(ixyz,iat,idir,dir)
         write(iw,9002) ixyz,dir,idir
         call PrintBeta(iw,dbdx(1,ixyz))
      end do
 9001 format(/,2X,A,'(',F10.6,';',F10.6,',',F10.6,')')
 9002 format(2X,'param.',I4,':',1('d',A1,'(',I3,')'))
      end
C
C*MODULE TDXPRP  *DECK printdBetadQ
      subroutine PrintdBetadQ(iw,wa,dbdq,w,nmodes,title)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision wa(nmodes),dbdq(27,nmodes),w(3)
      character*(*) title
C     Variables
      integer imode
      double precision Rdata(8),mu(3)
      character*25 buffer
C     Called subprogram
      logical GotDip
C     Body of the subprogram
      mu(1) = 0.0d+00
      write(iw,9001) title,w(1),w(2),w(3)
      call GetDipole(mu,.true.,GotDip)
      if(.not.GotDip) then
         call abrtx("PRINDBETADQ")
      end if
      do imode=1,nmodes
         call InitHyperRamanData(Rdata,dbdq(1,imode),wa(imode))
         write(buffer,9002) imode,Rdata(1)
         call PrintBetaX(iw,dbdq(1,imode),mu,w,buffer)
         write(iw,9003) Rdata(3),Rdata(4),Rdata(5),Rdata(6),Rdata(7),
     &                  Rdata(8)
C
      end do
 9001 format(/,2X,A,'(',F10.6,';',F10.6,',',F10.6,')')
 9002 format(2X,'mode ',I4, '(',F6.1,' cm^-1)')
 9003 format(/,2X,'hyper-Raman',
     &       /,3X,'Biii^2=',F12.6,' Bijj^2=',F12.6,' Bijk^2=',F12.6,
     &       /,3X,'Intensity=',F10.4,
     &       /,3X,'l-depolarization ratio=',F15.4,
     &       /,3X,'n-depolarization ratio=',F15.4)
      end
C
C*MODULE TDXPRP  *DECK Printd2MudX2
      subroutine Printd2MudX2(iw,d2mdx2,nxyz,title)
      implicit none
C     Arguments
      integer iw,nxyz
      double precision d2mdx2(3,*)
      character*(*) title
C     Variables
      integer ixyz,jxyz,ijxyz,iat,jat,idir,jdir
      character*1 dirj,diri
C     Body of the subprogram
      write(iw,9001) title
      ijxyz=0
      do jxyz=1,nxyz
         call GetAtomAndDir(jxyz,jat,jdir,dirj)
         do ixyz=1,jxyz
            call GetAtomAndDir(ixyz,iat,idir,diri)
            ijxyz = ijxyz+1
            write(iw,9002) ijxyz,diri,iat,dirj,jat
            call PrintMu(iw,d2mdx2(1,ijxyz))
         end do
      end do
 9001 format(/,2X,A)
 9002 format(2X,'param.',I4,':',2('d',A1,'(',I3,')'))
      end
C
C*MODULE TDXPRP  *DECK Printd2MudQ2
      subroutine Printd2MudQ2(iw,wa,d2mdq2,nmodes,title)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision wa(nmodes),d2mdq2(3,*)
      character*(*) title
C     Variables
      integer imode,jmode,ijmode
      double precision freqi,freqj
      double precision FACT
      parameter (FACT=219474.6313710d+00)
C     Body of the subprogram
      write(iw,9001) title
      ijmode=0
      do jmode=1,nmodes
         freqj = wa(jmode)*FACT
         do imode=1,jmode
            ijmode=ijmode+1
C
            freqi = wa(imode)*FACT
            write(iw,9002) imode,freqi,jmode,freqj
            call PrintMu(iw,d2mdq2(1,ijmode))
         end do
      end do
 9001 format(/,2X,A)
 9002 format(2X,2('dQ(',I4, ')(w=',F6.1,' cm^-1)'))
      end
C
C*MODULE TDXPRP  *DECK Printd2AlphadX2
      subroutine Printd2AlphadX2(iw,d2adx2,w,nxyz,title)
      implicit none
C     Arguments
      integer iw,nxyz
      double precision d2adx2(9,*),w(2)
      character*(*) title
C     Variables
      integer ixyz,jxyz,ijxyz,iat,jat,idir,jdir
      character*1 dirj,diri
C     Body of the subprogram
      write(iw,9001) title,w(1),w(2)
      ijxyz=0
      do jxyz=1,nxyz
         call GetAtomAndDir(jxyz,jat,jdir,dirj)
         do ixyz=1,jxyz
            call GetAtomAndDir(ixyz,iat,idir,diri)
            ijxyz = ijxyz+1
            write(iw,9002) ijxyz,diri,iat,dirj,jat
            call PrintAlpha(iw,d2adx2(1,ijxyz))
         end do
      end do
 9001 format(/,2X,A,'(',F10.6,';',F10.6,')')
 9002 format(2X,'param.',I4,':',2('d',A1,'(',I3,')'))
      end
C
C*MODULE TDXPRP  *DECK Printd2AlphadQ2
      subroutine Printd2AlphadQ2(iw,wa,d2adq2,w,nmodes,title)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision wa(nmodes),d2adq2(9,*),w(2)
      character*(*) title
C     Variables
      integer imode,jmode,ijmode
      double precision freqi,freqj
      double precision FACT
      parameter (FACT=219474.6313710d+00)
      character*50 buffer
C     Body of the subprogram
      write(iw,9001) title,w(1),w(2)
      ijmode=0
      do jmode=1,nmodes
         freqj = wa(jmode)*FACT
         do imode=1,jmode
            ijmode=ijmode+1
C
            freqi = wa(imode)*FACT
            write(buffer,9002) imode,freqi,jmode,freqj
            call PrintAlphaX(iw,d2adq2(1,ijmode),w,buffer)
         end do
      end do
 9001 format(/,2X,A,'(',F10.6,';',F10.6,')')
 9002 format(2X,2('dQ(',I4, ')(w=',F6.1,' cm^-1)'))
      end
C
C*MODULE TDXPRP  *DECK Printd2BetadX2
      subroutine Printd2BetadX2(iw,d2bdx2,w,nxyz,title)
      implicit none
C     Arguments
      integer iw,nxyz
      double precision d2bdx2(27,*),w(3)
      character*(*) title
C     Variables
      integer ixyz,jxyz,ijxyz,iat,jat,idir,jdir
      character*1 dirj,diri
C     Body of the subprogram
      write(iw,9001) title,w(1),w(2),w(3)
      ijxyz=0
      do jxyz=1,nxyz
         call GetAtomAndDir(jxyz,jat,jdir,dirj)
         do ixyz=1,jxyz
            call GetAtomAndDir(ixyz,iat,idir,diri)
            ijxyz = ijxyz+1
            write(iw,9002) ijxyz,diri,iat,dirj,jat
            call PrintBeta(iw,d2bdx2(1,ijxyz))
         end do
      end do
 9001 format(/,2X,A,'(',F10.6,';',F10.6,',',F10.6,')')
 9002 format(2X,'param.',I4,':',2('d',A1,'(',I3,')'))
      end
C
C*MODULE TDXPRP  *DECK Printd2BetadQ2
      subroutine Printd2BetadQ2(iw,wa,d2bdq2,w,nmodes,title)
      implicit none
C     Arguments
      integer iw,nmodes
      double precision wa(nmodes),d2bdq2(27,*),w(3)
      character*(*) title
C     Variables
      integer imode,jmode,ijmode
      double precision freqi,freqj
      double precision FACT
      parameter (FACT=219474.6313710d+00)
      double precision mu(3)
      character*50 buffer
C     Called subprogram
      logical GotDip
C     Body of the subprogram
      mu(1) = 0.0d+00
      write(iw,9001) title,w(1),w(2),w(3)
      call GetDipole(mu,.true.,GotDip)
      if(.not.GotDip) then
         call abrtx("PRINTD2BETADQ2")
      end if
      ijmode=0
      do jmode=1,nmodes
         freqj = wa(jmode)*FACT
         do imode=1,jmode
            ijmode=ijmode+1
C
            freqi = wa(imode)*FACT
            write(buffer,9002) imode,freqi,jmode,freqj
            call PrintBetaX(iw,d2bdq2(1,ijmode),mu,w,buffer)
         end do
      end do
 9001 format(/,2X,A,'(',F10.6,';',F10.6,',',F10.6,')')
 9002 format(2X,2('dQ(',I4, ')(w=',F6.1,' cm^-1)'))
C--      end
C--c
C--C*MODULE TDXPRP  *DECK PrintdUpperAdX
C--      subroutine PrintdUpperAdX(iw,dupperadx,w,nxyz,title)
C--      implicit none
C--c     Arguments
C--      integer iw,nxyz
C--      double precision dupperadx(27,nxyz),w(2)
C--      character*(*) title
C--c     Variables
C--      integer ixyz,iat,idir
C--      character*1 dir
C--c     Body of the subprogram
C--      write(iw,9001) title,w(1),w(2)
C--      do ixyz=1,nxyz
C--         call GetAtomAndDir(ixyz,iat,idir,dir)
C--         write(iw,9002) ixyz,dir,iat
C--         call PrintUpperA(iw,dupperadx(1,ixyz))
C--      end do
C-- 9001 format(/,2X,A,'(',F10.6,';',F10.6,')')
C-- 9002 format(2X,'param.',I4,':',1('d',A1,'(',I3,')'))
C--      end
C--c
C--C*MODULE TDXPRP  *DECK PrintdUpperAdQ
C--      subroutine PrintdUpperAdQ(iw,wa,dupperadq,w,nmodes,title)
C--      implicit none
C--c     Arguments
C--      integer iw,nmodes
C--      double precision wa(nmodes),dupperadq(27,nmodes),w(2)
C--      character*(*) title
C--c     Variables and constants
C--      character*25 buffer
C--      integer imode
C--      double precision freqi
C--      double precision FACT
C--      parameter (FACT=219474.6313710d+00)
C--c     Body of the subprogram
C--      write(iw,9001) title,w(1),w(2)
C--      do imode=1,nmodes
C--         freqi = wa(imode)*FACT
C--         write(buffer,9002) imode,freqi
C--         call PrintUpperaX(iw,dupperadq(1,imode),w,buffer)
C--      end do
C--      return
C-- 9001 format(/,2X,A,'(',F10.6,';',F10.6,')')
C-- 9002 format(2X,'mode ',I4, '(',F6.1,' cm^-1)')
      return
      end
! Local Variables:
! mode:Fortran
! End:
