********************************************************************************
*                                                                              *
*  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 ### prepmesh ####
c
c     Establishes meshes for each subgrid and determines coefficients
c     of the interpolation polynomials needed in fill arrays.
c
c     Warning! At present at most 3 grids are allowed. 
c     See mesh for further detailes.
c
      subroutine prepmesh (cw_sor)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      integer*4 cw_sor	
      include 'commons8.inc'
      dimension cw_sor(*)

c     Array cw_sor is used to store index arrays needed to perform
c     sor relaxation on normal and extended grids. Two cases have to be
c     considered since all ngrids-1 subgrids can either be interior or 
c     last ones depending on the particular orbital or potential
c     function relaxed (outer diffused orbitals would need, say, a grid
c     consisting of the 3 subgrids while for the proper description of the 
c     the core, very much confined orbitals only one subgrid would sufice).

c     Division of the cw_sor array is as folows: 

c     ind1   - relaxation ordering on the extended grid
c     ind2   - relaxation ordering on the normal grid 
c     ind1ex - extrapolation ordering on the extended grid 
c     ind2ex - extrapolation ordering on the normal grid 
c     ind3   - number of points for relaxation and extrapolation
c     mxsize8==(nni+8)*(mxnmu+8)
c     mxnmu8=mxnmu+8

c     case 1 - all ngrids subgrids are considered as the last ones; 
c              in effect the last 4 columns of orbitals and potentials 
c	       are not relaxed two arrays of the indx1 type

c     case 2 - all (ngrids-1) subgrids are considered as the interior ones

c     cw_sor consists of
c     - ngrids     arrays of indx1 type (case 1);
c     - (ngrids-1) arrays of indx1 type (case 2);
c                  the length of each array is determined by ngsize 
c
c     - ngrids     arrays of indx2 type (case 1);
c     - (ngrids-1) arrays of indx2 type (case 2);
c		   the length of each array is determined by ngsize 
c
c     - ngrids	   arrays of indx6a (case 1)
c     - (ngrids-1) arrays of indx6a (case 2)
c     - ngrids	   arrays of indx6b (case 1)
c     - (ngrids-1) arrays of indx6b (case 2)
c     - ngrids	   arrays of indx7 (case 1)
c     - (ngrids-1) arrays of indx7 (case 2)
c  		   the length of each array is determined by nmu

c     ingr1(k,ig) and ingr2(k,ig),ig=1,ngrids, k=1,..,4, contain number of 
c     grids points to be relaxed in each subgrid (k=1) and number of points 
c     for which extrapolation is to be carried out (k=2 - ngrd6a, k=3 - ngrd6b, 
c     k=4 - ngrd7); ingr1 and ingr2 correspond to the extended and normal 
c     case, respectively.

c     i6b(1) - begining of ind1
c     i6b(2) - begining of ind2
c     i6b(3) - begining of index1 (indx6a)
c     i6b(4) - begining of index2 (indx6b)
c     i6b(5) - begining of index3 (indx7)

c     more space is reserved for the ordering arrays to avoid problems when
c     more than 3 subgrid would be used in the future

c     determine the largest subgrid length in mu variable

      mxsmu=0
      do ig=1,ngrids
         if (mxsmu.lt.nmu(ig)) mxsmu=nmu(ig)
      enddo
	
      is34=ngrids*mxsmu	  
      is5 =ngrids*nni
      mxnmu16=mxnmu+16
      ngrid16=nni*mxnmu16
      
      i6b(1)=1
      i6b(2)=i6b(1)+2*ngrid16
      i6b(3)=i6b(2)+2*ngrid16
      i6b(4)=i6b(3)+2*is34
      i6b(5)=i6b(4)+2*is5

c     i6b(4)=i6b(3)+is34
c     i6b(5)=i6b(4)+is5

c     determine ordering of grid points assuming the subgrids are the last
c     ones (grid points in the last 4 columns are excluded from relaxation
c     process

      imcase=1
      ibeg1=i6b(1)
      ibeg2=i6b(2)	
      ibeg3=i6b(3)		
      ibeg4=i6b(4)
      ibeg5=i6b(5)	
      imut=0		
      init=0	
      igs=0
      do ig=1,ngrids
         ib1=ibeg1+igs
         ib2=ibeg2+igs
         ib3=ibeg3+imut
         ib4=ibeg4+imut	  
         ib5=ibeg5+init	  	  
         iadext(ig)=ib1	  
         iadnor(ig)=ib2	  	  
         iadex1(ig)=ib3	  	  	  
         iadex2(ig)=ib4	  	  	  
         iadex3(ig)=ib5	  	  	  	  	  
         call mesh(imcase,ig,cw_sor(ib1),cw_sor(ib2),
     &        cw_sor(ib3),cw_sor(ib4),cw_sor(ib5)) 	
         imut=imut+nmu(ig)
         init=init+nni     
         igs=igs+ngsize(ig)
      enddo	  

c     determine ordering of grid points assuming the subgrids are 
c     the interior ones (only boundary points coressponding to 
c     mu=1, ni=1 and ni=nni are excluded from the relaxation process)

      imcase=2
      ibeg1=i6b(1)+ngrid16
      ibeg2=i6b(2)+ngrid16	
      ibeg3=i6b(3)+is34
      ibeg4=i6b(4)+is34			
      ibeg5=i6b(5)+is5				
      imut=0	
      init=0			
      igs=0
      do ig=1,ngrids-1
         ib1=ibeg1+igs
         ib2=ibeg2+igs
         ib3=ibeg3+imut
         ib4=ibeg4+imut	  
         ib5=ibeg5+init	  	  	  
         iadext(ig+ngrids)=ib1	  
         iadnor(ig+ngrids)=ib2	  	  	  
         iadex1(ig+ngrids)=ib3	  	  	  
         iadex2(ig+ngrids)=ib4	  	  	  
         iadex3(ig+ngrids)=ib5	  	  	  	  	  	  
         call mesh(imcase,ig,cw_sor(ib1),cw_sor(ib2),
     &        cw_sor(ib3),cw_sor(ib4),cw_sor(ib5)) 	     
         imut=imut+nmu(ig)
         init=init+nni          
         igs=igs+ngsize(ig)
      enddo	  

c     determine coefficients of interpolation polynomials

      if (ngrids.gt.1) call interpolq

      return 
      end