********************************************************************************
*                                                                              *
*  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 ### dointerp ###	
c
c	Performs interpolations of orbitals and potentials between grids. 
c       Orbitals and potentials are treated separately.
c
c work array is not needed any more

      subroutine dointerp (ic,nmuall_p,nmuall,fbefore,fafter) 
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)	
      real*8 vmuq
      include 'commons8.inc'
      include 'commons16.inc'
      dimension fbefore(nni_p,*),fafter(nni,*)

c  The previous grid containes nni_p point in \nu (\eta) direction and
c  nmuall_pt point in \mu (\xi) direction

      if ((nmuall.ne.nmuall_p.or.
     &     ngrids.ne.ngrids_p.or.
     &     rinf.ne.rinf_p).and.nni.eq.nni_p) then

c  prepare data for routines evaluating coefficients of the Lagrange
c  polynomial

         if (ic.eq.1) then

c  Interpolate in mu variable      

c  Interpolation of orbitals is done separately in two regions in mu 
c  variable. The inner region extends from mu=1 to the last maximum
c  and the outer (tail) region extends from that maximum to the practical 
c  infinity. The values in the latter are logarithmed before 
c  the interpolation is carried out.

c  Find the tail region: the maximum is sought along the line which is 
c  almost perpendicular to the intermolecular axis (to avoid picking up
c  zeros for ungerade orbitals)

c  The above scheme turned out to be unnecessary

            iord=iord_mu_orb
            kbeg=1
            kend=iord
            iord2=iord/2
            do imu=1,nmuall_p
               vmuq(imu)=vmu_p(imu)
            enddo
            
            nall_p=nmuall_p
            nall  =nmuall
            nstart_p  =1
            nstart    =1
            call dointerp_mu (nni_p,nstart_p,nall_p,nstart,nall,
     &           fbefore,fafter) 

         elseif (ic.eq.2) then
            iord=iord_mu_coul
            kbeg=1
            kend=iord
            iord2=iord/2
            do imu=1,nmuall_p
               vmuq(imu)=vmu_p(imu)
            enddo
            
            ninner_p=1
            ninner  =1
            nall_p=nmuall_p
            nall  =nmuall
            call dointerp_mu (nni_p,ninner_p,nall_p,ninner,nall,
     &           fbefore,fafter) 

         elseif (ic.eq.3) then
            iord=iord_mu_exch
            kbeg=1
            kend=iord
            iord2=iord/2
            do imu=1,nmuall_p
               vmuq(imu)=vmu_p(imu)
            enddo

            ninner_p=1
            ninner  =1
            nall_p=nmuall_p
            nall  =nmuall
            call dointerp_mu (nni_p,ninner_p,nall_p,ninner,nall,
     &           fbefore,fafter) 
         endif

c        Extrapolation (needed when R_infy is being increased seems to 
c        cause degradation of accuracy. That is why these values are 
c        taken equal (for each nu) to those corresponding to nmuall_p, 
c        i.e tp the values of the last column of fbefore

         imu_bext=nmuall
         do imu=1,nmuall
            if (vmu(imu).ge.vmu_p(nmuall_p)) then
               imu_bext=imu
               goto 1234
            endif
         enddo
01234    continue
         do in=1,nni_p
            do im=imu_bext,nmuall
               fafter(in,im)=fbefore(in,nmuall_p)
            enddo
         enddo
         
      elseif (nni.ne.nni_p.and.
     &        (nmuall.eq.nmuall_p.and.
     &	           ngrids.eq.ngrids_p.and.rinf.eq.rinf_p)) then

c        interpolate in nu variable      

         if (nni_p.gt.mxnmu) then
            print *,"dointerp: error! array vmuq is too short"
            stop
         endif

         do i=1,nni_p
            vmuq(i)=vni_p(i)
         enddo

         if (ic.eq.1) then
            iord=iord_nu_orb
            kbeg=1
            kend=iord
            iord2=iord/2
c            stop "dointerp: error check array dimension"
            call dointerp_nu (nmuall_p,fbefore,fafter) 
         elseif (ic.eq.2) then
            iord=iord_nu_coul
            kbeg=1
            kend=iord
            iord2=iord/2
            call dointerp_nu (nmuall_p,fbefore,fafter) 
         elseif (ic.eq.3) then
            iord=iord_nu_exch
            kbeg=1
            kend=iord
            iord2=iord/2
            call dointerp_nu (nmuall_p,fbefore,fafter) 
         endif
      else
         write(*,'(/,1x,"dointerp: error! cannot perform interpolation")')
	 stop "dointerp"   
      endif
      return
      end