! $Id$
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM integer, parameter :: nghost = 3
!
!***************************************************************
module Deriv
!
! Experimental module with explicit difference formulae for nonequidistant grids
! not yet fully tested and documented
!
! 26-Mar-12/MR: derived from deriv.f90
!
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages, only: fatal_error, warning, not_implemented
!
  implicit none
!
  private
!
  include 'deriv.h'
!
!  interface deri_3d                  !not working, why?
!    module procedure deri_3d_ind
!    module procedure deri_3d_inds
!  endinterface
!
  contains
!
!***********************************************************************
    subroutine initialize_deriv()
!
!  Initialize stencil coefficients
!
      integer :: i

      do i=1,6

        if (lspherical_coords) then                         !-> grid.f90
          r1i  (:,i) = r1_mn**i
          sth1i(:,i) = sin1th**i
        endif
        if (lcylindrical_coords) r1i(:,i) = rcyl_mn1**i

      enddo

      der_coef(:,1) = (/  0.   ,   3./4., -3./20.,  1./60./)
      der_coef(:,3) = (/  0.   , -13./8.,  1.    , -1./8. /)
      der_coef(:,4) = (/ 28./3., -13./2.,  2.    , -1./6. /)
      der_coef(:,5) = (/  0.   ,   2.5  , -2.    ,  0.5   /)
      der_coef(:,6) = (/-20.   ,  15.   , -6.    ,  1.    /)

      select case (der2_type)
!
      case ('standard')
        der_coef(:,2) = (/ -49./18., 3./2.,-3./20.,1./90. /)
!
      case ('tuned1')
        der_coef(:,2) = (/ -0.75, 0.34375, 0.125, -0.09375 /)  !tbc
!
      case default
        call fatal_error('initialize_deriv',"der2_type doesn't exist")
!
      endselect

      if ( lequidist(1) ) then

        if (.not.allocated(coeffsx)) allocate( coeffsx(-3:3,8,2) )

        do i=1,6
          coeffsx( 0:3,i,1) = der_coef(:,i)*dx_1(1)**i
          coeffsx( 0:3,i,2) = der_coef(:,i)
        enddo

      else

        if (.not.allocated(coeffsx)) allocate( coeffsx(-3:3,8,nx) )
        if (.not.allocated(coeffsx_1)) allocate( coeffsx_1(-2:2,2,2) )

        do i=l1,l2
          call calc_coeffs( x(i-nghost+1:i+nghost)-x(i-nghost:i+nghost-1), coeffsx(-3,1,i-l1+1) )   !use dx_1 ?
        enddo
!        call test( x, dx_1, nx, coeffsx )

        call calc_coeffs_1( x(l1-1:l1+2)-x(l1-2:l1+1), coeffsx_1(-2,1,1) )
        call calc_coeffs_1( x(l2-1:l2+2)-x(l2-2:l2+1), coeffsx_1(-2,1,2) )

      endif

      if ( lequidist(2) ) then

        if (.not.allocated(coeffsy)) allocate( coeffsy(-3:3,8,2) )
        do i=1,6
          coeffsy( 0:3,i,1) = der_coef(:,i)*dy_1(1)**i
          coeffsy( 0:3,i,2) = der_coef(:,i)
        enddo

      else

        if (.not.allocated(coeffsy)) allocate( coeffsy(-3:3,8,ny) )
        if (.not.allocated(coeffsy_1)) allocate( coeffsy_1(-2:2,2,2) )

        do i=m1,m2
          call calc_coeffs( y(i-nghost+1:i+nghost)-y(i-nghost:i+nghost-1), coeffsy(-3,1,i-m1+1) )
        enddo

        call calc_coeffs_1( y(m1-1:m1+2)-y(m1-2:m1+1), coeffsy_1(-2,1,1) )
        call calc_coeffs_1( y(m2-1:m2+2)-y(m2-2:m2+1), coeffsy_1(-2,1,2) )

      endif

      if ( lequidist(3) ) then

        if (.not.allocated(coeffsz)) allocate( coeffsz(-3:3,8,2) )

        do i=1,6
          coeffsz( 0:3,i,1) = der_coef(:,i)*dz_1(1)**i
          coeffsz( 0:3,i,2) = der_coef(:,i)
        enddo

      else

        if (.not.allocated(coeffsz)) allocate( coeffsz(-3:3,8,nz) )
        if (.not.allocated(coeffsz_1)) allocate( coeffsz_1(-2:2,2,2) )

        do i=n1,n2
          call calc_coeffs( z(i-nghost+1:i+nghost)-z(i-nghost:i+nghost-1), coeffsz(-3,1,i-n1+1) )
        enddo
!        call test( z, dz_1, nz, coeffsz )

        call calc_coeffs_1( z(n1-1:n1+2)-z(n1-2:n1+1), coeffsz_1(-2,1,1) )
        call calc_coeffs_1( z(n2-1:n2+2)-z(n2-2:n2+1), coeffsz_1(-2,1,2) )

!        call test_1( z, dz_1, nz, coeffsz_1 )
      endif
!
    endsubroutine initialize_deriv
!***********************************************************************
    subroutine finalize_deriv()

      deallocate( coeffsx, coeffsy, coeffsz )
      if (allocated(coeffsx_1)) deallocate(coeffsx_1)
      if (allocated(coeffsy_1)) deallocate(coeffsy_1)
      if (allocated(coeffsz_1)) deallocate(coeffsz_1)

    endsubroutine finalize_deriv
!***********************************************************************
    subroutine der_main(f,k,df,j,ignoredx)
!
!  calculates derivative df_k/dx_j
!  accurate to 6th order, explicit, periodic (??)
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: df
      logical, intent(in), optional :: ignoredx
      integer :: j,k
!
      intent(in)  :: f,k,j
      intent(out) :: df
!
      if (present(ignoredx)) call fatal_error('der_main', 'optional argument ignoredx is not implemented. ')
!
      call deri_3d( f(1,1,1,k), df, 1, j )
!
    endsubroutine der_main
!***********************************************************************
    subroutine der_other(f,df,j)
!
!  Along one pencil in NON f variable
!  calculate derivative of a scalar, get scalar
!  accurate to 6th order, explicit, periodic
!
!  26-nov-02/tony: coded - duplicate der_main but without k subscript
!                          then overload the der interface.
!  25-jun-04/tobi+wolf: adapted for non-equidistant grids
!  21-feb-07/axel: added 1/r and 1/pomega factors for non-coord basis
!
      real, dimension (mx,my,mz) :: f
      real, dimension (nx) :: df
      integer :: j
!
      intent(in)  :: f,j
      intent(out) :: df
!
!debug      if (loptimise_ders) der_call_count(1,icount_der_other,j,1) = &
!debug                          der_call_count(1,icount_der_other,j,1) + 1
!
      call deri_3d(f,df,1,j)
!
    endsubroutine der_other
!***********************************************************************
    subroutine deri_3d(f,df,i,j,lignored,lnometric)
!
!  Along one pencil in NON f variable
!  calculate derivative of a scalar, get scalar
!  accurate to 6th order, explicit, periodic
!
!  26-nov-02/tony: coded - duplicate der_main but without k subscript
!                          then overload the der interface.
!  25-jun-04/tobi+wolf: adapted for non-equidistant grids
!  21-feb-07/axel: added 1/r and 1/pomega factors for non-cart basis
!
      real, dimension (mx,my,mz)          :: f
      real, dimension (nx)                :: df
      integer                             :: i,j
      logical,                   optional :: lignored, lnometric
!
      intent(in)  :: f,i,j,lignored,lnometric
      intent(out) :: df
!
      integer :: mm, nn, is
      logical :: lmetric
