********************************************************************************
*                                                                              *
*  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 ### grid_old ###
c
c     Initializes grid info to enable interpolation.
c
      subroutine grid_old
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'

      ibmu_p(1)=1
      iemu_p(1)=nmu_p(1)
      mxnmu_p=nmu_p(1)
      ioffs_p(1)=0
      do ig=2,ngrids_p
         mxnmu_p=mxnmu_p+nmu_p(ig)-1
         ibmu_p(ig)=iemu_p(ig-1)	 	 
         iemu_p(ig)=mxnmu_p
         ioffs_p(ig)=nni_p*(iemu_p(ig-1)-1)
      enddo

c     determine the size of each grid
      
      do ig=1,ngrids_p
         ngsize_p(ig)=nni_p*nmu_p(ig)
      enddo

c     total no of grid points 

      mxsize_p=nni_p*mxnmu_p

c     hni - step in ni variable
c     hmu - step in mu variable

c     determine step size in ni variable 

      hni_p=pii/dble(nni_p-1)

c     determine step size in mu variable for each grid

      if (ngrids_p.eq.1.and.rgrid_p(1).ge.0.5d0) then
         xmi0=2.d0*rgrid_p(1)/r_p
         xmi0=log(xmi0+sqrt(xmi0*xmi0-1.d0))
         hmu_p(1)=xmi0/dble(mxnmu_p-1)
      else
         do ig=1,ngrids_p
            if(ig.eq.1) then 
               hmu_p(ig)=rgrid_p(ig)
            else
               hmu_p(ig)=hmu_p(1)*rgrid_p(ig)
            endif
         enddo
      endif

c     initialize mu and xi arryas

      vmu_p(1) =0.d0
      vxi_p(1) =cosh(vmu_p(1))
      vxisq_p(1)=vxi_p(1)*vxi_p(1)	    	    
      vxi1_p(1)=sqrt(vxi_p(1)*vxi_p(1)-1.d0)
      vxi2_p(1)=0.0d0

      ib=2
      
      do ig=1,ngrids_p
         ie=iemu_p(ig)
         do i=ib,ie
            vmu_p(i)=vmu_p(i-1)+hmu_p(ig)
            vxi_p(i)=cosh(vmu_p(i))
            vxisq_p(i)=vxi_p(i)*vxi_p(i)	    	    

c           sqrt(vxi*vxi-1)  and  ssf*vxi/sqrt(vxi*vxi-1)
            
            vxi1_p(i)=sqrt(vxi_p(i)*vxi_p(i)-1.d0)
            vxi2_p(i)=vxi_p(i)/vxi1_p(i)
c           enddo over i
         enddo 
         ib=ie+1
c        enddo over ig
      enddo

c     determine rinf
      
      write(iout6,*) 
      write(iout6,*) 'about to interpolate: previous grid info:'
      write(iout6,*) 'subgrid  nni    nmu   hni  hmu      rb '
      rinf_p=r_p*vxi_p(iemu_p(ngrids_p))/2.d0
      do ig=1,ngrids_p
         rinfig_p=r_p*vxi_p(iemu_p(ig))/2.d0
         write(iout6,1000) ig,nni_p,nmu_p(ig),hni_p,hmu_p(ig),rinfig_p
      enddo
01000 format(4x,i2,3x,i4,2x,i4,2f8.5,f8.3)

c    initialize ni and eta arrays

      do i=1,nni_p
         vni_p(i)=dble((i-1))*hni_p
         veta_p(i)=cos(vni_p(i))
         vetasq_p(i)=veta_p(i)*veta_p(i)		

c        sqrt(1-veta*veta) and  veta/sqrt(1-veta*veta)
c        
         veta1_P(i)=sqrt(1.d0-veta_p(i)*veta_p(i))
         if (veta1_p(i).lt.precis) then
            veta2_p(i)=0.0
         else
            veta2_p(i)=veta_p(i)/veta1_p(i)
         endif
      enddo

      return
      end