subroutine diff(h,f,df,dir,m,d,nl) 
implicit none 
integer                       :: nl
real*8                        :: h
real*8,dimension(:,:,:)    :: f, df
real*8,dimension(:)        :: d
integer,dimension(:)       :: m
character                     :: dir 

interface 
   subroutine diff_int(h,f,df,dir)
      implicit none
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      character                     :: dir
   end subroutine

   subroutine diff_ob(h,f,df,dir)
      implicit none
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      character                     :: dir
   end subroutine

   subroutine diff_ib(h,f,df,dir,m,d,nl)
      implicit none
      integer                       :: x, ny, nz, nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      character                     :: dir
   end subroutine
end interface

 
   call diff_int(h,f,df,dir) 
   call diff_ob(h,f,df,dir) 
 
   if (nl.gt.1) call diff_ib(h,f,df,dir,m,d,nl) 
 
return 
end 
 
subroutine diff_int(h,f,df,dir) 
implicit none 
real*8                        :: h
real*8,dimension(:,:,:)    :: f, df
character                     :: dir 
interface
   subroutine diffx_int(h,f,df)
      implicit none
      real*8                        :: h, h2
      real*8,dimension(:,:,:)    :: f, df
      character                     :: dir
   end subroutine

   subroutine diffy_int(h,f,df)
      implicit none
      real*8                        :: h, h2
      real*8,dimension(:,:,:)    :: f, df
      character                     :: dir
   end subroutine

   subroutine diffz_int(h,f,df)
      implicit none
      real*8                        :: h, h2
      real*8,dimension(:,:,:)    :: f, df
      character                     :: dir
   end subroutine
end interface

   if (dir.eq.'x') then 
      call diffx_int(h,f,df) 
   else if (dir.eq.'y') then 
      call diffy_int(h,f,df) 
   else if (dir.eq.'z') then 
      call diffz_int(h,f,df) 
   else 
      call err('diff_int','wrong derivative')    
   end if 
 
   return 
   end 
 
subroutine diff_ob(h,f,df,dir) 
implicit none
real*8                        :: h
real*8,dimension(:,:,:)    :: f, df
character                     :: dir

interface
   subroutine diffx_ob(h,f,df)
      implicit none
      real*8,dimension(:,:,:)     :: f,df
      real*8                      :: h2,hi,h
      integer,dimension(3)        :: shp
      integer                     :: nx,ny,nz
      integer                       :: i,j,k
	end subroutine

   subroutine diffy_ob(h,f,df)
      implicit none
      real*8,dimension(:,:,:)     :: f,df
      real*8                      :: h2,hi,h
      integer,dimension(3)          :: shp
      integer                       :: nx,ny,nz
      integer                       :: i,j,k
	end subroutine

   subroutine diffz_ob(h,f,df)
      implicit none
      real*8,dimension(:,:,:)     :: f,df
      real*8                      :: h2,hi,h
      integer,dimension(3)          :: shp
      integer                       :: nx,ny,nz
      integer                       :: i,j,k
	end subroutine
end interface

   if (dir.eq.'x') then 
      call diffx_ob(h,f,df) 
   else if (dir.eq.'y') then 
      call diffy_ob(h,f,df) 
   else if (dir.eq.'z') then 
      call diffz_ob(h,f,df) 
   else 
      call err('diff_ob','wrong derivative')    
   end if 
 
   return 
   end 
 
subroutine diff_ib(h,f,df,dir,m,d,nl) 
implicit none
integer                       :: x, ny, nz, nl
real*8                        :: h
real*8,dimension(:,:,:)    :: f, df
real*8,dimension(:)        :: d
integer,dimension(:)       :: m
character                     :: dir

interface
   subroutine diffy_ib(h,f,df,m,d,nl)
      implicit none
      integer                      :: nl
      real*8,dimension(:,:,:)      :: f,df
      real*8                       :: h,h2
      real*8,dimension(:)          :: d
      integer,dimension(:)         :: m
      integer                      :: i,j,k,l
      real*8                       :: d1,d2,d3,da
      integer,dimension(3)          :: shp
      integer                       :: nx,ny,nz
	end subroutine

   subroutine diffx_ib(h,f,df,m,d,nl)
      implicit none
      integer                      :: nl
      real*8,dimension(:,:,:)      :: f,df
      real*8                       :: h,h2
      real*8,dimension(:)          :: d
      integer,dimension(:)         :: m
      integer                      :: i,j,k,l
      real*8                       :: d1,d2,d3,da
      integer,dimension(3)          :: shp
      integer                       :: nx,ny,nz
	end subroutine

   subroutine diffz_ib(h,f,df,m,d,nl)
      implicit none
      integer                      :: nl
      real*8,dimension(:,:,:)      :: f,df
      real*8                       :: h,h2
      real*8,dimension(:)          :: d
      integer,dimension(:)         :: m
      integer                      :: i,j,k,l
      real*8                       :: d1,d2,d3,da
      integer,dimension(3)          :: shp
      integer                       :: nx,ny,nz
	end subroutine
