********************************************************************************
*                                                                              *
*  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_nu ###	
c
c     Performs interpolations of functions in ni variable.
c
      subroutine dointerp_nu (nmuall,fbefore,fafter) 
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)	
      real*8 vmuq,coeffq(9),coeffq2(9,9),xni,vpoly1q
      include 'commons8.inc'
      include 'commons16.inc'
      dimension fbefore(nni_p,*),fafter(nni,*)

c interpolation for the first iord2 points in ni variable

      do ini=1,nni
         if(vni(ini).ge.vni_p(iord2)) goto 115
      enddo
00115 continue
      nni_first=ini 
      
      do k=1,kend
         call lpcoeffq(iord2+1,k,coeffq)
         do i=1,kend
            coeffq2(i,k)=coeffq(i)
         enddo
      enddo
      
      do ini=1,nni_first
         xni=vni(ini)            
         do imu=1,nmuall
            fafter(ini,imu)=0.d0
            do k=1,kend
               fafter(ini,imu)=fafter(ini,imu)+
     &              fbefore(k,imu)
     &              *vpoly1q(xni,coeffq2(1,k))
            enddo
            if (ini.eq.1.and.(imu.eq.1.or.imu.eq.nmuall)) then
               fafter(ini,imu)=fbefore(ini,imu)
            endif
         enddo
      enddo

c  Interpolation for the last iord2 points in ni variable.
c  Determine the location of the tail region in the new grid

      do ini=1,nni
         if(vni(ini).ge.vni_p(nni_p-iord2+1)) goto 120
      enddo
00120 continue
      if (ini.gt.nni) then
         nni_last=nni
      else
         nni_last=ini
      endif
      
      do k=1,kend
         call lpcoeffq(nni_p-iord2,k,coeffq)
         do i=1,kend
            coeffq2(i,k)=coeffq(i)
         enddo
      enddo
      
      do ini=nni_last-iord2+1,nni
         xni=vni(ini)            
         do imu=1,nmuall
            fafter(ini,imu)=0.d0
            do k=1,kend
               fafter(ini,imu)=fafter(ini,imu)+
     &              fbefore(nni_p-iord+k,imu)
     &              *vpoly1q(xni,coeffq2(1,k))
            enddo
            if (ini.eq.nni.and.(imu.eq.1.or.imu.eq.nmuall)) then
               fafter(nni,imu)=fbefore(nni_p,imu)
            endif
         enddo
      enddo
      
c     interpolation for the inner points in this region

      do ini=nni_first+1,nni_last-iord2
         xni=vni(ini)            
         do i=1,nni_p
            if(vni_p(i).ge.xni) goto 130
         enddo
00130    continue
         if (i.gt.nni_p) then
            ini_p=nni_p
         else
            ini_p=i
         endif
         
         do k=1,kend
            call lpcoeffq(ini_p,k,coeffq)
            do i=1,kend
               coeffq2(i,k)=coeffq(i)
            enddo
         enddo
         
         do imu=1,nmuall
            fafter(ini,imu)=0.d0
            do k=1,kend
               fafter(ini,imu)=fafter(ini,imu)+
     &              fbefore(ini_p-iord2-1+k,imu)
     &              *vpoly1q(xni,coeffq2(1,k))
            enddo
         enddo
      enddo
      return
      end