!
!debug      if (loptimise_ders) der_call_count(1,icount_der_other,j,1) = &
!debug                          der_call_count(1,icount_der_other,j,1) + 1
!
      is = 1
      if ( present(lignored) ) then
        if (lignored) is=2
      endif

      if (present(lnometric)) then
        lmetric = .not.lnometric
      else
        lmetric = .true.
      endif

      if (j==1) then
        if (nxgrid/=1) then
          if ( lequidist(j) ) then
            if ( mod(i,2) == 0 ) then
              df =  coeffsx(0,i,is)* f(l1  :l2  ,m,n)                   &
                  + coeffsx(1,i,is)*(f(l1+1:l2+1,m,n)+f(l1-1:l2-1,m,n)) &
                  + coeffsx(2,i,is)*(f(l1+2:l2+2,m,n)+f(l1-2:l2-2,m,n)) &
                  + coeffsx(3,i,is)*(f(l1+3:l2+3,m,n)+f(l1-3:l2-3,m,n))
            else
              df =  coeffsx(1,i,is)*(f(l1+1:l2+1,m,n)-f(l1-1:l2-1,m,n)) &
                  + coeffsx(2,i,is)*(f(l1+2:l2+2,m,n)-f(l1-2:l2-2,m,n)) &
                  + coeffsx(3,i,is)*(f(l1+3:l2+3,m,n)-f(l1-3:l2-3,m,n))
            endif
          else
            df =  coeffsx(-3,i,:)*f(l1-3:l2-3,m,n) &
                + coeffsx(-2,i,:)*f(l1-2:l2-2,m,n) &
                + coeffsx(-1,i,:)*f(l1-1:l2-1,m,n) &
                + coeffsx( 0,i,:)*f(l1  :l2  ,m,n) &
                + coeffsx( 1,i,:)*f(l1+1:l2+1,m,n) &
                + coeffsx( 2,i,:)*f(l1+2:l2+2,m,n) &
                + coeffsx( 3,i,:)*f(l1+3:l2+3,m,n)
             !!lignore???
          endif
        else
          df=0.
          if (ip<=5) print*, 'deri_3d: Degenerate case in x-direction'
        endif
      elseif (j==2) then
        if (nygrid/=1) then
          if ( lequidist(j) ) then
            if ( mod(i,2) == 0 ) then
              df =  coeffsy(0,i,is)* f(l1:l2,m  ,n)                 &
                  + coeffsy(1,i,is)*(f(l1:l2,m+1,n)+f(l1:l2,m-1,n)) &
                  + coeffsy(2,i,is)*(f(l1:l2,m+2,n)+f(l1:l2,m-2,n)) &
                  + coeffsy(3,i,is)*(f(l1:l2,m+3,n)+f(l1:l2,m-3,n))
            else
              df =  coeffsy(1,i,is)*(f(l1:l2,m+1,n)-f(l1:l2,m-1,n)) &
                  + coeffsy(2,i,is)*(f(l1:l2,m+2,n)-f(l1:l2,m-2,n)) &
                  + coeffsy(3,i,is)*(f(l1:l2,m+3,n)-f(l1:l2,m-3,n))
            endif
          else
            mm = m-m1+1
            df =  coeffsy(-3,i,mm)*f(l1:l2,m-3,n) &
                + coeffsy(-2,i,mm)*f(l1:l2,m-2,n) &
                + coeffsy(-1,i,mm)*f(l1:l2,m-1,n) &
                + coeffsy( 0,i,mm)*f(l1:l2,m  ,n) &
                + coeffsy( 1,i,mm)*f(l1:l2,m+1,n) &
                + coeffsy( 2,i,mm)*f(l1:l2,m+2,n) &
                + coeffsy( 3,i,mm)*f(l1:l2,m+3,n)
             !!lignore???
          endif
          if (lmetric.and.(lspherical_coords .or. lcylindrical_coords)) df = df*r1i(:,i)
        else
          df=0.
          if (ip<=5) print*, 'deri_3d: Degenerate case in y-direction'
        endif
      elseif (j==3) then
        if (nzgrid/=1) then
          if ( lequidist(j) ) then
            if ( mod(i,2) == 0 ) then
              df =  coeffsz(0,i,is)* f(l1:l2,m,n  )                 &
                  + coeffsz(1,i,is)*(f(l1:l2,m,n+1)+f(l1:l2,m,n-1)) &
                  + coeffsz(2,i,is)*(f(l1:l2,m,n+2)+f(l1:l2,m,n-2)) &
                  + coeffsz(3,i,is)*(f(l1:l2,m,n+3)+f(l1:l2,m,n-3))
            else
              df =  coeffsz(1,i,is)*(f(l1:l2,m,n+1)-f(l1:l2,m,n-1)) &
                  + coeffsz(2,i,is)*(f(l1:l2,m,n+2)-f(l1:l2,m,n-2)) &
                  + coeffsz(3,i,is)*(f(l1:l2,m,n+3)-f(l1:l2,m,n-3))
            endif
          else
            nn = n-n1+1
            df =  coeffsz(-3,i,nn)*f(l1:l2,m,n-3) &
                + coeffsz(-2,i,nn)*f(l1:l2,m,n-2) &
                + coeffsz(-1,i,nn)*f(l1:l2,m,n-1) &
                + coeffsz( 0,i,nn)*f(l1:l2,m,n  ) &
                + coeffsz( 1,i,nn)*f(l1:l2,m,n+1) &
                + coeffsz( 2,i,nn)*f(l1:l2,m,n+2) &
                + coeffsz( 3,i,nn)*f(l1:l2,m,n+3)
            !!lignore??
          endif
          if (lmetric.and.lspherical_coords) df=df*r1i(:,i)*sth1i(m,i)
        else
          df=0.
          if (ip<=5) print*, 'deri_3d: Degenerate case in z-direction'
        endif
      endif
!
    endsubroutine deri_3d
!***********************************************************************
    subroutine deri_3d_inds(f,df,inds,j,lnometric)
!
! only for first derivatives, i.e., inds(i) = 1,7,8 !!!
!
      real, dimension (mx,my,mz)          :: f
      real, dimension (nx)                :: df
      integer, dimension(nx)              :: inds
      integer                             :: j
      logical,                   optional :: lnometric
!
      intent(in)  :: f,j,inds,lnometric
      intent(out) :: df
!
      integer :: mm, nn, ii, li, ini
      logical :: lmetric

      if (present(lnometric)) then
        lmetric = .not.lnometric
      else
        lmetric = .true.
      endif

      if (j==1) then
        if (nxgrid/=1) then
          do ii=1,nx
            li = l1+ii-1
            ini = inds(ii)
            df =  coeffsx(-3,ini,ii)*f(li-3,m,n) &
                + coeffsx(-2,ini,ii)*f(li-2,m,n) &
                + coeffsx(-1,ini,ii)*f(li-1,m,n) &
                + coeffsx( 0,ini,ii)*f(li,  m,n) &
                + coeffsx( 1,ini,ii)*f(li+1,m,n) &
                + coeffsx( 2,ini,ii)*f(li+2,m,n) &
                + coeffsx( 3,ini,ii)*f(li+3,m,n)
          enddo
        else
          df=0.
          if (ip<=5) print*, 'deri_3d: Degenerate case in x-direction'
        endif
      elseif (j==2) then
        if (nygrid/=1) then
          mm = m-m1+1
          df =  coeffsy(-3,inds,mm)*f(l1:l2,m-3,n) &
              + coeffsy(-2,inds,mm)*f(l1:l2,m-2,n) &
              + coeffsy(-1,inds,mm)*f(l1:l2,m-1,n) &
              + coeffsy( 0,inds,mm)*f(l1:l2,m  ,n) &
              + coeffsy( 1,inds,mm)*f(l1:l2,m+1,n) &
              + coeffsy( 2,inds,mm)*f(l1:l2,m+2,n) &
              + coeffsy( 3,inds,mm)*f(l1:l2,m+3,n)
          if (lmetric.and.(lspherical_coords .or. lcylindrical_coords)) df = df*r1i(:,1)
        else
          df=0.
          if (ip<=5) print*, 'deri_3d: Degenerate case in y-direction'
        endif
      elseif (j==3) then
        if (nzgrid/=1) then
          nn = n-n1+1
          df =  coeffsz(-3,inds,nn)*f(l1:l2,m,n-3) &
              + coeffsz(-2,inds,nn)*f(l1:l2,m,n-2) &
              + coeffsz(-1,inds,nn)*f(l1:l2,m,n-1) &
              + coeffsz( 0,inds,nn)*f(l1:l2,m,n  ) &
              + coeffsz( 1,inds,nn)*f(l1:l2,m,n+1) &
              + coeffsz( 2,inds,nn)*f(l1:l2,m,n+2) &
              + coeffsz( 3,inds,nn)*f(l1:l2,m,n+3)
          if (lmetric.and.lspherical_coords) df=df*r1i(:,1)*sth1i(m,1)
        else
          df=0.
          if (ip<=5) print*, 'deri_3d: Degenerate case in z-direction'
        endif
      endif
!
    endsubroutine deri_3d_inds
!***********************************************************************
    subroutine deri(f,df,i1,i2,i,coefs,lequi)
!
!  i-th derivative operating on a 1-D array
!
!   9-feb-07/axel: adapted from der_main; note that f is not the f array!
!
      real, dimension(*)       , intent(in)  :: f
      real, dimension(-3:3,8,*), intent(in)  :: coefs
      real, dimension(*)       , intent(out) :: df
      integer                  , intent(in)  :: i1,i2,i
      logical                  , intent(in)  :: lequi
!
      integer :: n1

      n1 = i2-i1+1

      if ( n1/=1 ) then
        if ( lequi ) then

          if ( mod(i,2) == 0 ) then
            df(1:n1) =  coefs(0,i,1)* f(i1  :i2  )               &
                      + coefs(1,i,1)*(f(i1+1:i2+1)+f(i1-1:i2-1)) &
                      + coefs(2,i,1)*(f(i1+2:i2+2)+f(i1-2:i2-2)) &
                      + coefs(3,i,1)*(f(i1+3:i2+3)+f(i1-3:i2-3))
          else
            df(1:n1) =  coefs(1,i,1)*(f(i1+1:i2+1)-f(i1-1:i2-1)) &
                      + coefs(2,i,1)*(f(i1+2:i2+2)-f(i1-2:i2-2)) &
                      + coefs(3,i,1)*(f(i1+3:i2+3)-f(i1-3:i2-3))
          endif

        else
          df(1:n1) =  coefs(-3,i,1:n1)*f(i1-3:i2-3) &
                    + coefs(-2,i,1:n1)*f(i1-2:i2-2) &
                    + coefs(-1,i,1:n1)*f(i1-1:i2-1) &
                    + coefs( 0,i,1:n1)*f(i1  :i2  ) &
                    + coefs( 1,i,1:n1)*f(i1+1:i2+1) &
                    + coefs( 2,i,1:n1)*f(i1+2:i2+2) &
                    + coefs( 3,i,1:n1)*f(i1+3:i2+3)
        endif

      else
        df(1:n1) = 0.
        if (ip<=5) print*, 'deri: Degenerate case!'
      endif
!
    endsubroutine deri
!***********************************************************************
    subroutine deri_2d(f,df,j,i,coefs,lequi)
!
!  i-th derivative operating on a 2-D array
!
!   9-feb-07/axel: adapted from der_main; note that f is not the f array!
!
      real, dimension (:,:)     , intent(in)  :: f
      real, dimension (-3:3,8,*), intent(in)  :: coefs
      real, dimension (*)       , intent(out) :: df
      integer                   , intent(in)  :: i,j
      logical                   , intent(in)  :: lequi
!
      integer :: n1

      n1 = size(f,1)

      if ( n1/=1 ) then

        if ( lequi ) then

          if ( mod(i,2) == 0 ) then
            df(1:n1) =  coefs(0,i,1)* f(:,j  )           &
                      + coefs(1,i,1)*(f(:,j+1)+f(:,j-1)) &
                      + coefs(2,i,1)*(f(:,j+2)+f(:,j-2)) &
                      + coefs(3,i,1)*(f(:,j+3)+f(:,j-3))
          else
            df(1:n1) =  coefs(1,i,1)*(f(:,j+1)-f(:,j-1)) &
                      + coefs(2,i,1)*(f(:,j+2)-f(:,j-2)) &
                      + coefs(3,i,1)*(f(:,j+3)-f(:,j-3))
          endif

        else
          df(1:n1) =  coefs(-3,i,1:n1)*f(:,j-3) &
                    + coefs(-2,i,1:n1)*f(:,j-2) &
                    + coefs(-1,i,1:n1)*f(:,j-1) &
                    + coefs( 0,i,1:n1)*f(:,j  ) &
                    + coefs( 1,i,1:n1)*f(:,j+1) &
                    + coefs( 2,i,1:n1)*f(:,j+2) &
                    + coefs( 3,i,1:n1)*f(:,j+3)
        endif

      else
        df(1:n1) = 0.
        if (ip<=5) print*, 'deri_2d: Degenerate case!'
      endif