end interface

   if (dir.eq.'x') then 
      call diffx_ib(h,f,df,m,d,nl) 
   else if (dir.eq.'y') then 
      call diffy_ib(h,f,df,m,d,nl) 
   else if (dir.eq.'z') then 
      call diffz_ib(h,f,df,m,d,nl) 
   else 
      call err('diff_ib','wrong derivative') 
   end if 
 
   return 
   end 
 
subroutine diffx_int(h,f,df) 
implicit none
real*8                        :: h, h2
real*8,dimension(:,:,:)    :: f, df
character                     :: dir

   h2 = 0.5d0 / h 
   df = h2 * ( cshift(f,1,1) - cshift(f,-1,1) )
 
return 
end 
 
subroutine diffy_int(h,f,df) 
implicit none
real*8                        :: h, h2
real*8,dimension(:,:,:)    :: f, df
character                     :: dir

      h2 = 0.5d0 / h 
      df = h2 * ( cshift(f,1,2) - cshift(f,-1,2) )

 
return 
end 
 
subroutine diffz_int(h,f,df) 
implicit none
real*8                        :: h, h2
real*8,dimension(:,:,:)    :: f, df
character                     :: dir

   h2 = 0.5d0 / h 
   df = h2 * ( cshift(f,1,3) - cshift(f,-1,3) )
 
return 
end 
 
 
subroutine diffx_ob(h,f,df) 
implicit none 
real*8,dimension(:,:,:)     :: f,df 
real*8                      :: h2,hi,h


integer,dimension(3)        :: shp
integer                     :: nx,ny,nz
integer                       :: i,j,k

   shp = shape(f)

   nx = shp(1)
   ny = shp(2)
   nz = shp(3)
 
   hi = 1.0d0/h 
   h2 = 0.5d0 / h 

   do k = 1 ,nz
   do j = 1 ,ny
      df(1,j,k) = h2 * ( -3.0d0*f(1,j,k)+ 4.0d0*f(1+1,j,k) - f(1+2,j,k)) 
      df(nx,j,k) = h2 * (3.0d0*f(nx,j,k)- 4.0d0*f(nx-1,j,k) + f(nx-2,j,k)) 
   end do
   end do
 
   do k = 1 , nz
   do i = 2 ,nx - 1
      df(i,1,k) = h2 *(f(i+1,1,k)-f(i-1,1,k) ) 
      df(i,ny,k) = h2 *(f(i+1,ny,k)-f(i-1,ny,k) ) 
   end do
   end do
 
   do j= 1 ,ny
   do i = 2 ,nx-1
      df(i,j,1) = h2 *(f(i+1,j,1)-f(i-1,j,1) ) 
      df(i,j,nz) = h2 *(f(i+1,j,nz)-f(i-1,j,nz) ) 
   end do
   end do
 
return 
end
 
subroutine diffy_ob(h,f,df)
implicit none
real*8,dimension(:,:,:)     :: f,df
real*8                      :: h2,hi,h

integer,dimension(3)          :: shp
integer                       :: nx,ny,nz
integer                       :: i,j,k

   shp = shape(f)

   nx = shp(1)
   ny = shp(2)
   nz = shp(3)

   hi = 1.0d0/h
   h2 = 0.5d0 / h

   do k = 1 , nz
   do i = 1 , nx
      df(i,1,k) = h2 * ( -3.0d0*f(i,1,k) + 4.0d0*f(i,1+1,k) - f(i,1+2,k)) 
      df(i,ny,k) = h2 * (3.0d0*f(i,ny,k)- 4.0d0*f(i,ny-1,k) + f(i,ny-2,k)) 
   end do
   end do
 
   do k = 1 , nz
   do j = 2 , ny-1
      df(1,j,k) = h2 *(f(1,j+1,k)-f(1,j-1,k) ) 
      df(nx,j,k) = h2 *(f(nx,j+1,k)-f(nx,j-1,k) ) 
   end do
   end do
 
   do i = 1 ,nx
   do j = 2 , ny -1
      df(i,j,1) = h2 *(f(i,j+1,1)-f(i,j-1,1) ) 
      df(i,j,nz) = h2 *(f(i,j+1,nz)-f(i,j-1,nz) ) 
   end do
   end do
 
return 
end 
 
subroutine diffz_ob(h,f,df)
implicit none
real*8,dimension(:,:,:)     :: f,df
real*8                      :: h2,hi,h

