********************************************************************************
*                                                                              *
*  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 ### fill3 ###
c
c     Immerses FUN array into WORK and adds boundary values: 
c     the interior subgrid
c
      subroutine fill3 (nni,nmi,fun,work)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)

      common/interp/ cint2(9,4),cint3l(9,4),cint3r(9,4),cint4(9,4),
     &		   iadint2(4),iadint3l(4),iadint3r(4),iadint4(4)
      common/sorwk/ omega,omega1,dmu2t(4),dmu1t(4),icase,ifill,
     &              isym,muoffs

      dimension fun(nni,*),work(nni+8,nmi+8),fint(_MAXMU_,4)	

c     muoffs=iemu(ig-1)	
c     fill the interior of work array
c
      do i=1,nmi
         do j=1,nni
            work(j+4,i+4)=fun(j,muoffs+i)
         enddo
      enddo

c     values over i=1 bondary are determined from the interpolation formula
c     mu=2,...,4 from interpolation (coefficients from cint3l)

      do i=2,4
         call mxv(fun(1,iadint3l(i)),nni,cint3l(1,i),9,fint(1,i))
      enddo

      do i=2,4
         do j=1,nni
	    work(j+4,i)= fint(j,i)
         enddo
      enddo

c     values over i=nmi bondary are determined from the interpolation formula
c     mu=nmu+5....nmu+8 from interpolation (coefficients from cint3r)

      do i=1,4
         call mxv(fun(1,iadint3r(i)),nni,cint3r(1,i),9,fint(1,i))
      enddo

      do i=1,4
         do j=1,nni
	    work(j+4,nmi+4+i)= fint(j,i)
         enddo
      enddo

c     isym = 1 - even symmetry, isym =-1 - odd symmetry

      if (isym.eq.1) then

c ni=1...4
         do i=1,nmi
            do j=2,5
               work(6-j,i+4)= fun(j,muoffs+i)
            enddo
         enddo

c ni=ni+4...ni+8
	  do i=1,nmi
             jj=0
             do j=nni-4,nni-1
                jj=jj+1
                work(nni+9-jj,i+4)= fun(j,muoffs+i)
             enddo
          enddo
      else
	  do i=1,nmi
             do j=2,5
                work(6-j,i+4)=-fun(j,muoffs+i)
             enddo
          enddo

	  do i=1,nmi
             jj=0
             do j=nni-4,nni-1
                jj=jj+1
                work(nni+9-jj,i+4)=-fun(j,muoffs+i)
             enddo
          enddo
      endif

      return
      end