********************************************************************************
*                                                                              *
*  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 ### mesh ###
c
c     Initializes index arrays determining the ordering of mesh points 
c     during sor process
c    
c     type of ordering (itype==iorder):
c     itype=1 -- natural column-wise 
c     itype=2 -- middle
c     itype=3 -- natural row-wise 
c     itype=4 -- reverse row-wise

      subroutine mesh(imcase,ig,indx1,indx2,indx6a,indx6b,indx7)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension indx1(*),indx2(*),indx6a(*),indx6b(*),indx7(*)

c   initialize variables needed in sor and mcsor

      nni1=nni+8
      nni2=nni1+nni1
      nni3=nni2+nni1
      nni4=nni3+nni1
      nni5=nni4+nni1

      nni8=nni+8
      nmut=nmu(ig)
      nmu8=nmut+8

      if (imcase.eq.1) then
         iend1=nmut
         iend2=nmut-4	  
      elseif (imcase.eq.2) then
         iend1=nmut+4
         iend2=nmut	  
      endif

      itype=iorder(ig)
      goto (100,200,300,400,500,600,700,800),itype

00100 continue

c     natural column-wise ordering. bondary points are excluded
c     the upper limit of mu variable is nmu.

c     ?  values for points (j,nmu+1), (j,nmu+2) and (j,nmu+3) are calculated
c     ?  from the asymptotic exspansion (the indexes refer to the larger
c     ?  grid, see routine fill and refill).
c
c     region i
c
c     in the case of more than one grids the upper limit in not nmut-4
c     but nmut

      ngrd1=0
      do i=6,iend1
         ii=(i-1)*nni8
         do j=6,nni+3
            ngrd1=ngrd1+1
c           ires2(j,i)=ngrd1
            indx1(ngrd1)=ii+j
         enddo
      enddo

      ngrd2=0
      do i=2,iend2
         ii=(i-1)*nni
         do j=2,nni-1
            ngrd2=ngrd2+1
c      ires1(j,i)=ngrd2
c      iresni(ngrd2)=j
c      iresmu(ngrd2)=i
            indx2(ngrd2)=ii+j
         enddo
      enddo

01000 continue

c     region vi
c     bondary values from extrapolation
c     values for points (j,nmut-3),...,(j,nmut) are calculated from
c     asymptotic expansion (orbitals) or multipole expansion (potentials)

      i0=0
      do i=5,iend1
         i0=i0+1
         indx6a(i0)=(i-1)*nni8+5
      enddo
      ngrd6a=i0

      i0=0
      do i=5,iend1
         i0=i0+1
         indx6b(i0)=(i-1)*nni8+nni+4
      enddo
      ngrd6b=i0

      i0=0
      do j=6,nni+3
         i0=i0+1
         indx7(i0)=4*nni8+j
      enddo
      ngrd7=i0

      if (imcase.eq.1) then
         ingr1(1,ig)=ngrd1
         ingr1(2,ig)=ngrd6a
         ingr1(3,ig)=ngrd6b
         ingr1(4,ig)=ngrd7	  	  	  
      elseif (imcase.eq.2) then	  
         ingr2(1,ig)=ngrd1
         ingr2(2,ig)=ngrd6a
         ingr2(3,ig)=ngrd6b
         ingr2(4,ig)=ngrd7	  	  	  
      endif	  	

      if (idbg(41).ne.0) then
         write(*,*) 'indx1 ',(indx1 (i),i=1,ngrd1)
         write(*,*) 'indx2 ',(indx2 (i),i=1,ngrd2)
         write(*,*) 'indx6a',(indx6a(i),i=1,ngrd6a)
         write(*,*) 'indx6b',(indx6b(i),i=1,ngrd6b)
         write(*,*) 'indx7 ',(indx7 (i),i=1,ngrd7 )
      endif

      if (idbg(42).ne.0) then
         write(*,*) 'ngrd1,ngrd2,ngrd6a,ngrd6b,ngrd7  ',
     &        ngrd1,ngrd2,ngrd6a,ngrd6b,ngrd7
      endif
      return

00200 continue

c     the following sequence of mesh points is a modification of
c     Laaksonen's "middle" type of sweeps

      i0=0
      nmux=(nmu8-1)/2
      nnix=(nni8-1)/2

      do i=6,iend1
         if (i.le.nmux) then
            ii=nmux-i+6
         else
            ii=i
         endif
         i2=ii-4

         do j=6,nni+3
            if(j.le.nnix) then
               jj=nnix-j+6
            else
               jj=j
            endif
            j2=jj-4

            i0=i0+1
            indx2(i0)=(i2-1)*nni  +j2
            indx1(i0)=(ii-1)*nni8 +jj
         enddo
      enddo

      ngrd1=i0
      ngrd2=i0

c     bondary points like in the itype=1 case

      goto 1000

00300 continue

c     natural row-wise ordering
      
c     region i

      ngrd1=0
      do j=6,nni+3
         do i=6,iend1
            ii=(i-1)*nni8
            ngrd1=ngrd1+1
            indx1(ngrd1)=ii+j
         enddo
      enddo

      ngrd2=0
      do j=2,nni-1
         do i=2,iend1-4
            ii=(i-1)*nni
            ngrd2=ngrd2+1
            indx2(ngrd2)=ii+j
         enddo
      enddo

c     bondary points like in the itype=1 case

      goto 1000

00400 continue

c     reverse  ordering: mi=nmut-4 sweep in ni variable, mi=nmut-5 followed
c     by another sweep in ni etc.. it should allow the bondary conditions
c     to propagate quicker into the solution.

c     region i

      ngrd1=0
      do j=nni+3,6,-1
         do i=iend1,6,-1
            ii=(i-1)*nni8
            ngrd1=ngrd1+1
            indx1(ngrd1)=ii+j
         enddo
      enddo

      ngrd2=0
      do j=nni-1,2,-1
         do i=iend1-4,2,-1
            ii=(i-1)*nni
            ngrd2=ngrd2+1
            indx2(ngrd2)=ii+j
         enddo
      enddo

c   bondary points like in the itype=1 case

      goto 1000

00500 continue
00600 continue
00700 continue
00800 continue
      write(*,*) 'undefined type of ordering'
      stop 'mesh'

      end