program mcdv-divide
      character*80 s80
c
      write(6,600)
600   format(' Devides gaussian output with numerical differentiation',/,
     1       ' to individual geometries',/,
     2       '    GV.OUT - only one calculation (# .... #) possible.',/)
c
      open(9,file='GV.OUT',status='old')
2     read(9,900,end=88,err=88)s80
900   format(a80)

      subroutine rdd(f,d,nat,nroot,u0,iwr)
c     read transition dipole derivatives from excited state freq calc
      implicit none
      character*(*) f
      real*8 d(*),u(3),u0(*),step,sign,stepau
      integer*4 nat,ln,ia,xa,i,nroot,xar,iar,ibr,ib,ic,ii,ix,iwr
      character*80 s80
      open(9,file=f,status='old')
      ln=0
c     where it starts:
1     read(9,900,end=88,err=88)s80
900   format(a80)
      ln=ln+1
      if(s80(1:2).eq.' #')then
       ic=0
       do 3 i=1,len(s80)-3
       if(s80(i:i+4).eq.' freq'.or.s80(i:i+4).eq.' freq')ic=ic+1
3      if(s80(i:i+2).eq.' td'  .or.s80(i:i+2).eq.' TD'  )ic=ic+1
       if(ic.eq.2)then
        write(6,*)'TD starts at line ',ln
        goto 2
       endif
      endif
      goto 1
88    close(88)
      call report('TD freq not found')

2     ia=1
      xa=1
      ib=0
      ii=0
      step=0.001d0
      do 7 ix=1,9*nat
7     d(ix)=0.0d0
4     read(9,900,end=78,err=78)s80
      if(s80(2:14).eq.'Nuclear step=')then
       write(6,*)s80(1:23)
       read(s80(15:23),*)step
      endif
      if(s80(2:24).eq.'Re-enter D2Numr: IAtom=')then
       read(s80(25:27),*)iar
       read(s80(34:34),*)xar
       read(s80(42:43),*)ibr
       if(ia.ne.iar)call report('ia <> iar')
       if(xa.ne.xar)call report('xa <> xar')
       if(ib.ne.ibr)call report('ib <> ibr')
      endif
      if(s80(2:65).eq.'Ground to excited state transition electric dipol
     1e moments (Au):')then
       if(ii.eq.0)then
        write(6,*)'Zero point'
        do 5 i=1,nroot
5       read(9,*)
        read(9,*)u0(1),(u0(ix),ix=1,3)
       else
        ib=ib+1
        if(ib.gt.2)then
         ib=1
         xa=xa+1
         if(xa.gt.3)then
          xa=1
          ia=ia+1
         endif
        endif
        if(iwr.gt.1)write(6,*)' d ',ia,xa,ib,iar,xar,ibr
        do 6 i=1,nroot
6       read(9,*)
        read(9,*)u(1),(u(ix),ix=1,3)
        if(ib.eq.1)then
         sign=1.0d0
        else
         sign=-1.0d0
        endif
        do 8 ix=1,3
8       d(ix+3*(xa-1)+9*(ia-1))=d(ix+3*(xa-1)+9*(ia-1))+sign*u(ix)
       endif
       ii=ii+1
       if(ia.eq.nat.and.ix.eq.3.and.ib.eq.2)goto 78
      endif
      goto 4
78    close(9)
      stepau=step/0.529177d0
      do 9 ix=1,9*nat
9     d(ix)=d(ix)/stepau
      open(40,file='DE.TEN')
      write(40,*)'atom coord dipx        dipy        dipz'
      do 10 ia=1,nat
      do 10 xa=1,3
 10   write(40,400)ia,xa,(d(ix+3*(xa-1)+9*(nat-1)),ix=1,3)
400   format(i5,i2,3g12.4)
      close(40)

      return
      end
c     ==============================================================
      SUBROUTINE READS(N3,S,E,NQ)
      IMPLICIT none
      integer*4 N3,NQ,N,nat3,nat,i,J,ix
      real*8 S(N3,N3),E(*),CM,SFAC
      CM=219470.0d0
      SFAC=0.0234280d0
C     1 a.u.= 2.1947E5 cm^-1
      open(4,file='F.INP')
      read(4,*)NQ,nat3,nat
      N=NAT3
      do 1 i=1,NAT
1     read(4,*)
      read(4,*)
      DO 2 I=1,NAT
      DO 2 J=1,NQ
2     read(4,*)(s(3*(i-1)+ix,J),ix=1,2),
     1(s(3*(i-1)+ix,J),ix=1,3)
      read(4,*)
      READ(4,4000)(E(I),I=1,NQ)
4000  FORMAT(6F11.6)
      close(4)
      DO 3 I=1,NQ
3     E(I)=E(I)/CM
      DO 5 I=1,N
      DO 5 J=1,N
5     S(I,J)=S(I,J)*SFAC
      RETURN
      end
c     ==============================================================
      subroutine trafo(dd,ddi,nat,N)
      implicit none
      integer*4 nat,NQ,N,ix,iq,ia,xa
      real*8 dd(*),ddi(*)
      real*8,allocatable::e(:),s(:,:)
      allocate(e(3*nat),s(3*nat,3*nat))
      call READS(3*nat,s,e,NQ)
      if(N.ne.NQ)call report('N <> NQ')
      do 1 ix=1,3
      do 1 iq=1,N
      ddi(ix+3*(iq-1))=0.0d0
      do 1 ia=1,nat
      do 1 xa=1,3
c     s-opposite order of normal modes:
1     ddi(ix+3*(iq-1))=ddi(ix+3*(iq-1))+dd(ix+3*(xa-1)+9*(ia-1))
     1*s(xa+3*(ia-1),N-iq+1)
      write(6,*)'dipoles transformed to normal modes'
      return
      end
c     ==============================================================
      subroutine zm(M,N)
      implicit none
      integer*4 N,i,j
      real*8 M(N,N)
      do 1 i=1,N
      do 1 j=1,N
1     M(i,j)=0.0d0
      return
      end
c     ==============================================================
      subroutine zv(J,N)
      implicit none
      real*8 J(*)
      integer*4 N,i
      do 1 i=1,N
1     J(i)=0.0d0
      return
      end