integer,dimension(3)          :: shp
integer                       :: nx,ny,nz
integer                       :: i,j,k

   shp = shape(f)

   nx = shp(1)
   ny = shp(2)
   nz = shp(3)

   hi = 1.0d0/h
   h2 = 0.5d0 / h

   do i = 1,nx
   do j= 1 ,ny

      df(i,j,1) = h2 * ( -3.0d0*f(i,j,1)        &
                + 4.0d0*f(i,j,1+1) - f(i,j,1+2) ) 
      df(i,j,nz) = h2 * (  3.0d0*f(i,j,nz)      &
                - 4.0d0*f(i,j,nz-1) + f(i,j,nz-2) ) 
   end do
   end do
 
   do i =1 ,nx
   do k = 2 ,nz-1
      df(i,1,k) = h2 *(f(i,1,k+1)-f(i,1,k-1) ) 
      df(i,ny,k) = h2 *(f(i,ny,k+1)-f(i,ny,k-1) ) 
   end do
   end do
 
   do j = 1, ny
   do k = 2, nz-1
      df(1,j,k) = h2 *(f(1,j,k+1)-f(1,j,k-1) ) 
      df(nx,j,k) = h2 *(f(nx,j,k+1)-f(nx,j,k-1) ) 
   end do
   end do
 

   return 
   end 
!------------------------------------------------------------------------------ 
subroutine diffy_ib(h,f,df,m,d,nl) 
implicit none 
integer                      :: nl 
real*8,dimension(:,:,:)      :: f,df
real*8                       :: h,h2 
real*8,dimension(:)          :: d
integer,dimension(:)         :: m
integer                      :: i,j,k,l
real*8                       :: d1,d2,d3,da 

integer,dimension(3)          :: shp
integer                       :: nx,ny,nz

   shp = shape(f)

   nx = shp(1)
   ny = shp(2)
   nz = shp(3)

 
   h2 = 0.5d0 / h 
 
   do l = 1 , nl 
 
      i  = m(l) 
      j  = m(nl+l) 
      k  = m(2*nl+l) 
      d1 = d(l) 
      d2 = d(nl+l) 
      d3 = d(2*nl+l) 
      da = abs(d2) 

		if (d2.eq.1.0d0) then
         df(i,j,k) = h2*(-3.0d0*f(i,j,k)+4.0d0*f(i,j+1,k) - f(i,j+2,k))
		else if (d2.eq.-1.0d0) then
         df(i,j,k) = -h2*(-3.0d0*f(i,j,k)+4.0d0*f(i,j-1,k) - f(i,j-2,k))
		else
         df(i,j,k) = f(i,j+1,k)-f(i,j-1,k)
		end if
   end do 
 
return 
end 
!------------------------------------------------------------------------------ 
subroutine diffx_ib(h,f,df,m,d,nl)
implicit none
integer                      :: nl
real*8,dimension(:,:,:)      :: f,df
real*8                       :: h,h2
real*8,dimension(:)          :: d
integer,dimension(:)         :: m
integer                      :: i,j,k,l
real*8                       :: d1,d2,d3,da

integer,dimension(3)          :: shp
integer                       :: nx,ny,nz

   shp = shape(f)

   nx = shp(1)
   ny = shp(2)
   nz = shp(3)


   h2 = 0.5d0 / h

   do l = 1 , nl 
 
      i  = m(l) 
      j  = m(nl+l) 
      k  = m(2*nl+l) 
      d1 = d(l) 
      d2 = d(nl+l) 
      d3 = d(2*nl+l) 
      da = abs(d1) 
      df(i,j,k) = h2*(d1*(-3.0d0*f(i,j,k)+4.0d0*f(i+int(d1),j,k)     &
                -      f(i+2*int(d1),j,k)) +     &
                  (1.0d0-da)*(f(i+1,j,k)-f(i-1,j,k) ) ) 
 
   end do 
 
return 
end 
!------------------------------------------------------------------------------ 
subroutine diffz_ib(h,f,df,m,d,nl)
implicit none
integer                      :: nl
real*8,dimension(:,:,:)      :: f,df
real*8                       :: h,h2
real*8,dimension(:)          :: d
integer,dimension(:)         :: m
integer                      :: i,j,k,l
real*8                       :: d1,d2,d3,da

integer,dimension(3)          :: shp
integer                       :: nx,ny,nz

   shp = shape(f)

   nx = shp(1)
   ny = shp(2)
   nz = shp(3)


   h2 = 0.5d0 / h
 
   do l = 1 , nl 
 
      i  = m(l) 
      j  = m(nl+l) 
      k  = m(2*nl+l) 
      d1 = d(l) 
      d2 = d(nl+l) 
      d3 = d(2*nl+l) 
      da = abs(d3) 
      df(i,j,k) = h2*(d3*(-3.0d0*f(i,j,k)+4.0d0*f(i,j,k+int(d3))     &
                -      f(i,j,k+2*int(d3))) +     &
                (1.0d0-da)*(f(i,j,k+1)-f(i,j,k-1) ) ) 

   end do 
 
return 
end 
