subroutine expfil
     &   ( fn,ifilt_x,ifilt_y,ifilt_z,iord,
     &     fiL, fiR, fjL, fjR, fkL, fkR,
     &     piL, piR, pjL, pjR, pkL, pkR,
     &     wiL, wiR, wjL, wjR, wkL, wkR,
     &     ierror,error,ierr_coeff,err_coeff)
!
      include 'filter_include.h'
!
!
      call expfil_coeff(iord,disi,disb,ierr_coeff,err_coeff)
!
!  EXPLICIT FILTERS - for even orders, 2, 4, 6, 8, 10, 12
!
      ka = iord / 2
      kb = iord
      kap1 = ka + 1
!
      width = iord + 1
c     I think we should check that the patch width > stencil width, not the
c     window width > stencil width. JR 04/10/03
c      nx = wiR - wiL + 1
c      ny = wjR - wjL + 1
c      nz = wkR - wkL + 1

      nx = piR - piL + 1
      ny = pjR - pjL + 1
      nz = pkR - pkL + 1
!
      if(( nx .lt. width ) .and. ( nx .gt. 1 )) then
        ierror = -10
        error = 'Expfil: Stencil width wider than x-domain'
        return
      endif
      if(( ny .lt. width ) .and. ( ny .gt. 1 )) then
        ierror = -20
        error = 'Expfil: Stencil width wider than y-domain'
        return
      endif
      if(( nz .lt. width ) .and. ( nz .gt. 1 )) then
        ierror = -30
        error = 'Expfil: Stencil width wider than z-domain'
        return
      endif
!
      do m = wkL, wkR
        do j = wjL, wjR
          do i = wiL, wiR
            dfn(i,j,m) = 0.0
          end do
        end do
      end do
!
!  FILTERING IN X-DIRECTION
!  ========================
!
      if (( ifilt_x .eq. 1) .and. ( nx .ge. width )) then
!
!  Interior filter in the x direction
! 
          iL = wiL + ka
          iR = wiR - ka
!
        do k = - ka, ka
         kmag = 1 + abs ( k )
          do m = wkL, wkR
           do j = wjL, wjR
            do i = iL, iR
             dfn(i,j,m) = dfn(i,j,m) + disi(kmag) * fn(i+k,j,m)
            end do
           end do
          end do
        end do
!
!  Boundary points
! 
        do k1 = wiL,  wiL - 1 + ka
          k1a =        wiR - (k1 - wiL)
          k1ind =      k1 - wiL + 1
         do k2 = wiL, wiL - 1 + kb
          k2a =       wiR - (k2-wiL)
          k2ind =     k2 - wiL + 1
           do m = wkL, wkR
            do j = wjL, wjR
               dfn(k1 ,j,m) = dfn(k1 ,j,m) + disb(k1ind,k2ind) * 
     &                                       fn(k2 ,j,m)
               dfn(k1a,j,m) = dfn(k1a,j,m) + disb(k1ind,k2ind) * 
     &                                       fn(k2a,j,m)
            end do
           end do
          end do
        end do
!
!  Apply Filter
!
        do m = wkL, wkR
          do j = wjL, wjR
            do i = wiL, wiR
              fn(i,j,m) = fn(i,j,m) + dfn(i,j,m)
              dfn(i,j,m) = 0.0
            end do
          end do
        end do

      endif
!
!  FILTERING IN Y-DIRECTION
!  ========================
!
      if(( ifilt_y .eq. 1) .and. ( ny .ge. width )) then
!
!  Interior filter in the y direction
!
        jL = wjL + ka
        jR = wjR - ka
!
        do k = - ka, ka
          kmag = 1 + abs ( k )
          do j = jL, jR
            do m = wkL, wkR
              do i = wiL, wiR
                dfn(i,j,m) = dfn(i,j,m) + disi(kmag)*fn(i,j+k,m)
              end do
            end do
          end do
        end do
!
!  Boundary points
!
        do k1 = wjL,   wjL - 1 + ka
          k1a =        wjR - (k1 - wjL)
          k1ind =      k1 - wjL + 1
          do k2 = wjL, wjL - 1 + kb
            k2a =      wjR - (k2 - wjL)
            k2ind =    k2 - wjL + 1
            do m = wkL, wkR
              do i = wiL, wiR
                 dfn(i,k1 ,m) = dfn(i,k1 ,m) + disb(k1ind,k2ind) * 
     &                                         fn(i,k2 ,m)
                 dfn(i,k1a,m) = dfn(i,k1a,m) + disb(k1ind,k2ind) * 
     &                                         fn(i,k2a,m)
              end do
            end do
          end do
        end do
!
!  Apply Filter
! 
        do m = wkL, wkR
          do j = wjL, wjR
            do i = wiL, wiR
              fn(i,j,m) = fn(i,j,m) + dfn(i,j,m)
              dfn(i,j,m) = 0.0
            end do
          end do
        end do

      endif
!
!  FILTERING IN Z-DIRECTION
!  ========================
!
      if(( ifilt_z .eq. 1) .and. ( nz .ge. width )) then
!
!  Interior filter in the z direction
!
        kL = wkL + ka
        kR = wkR - ka
!
        do k = - ka, ka
          kmag = 1 + abs ( k )
          do m = kL, kR
            do j = wjL, wjR
              do i = wiL, wiR
                dfn(i,j,m) = dfn(i,j,m) + disi(kmag)*fn(i,j,m+k)
              end do
            end do
          end do
        end do
!
!  Boundary points
!
        do k1 = wkL,   wkL - 1 + ka
          k1a =        wkR - (k1 - wkL)
          k1ind =      k1 - wkL + 1
          do k2 = wkL, wkL - 1 + kb
            k2a =      wkR - (k2 - wkL)
            k2ind =    k2 - wkL + 1
            do j = wjL, wjR
              do i = wiL, wiR
                dfn(i,j,k1 ) = dfn(i,j,k1 ) + disb(k1ind,k2ind) * 
     &                                        fn(i,j,k2 )
                dfn(i,j,k1a) = dfn(i,j,k1a) + disb(k1ind,k2ind) *
     &                                        fn(i,j,k2a)
              end do
            end do
          end do
        end do
!
!  Apply Filter
!
        do m = wkL, wkR
          do j = wjL, wjR
            do i = wiL, wiR
              fn(i,j,m) = fn(i,j,m) + dfn(i,j,m)
              dfn(i,j,m) = 0.0
            end do
          end do
        end do

      endif
!
      return
      end