********************************************************************************
*                                                                              * 
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### rfun_int ###	
c
c     Reads functions from disk in unformatted form and interpolates 
c     to a new grid.
c
      subroutine rfun_int (norb_p,cw_orb,cw_coul,cw_exch,cw_sctch)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)	
      include 'commons8.inc'
      dimension cw_orb(*),cw_coul(*),cw_exch(*),cw_sctch(*)
      dimension i1b_p(60),i2b_p(60),i3b_p(1830),
     &     i1e_p(60),i2e_p(60),i3e_p(1830),
     &     i1si_p(60),i2si_p(60),i3si_p(1830),
     &     i1ng_p(60),i2ng_p(60),i3ng_p(1830),
     &     i1mu_p(60),i2mu_p(60),i3mu_p(1830)	

      read (iinp11,err=1000) i1b_p,i2b_p,i3b_p,i1e_p,i2e_p,i3e_p,
     &     i1si_p,i2si_p,i3si_p,i1ng_p,i2ng_p,i3ng_p,
     &     i1mu_p,i2mu_p,i3mu_p	
      
      ioffset=norb-norb_p
      ica=1
      write(iout6,1100) 
01100 format('... interpolating orbitals:',/6x,$)
01110 format(i4,$)
      
      do i=1,norb_p
         write(iout6,1110) i
         
         ipex=mgx(6,i)
         ipex=ipex-2*(ipex/2)
         if (ipex.eq.0) then
            isym= 1 
         else 
            isym=-1
         endif
         call reada(iinp11,i1si_p(i+ioffset),cw_sctch(i5b(1)),ierr)
         call dointerp     (ica,i1mu_p(i),i1mu(i),cw_sctch(i5b(1)),
     &        cw_orb(i1b(i+ioffset)))
         
         if (ierr.ne.0) then
            write(iout6,*) 'error detected when reading orbital',i
            stop 'inifun'
         endif				   	
      enddo	   

c     retrieve the extra data from the orbital input file

      read(iinp11,end=1010,err=1010) area
      read(iinp11,end=1010,err=1010) eng
      read(iinp11,end=1010,err=1010) engo
      read(iinp11,end=1010,err=1010) cmulti
      read(iinp11,end=1010,err=1010) excdi
      read(iinp11,end=1010,err=1010) excqu
      read(iinp11,end=1010,err=1010) excoc
      read(iinp11,end=1010,err=1010) exche
      read(iinp11,end=1010,err=1010) exc5
      read(iinp11,end=1010,err=1010) exc6
      read(iinp11,end=1010,err=1010) exc7
      read(iinp11,end=1010,err=1010) exc8

      write(iout6,*)
      rewind iinp11	
      
      ica=2
      isym=1
      write(iout6,1102) 
01102 format('... interpolating Coulomb potentials:',
     &       /6x,$)
      do i=1,norb_p
         write(iout6,1110) i
         call reada(iinp12,i2si_p(i+ioffset),cw_sctch(i5b(1)),ierr)
         call dointerp     (ica,i2mu_p(i),i2mu(i),cw_sctch(i5b(1)),
     &        cw_coul(i2b(i+ioffset)),cw_sctch(i5b(2)) )
         
         if (ierr.ne.0) then
            write(iout6,*) 'error detected when reading coulomb',
     &           ' potential',i
            stop 'inifun'
         endif				   	
      enddo	   
      write(iout6,*)
      rewind iinp12	
      
      if (imethod.eq.2) then
         write(*,*) 'inifun: disk file without extension has been ',
     &        'retrieved '
         write(*,*) '      '
         rewind iinp13
         return
      endif
      
      if (imethod.eq.1) then
         ica=3

c        to interpolate exchange potentials they have to be read in as a
c        single file
         if (iform.eq.0.or.iform.eq.2) then
            write(*,'("rfun_int: cannot interpolate exchange "
     &        "potentials when they are retrieved as separate files")')
            stop "rfun_int"
         endif

         write(iout6,1104) 
01104    format('... interpolating exchange ',
     &        'potentials for orbitals:',/6x,$)
         do iorb1=1,norb_p
            write(iout6,1110) iorb1 
            
            do iorb2=iorb1,norb_p
               k=iorb1+iorb2*(iorb2-1)/2
               
               idel=iabs(mgx(6,iorb1)-mgx(6,iorb2))
               if (iorb1.eq.iorb2) idel=2*mgx(6,iorb1)
               ipex=idel-2*(idel/2)
               if (ipex.eq.0) then
                  isym= 1  
               else
                  isym=-1
               endif
               
               if (iorb1.eq.iorb2.and.ll(iorb1).eq.0) goto 50
               call reada(iinp13,i3si_p(k),cw_sctch(i5b(1)),ierr)
               call dointerp     (ica,i3mu_p(k),i3mu(k),cw_sctch(i5b(1)),
     &              cw_exch(i3b(k)),cw_sctch(i5b(2)) )
               
               if (ierr.ne.0) then
                  write(iout6,*) 'error detected when reading exchange',
     &                 ' potential',iorb1,iorb2,k
                  stop 'rfun_int'
               endif				   	
               if (iorb1.eq.iorb2) goto 50
               if (ll(iorb1).eq.0.or.ll(iorb2).eq.0) goto 50
               call reada(iinp13,i3si_p(k),cw_sctch(i5b(1)),ierr)
               call dointerp     (ica,i3mu_p(k),i3mu(k),cw_sctch(i5b(1)),
     &              cw_exch(i3b(k)+i3si(k)),cw_sctch(i5b(2)) )
               if (ierr.ne.0) then
                  write(iout6,*) 'error detected when reading exchange',
     &                 ' potential',iorb1,iorb2,k
                  stop 'rfun_int'
               endif				   	
00050          continue
            enddo
         enddo
         write(iout6,*)
         write(iout6,*)
         rewind iinp13
      endif
      
      return
01000 continue
      write(*,*) '... error detected when retrieving the disk', 
     &     ' file ... '
      stop 'rfun_int'
01010 write(*,*) '... error has been detected when retrieving ',
     &     'extension of the orbital input file ...'
      write(*,*) '... disk file without extension retrieved ... '
      idump=1
      return
      end