!
    endsubroutine deri_2d
!***********************************************************************
    subroutine der_z(f,df)
!
!  first z derivative operating on a z-dependent 1-D array
!
!   9-feb-07/axel: adapted from der_main; note that f is not the f array!
!
      real, dimension (mz), intent(in)  :: f
      real, dimension (nz), intent(out) :: df

      call deri(f,df,n1,n2,1,coeffsz,lequidist(3))
!
    endsubroutine der_z
!***********************************************************************
    subroutine der2_z(f,df2)
!
!  z derivative operating on a z-dependent 1-D array
!
!   2-jan-10/axel: adapted from der_z and der_main
!
      real, dimension (mz), intent(in)  :: f
      real, dimension (nz), intent(out) :: df2
!
      call deri(f,df2,n1,n2,2,coeffsz,lequidist(3))
!
    endsubroutine der2_z
!***********************************************************************
    subroutine der_pencil(j,pencil,df)
!
!  Calculate first derivative of any x, y or z pencil.
!
!  01-nov-07/anders: adapted from der
!
      real, dimension (:) :: pencil,df
      integer :: j
!
      intent(in)  :: j, pencil
      intent(out) :: df
!
!  x-derivative
!
      if (j==1) then

        if (size(pencil)/=mx) then
          if (lroot) print*, 'der_pencil: pencil must be of size mx for x derivative'
          call fatal_error('der_pencil','')
        endif

        call deri(pencil,df,l1,l2,1,coeffsx,lequidist(j))

      else if (j==2) then
!
!  y-derivative
!
        if (size(pencil)/=my) then
          if (lroot) print*, 'der_pencil: pencil must be of size my for y derivative'
          call fatal_error('der_pencil','')
        endif

        call deri(pencil,df,m1,m2,1,coeffsy,lequidist(j))
        if (lspherical_coords.or.lcylindrical_coords) df = df*r1i(:,1)        !tbc

      else if (j==3) then
!
!  z-derivative
!
        if (size(pencil)/=mz) then
          if (lroot) print*, 'der_pencil: pencil must be of size mz for z derivative'
          call fatal_error('der_pencil','')
        endif

        call deri(pencil,df,n1,n2,1,coeffsz,lequidist(j))
        if (lspherical_coords) df=df*r1i(:,1)*sth1i(m,1)   !tbc

      else
        if (lroot) print*, 'der_pencil: no such direction j=', j
        call fatal_error('der_pencil','')
      endif
!
    endsubroutine der_pencil
!***********************************************************************
    subroutine distr_der(arr,idir,der,order)
!
!  Dummy
!
    real, dimension(:,:) :: arr, der
    integer :: idir
    integer, optional :: order

    call not_implemented('distr_der','deriv_alt')
    call keep_compiler_quiet(idir)
    call keep_compiler_quiet(arr)
    call keep_compiler_quiet(der)
    
    endsubroutine distr_der
!***********************************************************************
    subroutine der2_main(f,k,df2,j)
!
!  calculate 2nd derivative d^2f_k/dx_j^2
!  accurate to 6th order, explicit, periodic
!  replace cshifts by explicit construction -> x6.5 faster!
!
!   1-oct-97/axel: coded
!   1-apr-01/axel+wolf: pencil formulation
!  25-jun-04/tobi+wolf: adapted for non-equidistant grids
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx)               :: df2
      integer                            :: j,k
!
      intent(in)  :: f,k,j
      intent(out) :: df2
!
!debug      if (loptimise_ders) der_call_count(k,icount_der2,j,1) = & !DERCOUNT
!debug                          der_call_count(k,icount_der2,j,1) + 1 !DERCOUNT
!
      call deri_3d( f(1,1,1,k), df2, 2, j )
!
    endsubroutine der2_main
!***********************************************************************
    subroutine der2_other(f,df2,j)
!
!  calculate 2nd derivative d^2f/dx_j^2 (of scalar f)
!  accurate to 6th order, explicit, periodic
!  replace cshifts by explicit construction -> x6.5 faster!
!
!   1-oct-97/axel: coded
!   1-apr-01/axel+wolf: pencil formulation
!  25-jun-04/tobi+wolf: adapted for non-equidistant grids
!
      real, dimension (mx,my,mz) :: f
      real, dimension (nx)       :: df2
      integer                    :: j
!
      intent(in)  :: f,j
      intent(out) :: df2
!
!debug      if (loptimise_ders) der_call_count(k,icount_der2,j,1) = & !DERCOUNT
!debug                          der_call_count(k,icount_der2,j,1) + 1 !DERCOUNT
!
!
      call deri_3d(f,df2,2,j)
!
    endsubroutine der2_other
!***********************************************************************
    subroutine der2_pencil(j,pencil,df2)
!
!  Calculate 2nd derivative of any x, y or z pencil.
!
!  01-nov-07/anders: adapted from der2
!
! non-cartesian???
!
      real, dimension (:) :: pencil,df2
      integer :: j
!
      intent(in)  :: j, pencil
      intent(out) :: df2
!
!  x-derivative
!
      if (j==1) then
        if (size(pencil)/=mx) then
          if (lroot) print*, 'der2_pencil: pencil must be of size mx for x derivative'
          call fatal_error('der2_pencil','')
        endif

        call deri(pencil,df2,l1,l2,2,coeffsx,lequidist(j))

      else if (j==2) then
!
!  y-derivative
!
        if (size(pencil)/=my) then
          if (lroot) &
              print*, 'der2_pencil: pencil must be of size my for y derivative'
          call fatal_error('der2_pencil','')
        endif

        call deri(pencil,df2,m1,m2,2,coeffsy,lequidist(j))
        if (lspherical_coords.or.lcylindrical_coords) df2 = df2*r1i(:,2)        !tbc

      else if (j==3) then
!
!  z-derivative
!
        if (size(pencil)/=mz) then
          if (lroot) &
              print*, 'der2_pencil: pencil must be of size mz for z derivative'
          call fatal_error('der2_pencil','')
        endif

        call deri(pencil,df2,n1,n2,2,coeffsz,lequidist(j))
        if (lspherical_coords) df2=df2*r1i(:,2)*sth1i(m,2)   !tbc

      else
        if (lroot) print*, 'der2_pencil: no such direction j=', j
        call fatal_error('der2_pencil','')
      endif
!
    endsubroutine der2_pencil
!***********************************************************************
    subroutine der3(f,k,df,j,ignoredx)
!
!  Calculate 3rd derivative of a scalar, get scalar
!
!  10-feb-06/anders: adapted from der5
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx)               :: df
      integer                            :: j,k
      logical, optional                  :: ignoredx
!
      intent(in)  :: f,k,j,ignoredx
      intent(out) :: df
!
!debug      if (loptimise_ders) der_call_count(k,icount_der5,j,1) = & !DERCOUNT
!debug                          der_call_count(k,icount_der5,j,1) + 1 !DERCOUNT
!
      if (present(ignoredx)) then
        call deri_3d(f(1,1,1,k),df,3,j,ignoredx)
      else
        call deri_3d(f(1,1,1,k),df,3,j)
      endif
!
    endsubroutine der3
!***********************************************************************
    subroutine der4(f,k,df,j,ignoredx,upwind)
!
!  Calculate 4th derivative of a scalar, get scalar
!    Used for hyperdiffusion that affects small wave numbers as little as
!  possible (useful for density).
!    The optional flag IGNOREDX is useful for numerical purposes, where
!  you want to affect the Nyquist scale in each direction, independent of
!  the ratios dx:dy:dz.
!
!   8-jul-02/wolf: coded
!   9-dec-03/nils: adapted from der6
!  10-feb-06/anders: corrected sign and factor
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx)               :: df
      integer                            :: j,k
      logical, optional                  :: ignoredx,upwind
!
      intent(in)  :: f,k,j,ignoredx,upwind
      intent(out) :: df
!
!debug      if (loptimise_ders) der_call_count(k,icount_der4,j,1) = & !DERCOUNT
!debug                          der_call_count(k,icount_der4,j,1) + 1 !DERCOUNT
!
      if (present(upwind)) &
        call warning('der4','upwinding not implemented')
!
      if (present(ignoredx)) then
        call deri_3d(f(1,1,1,k),df,4,j,ignoredx)
      else
        call deri_3d(f(1,1,1,k),df,4,j)
      endif
!
    endsubroutine der4
!***********************************************************************
    subroutine der5(f,k,df,j,ignoredx)
!
!  Calculate 5th derivative of a scalar, get scalar
!    Used for hyperdiffusion that affects small wave numbers as little as
!  possible (useful for density).
!    The optional flag IGNOREDX is useful for numerical purposes, where
!  you want to affect the Nyquist scale in each direction, independent of
!  the ratios dx:dy:dz.
!
!  29-oct-04/anders: adapted from der6
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx)               :: df
      integer                            :: j,k
      logical, optional                  :: ignoredx
!
      intent(in)  :: f,k,j,ignoredx
      intent(out) :: df
!
!debug      if (loptimise_ders) der_call_count(k,icount_der5,j,1) = & !DERCOUNT
!debug                          der_call_count(k,icount_der5,j,1) + 1 !DERCOUNT
!
      if (present(ignoredx)) then
        call deri_3d(f(1,1,1,k),df,5,j,ignoredx)
      else
        call deri_3d(f(1,1,1,k),df,5,j)
      endif
!
    endsubroutine der5
!***********************************************************************
    subroutine der6_main(f,k,df,j,ignoredx,upwind)
!
!  Calculate 6th derivative of a scalar, get scalar
!    Used for hyperdiffusion that affects small wave numbers as little as
!  possible (useful for density).
!    The optional flag IGNOREDX is useful for numerical purposes, where
!  you want to affect the Nyquist scale in each direction, independent of
!  the ratios dx:dy:dz.
!    The optional flag UPWIND is a variant thereof, which calculates
!  D^(6)*dx^5/60, which is the upwind correction of centered derivatives.
!
!   8-jul-02/wolf: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: df
      integer :: j,k
      logical, optional :: ignoredx,upwind
