********************************************************************************
*                                                                              *
*  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 ### difmu ###
c
c    This routine calculates 
c
c     (\frac{\partial^2}{\partial\mu^2} + 
c           b(\ni,\mu) \frac{\partial}{\partial \mu}) f(\ni,\mu)
c
c    Function f has been imersed in the array f(nni+8,nmu+8) in order
c    to calculated derivatives in all the grid points.
c    Originally the routine was used for a single grid of constatnt step
c    size hmu. Accordingly dmu array contained the first- and second-order
c    derivative coefficients (taken from the 8th-order Sterling
c    interpolation formula) multiplied by the b array. 
c
c    To make the routine work in the multigrid case (ngrids.ne.1) 
c    the values of dmu(k,imu) for 
c
c        imu=iemu(1)-3 ... iemu(1)+3
c        imu=iemu(2)-3 ... iemu(2)+3
c        .
c        imu=iemu(ngrids-1)-3 ... iemu(ngrids-1)+3
c
c    must be prepared with derivative coefficients which are based on
c    other interpolation formulae taking into account different grid
c    density to the left and right of the grid boundaries. See prepfix
c    for detailes.
c
c
c    This routine calculates 
c
c     {\partial^2 / \partial\mu^2 + 
c           b(\ni,\mu) \partial / \partial \mu} f(\ni,\mu)
c
c    Function f has been imersed in the array f(nni+8,nmu+8) in order
c    to calculated derivatives in all the grid points.
c    Originally the routine was used for a single grid of constatnt step
c    size hmu. Accordingly dmu array contained the first- and second-order
c    derivative coefficients (taken from the 8th-order Sterling
c    interpolation formula) multiplied by the B array. 
c
c    To make the routine work in the multigrid case (ngrids.ne.1!) 
c    the values of dmu(k,imu) for 
c
c        imu=iemu(1)-3 ... iemu(1)+3
c        imu=iemu(2)-3 ... iemu(2)+3
c        .
c        imu=iemu(ngrids-1)-3 ... iemu(ngrids-1)+3
c
c    must be prepared with derivative coefficients which are based on
c    other interpolation formula taking into account different grid
c    density to the left and right of the grid boundaries. See prepfix
c    for detailes.
c
      subroutine difmu (n,f,fd)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension f(nni+8,*),fd(nni+8,*)	

c     The following loop runs now till j=mxnmu so that the derivatives are
c     also defined in the tail region, i.e. for (i,mxnmu-3),...,(i,mxnmu)
c     To this end the fill routine provides extra values for
c     (i,mxnmu+1),...,(i,mxnmu+4) points.

      do j=1,n
         call mxv(f(1,j),nni+8,dmu(1,j),9,fd(1,j+4))
      enddo
      return
      end