PROGRAM PESFC
c     spectra generation
      IMPLICIT none
      character*10 s10
      integer*4 ia,na,ib,nb,istart,np,ie,ip,is,jp
      real*8 qmin,qmax,dq,o1,o2,fc
      real*8,allocatable::psi(:),fi(:),psie(:)
      logical lex
      write(6,6000)
6000  format(
     1' Franc-Condon factors',/,
     1'      Input : ',/,
     3'              N.PRN  wavefunctions in folder A, N=1,2,... ',/,
     3'              N.PRN  wavefunctions in folder B, N=1,2,... ',/,/,
     3'      Output: FC.TAB                                  ',/)
      ia=1
3     write(s10,100)ia
100   format(i10)
      do 1 istart=1,len(s10)
1     if(s10(istart:istart).ne.' ')goto 2
2     inquire(file='A/'//s10(istart:len(s10))//'.PRN',exist=lex)
      if(lex)then
       ia=ia+1
       goto 3       
      endif
      na=ia-1
      write(6,*)na,' A-functions'
      ib=1
4     write(s10,100)ib
      do 5 istart=1,len(s10)
5     if(s10(istart:istart).ne.' ')goto 6
6     inquire(file='B/'//s10(istart:len(s10))//'.PRN',exist=lex)
      if(lex)then
       ib=ib+1
       goto 4       
      endif
      nb=ib-1
      write(6,*)nb,' B-functions'
      if(na.eq.0.or.nb.eq.0)call report('not enough functions')

      np=0
      open(10,file='A/1.PRN')
7     read(10,*,end=901,err=901)qmax
      np=np+1
      if(np.eq.1)qmin=qmax
      goto 7
901   close(10)
      if(np.lt.2)call report('not enough points')
      dq=(qmax-qmin)/dble(np-1)
      write(6,6001)np,qmin,qmax,dq
6001  format(i10,' points ',/,
     1       10x,' Qmin: ',f12.5,' Qmax: ',f12.5,' dQ: ',f12.5)


      allocate(psi(np),fi(np),psie(np))
      write(6,*)'A norms:'
      do 9 ip=1,na
      call getpsi('A',psi,np,ip)
9     call testpsi(np,psi,dq)
      write(6,*)
      write(6,*)'B norms:'
      do 10 ip=1,nb
      call getpsi('B',psi,np,ip)
10    call testpsi(np,psi,dq)
      write(6,*)

      write(6,*)
      write(6,*)'AB overlaps:'
      do 8 ip=1,na
      call getpsi('A',psi,np,ip)
      do 8 jp=1,nb
      write(6,609)ip,jp
609   format(2i4,$)
      call getpsi('B',fi,np,jp)
8     call psifi(np,psi,fi,dq,o1,1)

      write(6,*)'<1|B><B|Ap> double overlaps'
      ip=1
      call getpsi('A',psi,np,ip)
      do 11 ie=1,na
      call getpsi('A',psie,np,ie)
      fc=0.0d0
      do 12 is=1,nb
      call getpsi('B',fi,np,is)
      call psifi(np,psi,fi,dq,o1,0)
      call psifi(np,psie,fi,dq,o2,0)
12    fc=fc+o1*o2
11    write(6,608)ip,ie,fc
608   format(2i4,' total FC ',f12.6)
      end
c     ------------------------------------------------------------
      subroutine testpsi(np,psi,dq)
      implicit none
      integer*4 np,j
      real*8 psi(np),dq,a
      a=0.0d0
      do 2 j=1,np
2     a=a+psi(j)**2
      a=a*dq
      write(6,6000)a
6000  format(f15.10,$)
      return
      end
c     ------------------------------------------------------------
      subroutine psifi(np,psi,fi,dq,b,io)
      implicit none
      integer*4 np,j,io
      real*8 psi(np),dq,a,fi(*),b
      a=0.0d0
      do 2 j=1,np
2     a=a+psi(j)*fi(j)
      a=a*dq
      if(io.eq.1)write(6,6000)a
6000  format(f15.10)
      b=a
      return
      end
c     ------------------------------------------------------------
c     ------------------------------------------------------------
      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end
c     ------------------------------------------------------------
c     ------------------------------------------------------------
      subroutine getpsi(ab,psi,np,ip)
      implicit none
      integer*4 np,ip,is,ix
      real*8 psi(np)
      character*5 s5
      character*1 ab
      write(s5,990)ip
990   format(i5)
      do 65 is=1,5
65    if(s5(is:is).ne.' ')goto 651
651   open(81,file=ab//'/'//s5(is:5)//'.PRN')
      do 8 ix=1,np
8     read(81,*)psi(ix),psi(ix)
      close(81)
      return
      end
c     ------------------------------------------------------------