!
      intent(in)  :: f,k,j,ignoredx,upwind
      intent(out) :: df
!
!debug      if (loptimise_ders) der_call_count(k,icount_der6,j,1) = & !DERCOUNT
!debug                          der_call_count(k,icount_der6,j,1) + 1 !DERCOUNT
!
      if (present(ignoredx)) then
        if (present(upwind)) then
          call der6_other(f(1,1,1,k),df,j,ignoredx,upwind)
        else
          call der6_other(f(1,1,1,k),df,j,ignoredx)
        endif
      else
        if (present(upwind)) then
          call der6_other(f(1,1,1,k),df,j,UPWIND=upwind)
        else
          call der6_other(f(1,1,1,k),df,j)
        endif
      endif
!
    endsubroutine der6_main
!***********************************************************************
    subroutine der6_pencil(j,pencil,df6)
!
!  Calculate 6th derivative of any x, y or z pencil.
!
      real, dimension (:) :: pencil,df6
      integer :: j
!
      intent(in)  :: j, pencil
      intent(out) :: df6

      call not_implemented('der6:pencil','deriv_alt')
      call keep_compiler_quiet(j)
      call keep_compiler_quiet(pencil)
      call keep_compiler_quiet(df6)

    endsubroutine der6_pencil
!***********************************************************************
    subroutine der2_minmod(f,j,delfk,delfkp1,delfkm1,k)
!
! Calculates gradient of a scalar along the direction j but
! get for the derivatives at the point i-1,i,i+1
!
      intent(in) :: f,k,j
      intent(out) :: delfk,delfkp1,delfkm1
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: delfk,delfkp1,delfkm1,fac
      real, dimension (nx,-1:1) :: delf
      real, dimension (0:nx+1) :: delfx
      integer :: j,k
      integer :: i,ii,ix
      real :: tmp1,tmp2,tmp3
! x-component
      select case (k)
        case(1)
          do i=l1-1,l2+1
            ix=i-nghost
            tmp1= f(i,m,n,j)-f(i-1,m,n,j)
            tmp2=(f(i+1,m,n,j)-f(i-1,m,n,j))/4.
            tmp3= f(i+1,m,n,j)-f(i,m,n,j)
            delfx(ix) = minmod(tmp1,tmp2,tmp3)
          enddo
          do i=1,nx;do ii=-1,1
            delf(i,ii) = delfx(i+ii)
          enddo;enddo
          fac = dx_1(l1:l2)
! y-component
        case(2)
          do i=l1,l2
            ix=i-nghost
            do ii=-1,1
              tmp1=f(i,m+ii,n,j)-f(i,m+ii-1,n,j)
              tmp2=(f(i,m+ii+1,n,j)-f(i,m+ii-1,n,j))/4.
              tmp3=f(i,m+ii+1,n,j)-f(i,m+ii,n,j)
              delf(ix,ii) = minmod(tmp1,tmp2,tmp3)*dy_1(i)
            enddo
          enddo
          fac = dy_1(m)
          if (lspherical_coords.or.lcylindrical_coords) fac = fac*r1i(:,1)
! z-component
        case(3)
          do i=l1,l2
            ix=i-nghost
            do ii=-1,1
              tmp1=f(i,m,n+ii,j)-f(i,m,n+ii-1,j)
              tmp2=(f(i,m,n+ii+1,j)-f(i,m,n+ii-1,j))/4.
              tmp3=f(i,m,n+ii+1,j)-f(i,m,n+ii,j)
              delf(ix,ii) = minmod(tmp1,tmp2,tmp3)
            enddo
          enddo
          fac = dz_1(n)
          if (lspherical_coords) fac = fac*r1i(:,1)*sth1i(m,i)
        case default
          call fatal_error('deriv:der2_minmod','wrong component')
        endselect
        delfkm1 = delf(:,-1)*fac
        delfk   = delf(:, 0)*fac
        delfkp1 = delf(:, 1)*fac
!
    endsubroutine der2_minmod
!***********************************************************************
    real function minmod(a,b,c)
!
!  DOCUMENT ME!
!
      real :: a,b,c
!
      if (((a>0) .and. (b>0) .and. (c>0))) then
        minmod=max(a,b,c)
      elseif (((a<0) .and. (b<0) .and. (c<0))) then
        minmod=min(a,b,c)
      else
        minmod=0.0
      endif
!
    endfunction minmod
!***********************************************************************
    subroutine der6_other(f,df,j,ignoredx,upwind)
!
!  Calculate 6th derivative of a scalar, get scalar
!    Used for hyperdiffusion that affects small wave numbers as little as
!  possible (useful for density).
!    The optional flag IGNOREDX is useful for numerical purposes, where
!  you want to affect the Nyquist scale in each direction, independent of
!  the ratios dx:dy:dz.
!    The optional flag UPWIND is a variant thereof, which calculates
!  D^(6)*dx^5/60, which is the upwind correction of centered derivatives.
!
!   8-jul-02/wolf: coded
!
      real, dimension (mx,my,mz)         :: f
      real, dimension (nx)               :: df
      integer                            :: j
      logical,                  optional :: ignoredx,upwind
!
      intent(in)  :: f,j,ignoredx,upwind
      intent(out) :: df
!
      logical :: igndx,upwnd
!
!debug      if (loptimise_ders) der_call_count(k,icount_der6,j,1) = & !DERCOUNT
!debug                          der_call_count(k,icount_der6,j,1) + 1 !DERCOUNT
!
      if (present(upwind)) then
        upwnd = upwind
        if ( upwnd .and. .not.lequidist(j) ) &
          call fatal_error('der6_other','Not to be used '// &    !tbc
          'for upwinding when grid is non-equidistant')
      else
        upwnd = .false.
        !if (.not.lcartesian_coords) &
        !     call fatal_error('der6_other','in non-cartesian coordinates '// &    !tbc
        !     'just works if upwinding is used')
      endif
!
      if (present(ignoredx)) then
        igndx = ignoredx
        upwnd = .false.   !tbc
      else
        igndx = upwnd
      endif

      call deri_3d(f,df,6,j,igndx,lnometric=upwnd)
!
      if ( upwnd .and. lequidist(j) ) then
        if (j==1) then
          df = df*(dx_1(1)/60.)           ! as dx is ignored for equidistant grids
        elseif (j==2) then
          df = df*(dy_1(1)/60.)
        elseif (j==3) then
          df = df*(dz_1(1)/60.)
        endif
      endif
!
    endsubroutine der6_other
!***********************************************************************
    subroutine derij_main(f,k,df,i,j)
!
!  calculate 2nd derivative with respect to two different directions
!  input: scalar, output: scalar
!  accurate to 6th order, explicit, periodic
!
!   8-sep-01/axel: coded
!  25-jun-04/tobi+wolf: adapted for non-equidistant grids
!  14-nov-06/wolf: implemented bidiagonal scheme
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: df
      integer :: i,j,k
!
      intent(in) :: f,k,i,j
      intent(out) :: df
!
!debug      if (loptimise_ders) der_call_count(k,icount_derij,i,j) = & !DERCOUNT
!debug                          der_call_count(k,icount_derij,i,j) + 1 !DERCOUNT
!
    call derij_other(f(1,1,1,k),df,i,j)
!
    endsubroutine derij_main
!***********************************************************************
    subroutine derij_other(f,df,i,j)
!
!  calculate 2nd derivative with respect to two different directions, i,j
!  input: scalar, output: scalar
!  accurate to 6th order, explicit, periodic
!   8-sep-01/axel: coded
!  25-jun-04/tobi+wolf: adapted for non-equidistant grids
!  14-nov-06/wolf: implemented bidiagonal scheme
!
      real, dimension (mx,my,mz) :: f
      real, dimension (nx)       :: df
      integer                    :: i,j
!
      intent(in)  :: f,i,j
      intent(out) :: df
!
      real, dimension (nx) :: fac
      real, dimension (nx,-3:3) :: work
      integer :: ii
      logical :: lbidiag
!
!debug      if (loptimise_ders) der_call_count(k,icount_derij,i,j) = & !DERCOUNT
!debug                          der_call_count(k,icount_derij,i,j) + 1 !DERCOUNT
!
      if ( i==j ) return     !warning??

      lbidiag = lbidiagonal_derij .and. lequidist(i) .and. lequidist(j)

      if (lbidiag) then
        !
        ! Use bidiagonal mixed-derivative operator, i.e.
        ! employ only the three neighbouring points on each of the four
        ! half-diagonals. This gives 6th-order mixed derivatives as the
        ! version below, but involves just 12 points instead of 36.
        !
        if ( i+j==3 ) then
          if (nxgrid/=1.and.nygrid/=1) then
            fac=(1./720.)*dx_1(1)*dy_1(1)
            df=fac*( 270.*( f(l1+1:l2+1,m+1,n)-f(l1-1:l2-1,m+1,n)  &
                           +f(l1-1:l2-1,m-1,n)-f(l1+1:l2+1,m-1,n)) &
                    - 27.*( f(l1+2:l2+2,m+2,n)-f(l1-2:l2-2,m+2,n)  &
                           +f(l1-2:l2-2,m-2,n)-f(l1+2:l2+2,m-2,n)) &
                    +  2.*( f(l1+3:l2+3,m+3,n)-f(l1-3:l2-3,m+3,n)  &
                           +f(l1-3:l2-3,m-3,n)-f(l1+3:l2+3,m-3,n))  )
          else
            df=0.
            if (ip<=5) print*, 'derij: Degenerate case in x- or y-direction'
          endif
        elseif ( i+j==5 ) then
          if (nygrid/=1.and.nzgrid/=1) then
            fac=(1./720.)*dy_1(1)*dz_1(1)
            df=fac*( 270.*( f(l1:l2,m+1,n+1)-f(l1:l2,m+1,n-1)  &
                           +f(l1:l2,m-1,n-1)-f(l1:l2,m-1,n+1)) &
                    - 27.*( f(l1:l2,m+2,n+2)-f(l1:l2,m+2,n-2)  &
                           +f(l1:l2,m-2,n-2)-f(l1:l2,m-2,n+2)) &
                    +  2.*( f(l1:l2,m+3,n+3)-f(l1:l2,m+3,n-3)  &
                           +f(l1:l2,m-3,n-3)-f(l1:l2,m-3,n+3))  )
          else
            df=0.
            if (ip<=5) print*, 'derij: Degenerate case in y- or z-direction'
          endif
        elseif ( i+j==4 ) then
          if (nzgrid/=1.and.nxgrid/=1) then
            fac=(1./720.)*dz_1(1)*dx_1(1)
            df=fac*( 270.*( f(l1+1:l2+1,m,n+1)-f(l1-1:l2-1,m,n+1)  &
                           +f(l1-1:l2-1,m,n-1)-f(l1+1:l2+1,m,n-1)) &
                    - 27.*( f(l1+2:l2+2,m,n+2)-f(l1-2:l2-2,m,n+2)  &
                           +f(l1-2:l2-2,m,n-2)-f(l1+2:l2+2,m,n-2)) &
                    +  2.*( f(l1+3:l2+3,m,n+3)-f(l1-3:l2-3,m,n+3)  &
                           +f(l1-3:l2-3,m,n-3)-f(l1+3:l2+3,m,n-3))  )
          else
            df=0.
            if (ip<=5) print*, 'derij: Degenerate case in x- or z-direction'
          endif
        endif
