********************************************************************************
*                                                                              *
*  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 ### fill ###
c
c     Immerses FUN array into WORK and adds boundary values.
      subroutine fill (nni,nmi,isym,fun,work)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      dimension fun(nni,nmi),work(nni+8,nmi+8)

c     fill the interior of work array

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

c     isym = 1 - even symmetry, isym =-1 - odd symmetry
c     values over i=nmi bondary are determined from the asymptotic expansion

      if (isym.eq.1.d0) then

c  mu=1...4	
         do i=2,5
            do j=1,nni
               work(j+4,6-i)= fun(j,i)
            enddo
         enddo

c        mu=nmu+5....nmu+8
c        the following code is necessary since the derivatives must be 
c        calculated up to i=nmi

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

c        ni=1...4
         do i=1,nmi
            do j=2,5
               work(6-j,i+4)= fun(j,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,i)
             enddo
          enddo
      else
         do i=2,5
            do j=1,nni
               work(j+4,6-i)=-fun(j,i)
            enddo
         enddo

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

         do i=1,nmi
            do j=2,5
               work(6-j,i+4)=-fun(j,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,i)
            enddo
         enddo
      endif

      return
      end