********************************************************************************
*                                                                              *
*  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 ### fill4 ###
c
c     Immerses FUN array into WORK and adds boundary values: 
c     the last subgrid
c
      subroutine fill4 (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(450,4)	

c     muoffs=iemu(ig-1)-1
c     fill the interior of work array

      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 cint4)

      do i=2,4
         call mxv(fun(1,iadint4(i)),nni,cint4(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     the following code is necessary since the derivatives must be calculated
c     up to i=nmi
      do i=1,4
         do j=1,nni
	    work(j+4,nmi+4+i)= fun(j,muoffs+nmi)
         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