!
      else                      ! not using bidiagonal mixed derivatives, necessary for non-uniform grid
        !
        ! This is the old, straight-forward scheme
        !
        if ( i+j==3 ) then
          if (nxgrid/=1.and.nygrid/=1) then

            do ii=-3,+3
              if ( ii/=0 .or. .not.lequidist(2) ) &
                call deri(f(1,m+ii,n),work(1,ii),l1,l2,1,coeffsx,lequidist(1))  ! x-derivative
            enddo

            call deri_2d(work,df,4,1,coeffsy,lequidist(2))                      ! y-derivative

          else
            df=0.
            if (ip<=5) print*, 'derij: Degenerate case in x- or y-direction'
          endif
        elseif ( i+j==5 ) then
          if (nygrid/=1.and.nzgrid/=1) then

            do ii=-3,+3
              if ( ii/=0 .or. .not.lequidist(3) ) &
                call deri_2d(f(l1:l2,:,n+ii),work(1,ii),m,1,coeffsy,lequidist(2))   ! y-derivative   performance?
            enddo

            call deri_2d(work,df,4,1,coeffsz,lequidist(3))                          ! z-derivative
          else
            df=0.
            if (ip<=5) print*, 'derij: Degenerate case in y- or z-direction'
          endif
        elseif ( i+j==4 ) then
          if (nzgrid/=1.and.nxgrid/=1) then

            do ii=-3,+3
              if ( ii/=0 .or. .not.lequidist(3) ) &
                call deri(f(1,m,n+ii),work(1,ii),l1,l2,1,coeffsx,lequidist(1))  ! x-derivative
            enddo

            call deri_2d(work,df,4,1,coeffsz,lequidist(3))                      ! z-derivative
          else
            df=0.
            if (ip<=5) print*, 'derij: Degenerate case in x- or z-direction'
          endif
        endif
!
      endif                     ! bidiagonal derij
!
!  Spherical polars. The comments about "minus extra terms" refer to the
!  presence of extra terms that are being evaluated later in gij_etc.
!
      if (lspherical_coords) then
        if ( i+j==3 ) df=df*r1i(:,1)                 !(minus extra terms)
        if ( i+j==4 ) df=df*r1i(:,1)*sth1i(m,1)      !(minus extra terms)
        if ( i+j==5 ) df=df*r1i(:,2)*sth1i(m,1)      !(minus extra terms)
      endif
!
      if (lcylindrical_coords) then
        if ( i+j==3 .or. i+j==5 ) df=df*r1i(:,1)
      endif
!
    endsubroutine derij_other
!***********************************************************************
    subroutine der5i1j(f,k,df,i,j)
!
!  Calculate 6th derivative with respect to two different directions.
!
!  05-dec-06/anders: adapted from derij
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: df
      integer :: i,j,k
!
      intent(in) :: f,k,i,j
      intent(out) :: df

      real, dimension(nx,-3:3) :: work
      integer :: ii

!debug      if (loptimise_ders) der_call_count(k,icount_derij,i,j) = & !DERCOUNT
!debug                          der_call_count(k,icount_derij,i,j) + 1 !DERCOUNT
!
      df=0.0
      if ((i==1.and.j==1)) then
        if (nxgrid/=1) then
          call der6(f,k,df,j)
        else
          df=0.
          if (ip<=5) print*, 'der5i1j: Degenerate case in x-direction'
        endif
      elseif ((i==2.and.j==2)) then
        if (nygrid/=1) then
          call der6(f,k,df,j)
        else
          df=0.
          if (ip<=5) print*, 'der5i1j: Degenerate case in y-direction'
        endif
      elseif ((i==3.and.j==3)) then
        if (nzgrid/=1) then
          call der6(f,k,df,j)
        else
          df=0.
          if (ip<=5) print*, 'der5i1j: Degenerate case in z-direction'
        endif
      else
        if ((i==1.and.j==2)) then                       ! 5th deriv in x, 1st in y
          if (nxgrid/=1.and.nygrid/=1) then

              do ii=-3,+3
                if ( ii/=0 .or. .not.lequidist(2) ) &
                  call deri(f(1,m+ii,n,k),work(1,ii),l1,l2,5,coeffsx,lequidist(1))      ! x-derivative
              enddo

              call deri_2d(work,df,4,1,coeffsy,lequidist(2))                            ! y-derivative
              if (lspherical_coords.or.lcylindrical_coords) df = df*r1i(:,1)            !(minus extra terms) tbc

          else
            df=0.
            if (ip<=5) print*, 'der5i1j: Degenerate case in x- or y-direction'
          endif
        elseif ((i==2.and.j==1)) then                      ! 5th deriv in y, 1st in x
          if (nygrid/=1.and.nxgrid/=1) then

              do ii=-3,+3
                if ( ii/=0 .or. .not.lequidist(2) ) &
                  call deri(f(1,m+ii,n,k),work(1,ii),l1,l2,1,coeffsx,lequidist(1))      ! x-derivative
              enddo

              call deri_2d(work,df,4,5,coeffsy,lequidist(2))                            ! y-derivative
              if (lspherical_coords.or.lcylindrical_coords) df = df*r1i(:,5)            !(minus extra terms) tbc

          else
            df=0.
            if (ip<=5) print*, 'der5i1j: Degenerate case in y- or x-direction'
          endif
        elseif ((i==1.and.j==3)) then                       ! 5th deriv in x, 1st in z
          if (nxgrid/=1.and.nzgrid/=1) then

              do ii=-3,+3
                if ( ii/=0 .or. .not.lequidist(3) ) &
                  call deri(f(1,m,n+ii,k),work(1,ii),l1,l2,5,coeffsx,lequidist(1))      ! x-derivative
              enddo

              call deri_2d(work,df,4,1,coeffsz,lequidist(3))                            ! z-derivative
              if (lspherical_coords) df = df*r1i(:,1)*sth1i(m,1)                        !(minus extra terms) tbc

          else
            df=0.
            if (ip<=5) print*, 'der5i1j: Degenerate case in x- or z-direction'
          endif
        elseif ((i==3.and.j==1)) then                       ! 5th deriv in z, 1st in x
          if (nzgrid/=1.and.nygrid/=1) then

              do ii=-3,+3
                if ( ii/=0 .or. .not.lequidist(3) ) &
                  call deri(f(1,m,n+ii,k),work(1,ii),l1,l2,1,coeffsx,lequidist(1))      ! x-derivative
              enddo

              call deri_2d(work,df,4,5,coeffsz,lequidist(3))                            ! z-derivative
              if (lspherical_coords) df = df*r1i(:,5)*sth1i(m,5)                        ! (minus extra terms) tbc

          else
            df=0.
            if (ip<=5) print*, 'der5i1j: Degenerate case in z- or x-direction'
          endif
        elseif ((i==2.and.j==3)) then                       ! 5th deriv in y, 1st in z
          if (nygrid/=1.and.nzgrid/=1) then

              do ii=-3,+3
                if ( ii/=0 .or. .not.lequidist(3) ) &
                  call deri_2d(f(l1:l2,:,n+ii,k),work(1,ii),m,5,coeffsy,lequidist(2))   ! y-derivative   performance?
              enddo

              call deri_2d(work,df,4,1,coeffsz,lequidist(3))                            ! z-derivative
              if (lspherical_coords  ) df = df*r1i(:,6)*sth1i(m,1)                      ! (minus extra terms) tbc
              if (lcylindrical_coords) df = df*r1i(:,5)                                 ! (minus extra terms) tbc

          else
            df=0.
            if (ip<=5) print*, 'der5i1j: Degenerate case in y- or z-direction'
          endif
        elseif ((i==3.and.j==2)) then                       ! 5th deriv in z, 1st in y
          if (nzgrid/=1.and.nygrid/=1) then

              do ii=-3,+3
                if ( ii/=0 .or. .not.lequidist(3) ) &
                  call deri_2d(f(l1:l2,:,n+ii,k),work(1,ii),m,1,coeffsy,lequidist(2))   ! y-derivative   performance?
              enddo

              call deri_2d(work,df,4,5,coeffsz,lequidist(3))                            ! z-derivative
              if (lspherical_coords  ) df = df*r1i(:,6)*sth1i(m,5)                      ! (minus extra terms) tbc
              if (lcylindrical_coords) df = df*r1i(:,1)                                 ! (minus extra terms) tbc

          else
            df=0.
            if (ip<=5) print*, 'der5i1j: Degenerate case in z- or y-direction'
          endif
        else
          print*, 'der5i1j: no such value for i,j=', i, j
          call fatal_error('der5i1j','')
        endif
!
      endif
!
    endsubroutine der5i1j
!***********************************************************************
    subroutine der_upwind1st(f,uu,k,df,j)
!
!  First order upwind derivative of variable
!
!  Useful for advecting non-logarithmic variables
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: uu
      real, dimension (nx) :: df
      integer :: j,k,l
!
      intent(in)  :: f,uu,k,j
      intent(out) :: df
!
!debug      if (loptimise_ders) der_call_count(k,icount_der_upwind1st,j,1) = & !DERCOUNT
!debug                          der_call_count(k,icount_der_upwind1st,j,1) + 1 !DERCOUNT
!
      if (lspherical_coords.or.lcylindrical_coords) &
          call fatal_error('der_upwind1st','NOT IMPLEMENTED for non-cartesian grid')
!
      if (j == 1) then
        if (nxgrid /= 1) then
          do l=1,nx
            if (uu(l,1) > 0.) then
              df(l) = (f(nghost+l  ,m,n,k) - f(nghost+l-1,m,n,k))*dx_1(nghost+l)
            else
              df(l) = (f(nghost+l+1,m,n,k) - f(nghost+l  ,m,n,k))*dx_1(nghost+l)
            endif
          enddo
        else
          df=0.
          if (ip<=5) print*, 'der_upwind1st: Degenerate case in x-direction'
        endif
      elseif (j == 2) then
        if (nygrid /= 1) then
          do l=1,nx
            if (uu(l,2) > 0.) then
              df(l) = (f(nghost+l,m  ,n,k) - f(nghost+l,m-1,n,k))*dy_1(m)
            else
              df(l) = (f(nghost+l,m+1,n,k) - f(nghost+l,m  ,n,k))*dy_1(m)
            endif
          enddo
        else
          df=0.
          if (ip<=5) print*, 'der_upwind1st: Degenerate case in y-direction'
        endif
      elseif (j == 3) then
        if (nzgrid /= 1) then
          do l=1,nx
            if (uu(l,3) > 0.) then
              df(l) = (f(nghost+l,m,n,k  ) - f(nghost+l,m,n-1,k))*dz_1(n)
            else
              df(l) = (f(nghost+l,m,n+1,k) - f(nghost+l,m,n,k  ))*dz_1(n)
            endif
          enddo
        else
          df=0.
          if (ip<=5) print*, 'der_upwind1st: Degenerate case in z-direction'
        endif
      endif
!
    endsubroutine der_upwind1st
!***********************************************************************
    subroutine der_onesided_4_slice_main(f,sgn,k,df,pos,j)
!
!   Calculate x/y/z-derivative on a yz/xz/xy-slice at gridpoint pos.
!   Uses a one-sided 4th order stencil.
!   sgn = +1 for forward difference, sgn = -1 for backwards difference.
!
!   Because of its original intended use in relation to solving
!   characteristic equations on boundaries (NSCBC), this sub should
!   return only PARTIAL derivatives, NOT COVARIANT. Applying the right
!   scaling factors and connection terms should instead be done when
!   solving the characteristic equations.
!
!   7-jul-08/arne: coded.
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (:,:) :: df
      real :: fac
      integer :: pos,k,sgn,j
!
      intent(in)  :: f,k,pos,sgn,j
      intent(out) :: df
!
      if (j==1) then
        if (nxgrid/=1) then
          fac=1./12.*dx_1(pos)
          df = fac*(-sgn*25*f(pos,m1:m2,n1:n2,k)&
                    +sgn*48*f(pos+sgn*1,m1:m2,n1:n2,k)&
                    -sgn*36*f(pos+sgn*2,m1:m2,n1:n2,k)&
                    +sgn*16*f(pos+sgn*3,m1:m2,n1:n2,k)&
                    -sgn*3 *f(pos+sgn*4,m1:m2,n1:n2,k))
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice: Degenerate case in x-directder_onesided_4_sliceion'
        endif
      elseif (j==2) then
        if (nygrid/=1) then
          fac=1./12.*dy_1(pos)
          df = fac*(-sgn*25*f(l1:l2,pos,n1:n2,k)&
                    +sgn*48*f(l1:l2,pos+sgn*1,n1:n2,k)&
                    -sgn*36*f(l1:l2,pos+sgn*2,n1:n2,k)&
                    +sgn*16*f(l1:l2,pos+sgn*3,n1:n2,k)&
                    -sgn*3 *f(l1:l2,pos+sgn*4,n1:n2,k))
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice: Degenerate case in y-direction'
        endif
      elseif (j==3) then
        if (nzgrid/=1) then
          fac=1./12.*dz_1(pos)
          df = fac*(-sgn*25*f(l1:l2,m1:m2,pos,k)&
                    +sgn*48*f(l1:l2,m1:m2,pos+sgn*1,k)&
                    -sgn*36*f(l1:l2,m1:m2,pos+sgn*2,k)&
                    +sgn*16*f(l1:l2,m1:m2,pos+sgn*3,k)&
                    -sgn*3 *f(l1:l2,m1:m2,pos+sgn*4,k))
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice: Degenerate case in z-direction'
        endif
      endif
    endsubroutine der_onesided_4_slice_main
!***********************************************************************
    subroutine der_onesided_4_slice_main_pt(f,sgn,k,df,lll,mmm,nnn,j)
!
!  made using der_onesided_4_slice_main. One sided derivative is calculated
!  at one point (lll,mmm,nnn)
!
!  15-oct-09/Natalia: coded.
!
      real, dimension (mx,my,mz,mfarray) :: f
      real  :: df
      real :: fac
      integer :: pos,lll,mmm,nnn,k,sgn,j
!
      intent(in)  :: f,k,lll,mmm,nnn,sgn,j
      intent(out) :: df
!
      if (j==1) then
       pos=lll
        if (nxgrid/=1) then
          fac=1./12.*dx_1(pos)
          df = fac*(-sgn*25*f(pos,mmm,nnn,k)&
                    +sgn*48*f(pos+sgn*1,mmm,nnn,k)&
                    -sgn*36*f(pos+sgn*2,mmm,nnn,k)&
                    +sgn*16*f(pos+sgn*3,mmm,nnn,k)&
                    -sgn*3 *f(pos+sgn*4,mmm,nnn,k))
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice: Degenerate case in x-directder_onesided_4_sliceion'
        endif
      elseif (j==2) then
       pos=mmm
        if (nygrid/=1) then
          fac=1./12.*dy_1(pos)
          df = fac*(-sgn*25*f(lll,pos,nnn,k)&
                    +sgn*48*f(lll,pos+sgn*1,nnn,k)&
                    -sgn*36*f(lll,pos+sgn*2,nnn,k)&
                    +sgn*16*f(lll,pos+sgn*3,nnn,k)&
                    -sgn*3 *f(lll,pos+sgn*4,nnn,k))
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice: Degenerate case in y-direction'
        endif
      elseif (j==3) then
       pos=nnn
        if (nzgrid/=1) then
          fac=1./12.*dz_1(pos)
          df = fac*(-sgn*25*f(lll,mmm,pos,k)&
                    +sgn*48*f(lll,mmm,pos+sgn*1,k)&
                    -sgn*36*f(lll,mmm,pos+sgn*2,k)&
                    +sgn*16*f(lll,mmm,pos+sgn*3,k)&
                    -sgn*3 *f(lll,mmm,pos+sgn*4,k))
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice: Degenerate case in z-direction'
        endif
      endif
    endsubroutine der_onesided_4_slice_main_pt
!***********************************************************************
    subroutine der_onesided_4_slice_other(f,sgn,df,pos,j)
!
!   Calculate x/y/z-derivative on a yz/xz/xy-slice at gridpoint pos.
!   Uses a one-sided 4th order stencil.
!   sgn = +1 for forward difference, sgn = -1 for backwards difference.
!
!   Because of its original intended use in relation to solving
!   characteristic equations on boundaries (NSCBC), this sub should
!   return only PARTIAL derivatives, NOT COVARIANT. Applying the right
!   scaling factors and connection terms should instead be done when
!   solving the characteristic equations.
!
!   7-jul-08/arne: coded.
!
      real, dimension (mx,my,mz) :: f
      real, dimension (:,:) :: df
      real :: fac
      integer :: pos,sgn,j
!
      intent(in)  :: f,pos,sgn,j
      intent(out) :: df
!
      if (j==1) then
        if (nxgrid/=1) then
          fac=1./12.*dx_1(pos)
          df = fac*(-sgn*25*f(pos,m1:m2,n1:n2)&
                    +sgn*48*f(pos+sgn*1,m1:m2,n1:n2)&
                    -sgn*36*f(pos+sgn*2,m1:m2,n1:n2)&
                    +sgn*16*f(pos+sgn*3,m1:m2,n1:n2)&
                    -sgn*3 *f(pos+sgn*4,m1:m2,n1:n2))
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice: Degenerate case in x-direction'
        endif
      elseif (j==2) then
        if (nygrid/=1) then
          fac=1./12.*dy_1(pos)
          df = fac*(-sgn*25*(f(l1:l2,pos,n1:n2)-f(l1:l2,pos+sgn*1,n1:n2))&
                    +sgn*23*(f(l1:l2,pos+sgn*1,n1:n2)-f(l1:l2,pos+sgn*2,n1:n2))&
                    -sgn*13*(f(l1:l2,pos+sgn*2,n1:n2)-f(l1:l2,pos+sgn*3,n1:n2))&
                    +sgn*3*(f(l1:l2,pos+sgn*3,n1:n2)-f(l1:l2,pos+sgn*4,n1:n2)))
!
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice: Degenerate case in y-direction'
        endif
      elseif (j==3) then
        if (nzgrid/=1) then
          fac=1./12.*dz_1(pos)
          df = fac*(-sgn*25*f(l1:l2,m1:m2,pos)&
                    +sgn*48*f(l1:l2,m1:m2,pos+sgn*1)&
                    -sgn*36*f(l1:l2,m1:m2,pos+sgn*2)&
                    +sgn*16*f(l1:l2,m1:m2,pos+sgn*3)&
                    -sgn*3 *f(l1:l2,m1:m2,pos+sgn*4))
!
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice: Degenerate case in z-direction'
        endif
      endif
!
    endsubroutine der_onesided_4_slice_other
!***********************************************************************
    subroutine der_onesided_4_slice_other_pt(f,sgn,df,lll,mmm,nnn,j)
!
!  One sided derivative is calculated
!  at one point (lll,mmm,nnn).
!
!  15-oct-09/Natalia: coded.
!  15-oct-09/axel: changed file name to shorter version
!
      real, dimension (mx,my,mz) :: f
      real :: df
      real :: fac
      integer :: pos,lll,mmm,nnn,sgn,j
!
      intent(in)  :: f,lll,mmm,nnn,sgn,j
      intent(out) :: df
!
      if (j==1) then
        pos=lll
        if (nxgrid/=1) then
          fac=1./12.*dx_1(pos)
          df = fac*(-sgn*25*f(pos,mmm,nnn)&
                    +sgn*48*f(pos+sgn*1,mmm,nnn)&
                    -sgn*36*f(pos+sgn*2,mmm,nnn)&
                    +sgn*16*f(pos+sgn*3,mmm,nnn)&
                    -sgn*3 *f(pos+sgn*4,mmm,nnn))
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice_other_pt: Degenerate case in x-direction'
        endif
      elseif (j==2) then
        pos=mmm
        if (nygrid/=1) then
          fac=1./12.*dy_1(pos)
          df = fac*(-sgn*25*(f(lll,pos,nnn)-f(lll,pos+sgn*1,nnn))&
                    +sgn*23*(f(lll,pos+sgn*1,nnn)-f(lll,pos+sgn*2,nnn))&
                    -sgn*13*(f(lll,pos+sgn*2,nnn)-f(lll,pos+sgn*3,nnn))&
                    +sgn* 3*(f(lll,pos+sgn*3,nnn)-f(lll,pos+sgn*4,nnn)))
!
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice_other_pt: Degenerate case in y-direction'
        endif
      elseif (j==3) then
        pos=nnn
        if (nzgrid/=1) then
          fac=1./12.*dz_1(pos)
          df = fac*(-sgn*25*f(lll,mmm,pos)&
                    +sgn*48*f(lll,mmm,pos+sgn*1)&
                    -sgn*36*f(lll,mmm,pos+sgn*2)&
                    +sgn*16*f(lll,mmm,pos+sgn*3)&
                    -sgn*3 *f(lll,mmm,pos+sgn*4))
!
        else
          df=0.
          if (ip<=5) print*, 'der_onesided_4_slice_other_pt: Degenerate case in z-direction'
        endif
      endif
    endsubroutine der_onesided_4_slice_other_pt
!***********************************************************************
   subroutine calc_coeffs( grid, coeffs )

   real, dimension(-2:3), intent(in)  :: grid
   real, dimension(-3:3,8), intent(out) :: coeffs

   real :: h1, h2, h3, hm1, hm2, hm3, h12, h23, h13, hm12, hm23, hm13, &
           h1m1, h1m12, h1m13, h12m1, h13m1, h12m12, h13m12, h12m13, h13m13

   h1    = grid(1);   h2    = grid(2); h3    = grid(3)
   hm1   = grid(0);   hm2   = grid(-1);hm3   = grid(-2)
   h12   = h1 + h2;   h23   = h2+h3;   h13   = h12+h3
   hm12  = hm1 + hm2; hm23  = hm2+hm3; hm13  = hm12+hm3
   h1m1  = h1 + hm1;  h1m12 = h1+hm12;      h1m13 = h1 + hm13
   h12m1 = h12 + hm1; h12m12 = h12 + hm12; h12m13 = h12 + hm13
   h13m1 = h13 + hm1; h13m12 = h13 + hm12; h13m13 = h13 + hm13

   coeffs(:,1) = (/ -h1*h12*h13*hm1*hm12/(hm3*hm13*hm23*h1m13*h12m13*h13m13), &
                     h1*h12*h13*hm1*hm13/(hm2*hm3*hm12*h1m12*h12m12*h13m12),  &
                    -h1*h12*h13*hm12*hm13/(hm1*hm2*hm23*h1m1*h12m1*h13m1),    &
                     0.,                                                      &
                     h12*h13*hm1*hm12*hm13/(h1*h2*h23*h1m1*h1m12*h1m13),      &
                    -h1*h13*hm1*hm12*hm13/(h2*h3*h12*h12m1*h12m12*h12m13),    &
                     h1*h12*hm1*hm12*hm13/(h3*h13*h23*h13m1*h13m12*h13m13)     /)

   coeffs(0,1) = -sum(coeffs(:,1))

   coeffs(:,2) = 2.*(/ (h12*h13*hm1*hm12 + h1*(h13*hm1*hm12 + h12*(hm1*hm12 - h13*(hm1 + hm12))))/       &
                       (hm3*hm13*hm23*h1m13*h12m13*h13m13),                                              &
!
                      -(h12*h13*hm1*hm13 + h1*(h13*hm1*hm13 + h12*(hm1*hm13 - h13*(hm1 + hm13))))/       &
                       (hm2*hm3*hm12*h1m12*h12m12*h13m12),                                               &
!
                       (h12*h13*hm12*hm13 + h1*(h13*hm12*hm13 + h12*(hm12*hm13 - h13*(hm12 + hm13))))/   &
                       (hm1*hm2*hm23*h1m1*h12m1*h13m1),                                                  &
!
                       0.,                                                                               &
!
                      -(h13*hm1*hm12*hm13 + h12*(hm1*hm12*hm13 - h13*(hm12*hm13 + hm1*(hm12 + hm13))))/  &
                       (h1*h2*h23*h1m1*h1m12*h1m13),                                                     &
!
                       (h13*hm1*hm12*hm13 + h1*(hm1*hm12*hm13 -  h13*(hm12*hm13 + hm1 * (hm12 + hm13))))/&
                       (h2*h3*h12*h12m1*h12m12*h12m13),                                                  &
!
                      -(h12*hm1*hm12*hm13 + h1*(hm1*hm12*hm13 - h12*(hm12*hm13 + hm1*(hm12 + hm13))))/   &
                       (h3*h13*h23*h13m1*h13m12*h13m13)                                                   /)

   coeffs(0,2) = -sum(coeffs(:,2))

   coeffs(:,3) = 6.*(/-(h13*hm1*hm12 + h12*(hm1*hm12 - h13*(hm1 + hm12)) +              &
                        h1*(hm1*hm12 - h13*(hm1 + hm12) + h12*(h13 - hm1 - hm12)))/     &
                       (hm3*hm13*hm23*h1m13*h12m13*h13m13),                             &
!
                       (h13*hm1*hm13 + h12*(hm1*hm13 - h13*(hm1 + hm13)) +              &
                        h1*(hm1*hm13 - h13*(hm1 + hm13) + h12*(h13 - hm1 - hm13)))/     &
                       (hm2*hm3*hm12*h1m12*h12m12*h13m12),                              &
!
                      -(h13*hm12*hm13 + h12*(hm12*hm13 - h13*(hm12 + hm13)) +           &
                        h1*(hm12*hm13 - h13*(hm12 + hm13) + h12*(h13 - hm12 - hm13)))/  &
                       (hm1*hm2*hm23*h1m1*h12m1*h13m1),                                 &
!
                       0.,                                                              &
!
                       (hm1*hm12*hm13 - h13*(hm12*hm13 + hm1*(hm12 + hm13)) -           &
                        h12*(hm12*hm13 + hm1*(hm12 + hm13) - h13*(hm1 + hm12 + hm13)))/ &
                       (h1*h2*h23*h1m1*h1m12*h1m13),                                    &
!
                      -(hm1*hm12*hm13 - h13*(hm12*hm13 + hm1*(hm12 + hm13)) -           &
                        h1*(hm12*hm13 + hm1*(hm12 + hm13) - h13*(hm1 + hm12 + hm13)))/  &
                       (h2*h3*h12*h12m1*h12m12*h12m13),                                 &
!
                       (hm1*hm12*hm13 - h12*(hm12*hm13 + hm1*(hm12 + hm13)) -           &
                        h1*(hm12*hm13 + hm1*(hm12 + hm13) - h12*(hm1 + hm12 + hm13)))/  &
                       (h3*h13*h23*h13m1*h13m12*h13m13)                                  /)

   coeffs(0,3) = -sum(coeffs(:,3))

   coeffs(:,4) = 24.*(/  (-h13*(hm1+hm12) + hm1*hm12  +(h1+h12)*(h13-hm1-hm12) + h1*h12)/ &
                         (hm3*hm13*hm23*h1m13*h12m13*h13m13), &
                        -(-h13*(hm1+hm13) + hm1*hm13  +(h1+h12)*(h13-hm1-hm13) + h1*h12)/ &
                         (hm2*hm3*hm12*h1m12*h12m12*h13m12),  &
                         (-h13*(hm12+hm13)+ hm12*hm13 +(h1+h12)*(h13-hm12-hm13)+ h1*h12)/ &
                         (hm1*hm2*hm23*h1m1*h12m1*h13m1),     &
                         0.,                                  &
                         (hm1*(hm12+hm13)+ hm12*hm13 -(h12+h13)*(hm1+hm12+hm13)+ h12*h13)/ &
                         (h1*h2*h23*h1m1*h1m12*h1m13),       &
                        -(hm1*(hm12+hm13)+ hm12*hm13 -(h1+h13)*(hm1+hm12+hm13) + h1*h13 )/ &
                         (h2*h3*h12*h12m1*h12m12*h12m13),    &
                         (hm1*(hm12+hm13)+ hm12*hm13 -(h1+h12)*(hm1+hm12+hm13) + h1*h12 )/ &
                         (h3*h13*h23*h13m1*h13m12*h13m13)     /)

   coeffs(0,4) = -sum(coeffs(:,4))

   coeffs(:,5) = 120.*(/-(h1 + h12 + h13 - hm1  - hm12)/(hm3*hm13*hm23*h1m13*h12m13*h13m13), &
                         (h1 + h12 + h13 - hm1  - hm13)/(hm2*hm3*hm12*h1m12*h12m12*h13m12),  &
                        -(h1 + h12 + h13 - hm12 - hm13)/(hm1*hm2*hm23*h1m1*h12m1*h13m1),     &
                          0.,                                                                &
                        -(h12 + h13 - hm1 - hm12 - hm13)/(h1*h2*h23*h1m1*h1m12*h1m13),       &
                         (h1  + h13 - hm1 - hm12 - hm13)/(h2*h3*h12*h12m1*h12m12*h12m13),    &
                        -(h1  + h12 - hm1 - hm12 - hm13)/(h3*h13*h23*h13m1*h13m12*h13m13)     /)

   coeffs(0,5) = -sum(coeffs(:,5))

   coeffs(:,6) = 720.*(/ 1./(hm3*hm13*hm23*h1m13*h12m13*h13m13),&
                        -1./(hm2*hm3*hm12*h1m12*h12m12*h13m12), &
                         1./(hm1*hm2*hm23*h1m1*h12m1*h13m1),    &
                         0.,                                    &
                         1./(h1*h2*h23*h1m1*h1m12*h1m13),       &
                        -1./(h2*h3*h12*h12m1*h12m12*h12m13),    &
                         1./(h3*h13*h23*h13m1*h13m12*h13m13)     /)

   coeffs(0,6) = -sum(coeffs(:,6))

!upwind 1st deriv, 5th order, u>0

   coeffs((/-3,-2,-1,1,2,3/),7) = (h1*h12*hm1*hm12*hm13/720.)*coeffs((/-3,-2,-1,1,2,3/),6)
   coeffs(0,7) = -sum(coeffs(:,7))

   coeffs(:,7) = (h1*h12*hm1*hm12*hm13/720.)*coeffs(:,6)

!upwind 1st deriv, 5th order, u<0, sign already inverted!   tbc

   coeffs((/-3,-2,-1,1,2,3/),8) = (h1*h12*h13*hm1*hm12/720.)*coeffs((/-3,-2,-1,1,2,3/),6)
   coeffs(0,8) = -sum(coeffs(:,8))

   coeffs(:,8) = (h1*h12*h13*hm1*hm12/720.)*coeffs(:,6)

   endsubroutine calc_coeffs
!*************************************************************************************************************
   subroutine calc_coeffs_1( grid, coeffs )

   real, dimension(-1:2), intent(in)  :: grid
   real, dimension(-2:2,2), intent(out) :: coeffs

   real :: h1, h2, hm1, hm2, h12, hm12, h1m1, h1m12, h12m1, h12m12

   h1    = grid(1);   h2    = grid(2)  ; hm1  = grid(0);   hm2   = grid(-1)
   h12   = h1 + h2;   hm12  = hm1 + hm2; h1m1 = h1 + hm1;  h1m12 = h1+hm12
   h12m1 = h12 + hm1; h12m12 = h12 + hm12

! 2nd order 1st derivative

   coeffs(:,1) = (/ 0., -h1/(hm1*h1m1), 0., hm1/(h1*h1m1), 0. /)
   coeffs(0,1) = -sum(coeffs(:,1))

! 4th order 1st dervative

   coeffs(:,2) = (/ h1*h12*hm1/(hm2*hm12*h1m12*h12m12), -h1*h12*hm12/(hm1*hm2*h1m1*h12m1), 0.,  &
                   +h12*hm1*hm12/(h1*h2*h1m1*h1m12),    -h1*hm1*hm12/(h2*h12*h12m1*h12m12)       /)
   coeffs(0,2) = -sum(coeffs(:,2))

   endsubroutine calc_coeffs_1
!*************************************************************************************************************
   real function derivative( func, i, coeffs, ldiff )

   integer                , intent(in)           :: i
   real, dimension(-2:*)  , intent(in)           :: func
   real, dimension(-3:3,*), intent(in)           :: coeffs
   logical,                 intent(in), optional :: ldiff

   real, dimension(-3:3) :: funcl

   if ( present(ldiff) ) then
     funcl = func(i-3:i+3) - func(i)
     derivative = sum(funcl*coeffs(:,i))
   else
     derivative = sum(func(i-3:i+3)*coeffs(:,i))
   endif

   endfunction derivative
!*************************************************************************************************************
   real function derivative_1( func, i, coeffs, ldiff )

   integer              , intent(in)           :: i
   real, dimension(-2:*), intent(in)           :: func
   real, dimension(-2:2), intent(in)           :: coeffs
   logical,               intent(in), optional :: ldiff

   real, dimension(-2:2) :: funcl

   if ( present(ldiff) ) then
     funcl = func(i-2:i+2) - func(i)
     derivative_1 = sum(funcl*coeffs)
   else
     derivative_1 = sum(func(i-2:i+2)*coeffs)
   endif

   endfunction derivative_1
!*************************************************************************************************************
   logical function heatflux_deriv_x( f, inh, fac, topbot )

     real, dimension(mx,my,mz,mfarray), intent(INOUT):: f
     real, dimension(my,mz)           , intent(IN)   :: inh
     real                             , intent(IN)   :: fac
     integer                          , intent(IN)   :: topbot

     real, dimension(:,:), allocatable :: work_yz
     integer                           :: i,ic,stat,ll, im, ia, ie, is, ima

     heatflux_deriv_x = .false.

     allocate(work_yz(my,mz),stat=stat)                                 !save this by working in f?
     if (stat>0) call fatal_error('heatflux_boundcond_x', &
         'Could not allocate memory for work_yz')

     if (topbot==1) then
       ll=l1; ia=1; ie=3; is=1
       ic=1;
     elseif (topbot==2) then
       ll=l2
       ic=nx; ia=-1; ie=-3; is=-1
     else
       call fatal_error('heatflux_boundcond_x', &
         'Illegal value for parameter "topbot"')
     endif

     do im=ia,ie,is

       ima = abs(im)
       work_yz = inh                                                      ! Fbot/(K*cs2)

       if (ima<3) then
         do i=-ima,ima                                                    ! dlnrho, 2nd+4th order
           if (ldensity_nolog) then
             work_yz = work_yz + coeffsx_1(i,ima,topbot)*log(f(ll+i,:,:,irho))
           else
             work_yz = work_yz + coeffsx_1(i,ima,topbot)*f(ll+i,:,:,ilnrho)
           endif
         enddo
       else
         do i=-ima,ima                                                    ! dlnrho, 6th order
           if (ldensity_nolog) then
             work_yz = work_yz + coeffsx(i,1,ic)*log(f(ll+i,:,:,irho))
           else
             work_yz = work_yz + coeffsx(i,1,ic)*f(ll+i,:,:,ilnrho)
           endif
         enddo
       endif

       work_yz = fac*work_yz

       if (ima<3) then

         do i=-im+is,im,is                                                ! dss, 2nd+4th order
           work_yz = work_yz + coeffsx_1(i,ima,topbot)*f(ll+i,:,:,iss)
         enddo

         f(ll-im,:,:,iss) = -work_yz/coeffsx_1(-im,ima,topbot)

       else

         do i=-im+is,im,is                                                ! dss, 6th order
           work_yz = work_yz + coeffsx(i,1,ic)*f(ll+i,:,:,iss)
         enddo

         f(ll-im,:,:,iss) = -work_yz/coeffsx(-im,1,ic)

       endif
     enddo

     heatflux_deriv_x = .true.

   endfunction heatflux_deriv_x
!*************************************************************************************************************
   subroutine test( coord, dc_1, nc, coefs )

   integer :: mc, nc, nc1, nc2
   real, dimension(nc+2*nghost) :: coord, dc_1
   real, dimension(-3:3,8,nc) :: coefs
   real, dimension(-2:2,2,2) :: coefs_1

   real, parameter, dimension(0:6) :: c=(/1.d0,1.1d0,1.2d0,0.d0,0.d0,0.d0,0.d0/)
   real, dimension(nc+2*nghost) :: func

   integer :: i,j
   logical :: t1

   t1=.false.
   goto 1
   entry test_1( coord, dc_1, nc, coefs_1 ) ! please use proper and clear interfaces [Bourdin.KIS]

   t1=.true.

 1 mc = nc+2*nghost

   do i=1,mc
     func(i) = c(6)
     do j=5,0,-1
       func(i) = func(i)*coord(i) + c(j)
     enddo
   enddo

   nc = mc-2*nghost; nc1 = nghost+1; nc2 = nc1+nc-1

   if (t1) then
     print*, derivative_1( func, 1, coefs_1(-2,1,1) )
     print*, c(1) + 2.e0*c(2)*coord(nc1) + 3.e0*c(3)*coord(nc1)**2 + 4.e0*c(4)*coord(nc1)**3 &
                  + 5.e0*c(5)*coord(nc1)**4 + 6.e0*c(6)*coord(nc1)**5

     print*, derivative_1( func, nc, coefs_1(:,1,2) )
     print*, c(1) + 2.e0*c(2)*coord(nc2) + 3.e0*c(3)*coord(nc2)**2 + 4.e0*c(4)*coord(nc2)**3 &
                  + 5.e0*c(5)*coord(nc2)**4 + 6.e0*c(6)*coord(nc2)**5
     stop
   endif

   !print*, (derivative( func, i, coefs(:,3,:) ), i=1,nc)
   !print*, '==='
   !print*, 6.e0*c(3) + 24.e0*c(4)*coord(nc1:nc2) + 60.e0*c(5)*coord(nc1:nc2)**2 + 120.e0*c(6)*coord(nc1:nc2)**3
   !print*, '==='
   !print*, (derivative( func, i, coefs(:,2,:) ), i=1,nc)
   !print*, '==='
   !print*, 2.e0*c(2) + 6.e0*c(3)*coord(nc1:nc2) + 12.e0*c(4)*coord(nc1:nc2)**2 + 20.e0*c(5)*coord(nc1:nc2)**3 + 30.e0*c(6)*coord(nc1:nc2)**4
   !print*, '==='
   print*, (derivative( func, i, coefs(:,1,:) ), i=1,nc)
   print*, '==='
   print*, c(1) + 2.e0*c(2)*coord(nc1:nc2) + 3.e0*c(3)*coord(nc1:nc2)**2 + 4.e0*c(4)*coord(nc1:nc2)**3 &
                + 5.e0*c(5)*coord(nc1:nc2)**4 + 6.e0*c(6)*coord(nc1:nc2)**5
   print*, '==='
   print*, ( dc_1(i)*(45.0*(func(i+1)-func(i-1))-9.0*(func(i+2)-func(i-2))+func(i+3)-func(i-3))/60., i=nc1,nc2 )


   stop
   endsubroutine test
!*************************************************************************************************************
   endmodule Deriv