! $Id$
!
!  This module contains useful subroutines.
!
!  Rules:
!
!    - Please do not put very specific subroutines here. If a subroutine
!      is only needed by a single module, put it directly in that module.
!
!    - Please DO NOT use large arrays or global arrays
!      [e.g. of size (mx,my,mz) or (nxgrid,nygrid,nzgrid)]
!
module Sub
!
  use Cdata
  use Cparam
  use Messages
!
  implicit none
!
  private
!
  public :: step,stepdown,der_stepdown
  public :: ylm, ylm_other
  public :: kronecker_delta
!
  public :: identify_bcs, parse_bc, parse_bc_rad, parse_bc_radg
  public :: inverse_parse_bc
!
  public :: poly
  public :: blob, vecout
  public :: cubic_step, cubic_der_step, quintic_step, quintic_der_step, erfunc
  public :: sine_step, interp1
  public :: hypergeometric2F1
  public :: gamma_function
  public :: one_minus_exp
!
  public :: get_nseed
  public :: get_where
!
  public :: grad, grad5, div, div_mn, curl, curli, curl_mn, curl_other
  public :: curl_horizontal
  public :: div_other
  public :: gij, g2ij, gij_etc
  public :: gijl_symmetric
  public :: symmetrise3x3_ut2lt
  public :: der_step
  public :: der6_step
  public :: u_dot_grad, h_dot_grad
  public :: u_dot_grad_alt
  public :: nou_dot_grad_scl
  public :: u_dot_grad_mat
  public :: del2, del2v, del2v_etc, del2fj,d2fi_dxj,del2fi_dxjk
  public :: del2m3x3_sym
  public :: del4v, del4, del2vi_etc
  public :: del6v, del6, del6_other, del6fj, del6fjv
  public :: gradf_upw1st, doupwind
  public :: matrix2linarray, linarray2matrix
!
  public :: dot, dot2, dot_mn, dot_mn_sv, dot_mn_sm, dot2_mn, dot_add, dot_sub, dot2fj
  public :: dot_mn_vm, dot_mn_vm_trans, div_mn_2tensor, trace_mn
  public :: dyadic2
  public :: cross, cross_mn, cross_mixed
  public :: sum_mn, max_mn
  public :: multsv, multsv_add, multsv_mn, multsv_mn_add
  public :: multvs, multvv_mat
  public :: multmm_sc
  public :: multm2, multm2_mn, multm2_sym, multm2_sym_mn
  public :: multmv, multmv_mn, multmv_transp
  public :: mult_matrix
!
  public :: read_line_from_file, remove_file, control_file_exists
  public :: noform
!
  public :: update_snaptime, read_snaptime
  public :: parse_shell
  public :: get_radial_distance, power_law
!
  public :: max_for_dt,unit_vector
!
  public :: write_dx_general, numeric_precision, wdim, rdim
  public :: write_prof, write_xprof, write_zprof, remove_zprof
  public :: read_zprof
!
  public :: tensor_diffusion_coef
!
  public :: smooth_kernel, despike
  public :: smooth, get_smooth_kernel
!
  public :: ludcmp, lubksb
  public :: bspline_basis, bspline_interpolation, bspline_precondition
!
  public :: gij_psi, gij_psi_etc
  public :: xlocation, ylocation, zlocation, location_in_proc, position
  public :: register_report_aux
  public :: fourier_single_mode
  public :: remove_mean,global_mean
  public :: insert
  public :: find_max_fvec, find_rms_fvec, find_xyrms_fvec
  public :: finalize_aver
  public :: fseek_pos, parallel_file_exists, parallel_count_lines, read_namelist
  public :: meanyz
  public :: calc_all_diff_fluxes
  public :: periodic_fold_back
  public :: lower_triangular_index
!
  interface poly                ! Overload the `poly' function
    module procedure poly_0
    module procedure poly_1
    module procedure poly_3
  endinterface
!
  interface grad                 ! Overload the `grad' function
    module procedure grad_main   ! grad of an 'mvar' variable
    module procedure grad_other  ! grad of another field (mx,my,mz)
  endinterface
!
  interface del2                 ! Overload the `del2' function
    module procedure del2_main
    module procedure del2_other
  endinterface
!
  interface cross
    module procedure cross_mn
    module procedure cross_0
    module procedure cross_mixed
  endinterface
!
  interface u_dot_grad
    module procedure u_dot_grad_scl
    module procedure u_dot_grad_vec
  endinterface
!
  interface u_dot_grad_alt
    module procedure u_dot_grad_scl_alt
    module procedure u_dot_grad_vec_alt
  endinterface
!
  interface h_dot_grad
    module procedure h_dot_grad_scl
    module procedure h_dot_grad_vec
  endinterface
!
  interface dot
    module procedure dot_mn_sv
    module procedure dot_mn
    module procedure dot_0
  endinterface
!
  interface dot2
    module procedure dot2_mn
    module procedure dot2_0
  endinterface
!
  interface dot_add
    module procedure dot_mn_add
  endinterface
!
  interface dot_sub
    module procedure dot_mn_sub
  endinterface
!
  interface multsv
    module procedure multsv_mn
  endinterface
!
  interface multsv_add
    module procedure multsv_add_mn
  endinterface
!
  interface multvs
    module procedure multvs_mn
  endinterface
!
  interface multvv_mat
    module procedure multvv_mat_mn
  endinterface
!
  interface multmm_sc
    module procedure multmm_sc_mn
  endinterface
!
  interface multm2
    module procedure multm2_mn
  endinterface
!
  interface multm2_sym
    module procedure multm2_sym_mn
  endinterface
!
  interface multmv_transp
    module procedure multmv_mn_transp
  endinterface
!
  interface multmv
    module procedure multmv_mn
  endinterface
!
  interface max_for_dt
    module procedure max_for_dt_nx_nx
    module procedure max_for_dt_1_nx
    module procedure max_for_dt_1_1_1_nx
  endinterface
!
  interface step
    module procedure step_scalar
    module procedure step_vector
  endinterface
!
  interface cubic_step
    module procedure cubic_step_pt
    module procedure cubic_step_mn
  endinterface
!
  interface cubic_der_step
    module procedure cubic_der_step_pt
    module procedure cubic_der_step_mn
  endinterface
!
  interface quintic_step
    module procedure quintic_step_pt
    module procedure quintic_step_mn
  endinterface
!
  interface quintic_der_step
    module procedure quintic_der_step_pt
    module procedure quintic_der_step_mn
  endinterface
!
  interface sine_step
    module procedure sine_step_pt
    module procedure sine_step_mn
  endinterface
!
  interface power_law
     module procedure power_law_pt
     module procedure power_law_mn
  endinterface
!
  interface insert                 ! Overload the 'insert' function
    module procedure insert_carray
    module procedure insert_carray_mult
    module procedure insert_rarray
  endinterface
!
  interface finalize_aver
    module procedure finalize_aver_1D
    module procedure finalize_aver_2D
    module procedure finalize_aver_3D
    module procedure finalize_aver_4D
  endinterface finalize_aver
!
  interface meanyz
     module procedure meanyz_s
     module procedure meanyz_v
  endinterface
!
!  extended intrinsic operators to do some scalar/vector pencil arithmetic
!  Tobi: Array valued functions do seem to be slower than subroutines,
!        hence commented out for the moment.
!
!  public :: operator(*),operator(+),operator(/),operator(-)
!
!  interface operator (*)
!    module procedure pencil_multiply1
!    module procedure pencil_multiply2
!  endinterface
!
!  interface operator (+)
!    module procedure pencil_add1
!    module procedure pencil_add2
!  endinterface
!
!  interface operator (/)
!    module procedure pencil_divide1
!    module procedure pencil_divide2
!  endinterface
!
!  interface operator (-)
!    module procedure pencil_subtract1
!    module procedure pencil_subtract2
!  endinterface
!
!ajwm Commented pending a C replacement
!  INTERFACE getenv
!    SUBROUTINE GETENV (VAR, VALUE)
!      CHARACTER(LEN=*) VAR, VALUE
!    endsubroutine
!  END INTERFACE
!
  real, dimension(7,7,7), parameter :: smth_kernel = reshape((/ &
 6.03438e-15,9.07894e-11,1.24384e-08,5.46411e-08,1.24384e-08,9.07894e-11,5.03438e-15,9.07894e-11,2.21580e-07,9.14337e-06,&
 2.69243e-05,9.14337e-06,2.21580e-07,9.07894e-11, 1.24384e-08, 9.14337e-06, 0.000183649, 0.000425400, 0.000183649, 9.14337e-06,&
 1.24384e-08,5.46411e-08,2.69243e-05,0.000425400, 0.000909623, 0.000425400, 2.69243e-05, 5.46411e-08, 1.24384e-08, 9.14337e-06,&
 0.000183649,0.000425400,0.000183649,9.14337e-06, 1.24384e-08, 9.07894e-11, 2.21580e-07, 9.14337e-06, 2.69243e-05, 9.14337e-06,&
 2.21580e-07,9.07894e-11,5.03438e-15,9.07894e-11, 1.24384e-08, 5.46411e-08, 1.24384e-08, 9.07894e-11, 5.03438e-15, 9.07894e-11,&
 2.21580e-07,9.14337e-06,2.69243e-05,9.14337e-06, 2.21580e-07, 9.07894e-11, 2.21580e-07, 7.31878e-05, 0.000909623, 0.00179548, &
 0.000909623,7.31878e-05,2.21580e-07,9.14337e-06, 0.000909623, 0.00550289, 0.00854438, 0.00550289, 0.000909623, 9.14337e-06,   &
 2.69243e-05, 0.00179548, 0.00854438,0.0122469, 0.00854438,      0.00179548, 2.69243e-05, 9.14337e-06, 0.000909623, 0.00550289,&
  0.00854438, 0.00550289,0.000909623,9.14337e-06, 2.21580e-07,  7.31878e-05, 0.000909623, 0.00179548, 0.000909623, 7.31878e-05,&
 2.21580e-07,9.07894e-11,2.21580e-07,9.14337e-06, 2.69243e-05, 9.14337e-06, 2.21580e-07, 9.07894e-11, 1.24384e-08, 9.14337e-06,&
 0.000183649,0.000425400,0.000183649,9.14337e-06, 1.24384e-08, 9.14337e-06, 0.000909623, 0.00550289, 0.00854438, 0.00550289,   &
 0.000909623,9.14337e-06,0.000183649,0.00550289, 0.0162043,    0.0197919, 0.0162043, 0.00550289, 0.000183649, 0.000425400,     &
  0.00854438,  0.0197919,  0.0223153,0.0197919, 0.00854438,        0.000425400, 0.000183649, 0.00550289, 0.0162043, 0.0197919, &
   0.0162043, 0.00550289,0.000183649,9.14337e-06, 0.000909623,  0.00550289, 0.00854438, 0.00550289, 0.000909623, 9.14337e-06,  &
 1.24384e-08,9.14337e-06,0.000183649,0.000425400, 0.000183649, 9.14337e-06, 1.24384e-08, 5.46411e-08, 2.69243e-05, 0.000425400,&
 0.000909623,0.000425400,2.69243e-05,5.46411e-08, 2.69243e-05, 0.00179548, 0.00854438, 0.0122469, 0.00854438, 0.00179548,      &
 2.69243e-05,0.000425400, 0.00854438,0.0197919, 0.0223153,      0.0197919, 0.00854438, 0.000425400, 0.000909623, 0.0122469,    &
   0.0223153,  0.0232260,  0.0223153,0.0122469, 0.000909623,       0.000425400, 0.00854438, 0.0197919, 0.0223153, 0.0197919,   &
  0.00854438,0.000425400,2.69243e-05,0.00179548, 0.00854438,   0.0122469, 0.00854438, 0.00179548, 2.69243e-05, 5.46411e-08,    &
 2.69243e-05,0.000425400,0.000909623,0.000425400, 2.69243e-05, 5.46411e-08, 1.24384e-08, 9.14337e-06, 0.000183649, 0.000425400,&
 0.000183649,9.14337e-06,1.24384e-08,9.14337e-06, 0.000909623, 0.00550289, 0.00854438, 0.00550289, 0.000909623, 9.14337e-06,   &
 0.000183649, 0.00550289,  0.0162043,0.0197919, 0.0162043,        0.00550289, 0.000183649, 0.000425400, 0.00854438, 0.0197919, &
   0.0223153,  0.0197919, 0.00854438,0.000425400, 0.000183649,    0.00550289, 0.0162043, 0.0197919, 0.0162043, 0.00550289,     &
 0.000183649,9.14337e-06,0.000909623,0.00550289, 0.00854438,   0.00550289, 0.000909623, 9.14337e-06, 1.24384e-08, 9.14337e-06, &
 0.000183649,0.000425400,0.000183649,9.14337e-06, 1.24384e-08, 9.07894e-11, 2.21580e-07, 9.14337e-06, 2.69243e-05, 9.14337e-06,&
 2.21580e-07,9.07894e-11,2.21580e-07,7.31878e-05, 0.000909623, 0.00179548, 0.000909623, 7.31878e-05, 2.21580e-07, 9.14337e-06, &
 0.000909623, 0.00550289, 0.00854438,0.00550289, 0.000909623,    9.14337e-06, 2.69243e-05, 0.00179548, 0.00854438, 0.0122469,  &
  0.00854438, 0.00179548,2.69243e-05,9.14337e-06, 0.000909623,  0.00550289, 0.00854438, 0.00550289, 0.000909623, 9.14337e-06,  &
 2.21580e-07,7.31878e-05,0.000909623,0.00179548, 0.000909623,  7.31878e-05, 2.21580e-07, 9.07894e-11, 2.21580e-07, 9.14337e-06,&
 2.69243e-05,9.14337e-06,2.21580e-07,9.07894e-11, 5.03438e-15, 9.07894e-11, 1.24384e-08, 5.46411e-08, 1.24384e-08, 9.07894e-11,&
 5.03438e-15,9.07894e-11,2.21580e-07,9.14337e-06, 2.69243e-05, 9.14337e-06, 2.21580e-07, 9.07894e-11, 1.24384e-08, 9.14337e-06,&
 0.000183649,0.000425400,0.000183649,9.14337e-06, 1.24384e-08, 5.46411e-08, 2.69243e-05, 0.000425400, 0.000909623, 0.000425400,&
 2.69243e-05,5.46411e-08,1.24384e-08,9.14337e-06, 0.000183649, 0.000425400, 0.000183649, 9.14337e-06, 1.24384e-08, 9.07894e-11,&
 2.21580e-07,9.14337e-06,2.69243e-05,9.14337e-06, 2.21580e-07, 9.07894e-11, 5.03438e-15, 9.07894e-11, 1.24384e-08, 5.46411e-08,&
 1.24384e-08,9.07894e-11,5.03438e-15 /), (/ 7, 7, 7 /))
!
  contains
!
!***********************************************************************
    subroutine max_mn(a,res)
!
!  successively calculate maximum of a, which is supplied at each call.
!  Start from scratch if lfirstpoint=.true.
!
!   1-apr-01/axel+wolf: coded
!
      real, dimension (nx), intent(in) :: a
      real, intent(inout) :: res
!
      if (lfirstpoint) then
        res=maxval(a)
      else
        res=max(res,maxval(a))
      endif
!
    endsubroutine max_mn
!***********************************************************************
    subroutine mean_mn(a,res)
!
!  Successively calculate mean of a, which is supplied at each call.
!  Start from zero if lfirstpoint=.true.
!
!  17-dec-01/wolf: coded
!  20-jun-07/dhruba:adapted for spherical polar coordinate system
!
      real, dimension (nx) :: a
      real :: res
      integer :: isum
!
      if (lfirstpoint) then
        if (lspherical_coords) then
          res = 0.
          do isum=l1,l2
            res = res+x(isum)*x(isum)*sinth(m)*a(isum)*a(isum)
          enddo
        else
          res=sum(a*1.D0)     ! sum at double precision to improve accuracy
        endif
      else
        if (lspherical_coords) then
          do isum=l1,l2
            res = res+x(isum)*x(isum)*sinth(m)*a(isum)*a(isum)
          enddo
        else
          res=res+sum(a*1.D0)
        endif
      endif
!
      if (lcylindrical_coords) &
          call fatal_error('mean_mn','not implemented for cylindrical')
!
    endsubroutine mean_mn
!***********************************************************************
    subroutine meanyz_s(f,iif,mean,lexp)
!
!  Calculates mean of variable at subscript iif of f over all x coordinate surfaces 
!  by proper integration.
!
!  7-jun-15/MR: coded
!
      use General, only: loptest

      real, dimension (mx,my,mz,mfarray) :: f
      integer :: iif
      real, dimension(nx) :: mean
      logical, optional :: lexp

      integer :: nn, ll, lll

      do ll=l1,l2
        lll=ll-l1+1
        mean(lll)=0.
        do nn=n1,n2
          if (loptest(lexp)) then
            mean(lll) = mean(lll) + dVol_z(nn)*sum(dVol_y(m1:m2)*exp(f(ll,m1:m2,nn,iif)))
          else
            mean(lll) = mean(lll) + dVol_z(nn)*sum(dVol_y(m1:m2)*f(ll,m1:m2,nn,iif))
          endif
        enddo
        mean(lll) = mean(lll)/Area_yz
      enddo

      call finalize_aver(nprocyz,23,mean)

    endsubroutine meanyz_s
!***********************************************************************
    subroutine meanyz_v(f,iif1,iif2,mean,lexp)
!
!  Calculates mean of variables at subscripts iif1 through iif2 of f over all x coordinate surfaces 
!  by proper integration.
!
!  7-jun-15/MR: coded
!
      use General, only: loptest

      real, dimension (mx,my,mz,mfarray) :: f
      integer :: iif1,iif2
      real, dimension(nx,iif2-iif1+1) :: mean
      logical, optional :: lexp

      integer :: nn, ll, lll, iif, iil

      do iif=iif1,iif2
        iil=iif-iif1+1
        do ll=l1,l2
          lll=ll-l1+1
          mean(lll,iil)=0.
          do nn=n1,n2
            if (loptest(lexp)) then
              mean(lll,iil) = mean(lll,iil) + dVol_z(nn)*sum(dVol_y(m1:m2)*exp(f(ll,m1:m2,nn,iif)))
            else
              mean(lll,iil) = mean(lll,iil) + dVol_z(nn)*sum(dVol_y(m1:m2)*f(ll,m1:m2,nn,iif))
            endif
          enddo
          mean(lll,iil) = mean(lll,iil)/Area_yz
        enddo
      enddo

      call finalize_aver(nprocyz,23,mean)

    endsubroutine meanyz_v
!***********************************************************************
    subroutine rms_mn(a,res)
!
!  Successively calculate rms of a, which is supplied at each call.
!  Start from zero if lfirstpoint=.true.
!
!   1-apr-01/axel+wolf: coded
!
      real, dimension (nx) :: a
      real :: res
      integer :: isum
!
      if (lfirstpoint) then
        if (lspherical_coords) then
          res = 0.
          do isum=l1,l2
            res = res+x(isum)*x(isum)*sinth(m)*a(isum)*a(isum)
          enddo
        else
          res=sum(a**2)
        endif
      else
        if (lspherical_coords) then
          do isum=l1,l2
            res = res+x(isum)*x(isum)*sinth(m)*a(isum)*a(isum)
          enddo
        else
          res=res+sum(a**2)
        endif
      endif
!
      if (lcylindrical_coords) &
          call fatal_error('rms_mn','not implemented for cylindrical')
!
    endsubroutine rms_mn
!***********************************************************************
    subroutine rms2_mn(a2,res)
!
!  Successively calculate rms of a, with a2=a^2 being supplied at each
!  call.
!  Start from zero if lfirstpoint=.true.
!
!   1-apr-01/axel+wolf: coded
!
      real, dimension (nx) :: a2
      real :: res
      integer :: isum
!
      if (lfirstpoint) then
        if (lspherical_coords) then
          res = 0.
          do isum=l1,l2
            res = res+x(isum)*x(isum)*sinth(m)*a2(isum)
          enddo
        else
          res=sum(a2)
        endif
      else
        if (lspherical_coords) then
          do isum=l1,l2
            res = res+x(isum)*x(isum)*sinth(m)*a2(isum)
          enddo
        else
          res=res+sum(a2)
        endif
      endif
!
      if (lcylindrical_coords) &
          call fatal_error('rms2_mn','not implemented for cylindrical')
!
    endsubroutine rms2_mn
!***********************************************************************
    subroutine sum_mn(a,res)
!
!  Successively calculate the sum over all points of a, which is supplied
!  at each call.
!  Start from zero if lfirstpoint=.true.
!
!   1-apr-01/axel+wolf: coded
!
      real, dimension (nx) :: a
      real :: res
      integer :: isum
!
      if (lfirstpoint) then
        if (lspherical_coords) then
          res = 0.
          do isum=l1,l2
            res = res+x(isum)*x(isum)*sinth(m)*a(isum)
          enddo
        else
          res=sum(a)
        endif
      else
        if (lspherical_coords) then
          do isum=l1,l2
            res = res+x(isum)*x(isum)*sinth(m)*a(isum)
          enddo
        else
          res=res+sum(a)
        endif
      endif
!
      if (lcylindrical_coords) &
          call fatal_error('sum_mn','not implemented for cylindrical')
!
    endsubroutine sum_mn
!***********************************************************************
    subroutine dot_mn(a,b,c,ladd)
!
!  Dot product, c=a.b, on pencil arrays
!
!   3-apr-01/axel+gitta: coded
!  24-jun-08/MR: ladd added for incremental work
!
      use General, only: loptest
!
      real, dimension (:,:) :: a,b
      real, dimension (:) :: c
!
      logical, optional :: ladd
!
      intent(in) :: a,b,ladd
      intent(inout) :: c
!
      integer :: i
      logical :: l0
!
      l0 = .not.loptest(ladd)

      do i=1,size(a,2)
        if (l0) then
          c=a(:,i)*b(:,i)
          l0=.false.
        else
          c=c+a(:,i)*b(:,i)
        endif
      enddo
!
    endsubroutine dot_mn
!***********************************************************************
    subroutine vec_dot_3tensor(a,b,c)
!
!  Dot product of a vector with 3 tensor,
!   c_ij = a_k b_ijk
!
!  28-aug-08/dhruba : coded
!
      real, dimension (nx,3) :: a
      real, dimension (nx,3,3) :: c
      real, dimension (nx,3,3,3) :: b
      integer :: i,j,k
!
      intent(in) :: a,b
      intent(out) :: c
!
      do i=1,3
        do j=1,3
          c(:,i,j) = 0.
          do k=1,3
            c(:,i,j) = c(:,i,j) + a(:,k)*b(:,i,j,k)
          enddo
        enddo
      enddo
!
    endsubroutine vec_dot_3tensor
!***********************************************************************
    subroutine contract_jk3(a,c)
!
!  Contracts the jk of a_ijk.
!
!  20-aug-08/dhruba: coded
!
      real, dimension (nx,3,3,3) :: a
      real, dimension (nx,3) :: c
      integer :: i,j,k
!
      intent(in) :: a
      intent(out) :: c
!
      c=0
      do i=1,3; do j=1,3; do k=1,3
        c(:,i)=c(:,i)+a(:,i,j,k)
      enddo; enddo; enddo
!
    endsubroutine contract_jk3
!***********************************************************************
    subroutine matrix2linarray(mm,aa)
!
! converts a 3X3 matrix to an array of length 9
!
!18-sep-15/dhruba: coded
!
      real,dimension(3,3),intent(in) :: mm
      real,dimension(9),intent(out) :: aa
!
      integer :: i,j,ij,k
      ij=0
      do i=1,3
        do j=1,3
          ij=ij+1
          aa(ij) = mm(i,j)
        enddo
      enddo
!
    endsubroutine matrix2linarray
!***********************************************************************
    subroutine linarray2matrix(aa,mm)
!
! converts a 3X3 matrix to an array of length 9
!
!18-sep-15/dhruba: coded
!
      real,dimension(9),intent(in) :: aa
      real,dimension(3,3),intent(out) :: mm
!
      integer :: i,j,ij,k
      ij=0
      do i=1,3
        do j=1,3
          ij=ij+1
          mm(i,j) = aa(ij)
       enddo
      enddo
!
    endsubroutine linarray2matrix
!***********************************************************************
    subroutine dot_mn_sv(a,b,c)
!
!  Dot product, c=a.b, between non-pencilized vector and pencil array.
!
!  10-oct-06/axel: coded
!
      real, dimension (3)    :: a
      real, dimension (nx,3) :: b
      real, dimension (nx)   :: c
!
      intent(in) :: a,b
      intent(out) :: c
!
      c=a(1)*b(:,1)+a(2)*b(:,2)+a(3)*b(:,3)
!
    endsubroutine dot_mn_sv
!***********************************************************************
    subroutine dot_mn_sm(a,b,c)
!
!  Dot product, c=a.b, between non-pencilized vector and pencil matrix.
!
!  10-oct-06/axel: coded
!
      real, dimension (3)      :: a
      real, dimension (nx,3,3) :: b
      real, dimension (nx,3)   :: c
      integer :: i
!
      intent(in) :: a,b
      intent(out) :: c
!
      do i=1,3
        c(:,i)=a(1)*b(:,i,1)+a(2)*b(:,i,2)+a(3)*b(:,i,3)
      enddo
!
    endsubroutine dot_mn_sm
!***********************************************************************
    subroutine dot_mn_vm(a,b,c)
!
!  Dot product, c=a.b, between pencil vector and pencil matrix.
!
!  10-oct-06/axel: coded
!
      real, dimension (nx,3)   :: a
      real, dimension (nx,3,3) :: b
      real, dimension (nx,3)   :: c
      integer :: i
!
      intent(in) :: a,b
      intent(out) :: c
!
      do i=1,3
        c(:,i)=a(:,1)*b(:,i,1)+a(:,2)*b(:,i,2)+a(:,3)*b(:,i,3)
      enddo
!
    endsubroutine dot_mn_vm
!***********************************************************************
    subroutine dot_mn_vm_trans(a,b,c)
!
!  Dot product, c=a.b, between pencil vector and pencil matrix.
!  I think the name of dot_mn_vm is not right and should have been transposed.
!
!  10-oct-06/axel: coded
!
      real, dimension (nx,3)      :: a
      real, dimension (nx,3,3) :: b
      real, dimension (nx,3)   :: c
      integer :: i
!
      intent(in) :: a,b
      intent(out) :: c
!
      do i=1,3
        c(:,i)=a(:,1)*b(:,1,i)+a(:,2)*b(:,2,i)+a(:,3)*b(:,3,i)
      enddo
!
    endsubroutine dot_mn_vm_trans
!***********************************************************************
    subroutine dot_0(a,b,c)
!
!  Dot product, c=a.b, of two simple 3-d arrays.
!
!  11-mar-04/wolf: coded
!
      real, dimension (:) :: a,b
      real :: c
!
      intent(in) :: a,b
      intent(out) :: c
!
      c = dot_product(a,b)
!
    endsubroutine dot_0
!***********************************************************************
    subroutine dot2_mn(a,b,fast_sqrt,precise_sqrt)
!
!  Dot product with itself, to calculate max and rms values of a vector.
!  FAST_SQRT is only correct for ~1e-18 < |a| < 1e18 (for single precision);
!  PRECISE_SQRT works for full range.
!
!  29-sep-97/axel: coded
!   1-apr-01/axel: adapted for cache-efficient sub-array formulation
!  25-jun-05/bing: added optional args for calculating |a|
!
      real, dimension (nx,3) :: a
      real, dimension (nx) :: b,a_max
      logical, optional :: fast_sqrt,precise_sqrt
      logical :: fast_sqrt1,precise_sqrt1
!
      intent(in) :: a,fast_sqrt,precise_sqrt
      intent(out) :: b
!
!     ifc treats these variables as SAVE so we need to reset
      if (present(fast_sqrt)) then
        fast_sqrt1=fast_sqrt
      else
        fast_sqrt1=.false.
      endif
!
      if (present(precise_sqrt)) then
        precise_sqrt1=precise_sqrt
      else
        precise_sqrt1=.false.
      endif
!
!  Rescale by factor a_max before taking sqrt.
!  In single precision this increases the dynamic range from 1e18 to 1e36.
!  To avoid division by zero when calculating a_max, we add tini.
!
      if (precise_sqrt1) then
        a_max=tini+maxval(abs(a),dim=2)
        b=(a(:,1)/a_max)**2+(a(:,2)/a_max)**2+(a(:,3)/a_max)**2
        b=a_max*sqrt(b)
      else
        b=a(:,1)**2+a(:,2)**2+a(:,3)**2
        if (fast_sqrt1) b=sqrt(b)
      endif
!
    endsubroutine dot2_mn
!***********************************************************************
    subroutine dot2_0(a,b)
!
!  Dot product, c=a.b, of two simple 3-d arrays.
!
!  11-mar-04/wolf: coded
!
      real, dimension (:) :: a
      real :: b
!
      intent(in) :: a
      intent(out) :: b
!
      b = dot_product(a,a)
!
    endsubroutine dot2_0
!***********************************************************************
    subroutine dot_mn_add(a,b,c)
!
!  Dot product, add to previous value.
!
!  11-nov-02/axel: adapted from dot_mn
!
      real, dimension (nx,3) :: a,b
      real, dimension (nx) :: c
!
      intent(in) :: a,b
      intent(inout) :: c
!
      c=c+a(:,1)*b(:,1)+a(:,2)*b(:,2)+a(:,3)*b(:,3)
!
    endsubroutine dot_mn_add
!***********************************************************************
    subroutine dot_mn_sub(a,b,c)
!
!  Dot product, subtract from previous value.
!
!  21-jul-03/axel: adapted from dot_mn_sub
!
      real, dimension (nx,3) :: a,b
      real, dimension (nx) :: c
!
      intent(in) :: a,b
      intent(inout) :: c
!
      c=c-(a(:,1)*b(:,1)+a(:,2)*b(:,2)+a(:,3)*b(:,3))
!
    endsubroutine dot_mn_sub
!***********************************************************************
    subroutine dot2fj(a,vec,b)
!
!  Dot product with itself, multiplied by anisotropic factor.
!
!  20-dec-12/wlad: adapted from dot2_mn
!
      real, dimension (nx,3) :: a
      real, dimension (3)    :: vec
      real, dimension (nx)   :: b
      integer :: j
!
      intent(in) :: a,vec
      intent(out) :: b
!
      b=0.
      do j=1,3
         b=b+vec(j)*a(:,j)**2
      enddo
!
    endsubroutine dot2fj
!***********************************************************************
    subroutine dyadic2(a,b)
!
!  Dyadic product with itself.
!
!  24-jan-09/axel: coded
!
      real, dimension (nx,3) :: a
      real, dimension (nx,3,3) :: b
!
      intent(in) :: a
      intent(out) :: b
!
!  diagonal components
!
      b(:,1,1)=a(:,1)**2
      b(:,2,2)=a(:,2)**2
      b(:,3,3)=a(:,3)**2
!
!  upper off-diagonal components
!
      b(:,1,2)=a(:,1)*a(:,2)
      b(:,1,3)=a(:,1)*a(:,3)
      b(:,2,3)=a(:,2)*a(:,3)
!
!  lower off-diagonal components
!
      b(:,2,1)=b(:,1,2)
      b(:,3,1)=b(:,1,3)
      b(:,3,2)=b(:,2,3)
!
    endsubroutine dyadic2
!***********************************************************************
    subroutine trace_mn(a,b)
!
!  Trace of a matrix.
!
!   3-apr-01/axel+gitta: coded
!
      real, dimension (nx,3,3) :: a
      real, dimension (nx) :: b
!
      intent(in) :: a
      intent(out) :: b
!
      b=a(:,1,1)+a(:,2,2)+a(:,3,3)
!
    endsubroutine trace_mn
!***********************************************************************
    subroutine multvv_mat_mn(a,b,c)
!
!  Vector multiplied with vector, gives matrix.
!
!   21-dec-01/nils: coded
!   16-jul-02/nils: adapted from pencil_mpi
!
      real, dimension (nx,3) :: a,b
      real, dimension (nx,3,3) :: c
      integer :: i,j
!
      do i=1,3; do j=1,3
        c(:,i,j)=a(:,j)*b(:,i)
      enddo; enddo
!
    endsubroutine multvv_mat_mn
!***********************************************************************
    subroutine multmm_sc_mn(a,b,c)
!
!  Matrix multiplied with matrix, gives scalar.
!
!   21-dec-01/nils: coded
!   16-jul-02/nils: adapted from pencil_mpi
!
      real, dimension (nx,3,3) :: a,b
      real, dimension (nx) :: c
      integer :: i,j
!
      c=0.0
      do i=1,3; do j=1,3
        c=c+a(:,i,j)*b(:,i,j)
      enddo; enddo
!
    endsubroutine multmm_sc_mn
!***********************************************************************
    subroutine mult_matrix(a,b,c)
!
!  Matrix multiplication of two pencil variables.
!
      real, dimension (nx,3,3) :: a,b
      real, dimension (nx,3,3) :: c
      integer :: ix
!
!  24-feb-11/dhruba: using the f90 command matmul
!
      do ix=1,nx
        c(ix,:,:)=matmul(a(ix,:,:),b(ix,:,:))
      enddo
!
    endsubroutine mult_matrix
!***********************************************************************
    subroutine multm2_mn(a,b)
!
!  Matrix squared, gives scalar.
!
!  11-nov-02/axel: adapted from multmm_sc_mn
!
      real, dimension (nx,3,3) :: a
      real, dimension (nx) :: b
      integer :: i,j
!
      b=0.0
      do i=1,3; do j=1,3
        b=b+a(:,i,j)**2
      enddo; enddo
!
    endsubroutine multm2_mn
!***********************************************************************
    subroutine multm2_sym_mn(a,b)
!
!  Symmetric matrix squared, gives scalar.
!
!  24-aug-2011/Bourdin.KIS: adapted from multm2_mn
!
      real, dimension (nx,3,3), intent(in) :: a
      real, dimension (nx), intent(out) :: b
!
      integer :: i, j
!
      b = a(:,1,1)**2
      do i = 2, 3
        b = b + a(:,i,i)**2
        do j = 1, i-1
          b = b + 2 * a(:,i,j)**2
        enddo
      enddo
!
    endsubroutine multm2_sym_mn
!***********************************************************************
    subroutine multmv_mn(a,b,c,ladd)
!
!  Matrix multiplied with vector, gives vector.
!
!  C_i = A_{i,j} B_j
!
!   3-apr-01/axel+gitta: coded
!  24-jun-08/MR: ladd added for incremental work
!
      use General, only: loptest
!
      real, dimension (nx,3,3) :: a
      real, dimension (nx,3) :: b,c
      real, dimension (nx) :: tmp
      integer :: i,j
      logical, optional :: ladd
!
      intent(in) :: a,b,ladd
      intent(out) :: c
!
      do i=1,3
!
        j=1
        tmp=a(:,i,j)*b(:,j)
        do j=2,3
          tmp=tmp+a(:,i,j)*b(:,j)
        enddo
!
        if (loptest(ladd)) then
          c(:,i)=c(:,i)+tmp
        else
          c(:,i)=tmp
        endif
!
      enddo
!
    endsubroutine multmv_mn
!***********************************************************************
    subroutine multmv_mn_transp(a,b,c,ladd)
!
!  Transposed matrix multiplied with vector, gives vector.
!  Could have called multvm_mn, but this may not be clear enough.
!
!  C_i = A_{j,i} B_j
!
!  21-jul-03/axel: adapted from multmv_mn
!  24-jun-08/MR: ladd added for incremental work
!
      use General, only: loptest
!
      real, dimension (nx,3,3) :: a
      real, dimension (nx,3) :: b,c
      real, dimension (nx) :: tmp
      integer :: i,j
      logical, optional :: ladd
!
      intent(in) :: a,b,ladd
      intent(inout) :: c
!
      do i=1,3
        j=1
        tmp=a(:,j,i)*b(:,j)
        do j=2,3
          tmp=tmp+a(:,j,i)*b(:,j)
        enddo
!
        if (loptest(ladd)) then
          c(:,i)=c(:,i)+tmp
        else
          c(:,i)=tmp
        endif
!
      enddo
!
    endsubroutine multmv_mn_transp
!***********************************************************************
    subroutine multsv_mn(a,b,c)
!
!  Vector multiplied with scalar, gives vector.
!
!  22-nov-01/nils erland: coded
!  10-oct-03/axel: a is now the scalar (now consistent with old routines)
!  24-jun-08/MR: ladd added for incremental work
!  28-feb-10/bing: removed ladd keyword, use multsv_mn_add instead
!
      intent(in) :: a,b
      intent(out) :: c
!
      real, dimension (nx,3) :: b,c
      real, dimension (nx) :: a
      integer :: i
!
      do i=1,3
        c(:,i)=a*b(:,i)
      enddo
!
    endsubroutine multsv_mn
!***********************************************************************
    subroutine multsv_mn_add(a,b,c)
!
!  Vector multiplied with scalar, gives vector.
!
!  22-nov-01/nils erland: coded
!  10-oct-03/axel: a is now the scalar (now consistent with old routines)
!  24-jun-08/MR: ladd added for incremental work
!
      intent(in) :: a,b
      intent(inout) :: c
!
      real, dimension (nx,3) :: b,c
      real, dimension (nx) :: a
      integer :: i
!
      do i=1,3
        c(:,i)=c(:,i)+a*b(:,i)
      enddo
!
    endsubroutine multsv_mn_add
!***********************************************************************
    subroutine multsv_add_mn(a,b,c,d)
!
!  Multiply scalar with a vector and add to another vector.
!
!  29-oct-97/axel: coded
!
      real, dimension (nx,3) :: a,c,d
      real, dimension (nx) :: b
      integer :: j
!
      intent(in) :: a,b,c
      intent(out) :: d
!
      do j=1,3
        d(:,j)=a(:,j)+b*c(:,j)
      enddo
!
    endsubroutine multsv_add_mn
!***********************************************************************
    subroutine multvs_mn(a,b,c)
!
!  Vector pencil multiplied with scalar pencil, gives vector pencil.
!
!  22-nov-01/nils erland: coded
!
      real, dimension (nx,3) :: a, c
      real, dimension (nx) :: b
      integer :: i
!
      do i=1,3
        c(:,i)=a(:,i)*b(:)
      enddo
!
    endsubroutine multvs_mn
!***********************************************************************
    subroutine cross_mn(a,b,c)
!
!  Cross product, c = a x b, for pencil variables.
!  Previously called crossp.
!
      real, dimension (nx,3) :: a,b,c
!
      intent(in) :: a,b
      intent(out) :: c
!
      c(:,1)=a(:,2)*b(:,3)-a(:,3)*b(:,2)
      c(:,2)=a(:,3)*b(:,1)-a(:,1)*b(:,3)
      c(:,3)=a(:,1)*b(:,2)-a(:,2)*b(:,1)
!
    endsubroutine cross_mn
!***********************************************************************
    subroutine cross_mixed(a,b,c)
!
!  Cross product, c = a x b, of a pencil and a non-pencil variable.
!
!  17-apr-2015/MR: coded
!
      real, dimension (nx,3) :: a,c
      real, dimension (3) :: b
!
      intent(in) :: a,b
      intent(out) :: c
!
      c(:,1)=a(:,2)*b(3)-a(:,3)*b(2)
      c(:,2)=a(:,3)*b(1)-a(:,1)*b(3)
      c(:,3)=a(:,1)*b(2)-a(:,2)*b(1)
!
    endsubroutine cross_mixed
!***********************************************************************
    subroutine cross_0(a,b,c)
!
!  Cross product, c = a x b, for simple 3-d vectors (independent of position).
!
      real, dimension (3) :: a,b,c
!
      intent(in) :: a,b
      intent(out) :: c
!
      c(1)=a(2)*b(3)-a(3)*b(2)
      c(2)=a(3)*b(1)-a(1)*b(3)
      c(3)=a(1)*b(2)-a(2)*b(1)
!
    endsubroutine cross_0
!***********************************************************************
    subroutine gij(f,k,g,nder)
!
!  Calculate gradient of a vector, return matrix.
!
!   3-apr-01/axel+gitta: coded
!
      use Deriv, only: der,der2,der3,der4,der5,der6
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3) :: g
      real, dimension (nx) :: tmp
      integer :: i,j,k,k1,nder
!
      intent(in) :: f,k
      intent(out) :: g
!
      k1=k-1
      do i=1,3; do j=1,3
        if (nder == 1) then
          call der(f,k1+i,tmp,j)
        elseif (nder == 2) then
          call der2(f,k1+i,tmp,j)
        elseif (nder == 3) then
          call der3(f,k1+i,tmp,j)
        elseif (nder == 4) then
          call der4(f,k1+i,tmp,j)
        elseif (nder == 5) then
          call der5(f,k1+i,tmp,j)
        elseif (nder == 6) then
          call der6(f,k1+i,tmp,j)
        endif
        g(:,i,j)=tmp
      enddo; enddo
!
    endsubroutine gij
!***********************************************************************
    subroutine gijl_symmetric(f,k,gijl)
!
!  Calculate gradient of a (symmetric) second rank matrix, return 3rd rank
!  matrix
!
!  18-aug-08/dhruba: coded
!
      use Deriv, only: der
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3,3) :: gijl
      real, dimension (nx,3,3,3) :: tmpg
      real, dimension (nx) :: tmp
      integer :: i,j,l,k,k1
!
      intent(in) :: f,k
      intent(out) :: gijl
!
      k1=k-1
      do i=1,3
        do j=i,3
          k1=k1+1
          do l=1,3
            call der(f,k1,tmp,l)
            tmpg(:,i,j,l) = tmp
          enddo
        enddo
      enddo
      do l=1,3
        call symmetrise3x3_ut2lt(tmpg(:,:,:,l))
      enddo
      gijl=tmpg
!
    endsubroutine gijl_symmetric
!***********************************************************************
    subroutine grad_main(f,k,g)
!
!  Calculate gradient of a scalar, get vector.
!
!  29-sep-97/axel: coded
!
      use Deriv, only: der
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: g
      integer :: k
!
      intent(in) :: f,k
      intent(out) :: g
!
      call der(f,k,g(:,1),1)
      call der(f,k,g(:,2),2)
      call der(f,k,g(:,3),3)
!
    endsubroutine grad_main
!***********************************************************************
    subroutine grad_other(f,g)
!
!  For non 'mvar' variable calculate gradient of a scalar, get vector
!
!  26-nov-02/tony: coded
!
      use Deriv, only: der
!
      real, dimension (mx,my,mz) :: f
      real, dimension (nx,3) :: g
      real, dimension (nx) :: tmp
!
      intent(in) :: f
      intent(out) :: g
!
!  Uses overloaded der routine.
!
      call der(f,tmp,1); g(:,1)=tmp
      call der(f,tmp,2); g(:,2)=tmp
      call der(f,tmp,3); g(:,3)=tmp
!
    endsubroutine grad_other
!***********************************************************************
    subroutine grad5(f,k,g)
!
!  Calculate 5th order gradient of a scalar, get vector.
!
!  03-jun-07/anders: adapted
!
      use Deriv, only: der5
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: g
      real, dimension (nx) :: tmp
      integer :: k
!
      intent(in) :: f,k
      intent(out) :: g
!
      call der5(f,k,tmp,1); g(:,1)=tmp
      call der5(f,k,tmp,2); g(:,2)=tmp
      call der5(f,k,tmp,3); g(:,3)=tmp
!
    endsubroutine grad5
!***********************************************************************
    subroutine div(f,k,g,ldiff_fluxes)
!
!  Calculate divergence of vector, get scalar.
!
!  13-dec-01/nils: coded
!  16-jul-02/nils: adapted from pencil_mpi
!  31-aug-07/wlad: adapted for cylindrical and spherical coords
!  28-sep-15/Joern+MR: adapted to use for slope-limited diffusive 
!                      flux given on a staggered grid. 
!
      use Deriv, only: der
      use General, only: loptest, ranges_dimensional
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: g
      integer :: k
      logical, optional :: ldiff_fluxes
!
      intent(in)  :: f,k,ldiff_fluxes
      intent(out) :: g
!
      integer :: k1,i
      integer, dimension(dimensionality) :: jrange
      real, dimension(nx) :: tmp
!
      k1=k-1
!
      if (loptest(ldiff_fluxes)) then
        call ranges_dimensional(jrange)
        g=0.
        do i=1,dimensionality
          call der_4th_stag(f,k1+i,tmp,jrange(i))
          g=g+tmp
        enddo
      else
        call der(f,k1+1,tmp,1)
        g=tmp
        call der(f,k1+2,tmp,2)
        g=g+tmp
        call der(f,k1+3,tmp,3)
        g=g+tmp
      endif
!
      if (lspherical_coords) &
        g=g+2.*r1_mn*f(l1:l2,m,n,k1+1)+r1_mn*cotth(m)*f(l1:l2,m,n,k1+2)
!
      if (lcylindrical_coords) &
        g=g+rcyl_mn1*f(l1:l2,m,n,k1+1)
!
    endsubroutine div
!***********************************************************************
    subroutine der_2nd(f,k,df,j)

      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      real, dimension(nx), intent(out) :: df
      integer, intent(in) :: j, k
!
      if (j==1) then
        if (nxgrid/=1) then
          df=(f(l1+1:l2+1,m,n,k)-f(l1-1:l2-1,m,n,k))/(2.*dx) 
        else
          df=0.
          if (ip<=5) print*, 'der_2nd: Degenerate case in x-direction'
        endif
      elseif (j==2) then
        if (nygrid/=1) then
          df=(f(l1:l2,m+1,n,k)-f(l1:l2,m-1,n,k))/(2.*dy)
        else
          df=0.
          if (ip<=5) print*, 'der_2nd: Degenerate case in y-direction'
        endif
      elseif (j==3) then
        if (nzgrid/=1) then
          df=(f(l1:l2,m,n+1,k)-f(l1:l2,m,n-1,k))/(2.*dz)
        else
          df=0.
          if (ip<=5) print*, 'der_2nd: Degenerate case in z-direction'
        endif
      endif

    endsubroutine der_2nd
!***********************************************************************
    subroutine der_4th_stag(f,k,df,j)
!
!  Calculates 1st order derivative by a 4th order difference scheme from
!  data given on a grid shifted by half a grid step w.r.t. the looked at point.
!  Only valid for equidistant grids!
!
!  30-sep-15/MR: coded
!   4-feb-16/MR: checked again
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      real, dimension(nx), intent(out) :: df
      integer, intent(in) :: j, k
!
      if (j==1) then
        if (nxgrid/=1) then
          df=( -    (f(l1+1:l2+1,m,n,k)-f(l1-2:l2-2,m,n,k))     &
               +27.*(f(l1  :l2  ,m,n,k)-f(l1-1:l2-1,m,n,k)) )/(24.*dx)
        else
          df=0.
          if (ip<=5) print*, 'der_4th_stag: Degenerate case in x-direction'
        endif
      elseif (j==2) then
        if (nygrid/=1) then
          df=( -    (f(l1:l2,m+1,n,k)-f(l1:l2,m-2,n,k))    &
               +27.*(f(l1:l2,m  ,n,k)-f(l1:l2,m-1,n,k)) )/(24.*dy) 
        else
          df=0.
          if (ip<=5) print*, 'der_4th_stag: Degenerate case in y-direction'
        endif
      elseif (j==3) then
        if (nzgrid/=1) then
          df=( -    (f(l1:l2,m,n+1,k)-f(l1:l2,m,n-2,k))    &
               +27.*(f(l1:l2,m,n  ,k)-f(l1:l2,m,n-1,k)) )/(24.*dz)
        else
          df=0.
          if (ip<=5) print*, 'der_4th_stag: Degenerate case in z-direction'
        endif
      endif

    endsubroutine der_4th_stag
!***********************************************************************
    subroutine div_other(f,g)
!
      use Deriv, only: der
!
      real, dimension (mx,my,mz,3) :: f
      real, dimension (nx) :: g, tmp
!
      call der(f(:,:,:,1),tmp,1)
      g=tmp
      call der(f(:,:,:,2),tmp,2)
      g=g+tmp
      call der(f(:,:,:,3),tmp,3)
      g=g+tmp
!
      if (lspherical_coords) then
        g=g+2.*r1_mn*f(l1:l2,m,n,1)+r1_mn*cotth(m)*f(l1:l2,m,n,2)
      endif
!
      if (lcylindrical_coords) then
        g=g+rcyl_mn1*f(l1:l2,m,n,1)
      endif
!
    endsubroutine div_other
!***********************************************************************
    subroutine div_mn(aij,b,a)
!
!  Calculate divergence from derivative matrix.
!
!  18-sep-04/axel: coded
!  21-feb-07/axel: corrected spherical coordinates
!  14-mar-07/wlad: added cylindrical coordinates
!
      real, dimension (nx,3,3) :: aij
      real, dimension (nx,3) :: a
      real, dimension (nx) :: b
!
      intent(in) :: aij,a
      intent(out) :: b
!
      b=aij(:,1,1)+aij(:,2,2)+aij(:,3,3)
!
!  Adjustments for other coordinate systems.
!
      if (lspherical_coords) then
        b=b+2.*r1_mn*a(:,1)+r1_mn*cotth(m)*a(:,2)
      endif
!
      if (lcylindrical_coords) then
        b=b+rcyl_mn1*a(:,1)
      endif
!
      if (lpipe_coords) then
        b=b+glnCrossSec*a(:,1)
      endif
!
    endsubroutine div_mn
!***********************************************************************
    subroutine div_mn_2tensor(aijk,bi)
!
!  Calculate divergence from derivative matrix.
!
!  07-aug-10/dhruba: coded
!
      real, dimension (nx,3,3,3) :: aijk
      real, dimension (nx,3) :: bi
!
      intent(in) :: aijk
      intent(out) :: bi
      integer :: i
!
      do i=1,3
        bi(:,i)=aijk(:,i,1,1)+aijk(:,i,2,2)+aijk(:,i,3,3)
      enddo
!
!  Adjustments for spherical coordinate system.
!
      if (lspherical_coords) then
        call fatal_error('div_mn_2tensor','not impelmented in sph-coordinate')
      endif
!
      if (lcylindrical_coords) then
        call fatal_error('div_mn_2tensor','not impelmented in cyl-coordinate')
      endif
!
    endsubroutine div_mn_2tensor
!***********************************************************************
    subroutine curl_mn(aij,b,a,lcovariant_derivative)
!
!  Calculate curl from derivative matrix.
!
!  21-jul-03/axel: coded
!  21-feb-07/axel: corrected spherical coordinates
!  14-mar-07/wlad: added cylindrical coordinates
!  16-jun-16/fred: added option to use covariant derivative bij
!
      use General, only: loptest
!
      real, dimension (nx,3,3), intent (in) :: aij
      real, dimension (nx,3), intent (in), optional :: a
      logical, intent (in), optional :: lcovariant_derivative
      real, dimension (nx,3), intent (out) :: b
      integer :: i1=1,i2=2,i3=3,i4=4,i5=5,i6=6,i7=7
!
      b(:,1)=aij(:,3,2)-aij(:,2,3)
      b(:,2)=aij(:,1,3)-aij(:,3,1)
      b(:,3)=aij(:,2,1)-aij(:,1,2)
!
!  Adjustments for spherical coordinate system.
!
      if (lspherical_coords) then
        if (.not.loptest(lcovariant_derivative)) then
          if (.not. present(a)) then
            call fatal_error('curl_mn','Need a for spherical curl')
          endif
          b(:,1)=b(:,1)+a(:,3)*r1_mn*cotth(m)
          b(:,2)=b(:,2)-a(:,3)*r1_mn
          b(:,3)=b(:,3)+a(:,2)*r1_mn
        endif
      endif
!
!  Adjustments for cylindrical coordinate system.
!  If we go all the way to the center, we need to put a regularity condition.
!  We do this here currently up to second order, and only for curl_mn.
!
      if (lcylindrical_coords.and.present(a)) then
        if (.not.loptest(lcovariant_derivative)) b(:,3)=b(:,3)+a(:,2)*rcyl_mn1
        if (rcyl_mn(1)==0.) b(i1,3)=(360.*b(i2,3)-450.*b(i3,3)+400.*b(i4,3) &
                                    -225.*b(i5,3)+72.*b(i6,3)-10.*b(i7,3))/147.
      endif
!
    endsubroutine curl_mn
!***********************************************************************
    subroutine curl_horizontal(f,k,g)
!
!  Calculate curl of a vector, whose z component is given.
!
!   8-oct-09/axel: adapted from
!
      use Deriv, only: der
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: g
      real, dimension (nx) :: tmp1,tmp2
      integer :: k
!
      intent(in) :: f,k
      intent(out) :: g
!
      call der(f,k,tmp1,2)
      g(:,1)=tmp1
!
      call der(f,k,tmp2,1)
      g(:,2)=-tmp2
!
!g(:,1)=0.
!g(:,2)=0.
      g(:,3)=0.
!
!  Adjustments for spherical corrdinate system.
!
      if (lspherical_coords) then
        g(:,1)=g(:,1)+f(l1:l2,m,n,k)*r1_mn*cotth(m)
        g(:,2)=g(:,2)-f(l1:l2,m,n,k)*r1_mn
      endif
!
!     if (lcylindrical_coords) then
!--     g(:,3)=g(:,3)+f(l1:l2,m,n,k1+2)*rcyl_mn1
!     endif
!
    endsubroutine curl_horizontal
!***********************************************************************
    subroutine curl(f, k, g, ignoredx)
!
!  Calculate curl of a vector, get vector.
!
!  12-sep-97/axel: coded
!  10-sep-01/axel: adapted for cache efficiency
!  11-sep-04/axel: began adding spherical coordinates
!  21-feb-07/axel: corrected spherical coordinates
!  14-mar-07/wlad: added cylindrical coordinates
!  20-sep-13/ccyang: added optional argument ignoredx
!
      use Deriv, only: der
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      real, dimension(nx,3), intent(out) :: g
      integer, intent(in) :: k
      logical, intent(in), optional :: ignoredx
!
      real, dimension(nx) :: tmp1, tmp2
      logical :: igdx
      integer :: k1
!
      if (present(ignoredx)) then
        igdx = ignoredx
      else
        igdx = .false.
      endif
!
      k1=k-1
!
      call der(f, k1+3, tmp1, 2, ignoredx=igdx)
      call der(f, k1+2, tmp2, 3, ignoredx=igdx)
      g(:,1)=tmp1-tmp2
!
      call der(f, k1+1, tmp1, 3, ignoredx=igdx)
      call der(f, k1+3, tmp2, 1, ignoredx=igdx)
      g(:,2)=tmp1-tmp2
!
      call der(f, k1+2, tmp1, 1, ignoredx=igdx)
      call der(f, k1+1, tmp2, 2, ignoredx=igdx)
      g(:,3)=tmp1-tmp2
!
!  Adjustments for spherical coordinate system.
!
      if (.not. igdx .and. lspherical_coords) then
        g(:,1)=g(:,1)+f(l1:l2,m,n,k1+3)*r1_mn*cotth(m)
        g(:,2)=g(:,2)-f(l1:l2,m,n,k1+3)*r1_mn
        g(:,3)=g(:,3)+f(l1:l2,m,n,k1+2)*r1_mn
      endif
!
      if (.not. igdx .and. lcylindrical_coords) then
        g(:,3)=g(:,3)+f(l1:l2,m,n,k1+2)*rcyl_mn1
      endif
!
    endsubroutine curl
!***********************************************************************
    subroutine curl_other(f,g)
!
!  Calculate curl of a non-mvar vector, get vector.
!
!  23-june-09/wlad: adapted from curl
!
      use Deriv, only: der
!
      real, dimension (mx,my,mz,3) :: f
      real, dimension (nx,3) :: g
      real, dimension (nx) :: tmp1,tmp2
!
      intent(in) :: f
      intent(out) :: g
!
      call der(f(:,:,:,3),tmp1,2)
      call der(f(:,:,:,2),tmp2,3)
      g(:,1)=tmp1-tmp2
!
      call der(f(:,:,:,1),tmp1,3)
      call der(f(:,:,:,3),tmp2,1)
      g(:,2)=tmp1-tmp2
!
      call der(f(:,:,:,2),tmp1,1)
      call der(f(:,:,:,1),tmp2,2)
      g(:,3)=tmp1-tmp2
!
!  Adjustments for spherical corrdinate system.
!
      if (lspherical_coords) then
        g(:,1)=g(:,1)+f(l1:l2,m,n,3)*r1_mn*cotth(m)
        g(:,2)=g(:,2)-f(l1:l2,m,n,3)*r1_mn
        g(:,3)=g(:,3)+f(l1:l2,m,n,2)*r1_mn
      endif
!
      if (lcylindrical_coords) then
        g(:,3)=g(:,3)+f(l1:l2,m,n,2)*rcyl_mn1
      endif
!
    endsubroutine curl_other
!***********************************************************************
    subroutine curli(f,k,g,i)
!
!  Calculate curl of a vector, get vector.
!
!  22-oct-02/axel+tarek: adapted from curl
!
      use Deriv, only: der
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: g
      real, dimension (nx) :: tmp1,tmp2
      integer :: k,k1,i
!
      intent(in) :: f,k,i
      intent(out) :: g
!
      k1=k-1
!
      select case (i)
!
!  1-component
!
      case (1)
        call der(f,k1+3,tmp1,2)
        call der(f,k1+2,tmp2,3)
!
!  corrections for spherical coordinates
!
        if (lspherical_coords) then
          g=tmp1-tmp2+f(l1:l2,m,n,k1+3)*r1_mn*cotth(m)
        else
          g=tmp1-tmp2
        endif
!
!  2-component
!
      case (2)
        call der(f,k1+1,tmp1,3)
        call der(f,k1+3,tmp2,1)
!
!  corrections for spherical coordinates
!
        if (lspherical_coords) then
          g=tmp1-tmp2-f(l1:l2,m,n,k1+3)*r1_mn
        else
          g=tmp1-tmp2
        endif
!
!  3-component
!
      case (3)
        call der(f,k1+2,tmp1,1)
        call der(f,k1+1,tmp2,2)
!
!  corrections for spherical and cylindrical coordinates
!
        if (lspherical_coords) then
          g=tmp1-tmp2+f(l1:l2,m,n,k1+2)*r1_mn
        elseif (lcylindrical_coords) then
          g=tmp1-tmp2+f(l1:l2,m,n,k1+2)*rcyl_mn1
        else
          g=tmp1-tmp2
        endif
!
      endselect
!
    endsubroutine curli
!***********************************************************************
    subroutine del2_main(f,k,del2f)
!
!  Calculate del2 of a scalar, get scalar.
!
!  12-sep-97/axel: coded
!   7-mar-07/wlad: added cylindrical coordinates
!
      use Deriv, only: der,der2
!
      intent(in) :: f,k
      intent(out) :: del2f
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: del2f,d2fdx,d2fdy,d2fdz,tmp
      integer :: k
!
      call der2(f,k,d2fdx,1)
      call der2(f,k,d2fdy,2)
      call der2(f,k,d2fdz,3)
      del2f=d2fdx+d2fdy+d2fdz
!
      if (lcylindrical_coords) then
        call der(f,k,tmp,1)
        del2f=del2f+tmp*rcyl_mn1
      endif
!
      if (lspherical_coords) then
        call der(f,k,tmp,1)
        del2f=del2f+2.*r1_mn*tmp
        call der(f,k,tmp,2)
        del2f=del2f+cotth(m)*r1_mn*tmp
      endif
!
    endsubroutine del2_main
!***********************************************************************
    subroutine del2_other(f,del2f)
!
!  Calculate del2 of a scalar, get scalar.
!   8-may-09/nils: adapted from del2
!
      use Deriv, only: der,der2
!
      intent(in) :: f
      intent(out) :: del2f
!
      real, dimension (mx,my,mz) :: f
      real, dimension (nx) :: del2f,d2fdx,d2fdy,d2fdz,tmp
!
      call der2(f,d2fdx,1)
      call der2(f,d2fdy,2)
      call der2(f,d2fdz,3)
      del2f=d2fdx+d2fdy+d2fdz
!
      if (lcylindrical_coords) then
        call der(f,tmp,1)
        del2f=del2f+tmp*rcyl_mn1
      endif
!
      if (lspherical_coords) then
        call der(f,tmp,1)
        del2f=del2f+2.*r1_mn*tmp
        call der(f,tmp,2)
        del2f=del2f+cotth(m)*r1_mn*tmp
      endif
!
    endsubroutine del2_other
!***********************************************************************
    subroutine del2v(f,k,del2f,fij,pff)
!
!  Calculate del2 of a vector, get vector.
!
!  28-oct-97/axel: coded
!  15-mar-07/wlad: added cylindrical coordinates
!
      use Deriv, only: der
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, optional, dimension(nx,3,3) :: fij
      real, optional, dimension(nx,3) :: pff
      real, dimension (nx,3) :: del2f
      real, dimension (nx) :: tmp
      integer :: i,k,k1
!
      intent(in) :: f,k
      intent(out) :: del2f
!
!  do the del2 diffusion operator
!
      k1=k-1
      do i=1,3
        call del2(f,k1+i,tmp)
        del2f(:,i)=tmp
      enddo
!
      if (lcylindrical_coords) then
        !del2 already contains the extra term 1/r*d(uk)/dt
        call der(f,k1+2,tmp,2)
        del2f(:,1)=del2f(:,1) -(2*tmp+f(l1:l2,m,n,k1+1))*rcyl_mn2
        call der(f,k1+1,tmp,2)
        del2f(:,2)=del2f(:,2) +(2*tmp-f(l1:l2,m,n,k1+2))*rcyl_mn2
      endif
!
      if (lspherical_coords) then
        if (.not. (present(fij) .and. present(pff))) then
           call fatal_error('del2v', &
               'Cannot do a spherical del2v without aij and aa')
        endif
!
! for r component (factors of line elements are taken care of inside p%uij
        del2f(:,1)= del2f(:,1)+&
            r1_mn*(2.*(fij(:,1,1)-fij(:,2,2)-fij(:,3,3) &
            -r1_mn*pff(:,1)-cotth(m)*r1_mn*pff(:,2) ) &
            +cotth(m)*fij(:,1,2) )
! for theta component
        del2f(:,2)=del2f(:,2)+&
            r1_mn*(2.*(fij(:,2,1)-cotth(m)*fij(:,3,3)&
            +fij(:,1,2) )&
            +cotth(m)*fij(:,2,2)-r1_mn*sin1th(m)*sin1th(m)*pff(:,2) )
! for phi component
        del2f(:,3)=del2f(:,3)+&
            r1_mn*(2.*(fij(:,3,1)+fij(:,1,3)&
            +cotth(m)*fij(:,2,3) ) &
            +cotth(m)*fij(:,3,2)-sin1th(m)*pff(:,3) )
      endif
!
    endsubroutine del2v
!***********************************************************************
    subroutine del2m3x3_sym(f,k,del2f)
!
!  Calculate del2 of a 3x3 symmetric matrix, get matrix
!  23-feb-11/dhruba: coded in a new manner
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3) :: del2f
      real, dimension (nx) :: tmp
      integer :: i,j,k,k1
!
      intent(in) :: f,k
      intent(out) :: del2f
!
!  simple del2 diffusion operator for each component
!
      k1=k-1
      do i=1,3
        do j=i,3
          k1=k1+1
          call del2(f,k1,tmp)
          del2f(:,i,j) = tmp
        enddo
      enddo
      call symmetrise3x3_ut2lt(del2f)
!
      if (lcylindrical_coords) then
        call fatal_error('del2m', &
            'del2m is not implemented in cylindrical coordiates')
      endif
!
      if (lspherical_coords) then
        call fatal_error('del2m', &
            'del2m is not implemented in spherical coordiates')
      endif
!
    endsubroutine del2m3x3_sym
!***********************************************************************
    subroutine del2fj(f,vec,k,del2f)
!
!  Calculate del2 of a scalar, get scalar, adding anisotropic factor.
!
!  20-dec-12/wlad: adapted from del2_main
!
      use Deriv, only: der,der2
!
      intent(in) :: f,k
      intent(out) :: del2f
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: del2f,d2fdx,d2fdy,d2fdz,tmp
      real, dimension (3) :: vec
      integer :: k
!
      call der2(f,k,d2fdx,1)
      call der2(f,k,d2fdy,2)
      call der2(f,k,d2fdz,3)
      del2f=vec(1)*d2fdx+vec(2)*d2fdy+vec(3)*d2fdz
!
      if (lcylindrical_coords.and.vec(1)/=0.) then
        call der(f,k,tmp,1)
        del2f=del2f+vec(1)*tmp*rcyl_mn1
      endif
!
      if (lspherical_coords) then
        if (vec(1)/=0.) then
          call der(f,k,tmp,1)
          del2f=del2f+vec(1)*2.*r1_mn*tmp
        endif
        if (vec(2)/=0.) then
          call der(f,k,tmp,2)
          del2f=del2f+vec(2)*cotth(m)*r1_mn*tmp
        endif
      endif
!
    endsubroutine del2fj
!***********************************************************************
    subroutine del2fi_dxjk(f,k,del2fkdxij)
!
!  Calculate \partial^2f/\partial x_j\partial x_k of a vector, get a 9 dimensional object
!
!
      use Deriv, only: der2,derij
!
      intent(in) :: f,k
      intent(out) :: del2fkdxij
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3,3) :: del2fkdxij
      real, dimension(nx) :: tmp
      integer :: k
      integer :: i,j,kincrement
!
!
      if (.not. lcartesian_coords) &
        call fatal_error('sub','del2fi_dxjk implemented only in cartesian')
!
      do kincrement=0,2
       do i=1,3; do j=i,3
         tmp=0.
         if (i.eq.j) then
           call der2(f,k+kincrement,tmp,i)
           del2fkdxij(:,k+kincrement,i,i)=tmp
         else
           call derij(f,k+kincrement,tmp,i,j)
           del2fkdxij(:,k+kincrement,i,j)=tmp
           del2fkdxij(:,k+kincrement,j,i)=tmp
         endif
       enddo;enddo
     enddo
!
   endsubroutine del2fi_dxjk
!***********************************************************************
    subroutine d2fi_dxj(f,k,d2fidxj)
!
!  Calculate d^2f_i/dx^2_j of a vector, get a six dimensional object
!
      use Deriv, only: derij
!
      intent(in) :: f,k
      intent(out) :: d2fidxj
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3) :: d2fidxj
      real, dimension(nx,3) :: tmp
      integer :: k
!
      if (.not. lcartesian_coords) &
        call fatal_error('sub','d2fidxj implemented only in cartesian')
!
      call d2f_dxj(f,k,tmp)
      d2fidxj(:,1,:) = tmp
      call d2f_dxj(f,k+1,tmp)
      d2fidxj(:,2,:) = tmp
      call d2f_dxj(f,k+2,tmp)
      d2fidxj(:,3,:) = tmp
!
    endsubroutine d2fi_dxj
!***********************************************************************
    subroutine d2f_dxj(f,k,d2fdxj)
!
!  Calculate d^2f/dx^2_j of a scalar, get a three dimensional object
!
      use Deriv, only: der,der2
!
      intent(in) :: f,k
      intent(out) :: d2fdxj
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: d2fdxj
      integer :: k
!
      if (.not. lcartesian_coords) &
        call fatal_error('sub','d2f_dxj implemented only in cartesian')
!
      call der2(f,k,d2fdxj(:,1),1)
      call der2(f,k,d2fdxj(:,2),2)
      call der2(f,k,d2fdxj(:,3),3)
!
    endsubroutine d2f_dxj
!***********************************************************************
    subroutine symmetrise3x3_ut2lt (matrix_ut3x3)
!
! sets the lower triangular values of a matrix to its upper
! triangular values. Does not touch the diagonal. Applies
! to 3x3 matrices (pencil) only.
!
!  23-dhruba-11/dhruba: coded
!
      real, dimension(nx,3,3) :: matrix_ut3x3
      integer :: i,j
!
      do i=1,3
        do j=1,i
          if (i/=j) &
            matrix_ut3x3(:,i,j) = matrix_ut3x3(:,j,i)
        enddo
      enddo
!
    endsubroutine symmetrise3x3_ut2lt
!***********************************************************************
    subroutine del2v_etc(f,k,del2,graddiv,curlcurl,gradcurl)
!
!  Calculates a number of second derivative expressions of a vector
!  outputs a number of different vector fields.
!  gradcurl is not the vector gradient.
!  Surprisingly, calling derij only if graddiv or curlcurl are present
!  does not speed up the code on Mephisto @ 32x32x64.
!
!  12-sep-01/axel: coded
!  15-mar-07/wlad: added cylindrical coordinates
!
      use Deriv, only: der,der2,derij
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3) :: fjji,fijj
      real, dimension (nx,3,3), optional :: gradcurl
      real, dimension (nx,3), optional :: del2,graddiv,curlcurl
      real, dimension (nx,3) ::  fjik
      real, dimension (nx) :: tmp
      integer :: i,j,k,k1
!
      intent(in) :: f,k
      intent(out) :: del2,graddiv,curlcurl,gradcurl
!
!  calculate f_{i,jj} and f_{j,ji}
!  AJ: graddiv needs diagonal elements from the first tmp (derij only sets
!      off-diagonal elements)
!
      k1=k-1
      do i=1,3
      do j=1,3
        if (present(del2) .or. present(curlcurl) .or. present(gradcurl) .or. &
            present(graddiv)) then
          call der2 (f,k1+i,tmp,  j); fijj(:,i,j)=tmp  ! f_{i,jj}
        endif
        if (present(graddiv) .or. present(curlcurl).or. present(gradcurl)) then
          call derij(f,k1+j,tmp,j,i); fjji(:,i,j)=tmp  ! f_{j,ji}
        endif
      enddo
      enddo
!
!  the diagonal terms have not been set in derij; do this now
!  ** They are automatically set above, because derij   **
!  ** doesn't overwrite the value of tmp for i=j!       **
!
!     do j=1,3
!       fjji(:,j,j)=fijj(:,j,j)
!     enddo
!
!  calculate f_{i,jk} for i /= j /= k
!
      if (present(gradcurl)) then
        call derij(f,k1+1,tmp,2,3)
        fjik(:,1)=tmp
        call derij(f,k1+2,tmp,1,3)
        fjik(:,2)=tmp
        call derij(f,k1+3,tmp,1,2)
        fjik(:,3)=tmp
        if (lcylindrical_coords.or.lspherical_coords) &
          call fatal_error('del2v_etc','gradcurl at del2v_etc not '// &
              'implemented for non-cartesian coordinates')
      endif
!
!  del2
!
      if (present(del2)) then
        do i=1,3
          del2(:,i)=fijj(:,i,1)+fijj(:,i,2)+fijj(:,i,3)
        enddo
        if (lcylindrical_coords) then
          !r-component
          call der(f,k1+2,tmp,2)
          del2(:,1)=del2(:,1) -(2*tmp+f(l1:l2,m,n,k1+1))*rcyl_mn2
          call der(f,k1+1,tmp,1)
          del2(:,1)=del2(:,1) + tmp*rcyl_mn1
          !phi-component
          call der(f,k1+1,tmp,2)
          del2(:,2)=del2(:,2) +(2*tmp-f(l1:l2,m,n,k1+2))*rcyl_mn2
          call der(f,k1+2,tmp,1)
          del2(:,2)=del2(:,2) + tmp*rcyl_mn1
          !z-component
          call der(f,k1+3,tmp,1)
          del2(:,3)=del2(:,3) + tmp*rcyl_mn1
        endif
        if (lspherical_coords) call fatal_error('del2v_etc', &
            'del2 at del2v_etc not implemented for spherical '// &
            'coordinates - use gij_etc')
      endif
!
      if (present(graddiv)) then
        do i=1,3
          graddiv(:,i)=fjji(:,i,1)+fjji(:,i,2)+fjji(:,i,3)
        enddo
        if (lcylindrical_coords) then
           call der(f,k1+1,tmp,1)
           graddiv(:,1)=graddiv(:,1)+tmp*rcyl_mn1 - f(l1:l2,m,n,k1+1)*rcyl_mn2
           call der(f,k1+1,tmp,2)
           graddiv(:,2)=graddiv(:,2)+tmp*rcyl_mn1
           call der(f,k1+1,tmp,3)
           graddiv(:,3)=graddiv(:,3)+tmp*rcyl_mn1
        endif
        if (lspherical_coords) call fatal_error('del2v_etc', &
            'graddiv is implemented in gij_etc for spherical '// &
            'coords - use gij_etc')
      endif
!
!  curlcurl
!
      if (present(curlcurl)) then
        curlcurl(:,1)=fjji(:,1,2)-fijj(:,1,2)+fjji(:,1,3)-fijj(:,1,3)
        curlcurl(:,2)=fjji(:,2,3)-fijj(:,2,3)+fjji(:,2,1)-fijj(:,2,1)
        curlcurl(:,3)=fjji(:,3,1)-fijj(:,3,1)+fjji(:,3,2)-fijj(:,3,2)
        if (lcylindrical_coords) call fatal_error('del2v_etc', &
            'curlcurl not implemented for non-cartesian coordinates')
        if (lspherical_coords) call fatal_error('del2v_etc', &
            'curlcurl not implemented for non-cartesian coordinates - '// &
            'use gij_etc')
      endif
!
!  gradcurl (as tensor)
!
      if (present(gradcurl)) then
         gradcurl(:,1,1) = fjik(:,3)   - fjik(:,2)
         gradcurl(:,1,2) = fjji(:,1,3) - fijj(:,3,1)
         gradcurl(:,1,3) = fijj(:,2,1) - fjji(:,1,2)
!
         gradcurl(:,2,1) = fijj(:,3,2) - fjji(:,2,3)
         gradcurl(:,2,2) = fjik(:,1)   - fjik(:,3)
         gradcurl(:,2,3) = fjji(:,2,1) - fijj(:,1,2)
!
         gradcurl(:,3,1) = fjji(:,3,2) - fijj(:,2,3)
         gradcurl(:,3,2) = fijj(:,1,3) - fjji(:,3,1)
         gradcurl(:,3,3) = fjik(:,2)   - fjik(:,1)
!
         if (lcylindrical_coords.or.lspherical_coords) &
             call fatal_error('del2v_etc','gradcurl not implemented '// &
             'for non-cartesian coordinates')
      endif
!
    endsubroutine del2v_etc
!***********************************************************************
    subroutine del2vi_etc(f,k,ii,del2,graddiv,curlcurl)
!
!  Calculates a number of second derivative expressions of a vector.
!  Outputs a number of different vector fields.
!  Surprisingly, calling derij only if graddiv or curlcurl are present
!  does not speed up the code on Mephisto @ 32x32x64.
!  Just do the ith component
!
!   7-feb-04/axel: adapted from del2v_etc
!
      use Deriv, only: der2,derij
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3) :: fjji,fijj
      real, dimension (nx), optional :: del2,graddiv,curlcurl
      real, dimension (nx) :: tmp
      integer :: i,j,k,k1,ii
!
      intent(in) :: f,k,ii
      intent(out) :: del2,graddiv,curlcurl
!
!  Do the del2 diffusion operator.
!
      k1=k-1
      do i=1,3
      do j=1,3
        call der2 (f,k1+i,tmp,  j); fijj(:,i,j)=tmp  ! f_{i,jj}
        call derij(f,k1+j,tmp,j,i); fjji(:,i,j)=tmp  ! f_{j,ji}
      enddo
      enddo
!
      if (present(del2)) then
        del2=fijj(:,ii,1)+fijj(:,ii,2)+fijj(:,ii,3)
      endif
!
      if (present(graddiv)) then
        graddiv=fjji(:,ii,1)+fjji(:,ii,2)+fjji(:,ii,3)
      endif
!
      if (present(curlcurl)) then
        select case (ii)
        case (1); curlcurl=fjji(:,1,2)-fijj(:,1,2)+fjji(:,1,3)-fijj(:,1,3)
        case (2); curlcurl=fjji(:,2,3)-fijj(:,2,3)+fjji(:,2,1)-fijj(:,2,1)
        case (3); curlcurl=fjji(:,3,1)-fijj(:,3,1)+fjji(:,3,2)-fijj(:,3,2)
        endselect
      endif
!
      if (lcylindrical_coords.or.lspherical_coords) &
          call fatal_error('del2vi_etc', &
          'not implemented for non-cartesian coordinates')
!
    endsubroutine del2vi_etc
!***********************************************************************
    subroutine del4v(f,k,del4f)
!
!  Calculate del4 of a vector, get vector.
!
!  09-dec-03/nils: adapted from del6v
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: del4f
      real, dimension (nx) :: tmp
      integer :: i,k,k1
!
      intent(in) :: f,k
      intent(out) :: del4f
!
!  Exit if this is requested for non-cartesian runs.
!
      if (lcylindrical_coords.or.lspherical_coords) &
          call fatal_error('del4v', &
          'not implemented for non-cartesian coordinates')
!
!  Do the del4 diffusion operator.
!
      k1=k-1
      do i=1,3
        call del4(f,k1+i,tmp)
        del4f(:,i)=tmp
      enddo
!
    endsubroutine del4v
!***********************************************************************
    subroutine del6v(f,k,del6f)
!
!  Calculate del6 of a vector, get vector.
!
!  28-oct-97/axel: coded
!  24-apr-03/nils: adapted from del2v
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: del6f
      real, dimension (nx) :: tmp
      integer :: i,k,k1
!
      intent(in) :: f,k
      intent(out) :: del6f
!
!  Do the del6 diffusion operator.
!
      k1=k-1
      do i=1,3
        call del6(f,k1+i,tmp)
        del6f(:,i)=tmp
      enddo
!
!  Exit if this is requested for non-cartesian runs.
!
      !!!!if (lcylindrical_coords.or.lspherical_coords) &
      !!!!    call fatal_error('del6v', &
      !!!!    'not implemented for non-cartesian coordinates')
!
    endsubroutine del6v
!***********************************************************************
    subroutine gij_etc(f,iref,aa,aij,bij,del2,graddiv,lcovariant_derivative)
!
!  Calculate B_i,j = eps_ikl A_l,jk and A_l,kk.
!
!  21-jul-03/axel: coded
!  26-jul-05/tobi: do not calculate both d^2 A/(dx dy) and d^2 A/(dy dx)
!  23-feb-07/axel: added spherical coordinates
!   7-mar-07/wlad: added cylindrical coordinates
!  29-aug-13/MR: made bij optional; added error messages for missing optional parameters
!  15-jun-16/fred: bij optional covariant derivative terms; spherical/cylindrical
!
      use Deriv, only: der2,derij
      use General, only: loptest
!
      real, dimension (mx,my,mz,mfarray), intent (in) :: f
      integer, intent (in) :: iref
      logical, intent (in), optional :: lcovariant_derivative
      real, dimension (nx,3), intent (in), optional :: aa
      real, dimension (nx,3,3), intent (in), optional :: aij
      real, dimension (nx,3,3), intent (out), optional :: bij
      real, dimension (nx,3), intent (out), optional :: del2,graddiv
!
!  Locally used variables.
!
      real, dimension (nx,3,3,3) :: d2A
      real, dimension (nx) :: tmp
      integer :: iref1,i,j
!
!  Reference point of argument.
!
      iref1=iref-1
!
!  Calculate all mixed and non-mixed second derivatives
!  of the vector potential (A_k,ij).
!
!  Do not calculate both d^2 A/(dx dy) and d^2 A/(dy dx).
!  (This wasn't spotted by me but by a guy from SGI...)
!  Note: for non-cartesian coordinates there are different correction terms,
!  see below.
!
      do i=1,3
        do j=1,3
          call der2(f,iref1+i,tmp,j); d2A(:,j,j,i)=tmp
        enddo
        call derij(f,iref1+i,tmp,2,3); d2A(:,2,3,i)=tmp; d2A(:,3,2,i)=tmp
        call derij(f,iref1+i,tmp,3,1); d2A(:,3,1,i)=tmp; d2A(:,1,3,i)=tmp
        call derij(f,iref1+i,tmp,1,2); d2A(:,1,2,i)=tmp; d2A(:,2,1,i)=tmp
      enddo
!
!  Corrections for spherical polars from swapping mixed derivatives:
!  Psi_{,theta^ r^} = Psi_{,r^ theta^} - Psi_{,\theta^}/r
!  Psi_{,phi^ r^} = Psi_{,r^ phi^} - Psi_{,\phi^}/r
!  Psi_{,phi^ theta^} = Psi_{,theta^ phi^} - Psi_{,\phi^}*r^{-1}*cot(theta)
!
      if (lspherical_coords) then
        if (.not.present(aij)) &
          call fatal_error('gij_etc', 'aij needed for spherical coordinates')
        do i=1,3
          d2A(:,2,1,i)=d2A(:,2,1,i)-aij(:,i,2)*r1_mn
          d2A(:,3,1,i)=d2A(:,3,1,i)-aij(:,i,3)*r1_mn
          d2A(:,3,2,i)=d2A(:,3,2,i)-aij(:,i,3)*r1_mn*cotth(m)
        enddo
      endif
!
!  For cylindrical, only
!  Psi_{,phi^ pom^} = Psi_{,pom^ phi^} - Psi_{,\phi^}/pom .
!
      if (lcylindrical_coords) then
        if (.not.present(aij)) &
          call fatal_error('gij_etc', 'aij needed for cylindrical coordinates')
        do i=1,3
          d2A(:,2,1,i)=d2A(:,2,1,i)-aij(:,i,2)*rcyl_mn1
        enddo
      endif
!
!  Calculate optionally b_i,j = eps_ikl A_l,kj,
!  del2_i = A_i,jj and graddiv_i = A_j,ji .
!
      if (present(bij)) then
!
        bij(:,1,:)=d2A(:,2,:,3)-d2A(:,3,:,2)
        bij(:,2,:)=d2A(:,3,:,1)-d2A(:,1,:,3)
        bij(:,3,:)=d2A(:,1,:,2)-d2A(:,2,:,1)
!
!  Corrections for spherical coordinates.
!
        if (lspherical_coords) then
          if (.not.present(aa)) &
            call fatal_error('gij_etc', 'aa needed for spherical coordinates')
          bij(:,3,2)=bij(:,3,2)+aij(:,2,2)*r1_mn
          bij(:,2,3)=bij(:,2,3)-aij(:,3,3)*r1_mn
          bij(:,1,3)=bij(:,1,3)+aij(:,3,3)*r1_mn*cotth(m)
          bij(:,3,1)=bij(:,3,1)+aij(:,2,1)*r1_mn         -aa(:,2)*r2_mn
          bij(:,2,1)=bij(:,2,1)-aij(:,3,1)*r1_mn         +aa(:,3)*r2_mn
          bij(:,1,2)=bij(:,1,2)+aij(:,3,2)*r1_mn*cotth(m)-aa(:,3)*r2_mn*sin2th(m)
          if (loptest(lcovariant_derivative)) then
            bij(:,1,1)=bij(:,1,1)+(aij(:,3,2)*r1_mn-aa(:,3)*r2_mn)*cotth(m)
            bij(:,1,2)=bij(:,1,2)+(aij(:,3,1)-aij(:,1,3))*r1_mn+aa(:,3)*r2_mn
            bij(:,1,3)=bij(:,1,3)+(aij(:,1,2)-aij(:,2,1))*r1_mn-aa(:,2)*r2_mn
            !bij(:,2,1)=bij(:,2,1)
            bij(:,2,2)=bij(:,2,2)-aij(:,2,3)*r1_mn         +aa(:,3)*r2_mn*cotth(m)
            bij(:,2,3)=bij(:,2,3)+(aij(:,1,2)*r1_mn-aij(:,2,1)*r1_mn-aa(:,2)*r2_mn)&
                                                                         *cotth(m)
            !bij(:,3,1)=bij(:,3,1)
            !bij(:,3,2)=bij(:,3,2)
            bij(:,3,3)=bij(:,3,3)+(aij(:,3,2)+(aij(:,1,3)-aij(:,3,1))*cotth(m))*r1_mn
          endif
        endif
!
!  Corrections for cylindrical coordinates.
!
        if (lcylindrical_coords) then
          if (.not.present(aa)) &
            call fatal_error('gij_etc', 'aa needed for cylindrical coordinates')
          bij(:,3,2)=bij(:,3,2)+ aij(:,2,2)*rcyl_mn1
!          !bij(:,3,1)=bij(:,3,1)+(aij(:,2,1)+aij(:,1,2))*rcyl_mn1-aa(:,2)*rcyl_mn2
!  FAG:Partial correction to -d2A(:,2,1,1) already applied above +aij(:,i,2)*rcyl_mn1
!
          bij(:,3,1)=bij(:,3,1)+aij(:,2,1)*rcyl_mn1-aa(:,2)*rcyl_mn2
          if (loptest(lcovariant_derivative)) then
            !bij(:,1,1)=bij(:,1,1)
            bij(:,1,2)=bij(:,1,2)+(aij(:,3,1)-aij(:,1,3))*rcyl_mn1
            !bij(:,1,3)=bij(:,1,3)
            !bij(:,2,1)=bij(:,2,1)
            bij(:,2,2)=bij(:,2,2)+(aij(:,3,2)-aij(:,2,3))*rcyl_mn1
            !bij(:,2,3)=bij(:,2,3)
            !bij(:,3,1)=bij(:,3,1)
            !bij(:,3,2)=bij(:,3,2)
            bij(:,3,3)=bij(:,3,3)+aij(:,2,3)*rcyl_mn1
          endif
        endif
      endif
!
!  Calculate del2 and graddiv, if requested.
!
      if (present(graddiv)) then
        graddiv(:,:)=d2A(:,1,:,1)+d2A(:,2,:,2)+d2A(:,3,:,3)
        if (lspherical_coords) then
          if (.not.present(aa)) &
            call fatal_error('gij_etc', 'aa needed for spherical coordinates')
          graddiv(:,1)=graddiv(:,1)+aij(:,1,1)*r1_mn*2+ &
             aij(:,2,1)*r1_mn*cotth(m)-aa(:,2)*r2_mn*cotth(m)-aa(:,1)*r2_mn*2
          graddiv(:,2)=graddiv(:,2)+aij(:,1,2)*r1_mn*2+ &
             aij(:,2,2)*r1_mn*cotth(m)-aa(:,2)*r2_mn*sin2th(m)
          graddiv(:,3)=graddiv(:,3)+aij(:,1,3)*r1_mn*2+ &
             aij(:,2,3)*r1_mn*cotth(m)
        endif
      endif
!
      if (present(del2)) then
        del2(:,:)=d2A(:,1,1,:)+d2A(:,2,2,:)+d2A(:,3,3,:)
        if (lspherical_coords) then
          if (.not.present(aa)) &
            call fatal_error('gij_etc', 'aa needed for spherical coordinates')
          del2(:,1)= del2(:,1)+&
            r1_mn*(2.*(aij(:,1,1)-aij(:,2,2)-aij(:,3,3)&
            -r1_mn*aa(:,1)-cotth(m)*r1_mn*aa(:,2) ) &
            +cotth(m)*aij(:,1,2) )
          del2(:,2)=del2(:,2)+&
            r1_mn*(2.*(aij(:,2,1)-cotth(m)*aij(:,3,3)&
            +aij(:,1,2) )&
            +cotth(m)*aij(:,2,2)-r1_mn*sin2th(m)*aa(:,2) )
          del2(:,3)=del2(:,3)+&
            r1_mn*(2.*(aij(:,3,1)+aij(:,1,3)&
            +cotth(m)*aij(:,2,3) ) &
            +cotth(m)*aij(:,3,2)-r1_mn*sin2th(m)*aa(:,3) )
        endif
        if (lcylindrical_coords) then
          if (.not.present(aa)) &
            call fatal_error('gij_etc', 'aa needed for cylindrical coordinates')
          del2(:,1)= del2(:,1)+&
            rcyl_mn1*(aij(:,1,1)-2*aij(:,2,2))-rcyl_mn2*aa(:,1)
          del2(:,2)=del2(:,2)+&               
            rcyl_mn1*(aij(:,2,1)-2*aij(:,1,2))-rcyl_mn2*aa(:,2)
        endif
      endif
!
    endsubroutine gij_etc
!***********************************************************************
    subroutine g2ij(f,k,g)
!
!  Calculates the Hessian, i.e. all second derivatives of a scalar.
!
!  11-jul-02/axel: coded
!
      use Deriv, only: der2,derij
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3) :: g
      real, dimension (nx) :: tmp
      integer :: i,j,k
!
      intent(in) :: f,k
      intent(out) :: g
!
!  Run though all 9 possibilities, treat diagonals separately.
!
      do j=1,3
        call der2 (f,k,tmp,j); g(:,j,j)=tmp
        do i=j+1,3
          call derij(f,k,tmp,i,j); g(:,i,j)=tmp; g(:,j,i)=tmp
        enddo
      enddo
!
    endsubroutine g2ij
!***********************************************************************
    subroutine del4(f,k,del4f,ignoredx)
!
!  Calculate del4 (defined here as d^4/dx^4 + d^4/dy^4 + d^4/dz^4, rather
!  than del2^3) of a scalar for hyperdiffusion.
!
!   8-jul-02/wolf: coded
!   9-dec-03/nils: adapted from del6
!
      use Deriv, only: der4
!
      intent(in) :: f,k,ignoredx
      intent(out) :: del4f
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: del4f,d4fdx,d4fdy,d4fdz
      integer :: k
      logical, optional :: ignoredx
      logical :: ignore_dx
!
      if (present(ignoredx)) then
        ignore_dx = ignoredx
      else
        ignore_dx = .false.
!
!  Exit if this is requested for non-cartesian runs.
!
        if (lcylindrical_coords.or.lspherical_coords) &
            call fatal_error('del4', &
            'not implemented for non-cartesian coordinates')
      endif
!
      call der4(f,k,d4fdx,1,ignore_dx)
      call der4(f,k,d4fdy,2,ignore_dx)
      call der4(f,k,d4fdz,3,ignore_dx)
      del4f = d4fdx + d4fdy + d4fdz
!
    endsubroutine del4
!***********************************************************************
    subroutine del6(f,k,del6f,ignoredx)
!
!  Calculate del6 (defined here as d^6/dx^6 + d^6/dy^6 + d^6/dz^6, rather
!  than del2^3) of a scalar for hyperdiffusion. Using INGOREDX
!  calculates something similar to del6, but ignoring the steps dx, dy, dz.
!  Useful for Nyquist filtering, where you just want to remove the
!  Nyquist frequency fully, while retaining the amplitude in small wave
!  numbers.
!
!  8-jul-02/wolf: coded
!  22-jul-11/bing: added ignoredx
!
      use Deriv, only: der6
!
      intent(in) :: f,k,ignoredx
      intent(out) :: del6f
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: del6f,d6fdx,d6fdy,d6fdz
      integer :: k
      logical, optional :: ignoredx
      logical :: ignore_dx
!
      if (present(ignoredx)) then
        ignore_dx = ignoredx
      else
        ignore_dx = .false.
!
!  Exit if this is requested for lspherical_coords run.
!
!!!!        if (lspherical_coords.or.lcylindrical_coords) &
!!!!            call fatal_error('del6', &
!!!!            'not implemented for non-cartesian coordinates')
      endif
!
      call der6(f,k,d6fdx,1,ignore_dx)
      call der6(f,k,d6fdy,2,ignore_dx)
      call der6(f,k,d6fdz,3,ignore_dx)
!
      del6f = d6fdx + d6fdy + d6fdz
!
    endsubroutine del6
!***********************************************************************
    subroutine del6_other(f,del6f)
!
!  Calculate del6 (defined here as d^6/dx^6 + d^6/dy^6 + d^6/dz^6, rather
!  than del2^3) of a scalar for hyperdiffusion.
!
!  13-jun-05/anders: adapted from del6
!
      use Deriv, only: der6_other
!
      intent(in) :: f
      intent(out) :: del6f
!
      real, dimension (mx,my,mz) :: f
      real, dimension (nx) :: del6f,d6fdx,d6fdy,d6fdz
!
      call der6_other(f,d6fdx,1)
      call der6_other(f,d6fdy,2)
      call der6_other(f,d6fdz,3)
      del6f = d6fdx + d6fdy + d6fdz
!
!  Exit if this is requested for non-cartesian runs.
!
      if (lcylindrical_coords.or.lspherical_coords) &
          call fatal_error('del6_other', &
          'not implemented for non-cartesian coordinates')
!
    endsubroutine del6_other
!***********************************************************************
    subroutine del6fj(f,vec,k,del6f)
!
!  Calculates fj*del6 (defined here as (vecx*d^6/dx^6 + vecy*d^6/dy^6 + vecz*d^6/dz^6)f )
!  needed for hyperdissipation of a scalar (diffrho) with non-cubic cells where
!  the coefficient depends on resolution. Returns scalar.
!
!  30-oct-06/wlad: adapted from del6
!
      use Deriv, only: der6
!
      intent(in) :: f,k
      intent(out) :: del6f
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: del6f,d6fdx,d6fdy,d6fdz
      real, dimension (3) :: vec
      integer :: k
!
      call der6(f,k,d6fdx,1)
      call der6(f,k,d6fdy,2)
      call der6(f,k,d6fdz,3)
      del6f = vec(1)*d6fdx + vec(2)*d6fdy + vec(3)*d6fdz
!
!  Exit if this is requested for non-cartesian runs.
!
      if (lcylindrical_coords.or.lspherical_coords) &
          call fatal_error('del6fj', &
          'not implemented for non-cartesian coordinates')
!
    endsubroutine del6fj
!***********************************************************************
    subroutine del6fjv(f,vec,k,del6f)
!
!  Calculates fj*del6 (defined here as fx*d^6/dx^6 + fy*d^6/dy^6 + fz*d^6/dz^6)
!  needed for hyperdissipation of vectors (visc, res) with non-cubic cells
!  where the coefficient depends on resolution. Returns vector.
!
!  30-oct-06/wlad: adapted from del6v
!
      intent(in) :: f,k
      intent(out) :: del6f
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: del6f
      real, dimension (nx) :: tmp
      integer :: i,k,k1
      real, dimension (3) :: vec
!
      k1=k-1
      do i=1,3
        call del6fj(f,vec,k1+i,tmp)
        del6f(:,i)=tmp
      enddo
!
!  Exit if this is requested for non-cartesian runs.
!
      if (lcylindrical_coords.or.lspherical_coords) &
          call fatal_error('del2fjv', &
          'not implemented for non-cartesian coordinates')
!
    endsubroutine del6fjv
!***********************************************************************
    subroutine u_dot_grad_vec(f,k,gradf,uu,ugradf,upwind,ladd)
!
!  u.gradu
!
!  21-feb-07/axel+dhruba: added spherical coordinates
!   7-mar-07/wlad: added cylindrical coordinates
!  24-jun-08/MR: ladd added for incremental work
!
      use General, only: loptest
!
      intent(in) :: f,k,gradf,uu,upwind
      intent(out) :: ugradf
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3) :: gradf
      real, dimension (nx,3) :: uu,ff,ugradf,grad_f_tmp
      real, dimension (nx) :: tmp
      integer :: j,k
      logical, optional :: upwind,ladd
!
      if (k<1 .or. k>mfarray) then
        call fatal_error('u_dot_grad_vec','variable index is out of bounds')
        return
      endif
!
      do j=1,3
!
        grad_f_tmp = gradf(:,j,:)
        call u_dot_grad_scl(f,k+j-1,grad_f_tmp,uu,tmp,UPWIND=loptest(upwind))
        if (loptest(ladd)) then
          ugradf(:,j)=ugradf(:,j)+tmp
        else
          ugradf(:,j)=tmp
        endif
!
      enddo
!
!  Adjustments for spherical coordinate system.
!  The following now works for general u.gradA.
!
      if (lspherical_coords) then
        ff=f(l1:l2,m,n,k:k+2)
        ugradf(:,1)=ugradf(:,1)-r1_mn*(uu(:,2)*ff(:,2)+uu(:,3)*ff(:,3))
        ugradf(:,2)=ugradf(:,2)+r1_mn*(uu(:,2)*ff(:,1)-uu(:,3)*ff(:,3)*cotth(m))
        ugradf(:,3)=ugradf(:,3)+r1_mn*(uu(:,3)*ff(:,1)+uu(:,3)*ff(:,2)*cotth(m))
      endif
!
!  The following now works for general u.gradA.
!
      if (lcylindrical_coords) then
        ff=f(l1:l2,m,n,k:k+2)
        ugradf(:,1)=ugradf(:,1)-rcyl_mn1*(uu(:,2)*ff(:,2))
        ugradf(:,2)=ugradf(:,2)+rcyl_mn1*(uu(:,1)*ff(:,2))
      endif
!
    endsubroutine u_dot_grad_vec
!***********************************************************************
    subroutine u_dot_grad_vec_alt(f,k,gradf,uu,ugradf,iadvec,ladd)
!
!  u.gradu
!
!  21-feb-07/axel+dhruba: added spherical coordinates
!   7-mar-07/wlad: added cylindrical coordinates
!  24-jun-08/MR: ladd added for incremental work
!
      use General, only: loptest
!
      intent(in) :: f,k,gradf,uu,iadvec
      intent(out) :: ugradf
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3) :: gradf
      real, dimension (nx,3) :: uu,ff,ugradf,grad_f_tmp
      real, dimension (nx) :: tmp
      integer :: j,k,iadvec
      logical, optional :: ladd
      logical :: ladd1
!
      if (k<1 .or. k>mfarray) then
        call fatal_error('u_dot_grad_vec','variable index is out of bounds')
        return
      endif
!
      ladd1=loptest(ladd)
!
      do j=1,3
        grad_f_tmp = gradf(:,j,:)
        call u_dot_grad_scl_alt(f,k+j-1,grad_f_tmp,uu,tmp,iadvec)
        if (ladd1) then
          ugradf(:,j)=ugradf(:,j)+tmp
        else
          ugradf(:,j)=tmp
        endif
!
      enddo
!
!  Adjustments for spherical coordinate system.
!  The following now works for general u.gradA.
!
      if (lspherical_coords) then
        ff=f(l1:l2,m,n,k:k+2)
        ugradf(:,1)=ugradf(:,1)-r1_mn*(uu(:,2)*ff(:,2)+uu(:,3)*ff(:,3))
        ugradf(:,2)=ugradf(:,2)+r1_mn*(uu(:,2)*ff(:,1)-uu(:,3)*ff(:,3)*cotth(m))
        ugradf(:,3)=ugradf(:,3)+r1_mn*(uu(:,3)*ff(:,1)+uu(:,3)*ff(:,2)*cotth(m))
      endif
!
!  The following now works for general u.gradA.
!
      if (lcylindrical_coords) then
        ff=f(l1:l2,m,n,k:k+2)
        ugradf(:,1)=ugradf(:,1)-rcyl_mn1*(uu(:,2)*ff(:,2))
        ugradf(:,2)=ugradf(:,2)+rcyl_mn1*(uu(:,1)*ff(:,2))
      endif
!
    endsubroutine u_dot_grad_vec_alt
!***********************************************************************
    subroutine u_dot_grad_mat(f,k,gradM,uu,ugradM,upwind)
!
!  Computes  u.grad(M)
!  where M is a second rank matrix.
!
!  07-aug-10/dhruba: coded
! 24-nov-11/dhruba: added upwinding
! 26-mar-12/MR: doupwind introduced
!
      intent(in) :: gradM,f,k
      intent(out) :: ugradM
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3,3,3) :: gradM
      real,dimension(nx,3) :: uu
      real, dimension (nx,3,3) :: ugradM
      integer :: k
      logical, optional :: upwind
      integer :: ipi,ipj,ipk
!
      if (k<1 .or. k>mfarray) then
        call fatal_error('u_dot_grad_mat','variable index is out of bounds')
        return
      endif
!
      call vec_dot_3tensor(uu,gradM,ugradM)
!
!  Test if Upwind is used.
!
      if (present(upwind)) then
        if (upwind) then
!
! The same operation needs to be done to each element
! of the matrix ugradM. We assume below that this matrix
! is symmetric. Otherwise the following should be rewritten.
!
          ipk=0
          do ipi=1,3
            do ipj=ipi,3
              call doupwind(f,k+ipk,uu,ugradM(1,ipi,ipj))
              ipk=ipk+1
            enddo
          enddo
          call symmetrise3x3_ut2lt(ugradM)
        endif
      endif
!
!  Spherical and cylindrical coordinates are not
!  implemented for this subroutine.
!
      if (lspherical_coords) then
        call fatal_error('u_dot_grad_mat','not implemented in sph-coordinates')
      endif
!
      if (lcylindrical_coords) then
        call fatal_error('u_dot_grad_mat','not implemented in cyl-coordinates')
      endif
!
    endsubroutine u_dot_grad_mat
!***********************************************************************
    subroutine u_dot_grad_scl(f,k,gradf,uu,ugradf,upwind,ladd)
!
!  Do advection-type term u.grad f_k.
!  Assumes gradf to be known, but takes f and k as arguments to be able
!  to calculate upwind correction.
!
!  28-Aug-2007/dintrans: attempt of upwinding in cylindrical coordinates
!  29-Aug-2007/dhruba: attempt of upwinding in spherical coordinates.
!  28-Sep-2009/MR: ladd added for incremental work
!  26-mar-12/MR: doupwind introduced
!
      use General, only: loptest
!
      intent(in) :: f,k,gradf,uu,upwind,ladd
      intent(out) :: ugradf
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: uu,gradf
      real, dimension (nx) :: ugradf
      integer :: k
      logical, optional :: upwind, ladd
!
      if (k<1 .or. k>mfarray) then
        call fatal_error('u_dot_grad_scl','variable index is out of bounds')
        return
      endif
!
      call dot_mn(uu,gradf,ugradf,loptest(ladd))
!
!  Upwind correction
!
      if (present(upwind)) then
        if (upwind) call doupwind(f,k,uu,ugradf)
      endif
!
    endsubroutine u_dot_grad_scl
!***********************************************************************
    subroutine u_dot_grad_scl_alt(f,k,gradf,uu,ugradf,iadvec,ladd)
!
!  Do advection-type term u.grad f_k.
!  Assumes gradf to be known, but takes f and k as arguments to be able
!  to calculate upwind correction.
!
!  28-Aug-2007/dintrans: attempt of upwinding in cylindrical coordinates
!  29-Aug-2007/dhruba: attempt of upwinding in spherical coordinates.
!  28-Sep-2009/MR: ladd added for incremental work
!  22-Jun-2011/dhruba: made this alternative version which also incorporated the
! kurganov-tadmore scheme.
!  26-mar-12/MR: doupwind introduced
!
      use General, only: loptest
!
      intent(in) :: f,k,gradf,uu,iadvec,ladd
      intent(out) :: ugradf
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: uu,gradf
      real, dimension (nx) :: ugradf, udelf
      logical, optional :: ladd
      integer :: k,iadvec
      logical :: ladd1
!
! iadvec=0 normal scheme
! iadvec=1 upwinding
! iadvec=2 Kurganov-Tadmor (KT)
!
      ladd1=loptest(ladd)
!
      select case (iadvec)
      case (0)
        call dot_mn(uu,gradf,ugradf,ladd1)
      case (1)
        call dot_mn(uu,gradf,ugradf,ladd1)
        call doupwind(f,k,uu,ugradf)
      case (2)
!
! x, y and z directions respectively
!
        call u_grad_kurganov_tadmore(f,k,udelf,1)
        ugradf=udelf
        call u_grad_kurganov_tadmore(f,k,udelf,2)
        ugradf=ugradf+udelf
        call u_grad_kurganov_tadmore(f,k,udelf,3)
        ugradf=ugradf+udelf
!
      case default
        if (lroot) print*, 'sub:u_dot_grad: iadvec must be 0,1 or 2, iadvec=', &
              iadvec
          call fatal_error('u_dot_grad_scl','')
        endselect
!
    endsubroutine u_dot_grad_scl_alt
!***********************************************************************
    subroutine u_grad_kurganov_tadmore(f,k,udelf,j)
!
      use Deriv, only: der2_minmod
!
      intent(in) :: f,k,j
      intent(out) :: udelf
!
      real, dimension (mx,my,mz,mfarray) :: f
      real,dimension(nx) :: udelf
      real, dimension (nx) :: vel,velpj,velmj,amhalf,aphalf,delfj,delfjp1,delfjm1
      integer :: k,j
      integer :: ix,iix
!
      vel=f(l1:l2,m,n,iuu+j-1)
      call der2_minmod(f,k,delfj,delfjp1,delfjm1,j)
      select case(j)
      case(1)
        do ix=l1,l2
          iix=ix-nghost
          velpj(iix) = f(ix+1,m,n,iux)
          velmj(iix) = f(ix-1,m,n,iux)
        enddo
      case(2)
        velpj = f(l1:l2,m+1,n,iuy)
        velmj = f(l1:l2,m-1,n,iuy)
      case(3)
        velpj = f(l1:l2,m,n+1,iuz)
        velmj = f(l1:l2,m,n-1,iuz)
      case default
        call fatal_error('sub:u_grad_kurganov_tadmore:','wrong component')
      endselect
      aphalf = abs(vel+velpj)/2.
      amhalf = abs(vel+velmj)/2.
      udelf = (aphalf-vel)*(delfjp1+delfj)/2. + &
          (amhalf-vel)*(delfjm1+delfj)/2.
!
    endsubroutine u_grad_kurganov_tadmore
!***********************************************************************
    subroutine nou_dot_grad_scl(gradf,uu,ugradf,del6u,upwind,ladd)
!
!  Do advection-type term u.grad f_k.
!  Assumes gradf to be known, but takes f and k as arguments to be able
!  to calculate upwind correction
!
!  28-Aug-2007/dintrans: attempt of upwinding in cylindrical coordinates
!  29-Aug-2007/dhruba: attempt of upwinding in spherical coordinates.
!  28-Sep-2009/MR: ladd added for incremental work
!
      use General, only: loptest
!
      intent(in) :: gradf,uu,upwind,ladd
      intent(out) :: ugradf
!
      real, dimension (nx,3) :: uu,gradf,del6u
      real, dimension (nx) :: ugradf
      logical, optional :: upwind,ladd
!
      call dot_mn(uu,gradf,ugradf,loptest(ladd))
!
!  Upwind correction (currently just for z-direction).
!
      if (present(upwind)) then; if (upwind) then
!
!  x-direction.
!
        ugradf=ugradf-abs(uu(:,1))*del6u(:,1)
!
!  y-direction.
!
        if (lcartesian_coords) then
          ugradf=ugradf-abs(uu(:,2))*del6u(:,2)
        else
          if (lcylindrical_coords) &
             ugradf=ugradf-rcyl_mn1*abs(uu(:,2))*del6u(:,2)
          if (lspherical_coords) &
             ugradf=ugradf-r1_mn*abs(uu(:,2))*del6u(:,2)
        endif
!
!  z-direction.
!
        if ((lcartesian_coords).or.(lcylindrical_coords)) then
          ugradf=ugradf-abs(uu(:,3))*del6u(:,3)
        else
          if (lspherical_coords) &
             ugradf=ugradf-r1_mn*sin1th(m)*abs(uu(:,3))*del6u(:,3)
        endif
!
      endif; endif
!
    endsubroutine nou_dot_grad_scl
!***********************************************************************
    subroutine h_dot_grad_vec(hh,gradf,ff,hgradf)
!
!  h.gradf for vectors h and f.
!
!  23-mar-08/axel: adapted from u_dot_grad_vec
!
      intent(in) :: hh,gradf,ff
      intent(out) :: hgradf
!
      real, dimension (nx,3,3) :: gradf
      real, dimension (nx,3) :: hh,ff,hgradf
      real, dimension (nx) :: tmp
      integer :: j
!
!  Dot product for each of the three components of gradf .
!
      do j=1,3
        call h_dot_grad_scl(hh,gradf(:,j,:),tmp)
        hgradf(:,j)=tmp
      enddo
!
!  Adjustments for spherical coordinate system.
!  The following now works for general u.gradA .
!
      if (lspherical_coords) then
        hgradf(:,1)=hgradf(:,1)-r1_mn*(hh(:,2)*ff(:,2)+hh(:,3)*ff(:,3))
        hgradf(:,2)=hgradf(:,2)+r1_mn*(hh(:,2)*ff(:,1)-hh(:,3)*ff(:,3)*cotth(m))
        hgradf(:,3)=hgradf(:,3)+r1_mn*(hh(:,3)*ff(:,1)+hh(:,3)*ff(:,2)*cotth(m))
      endif
!
!  The following now works for general u.gradA .
!
      if (lcylindrical_coords) then
        hgradf(:,1)=hgradf(:,1)-rcyl_mn1*(hh(:,2)*ff(:,2))
        hgradf(:,2)=hgradf(:,2)+rcyl_mn1*(hh(:,1)*ff(:,2))
      endif
!
    endsubroutine h_dot_grad_vec
!***********************************************************************
    subroutine h_dot_grad_scl(hh,gradf,hgradf)
!
!  Do advection-type term h.grad f_k, but h is not taken from f array.
!
!  23-mar-08/axel: adapted from u_dot_grad_scl
!
      intent(in) :: hh,gradf
      intent(out) :: hgradf
!
      real, dimension (nx,3) :: hh,gradf
      real, dimension (nx) :: hgradf
!
      call dot_mn(hh,gradf,hgradf)
!
    endsubroutine h_dot_grad_scl
!***********************************************************************
    subroutine gradf_upw1st(f,uu,k,gradf)
!
!  Do advection-type term u.grad f_k for upwind 1st order der scheme.
!
      use Deriv, only: der_upwind1st
!
      intent(in) :: f,uu
      intent(out) :: gradf
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: uu,gradf
      integer :: j,k
!
      do j=1,3
        call der_upwind1st(f,uu,k,gradf(:,j),j)
      enddo
!
    endsubroutine gradf_upw1st
!***********************************************************************
    character function numeric_precision()
!
!  Return 'S' if running in single, 'D' if running in double precision.
!
!  12-jul-06/wolf: extracted from wdim()
!
      integer :: real_prec
!
      real_prec = precision(1.)
      if (real_prec==6 .or. real_prec==7) then
        numeric_precision = 'S'
      elseif (real_prec == 15) then
        numeric_precision = 'D'
      else
        print*, 'WARNING: encountered unknown precision ', real_prec
        numeric_precision = '?'
      endif
!
    endfunction numeric_precision
!***********************************************************************
    subroutine wdim(file,mxout,myout,mzout,lglobal)
!
!  Write dimension to file.
!
!   8-sep-01/axel: adapted to take myout,mzout
!  30-sep-14/MR  : added parameter lglobal to be able to output
!                  dimensions different from mx,my,mz both locally and globally
!
      use General, only: loptest
!
      character (len=*) :: file
      character         :: prec
      integer, optional :: mxout,myout,mzout
      logical, optional :: lglobal
      integer           :: mxout1,myout1,mzout1,iprocz_slowest=0
!
!  Determine whether mxout=mx (as on each processor),
!  or whether mxout is different (e.g. when writing out full array).
!
      if (present(mzout)) then
        mxout1=mxout
        myout1=myout
        mzout1=mzout
      elseif (lmonolithic_io) then
        mxout1=nxgrid+2*nghost
        myout1=nygrid+2*nghost
        mzout1=nzgrid+2*nghost
      else
        mxout1=mx
        myout1=my
        mzout1=mz
      endif
!
!  Only root writes allprocs/dim.dat (with io_mpio.f90),
!  but everybody writes to their procN/dim.dat (with io_dist.f90).
!
      if (lroot .or. .not. lmonolithic_io) then
        open(1,file=file)
        write(1,'(3i7,3i5)') mxout1,myout1,mzout1,mvar,maux,mglobal
!
!  Check for double precision.
!
        prec = numeric_precision()
        write(1,'(a)') prec
!
!  Write number of ghost cells (could be different in x, y and z).
!
        write(1,'(3i5)') nghost, nghost, nghost
        if (loptest(lglobal)) then
          if (lprocz_slowest) iprocz_slowest=1
          write(1,'(4i5)') nprocx, nprocy, nprocz, iprocz_slowest
        else
          write(1,'(3i5)') ipx, ipy, ipz
        endif
!
        close(1)
      endif
!
    endsubroutine wdim
!***********************************************************************
    subroutine rdim(file,mx_in,my_in,mz_in,mvar_in,maux_in,mglobal_in,&
        prec_in,nghost_in,ipx_in, ipy_in, ipz_in)
!
!  Read dimension from file.
!
!  15-sep-09/nils: adapted from rdim
!
      character (len=*) :: file
      character         :: prec_in
      integer           :: mx_in,my_in,mz_in
      integer           :: mvar_in,maux_in,mglobal_in,nghost_in
      integer           :: ipx_in, ipy_in, ipz_in
!
!  Every processor writes to their procN/dim.dat (with io_dist.f90).
!
      open(124,file=file,FORM='formatted')
      read(124,*) mx_in,my_in,mz_in,mvar_in,maux_in,mglobal_in
      read(124,*) prec_in
      read(124,*) nghost_in, nghost_in, nghost_in
      read(124,*) ipx_in, ipy_in, ipz_in
!
      close(124)
!
    endsubroutine rdim
!***********************************************************************
    subroutine read_snaptime(file,tout,nout,dtout,t)
!
!  Read in output time for next snapshot (or similar) from control file.
!
!  30-sep-97/axel: coded
!  24-aug-99/axel: allow for logarithmic spacing
!   9-sep-01/axel: adapted for MPI
!  10-sep-15/MR  : tout set to t if file is missing and dtout>0
!
      use Mpicomm, only: mpibcast_real, MPI_COMM_WORLD
!
      character (len=*), intent(in) :: file
      real, intent(out) :: tout
      integer, intent(out) :: nout
      real, intent(in) :: dtout
      double precision, intent(in) :: t
!
      integer, parameter :: lun = 31
      logical :: exist
      integer, parameter :: nbcast_array=2
      real, dimension(nbcast_array) :: bcast_array
      double precision :: t0
!
      if (lroot) then
!
!  Depending on whether or not file exists, we need to
!  either read or write tout and nout from or to the file.
!
        inquire(FILE=trim(file),EXIST=exist)
        open(lun,FILE=trim(file))
        if (exist) then
          read(lun,*) tout,nout
        else
!
!  Special treatment when dtout is negative.
!  Now tout and nout refer to the next snapshopt to be written.
!
          settout: if (dtout < 0.0) then
            tout = log10(t)
          elseif (dtout /= 0.0) then settout
            !  make sure the tout is a good time
            t0 = max(t - dt, 0.0D0)
            tout = t0 + (dble(dtout) - modulo(t0, dble(dtout)))
          else settout
            call warning("read_snaptime", "Writing snapshot every time step. ")
            tout = 0.0
          endif settout
          nout=1
          write(lun,*) tout,nout
        endif
        close(lun)
!
      endif
!
!  Broadcast tout and nout in one go.
!
      if (lroot) then
        bcast_array(1) = tout
        bcast_array(2) = nout
      endif
!
      call mpibcast_real(bcast_array,nbcast_array,comm=MPI_COMM_WORLD)
      tout = bcast_array(1)
      nout = bcast_array(2)
!
    endsubroutine read_snaptime
!***********************************************************************
    subroutine update_snaptime(file,tout,nout,dtout,t,lout,ch,nowrite)
!
!  Check whether we need to write snapshot; if so, update the snapshot
!  file (e.g. tsnap.dat). Done by all processors.
!
!  30-sep-97/axel: coded
!  24-aug-99/axel: allow for logarithmic spacing
!  27-jul-15/MR  : try to fix a strange behavior with gfortran:
!                  when crashing, a big number of unmotivated snapshots
!                  is output -> test of NaN in t
!
      use General, only: itoa, notanumber_0d
!
      character (len=*), intent(in) :: file
      real, intent(inout) :: tout
      integer, intent(inout) :: nout
      real, intent(in) :: dtout
      double precision, intent(in) :: t
      logical, intent(out) :: lout
      logical, intent(in), optional :: nowrite
      character (len=intlen), intent(out), optional :: ch
!
      integer, parameter :: lun = 31
      logical :: lwrite
      real :: t_sp   ! t in single precision for backwards compatibility
      logical, save :: lfirstcall=.true.
      real, save :: deltat_threshold
!
      if (notanumber_0d(t)) then
        lout=.false.
        return
      endif
!
!  Use t_sp as a shorthand for either t or lg(t).
!
      if (dtout<0.0) then
        t_sp=log10(t)
      else
        t_sp=t
      endif
!
!  Check if no writing tout is requested.
!
      lwrite = .true.
      if (present(nowrite)) lwrite = .not. nowrite
!
!  Generate a running file number, if requested.
!
      if (present (ch)) ch = itoa (nout)
!
!  Mark lout=.true. when time has exceeded the value of tout do while loop to
!  make sure tt is always larger than tout.
!  (otherwise slices are written just to catch up with tt.)
!
!  WL: Add possibility that there should be a small threshold in this
!      comparison. Needed for outputing at the exact tsnap, otherwise
!      a difference between tsp and tout to machine precision can be
!      interpreted as stating that the output is to be done at the next,
!      not the current, timestep.
!
      if (lfirstcall) then
        if (.not.loutput_varn_at_exact_tsnap) then
          deltat_threshold=0.0
        else
          deltat_threshold=dtmin
        endif
        lfirstcall=.false.
      endif
!
      if (    (t_sp       >= tout)             .or. &
          (abs(t_sp-tout) <  deltat_threshold)) then
        tout=tout+abs(dtout)
        nout=nout+1
        lout=.true.
!
!  Write corresponding value of tout to file to make sure we have it, in case
!  the code craches. If the disk is full, however, we need to reset the values
!  manually.
!
        writenext: if (lroot .and. lwrite) then
          open(lun,FILE=trim(file))
          write(lun,*) tout,nout
          write(lun,*) 'This file is written automatically (routine'
          write(lun,*) 'check_snaptime in sub.f90). The values above give'
          write(lun,*) 'time and number of the *next* snapshot. These values'
          write(lun,*) 'are only read once in the beginning. You may adapt'
          write(lun,*) 'them by hand (eg after a crash).'
          close(lun)
        endif writenext
      else
        lout=.false.
      endif
!
    endsubroutine update_snaptime
!***********************************************************************
    subroutine vecout(lun,file,vv,thresh,nvec)
!
!  Write vectors to disc if their length exceeds thresh.
!
!  22-jul-03/axel: coded
!
      character (len=*) :: file
      real, dimension(nx,3) :: vv
      real, dimension(nx) :: v2
      real :: thresh,thresh2,dummy=0.
      integer :: l,lun,nvec
      real :: t_sp   ! t in single precision for backwards compatibility
!
      t_sp = t
!
!  Return if thresh=0 (default).
!
      if (thresh==0.) return
!
!  Open files when first data point.
!
      if (lfirstpoint) then
        open(lun,FILE=trim(file)//'.dat',form='unformatted',position='append')
        write(lun) 0,0,0,t_sp,dummy,dummy  !(marking first line)
        nvec=0
      endif
!
!  Write data.
!
      thresh2=thresh**2
      v2=vv(:,1)**2+vv(:,2)**2+vv(:,3)**2
      do l=1,nx
        if (v2(l)>=thresh2) then
          write(lun) l,m-nghost,n-nghost,vv(l,:)
          nvec=nvec+1
        endif
      enddo
!
!  Close file, and write number of vectors to a separate file.
!
      if (llastpoint) then
        close(lun)
        open(lun,FILE=trim(file)//'.num',position='append')
        write(lun,*) t_sp,nvec
        close(lun)
      endif
!
    endsubroutine vecout
!***********************************************************************
    subroutine despike(f,j,retval,factor)
!
!  Remove large spikes from
!
!  14-aug-06/tony: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension(nx) :: retval
      real, dimension (mx) :: tmp_penc
      real, dimension (mx) :: meanf
      real :: factor
      real, parameter :: t1 = 1./26.
      real, parameter :: t2 = 0.70710678/26.
      real, parameter :: t3 = 0.57735027/26.
      real, parameter :: t4 = 0.
      real, parameter, dimension (-1:1,-1:1,-1:1) :: interp3D = reshape(&
            (/ t3, t2, t3, &
               t2, t1, t2, &
               t3, t2, t3, &
               t2, t1, t2, &
               t1, t4, t1, &
               t2, t1, t2, &
               t3, t2, t3, &
               t2, t1, t2, &
               t3, t2, t3 /),&
            (/ 3,3,3 /))
      integer :: ii,jj,kk
      integer :: j
!
      meanf=0.
      if ((nxgrid/=1).and.(nygrid/=1).and.(nzgrid/=1)) then
        tmp_penc=f(:,m,n,j)
        do kk=-1,1
        do jj=-1,1
        do ii=-1,1
          if (ii/=0.or.jj/=0.or.kk/=0) &
          meanf(3:mx-2)=meanf(3:mx-2)+interp3D(ii,jj,kk)*(f(3+ii:mx-2+ii,m+jj,n+kk,j)-tmp_penc(3:mx-2))
        enddo
        enddo
        enddo
      else
        call fatal_error('shock_max3_pencil_interp', &
            'Tony got lazy and only implemented the 3D case')
      endif
!
!      factor1=1./factor
      retval=max(meanf(l1:l2)*factor,f(l1:l2,m,n,j))
!      retval=max(meanf(l1:l2)*factor,retval)
!
    endsubroutine despike
!***********************************************************************
    subroutine smooth_kernel(f,j,smth)
!
!  Smooth scalar field FF using predefined constant gaussian like kernel.
!
!  20-jul-06/tony: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension(nx) :: smth
      integer :: j,l
!
      do l=l1,l2
        smth(l-l1+1)=sum(smth_kernel*f(l-3:l+3,m-3:m+3,n-3:n+3,j))
      enddo
!
    endsubroutine smooth_kernel
!***********************************************************************
    subroutine smooth(f, ivar)
!
!  Smoothes the f-variable ivar with a polynomial kernel.  It assumes
!  that the boundary conditions for ivar have been applied, and the
!  ghost cells are not treated upon return.
!
!  23-jan-14/ccyang: coded.
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
      integer, intent(in) :: ivar
!
      real, dimension(-nghost:nghost,-nghost:nghost,-nghost:nghost) :: kernel = 0.0
      real, dimension(nx,ny,nz) :: work
      logical :: lfirstcall = .true.
      integer :: i1 = 0, i2 = 0, j1 = 0, j2 = 0, k1 = 0, k2 = 0
      integer :: i, j, k
!
!  Initialization at first call.
!
      init: if (lfirstcall) then
        call get_smooth_kernel(kernel)
        xdim: if (nxgrid > 1) then
          i1 = -nghost
          i2 = nghost
        endif xdim
        ydim: if (nygrid > 1) then
          j1 = -nghost
          j2 = nghost
        endif ydim
        zdim: if (nzgrid > 1) then
          k1 = -nghost
          k2 = nghost
        endif zdim
        lfirstcall = .false.
      endif init
!
!  Smooth the field (assuming boundary conditions for ivar has been applied).
!
      work = 0.0
      do k = k1, k2
        do j = j1, j2
          do i = i1, i2
            work = work + kernel(i,j,k) * f(l1+i:l2+i,m1+j:m2+j,n1+k:n2+k,ivar)
          enddo
        enddo
      enddo
      f(l1:l2,m1:m2,n1:n2,ivar) = work
!
    endsubroutine smooth
!***********************************************************************
    subroutine get_smooth_kernel(kernel)
!
!  Gets the smooth kernel used by subroutine smooth.
!
!  15-feb-14/ccyang: coded
!
      real, dimension(-nghost:nghost,-nghost:nghost,-nghost:nghost), intent(out) :: kernel
!
      integer, dimension(-nghost:nghost) :: b
      real, dimension(-nghost:nghost) :: weight
      integer :: k
!
      call binomial(b)
      weight = real(b)
      k = 2 * nghost + 1
      kernel = 0.0
      kernel(0,0,0) = 1.0
      if (nxgrid > 1) kernel(:,0,0) = kernel(0,0,0) * weight
      if (nygrid > 1) kernel(:,:,0) = spread(kernel(:,0,0), 2, k) * spread(weight, 1, k)
      if (nzgrid > 1) kernel = spread(kernel(:,:,0), 3, k) * spread(spread(weight, 1, k), 1, k)
      kernel = kernel / sum(kernel)
!
    endsubroutine get_smooth_kernel
!***********************************************************************
    subroutine binomial(b)
!
!  Finds the binomial coefficients.
!
!  23-jan-14/ccyang: coded.
!
      integer, dimension(:), intent(out) :: b
!
      integer :: n, i
!
      n = size(b) - 1
!
      b(1) = 1
      do i = 1, n - 1
        b(i+1) = b(i) * (n - i + 1) / i
      enddo
      b(n+1) = 1
!
    endsubroutine binomial
!***********************************************************************
    subroutine identify_bcs(varname_input,idx)
!
!  Print boundary conditions for scalar field.
!
!  19-jul-02/wolf: coded
!  29-may-04/axel: allowed variable name to be 8 chars long
!
      character (len=*) :: varname_input
      integer :: idx
!
      write(*,'(A,A10,",  x: <",A8,">, y: <",A8,">,  z: <",A8,">")') &
           'Bcs for ', varname_input, &
           trim(bcx(idx)), trim(bcy(idx)), trim(bcz(idx))
!
    endsubroutine identify_bcs
!***********************************************************************
    function noform(cname,lcomplex)
!
!  Given a string of the form `name(format)',
!  returns the name without format, fills empty space
!  of correct length (depending on format) with dashes.

!  For output as legend.dat and first line of time_series.dat.
!
!  22-jun-02/axel: coded
!  20-aug-13/MR: optional parameter lcomplex added
!  26-aug-13/MR: unnecessary p descriptors removed from cform
!
      use Cparam, only: max_col_width
      use General, only: loptest
!
      character (len=*) :: cname
      logical, optional :: lcomplex
      character (len=max_col_width) :: noform,cform,cnumber,dashes
      integer :: index_e,index_f,index_g,index_i,index_d,index_r,index1,index2
      integer :: iform0,iform1,iform2,length,number,number1,number2,io_code
!
      intent(in)  :: cname
!
!  Fill DASHES with, well, dashes.
!
      dashes = repeat('-', max_col_width)
!
!  Find position of left bracket to isolate format, cform.
!
      iform0=index(cname,' ')
      iform1=index(cname,'(')
      iform2=index(cname,')')
!
!  Set format; use default if not given.
!  Here we keep the parenthesis in cform.
!
      if (iform1>0) then
        cform=cname(iform1:iform2)
        length=iform1-1
      else
        cform='(e10.2)'
        length=iform0-1
      endif
!
!  Find length of formatted expression, examples: f10.2, e10.3, g12.1 .
!  index_1 is the position of the format type (f,e,g), and
!  index_d is the position of the dot.
!
      index_e=scan(cform,'eE')
      index_f=scan(cform,'fF')
      index_g=scan(cform,'gG')
      index_i=scan(cform,'iI')
      index_d=index(cform,'.')
      index_r=index(cform,')')
      index1=max(index_e,index_f,index_g,index_i)
      index2=index_d; if (index_d==0) index2=index_r
!
!  Calculate the length of the format and assemble expression for legend.
!
      cnumber=cform(index1+1:index2-1)
      read(cnumber,'(i4)',iostat=io_code) number
      if (io_code /= 0) then
!
!  Error while reading or end of file.
!
        print*,'noform: formatting problem'
        print*,'problematic cnumber= <',cnumber,'>'
        number=10
      endif
      if (loptest(lcomplex)) number = 2*number+4
      number1=max(0,(number-length)/2)
      number2=max(1,number-length-number1) ! at least one separating dash
!
!  Sanity check.
!
      if (number1+length+number2 > max_col_width) then
        call error("noform", &
                   "Increase max_col_width or sanitize print.in{,.double}")
      endif
!
      noform=dashes(1:number1)//cname(1:length)//dashes(1:number2)
!
    endfunction noform
!***********************************************************************
    function levi_civita(i,j,k)
!
!  Totally antisymmetric tensor.
!
!  20-jul-03/axel: coded
!
      real :: levi_civita
      integer :: i,j,k
!
      if ( &
        (i==1 .and. j==2 .and. k==3) .or. &
        (i==2 .and. j==3 .and. k==1) .or. &
        (i==3 .and. j==1 .and. k==2) ) then
        levi_civita=1.
      elseif ( &
        (i==3 .and. j==2 .and. k==1) .or. &
        (i==1 .and. j==3 .and. k==2) .or. &
        (i==2 .and. j==1 .and. k==3) ) then
        levi_civita=-1.
      else
        levi_civita=0.
      endif
!
    endfunction levi_civita
!***********************************************************************
    function kronecker_delta(i,j)
!
!  \delta_{ij} = 1 if i==j, 0 otherwise
!
!  28-oct-11/dhruba: coded
!
      real :: kronecker_delta
      integer :: i,j
!
      if (i==j)  then
        kronecker_delta = 1.
      else
        kronecker_delta=0.
      endif
!
    endfunction
!***********************************************************************
    function poly_1(coef, x)
!
!  Horner's scheme for polynomial evaluation.
!  Version for 1d array.
!
!  17-jan-02/wolf: coded
!
      real, dimension(:) :: coef
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: poly_1
      integer :: Ncoef,i
!
      Ncoef = size(coef,1)
!
      poly_1 = coef(Ncoef)
      do i=Ncoef-1,1,-1
        poly_1 = poly_1*x+coef(i)
      enddo
!
    endfunction poly_1
!***********************************************************************
    function poly_0(coef, x)
!
!  Horner's scheme for polynomial evaluation.
!  Version for scalar.
!
!  17-jan-02/wolf: coded
!
      real, dimension(:) :: coef
      real :: x
      real :: poly_0
      integer :: Ncoef,i
!
      Ncoef = size(coef,1)
!
      poly_0 = coef(Ncoef)
      do i=Ncoef-1,1,-1
        poly_0 = poly_0*x+coef(i)
      enddo
!
    endfunction poly_0
!***********************************************************************
    function poly_3(coef, x)
!
!  Horner's scheme for polynomial evaluation.
!  Version for 3d array.
!
!  17-jan-02/wolf: coded
!
      real, dimension(:) :: coef
      real, dimension(:,:,:) :: x
      real, dimension(size(x,1),size(x,2),size(x,3)) :: poly_3
      integer :: Ncoef,i
!
      Ncoef = size(coef,1)
!
      poly_3 = coef(Ncoef)
      do i=Ncoef-1,1,-1
        poly_3 = poly_3*x+coef(i)
      enddo
!
    endfunction poly_3
!***********************************************************************
    subroutine lower_triangular_index(ij,i1,j1)
      integer,intent(out)::ij
      integer,intent(in) :: i1,j1
      integer :: ii,jj
      ii=i1;jj=j1
      if (i1.lt.j1) then
        ii=j1
        jj=i1
      endif
      ij=ii*(ii-1)/2 + jj
    endsubroutine lower_triangular_index
!***********************************************************************
    recursive function ylm(ell,emm,der) result (sph_har)
!
!  Spherical harmonic
!
!   24-nov-14/dhruba: copied from step
!
      real :: sph_har
      real :: theta,phi
      integer :: ell,emm
      real, optional :: der
!
      real :: cos2p,cost,sint,cosp,sinp
      integer :: aemm
      logical :: lother
!
! the one over pi, cosines and sines below may be pre-calculated
!
      lother=.false.
      cost=costh(m); sint=sinth(m); cosp=cosph(n); sinp=sinph(n)
      goto 1

      entry ylm_other(theta,phi,ell,emm,der) result (sph_har)

      lother=.true.
      cost=cos(theta); sint=sin(theta); cosp=cos(phi); sinp=sin(phi)
      
 1    aemm=iabs(emm)
      select case (ell)
          case (0)
            sph_har=(0.5)*sqrt(1./pi)
          case (1)
            select case(aemm)
              case (0)
                sph_har=(0.5)*sqrt(3./pi)*cost
              case (1) 
                sph_har=(0.5)*sqrt(3./(2*pi))*sint*cosp
                if (emm<0) sph_har = -sph_har       ! Condon-Shortley phase
              case default
                call fatal_error('sub:ylm','l=1 wrong m ')
              endselect
         case (2)
            if (aemm==2) cos2p=2*cosp*cosp-1
            select case(aemm)
              case (0)
                sph_har=(0.25)*sqrt(5./pi)*(3.*cost*cost-1.)
              case (1) 
                sph_har=-(0.5)*sqrt(15./(2*pi))*sint*cost*cosp
                if (emm<0) sph_har = -sph_har       ! Condon-Shortley phase
              case (2)
                sph_har=(0.25)*sqrt(15./(2*pi))*sint*sint*cos2p
              case default
                call fatal_error('sub:ylm','l=2 wrong m ')
              endselect
          case default
            call fatal_error('sub:ylm','your ylm is not implemented')
      endselect

      if (present(der)) then
        if (lother) then
          der = ( ell*cost*sph_har - sqrt((2.*ell+1.)*(ell-aemm)*(ell+aemm)/(2.*ell-1.))* &
                  ylm_other(theta,phi,ell-1,emm) )/sint
        else
          der = ( ell*cost*sph_har - sqrt((2.*ell+1.)*(ell-aemm)*(ell+aemm)/(2.*ell-1.))* &
                  ylm(ell-1,emm) )/sint
        endif
      endif
!
    endfunction ylm
!***********************************************************************
    function step_scalar(x,x0,width)
!
!  Smooth unit step function centred at x0; implemented as tanh profile.
!
!   5-sep-08/dhruba: copied from step
!   9-nov-10/axel: no need to have the tini here
!
      real :: x
      real :: step_scalar
      real :: x0,width
!
!  check whether width is finite.
!
      if (width==0) call fatal_error('step_scalar','width must not be zero')
      step_scalar = 0.5*(1+tanh((x-x0)/width))
!
    endfunction step_scalar
!***********************************************************************
    function step_vector(x,x0,width)
!
!  Smooth unit step function centred at x0; implemented as tanh profile
!
!  23-jan-02/wolf: coded
!
      real, dimension(:) :: x
      real, dimension(size(x)) :: step_vector
      real :: x0,width
!
      step_vector = 0.5*(1+tanh((x-x0)/(width+tini)))
!
    endfunction step_vector
!***********************************************************************
    function der_step(x,x0,width)
!
!  Derivative of smooth unit STEP() function given above (i.e. a bump profile).
!  Adapt this if you change the STEP() profile, or you will run into
!  inconsistenies.
!  MR: perhaps to be merged with step
!
!  23-jan-02/wolf: coded
!
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: der_step,arg
      real :: x0,width
!
!  Some argument gymnastics to avoid `floating overflow' for large
!  arguments.
!
      arg = abs((x-x0)/(width+tini))
      arg = min(arg,8.)         ! cosh^2(8) = 3e+27
      der_step = 0.5/(width*cosh(arg)**2)
!
      endfunction der_step
!***********************************************************************
    function der6_step(x,x0,width)
!
!  6th order derivative of smooth unit STEP() function given
!  above (i.e. a bump profile).
!  Adapt this if you change the STEP() profile, or you will run into
!  inconsistenies.
!
!  08-dec-09/dhruba: aped from der_step
!
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: der6_step,arg,sechx,tanhx
      real :: x0,width
!
!  Some argument gymnastics to avoid `floating overflow' for large
!  arguments.
!
      arg = abs((x-x0)/(width+tini))
      tanhx=tanh(arg)
      arg = min(arg,8.)         ! cosh^2(8) = 3e+27
      sechx=1./cosh(arg)
      der6_step = (1./(2*width**6))*(&
             -272.0*(sechx**6)*tanhx+416.0*(sechx**4)*(tanhx**3) &
             -32.0*(sechx**2)*(tanhx**5) )
!
      endfunction der6_step
!***********************************************************************
    function stepdown(x,x0,width)
!
!  Smooth unit step function centred at x0; implemented as tanh profile
!
!  23-jan-02/wolf: coded
!
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: stepdown
      real :: x0,width
!
      stepdown = -0.5*(1+tanh((x-x0)/(width+tini)))
!
    endfunction stepdown
!***********************************************************************
    function der_stepdown(x,x0,width)
!
!  Derivative of smooth unit STEPDOWN() function given above (i.e. a bump profile).
!  Adapt this if you change the STEP() profile, or you will run into
!  inconsistenies.
!
!  27-mar-10/dhruba: aped from der_step
!
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: der_stepdown,arg
      real :: x0,width
!
!  Some argument gymnastics to avoid `floating overflow' for large
!  arguments.
!
      arg = abs((x-x0)/(width+tini))
      arg = min(arg,8.)         ! cosh^2(8) = 3e+27
      der_stepdown = -0.5/(width*cosh(arg)**2)
!
    endfunction der_stepdown
!***********************************************************************
    function cubic_step_pt(x,x0,width,shift)
!
!  Smooth unit step function with cubic (smooth) transition over [x0-w,x0+w].
!  Optional argument SHIFT shifts center:
!  for shift=1. the interval is [x0    ,x0+2*w],
!  for shift=-1. it is          [x0-2*w,x0    ].
!  This is to make sure the interior region is not affected.
!  Maximum slope is 3/2=1.5 times that of a linear profile.
!
!  This version is for scalar args.
!
!  18-apr-04/wolf: coded
!
      real :: x
      real :: cubic_step_pt,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      xi = (x-x0)/(width+tini) - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      cubic_step_pt = 0.5 + xi*(0.75-xi**2*0.25)
!
    endfunction cubic_step_pt
!***********************************************************************
    function cubic_step_mn(x,x0,width,shift)
!
!  Smooth unit step function with cubic (smooth) transition over [x0-w,x0+w].
!  Version for 1d arg (in particular pencils).
!
!  18-apr-04/wolf: coded
!
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: cubic_step_mn,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      xi = (x-x0)/(width+tini) - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      cubic_step_mn = 0.5 + xi*(0.75-xi**2*0.25)
!
    endfunction cubic_step_mn
!***********************************************************************
    function cubic_der_step_pt(x,x0,width,shift)
!
!  Derivative of smooth unit step function, localized to [x0-w,x0+w].
!  This version is for scalar args.
!
!  12-jul-05/axel: adapted from cubic_step_pt
!
      real :: x
      real :: cubic_der_step_pt,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift,width1
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      width1 = 1./(width+tini)
      xi = (x-x0)*width1 - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      cubic_der_step_pt = (0.75-xi**2*0.75) * width1
!
    endfunction cubic_der_step_pt
!***********************************************************************
    function cubic_der_step_mn(x,x0,width,shift)
!
!  Derivative of smooth unit step function, localized to [x0-w,x0+w].
!  Version for 1d arg (in particular pencils).
!
!  12-jul-05/axel: adapted from cubic_step_mn
!
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: cubic_der_step_mn,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift,width1
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      width1 = 1./(width+tini)
      xi = (x-x0)*width1 - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      cubic_der_step_mn = (0.75-xi**2*0.75) * width1
!
    endfunction cubic_der_step_mn
!***********************************************************************
    function quintic_step_pt(x,x0,width,shift)
!
!  Smooth unit step function with quintic (smooth) transition over [x0-w,x0+w].
!  Optional argument SHIFT shifts center:
!  for shift=1. the interval is [x0    ,x0+2*w],
!  for shift=-1. it is          [x0-2*w,x0    ].
!  Maximum slope is 15/8=1.875 times that of a linear profile.
!
!  This version is for scalar args.
!
!  09-aug-05/wolf: coded
!
      real :: x
      real :: quintic_step_pt,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      xi = (x-x0)/(width+tini) - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      quintic_step_pt = 0.5 + xi*(0.9375 + xi**2*(-0.625 + xi**2*0.1875))
!
    endfunction quintic_step_pt
!***********************************************************************
    function quintic_step_mn(x,x0,width,shift)
!
!  Smooth unit step function with quintic (smooth) transition over [x0-w,x0+w].
!
!  Version for 1d arg (in particular pencils).
!
!  09-aug-05/wolf: coded
!
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: quintic_step_mn,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      xi = (x-x0)/(width+tini) - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      quintic_step_mn = 0.5 + xi*(0.9375 + xi**2*(-0.625 + xi**2*0.1875))
!
    endfunction quintic_step_mn
!***********************************************************************
    function quintic_der_step_pt(x,x0,width,shift)
!
!  Derivative of smooth unit step function, localized to [x0-w,x0+w].
!
!  This version is for scalar args.
!
!  09-aug-05/wolf: coded
!
      real :: x
      real :: quintic_der_step_pt,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift,width1
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      width1 = 1./(width+tini)
      xi = (x-x0)*width1 - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      quintic_der_step_pt = (0.9375 + xi**2*(-1.875 + xi**2*0.9375)) &
                              * width1
!
    endfunction quintic_der_step_pt
!***********************************************************************
    function quintic_der_step_mn(x,x0,width,shift)
!
!  Derivative of smooth unit step function, localized to [x0-w,x0+w].
!
!  Version for 1d arg (in particular pencils).
!
!  09-aug-05/wolf: coded
!
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: quintic_der_step_mn,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift,width1
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      width1 = 1./(width+tini)
      xi = (x-x0)*width1 - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      quintic_der_step_mn = (0.9375 + xi**2*(-1.875 + xi**2*0.9375)) &
                              * width1
!
    endfunction quintic_der_step_mn
!***********************************************************************
    function sine_step_pt(x,x0,width,shift)
!
!  Smooth unit step function with sine (smooth) transition over [x0-w,x0+w].
!  Optional argument SHIFT shifts center:
!  for shift=1. the interval is [x0    ,x0+2*w],
!  for shift=-1. it is          [x0-2*w,x0    ].
!  Maximum slope is 15/8=1.875 times that of a linear profile.
!
!  This version is for scalar args.
!
!  13-jun-06/tobi: Adapted from cubic_step
!
      real :: x
      real :: sine_step_pt,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      xi = (x-x0)/(width+tini) - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      sine_step_pt = 0.5*(1+sin(0.5*pi*xi))
!
    endfunction sine_step_pt
!***********************************************************************
    function sine_step_mn(x,x0,width,shift)
!
!  Smooth unit step function with sine (smooth) transition over [x0-w,x0+w].
!
!  Version for 1d arg (in particular pencils).
!
!  13-jun-06/tobi: Adapted from cubic_step
!
      real, dimension(:) :: x
      real, dimension(size(x,1)) :: sine_step_mn,xi
      real :: x0,width
      real, optional :: shift
      real :: relshift
!
      if (present(shift)) then; relshift=shift; else; relshift=0.0; endif
      xi = (x-x0)/(width+tini) - relshift
      xi = max(xi,-1.0)
      xi = min(xi, 1.0)
      sine_step_mn = 0.5*(1+sin(0.5*pi*xi))
!
    endfunction sine_step_mn
!***********************************************************************
    subroutine nan_inform(f,msg,region,int1,int2,int3,int4,lstop)
!
!  Check input array (f or df) for NaN, -Inf, Inf, and output location in
!  array.
!
!  30-apr-04/anders: coded
!  12-jun-04/anders: region or intervals supplied in call
!
      use General, only: notanumber
!
      real, dimension(:,:,:,:) :: f
      character (len=*) :: msg
      integer :: a,b,c,d,a1=1,a2=mx,b1=1,b2=my,c1=1,c2=mz,d1=1,d2=1
      integer, dimension(2), optional :: int1,int2,int3,int4
      character (len=*), optional :: region
      logical, optional :: lstop
!
!  Must set d2 according to whether f or df is considered.
!
      d2 = size(f,4)
!
!  Set intervals for different predescribed regions.
!
      if (present(region)) then
!
        select case (region)
        case ('f_array')
        case ('pencil')
          b1=m
          b2=m
          c1=n
          c2=n
        case ('default')
          call fatal_error('nan_inform','No such region')
        endselect
!
      endif
!
!  Overwrite with supplied intervals.
!
      if (present(int1)) then  ! x
        a1=int1(1)
        a2=int1(2)
      endif
!
      if (present(int2)) then  ! y
        b1=int2(1)
        b2=int2(2)
      endif
!
      if (present(int3)) then  ! z
        c1=int3(1)
        c2=int3(2)
      endif
!
      if (present(int4)) then  ! variable
        d1=int4(1)
        d2=int4(2)
      endif
!
!  Look for NaN and inf in resulting interval.
!
      do a=a1,a2; do b=b1,b2; do c=c1,c2; do d=d1,d2
        if (notanumber(f(a,b,c,d))) then
          print*,'nan_inform: NaN with message "', msg, &
              '" encountered in the variable ', varname(d)
          print*,'nan_inform: ', varname(d), ' = ', f(a,b,c,d)
          print*,'nan_inform: t, it, itsub   = ', t, it, itsub
          print*,'nan_inform: l, m, n, iproc = ', a, b, c, iproc_world
          print*,'----------------------------'
          if (present(lstop)) then
            if (lstop) call fatal_error('nan_stop','')
          endif
        endif
      enddo; enddo; enddo; enddo
!
    endsubroutine nan_inform
!***********************************************************************
    subroutine parse_bc(bc,bc12)
!
!  Parse boundary conditions, which may be in the form `a' (applies to
!  both `lower' and `upper' boundary) or `a:s' (use `a' for lower,
!  `s' for upper boundary.
!
!  24-jan-02/wolf: coded
!
      character (len=2*bclen+1), dimension(mcom) :: bc
      character (len=bclen), dimension(mcom,2) :: bc12
      integer :: j,isep
!
      intent(in) :: bc
      intent(out) :: bc12
!
      do j=1,mcom
        if (bc(j) == '') then ! will probably never happen due to default='p'
          if (lroot) print*, 'Empty boundary condition No. ', &
              j, 'in (x, y, or z)'
          call fatal_error('parse_bc','')
        endif
        isep = index(bc(j),':')
        if (isep > 0) then
          bc12(j,1) = bc(j)(1:isep-1)
          bc12(j,2) = bc(j)(isep+1:)
        else
          bc12(j,:) = bc(j)(1:bclen)
        endif
      enddo
!
    endsubroutine parse_bc
!***********************************************************************
    subroutine inverse_parse_bc(bc,bc12)
!
!  27-jun-12/joern+dhruba: coded
!
      character (len=2*bclen+1), dimension(mcom) :: bc
      character (len=bclen), dimension(mcom,2) :: bc12
      integer :: j
!
      intent(out) :: bc
      intent(in) :: bc12
!
      do j=1,mcom
        if ((bc12(j,1) == '') .or. (bc12(j,2) == '')) then
! will probably never happen due to default='p'
          if (lroot) print*, 'Empty boundary condition No. ', &
              j, 'in (x, y, or z)'
          call fatal_error('inverse_parse_bc','')
        endif
        bc(j)(1:bclen) = bc12(j,1)
        bc(j)(bclen+1:bclen+1) = ':'
        bc(j)(bclen+2:2*bclen+1)=bc12(j,2)
      enddo
!
    endsubroutine inverse_parse_bc
!***********************************************************************
    subroutine parse_bc_rad(bc,bc1,bc2)
!
!  Parse boundary conditions, which may be in the form `a' (applies to
!  both `lower' and `upper' boundary) or `a:s' (use `a' for lower,
!  `s' for upper boundary.
!
!   6-jul-03/axel: adapted from parse_bc
!
      character (len=2*bclen+1), dimension(3) :: bc
      character (len=bclen), dimension(3) :: bc1,bc2
      integer :: j,isep
!
      intent(in) :: bc
      intent(out) :: bc1,bc2
!
      do j=1,3
        if (bc(j) == '') then ! will probably never happen due to default='p'
          if (lroot) print*, 'Empty boundary condition No. ', &
                 j, 'in (x, y, or z)'
          call fatal_error('parse_bc','')
        endif
        isep = index(bc(j),':')
        if (isep > 0) then
          bc1(j) = bc(j)(1:isep-1)
          bc2(j) = bc(j)(isep+1:)
        else
          bc1(j) = bc(j)(1:bclen)
          bc2(j) = bc(j)(1:bclen)
        endif
      enddo
!
    endsubroutine parse_bc_rad
!***********************************************************************
    subroutine parse_bc_radg(bc,bc1,bc2)
!
!  Parse boundary conditions, which may be in the form `a' (applies to
!  both `lower' and `upper' boundary) or `a:s' (use `a' for lower,
!  `s' for upper boundary.
!
!   6-jul-03/axel: adapted from parse_bc
!
      character (len=2*bclen+1) :: bc
      character (len=bclen) :: bc1,bc2
      integer :: isep
!
      intent(in) :: bc
      intent(out) :: bc1,bc2
!
      if (bc == '') then
        if (lroot) print*, 'Empty boundary condition in (x, y, or z)'
        call fatal_error('parse_bc_radg','')
      endif
      isep = index(bc,':')
      if (isep > 0) then
        bc1 = bc(1:isep-1)
        bc2 = bc(isep+1:)
      else
        bc1 = bc(1:bclen)
        bc2 = bc(1:bclen)
      endif
!
    endsubroutine parse_bc_radg
!***********************************************************************
    subroutine parse_shell(strin,strout)
!
!  Parse string replacing all $XXXX sequences with appropriate
!  values from the environment. Return the parsed result in strout.
!
      use General, only: safe_character_assign
!
      character (len=*) :: strin, strout
      character (len=255) :: envname, chunk !, envvalue
      character (len=1) :: chr
      character (len=64), parameter :: envnamechars = 'abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_-'
      integer :: inptr, inlen, envstart, nameptr
!
      intent(in)    :: strin
      intent(inout)   :: strout
!
      inptr=1
      inlen=len(trim(strin))
      strout=''
!
dlrloop:do
        envstart =index(strin(inptr:inlen),'$')
        if (envstart <= 0) exit;
        chunk = trim(strin(inptr:envstart-1))
        if (envstart > inptr) call safe_character_assign(strout,trim(strout)//trim(chunk))
        inptr = envstart + 1;
        if (inptr > inlen) exit dlrloop
!
        nameptr = inptr
nameloop: do
          chr = trim(strin(nameptr:nameptr))
          if (index(envnamechars,chr) > 0) then
            nameptr=nameptr+1
          else
            exit nameloop
          endif
!
          if (nameptr > inlen) exit nameloop
        enddo nameloop
        if ((nameptr-1) >= inptr) then
         envname=trim(strin(inptr:nameptr-1))
! Commented pending a C replacement
!         call getenv(trim(envname),envvalue)
!         call safe_character_assign(strout,trim(strout)//trim(envvalue))
        endif
!
        inptr=nameptr
        if (inptr>inlen) exit dlrloop
!
      enddo dlrloop
!
      if (inptr <= inlen) then
        chunk = trim(strin(inptr:inlen))
        call safe_character_assign(strout,trim(strout)//trim(chunk))
      endif
!
    endsubroutine parse_shell
!***********************************************************************
    subroutine remove_file(fname)
!
!  Remove a file.
!
!  5-mar-02/wolf: coded
!
      use General, only: file_exists
!
      character (len=*), intent(in) :: fname
!
      logical :: removed
!
      removed = file_exists(fname,DELETE=.true.)
      if (removed .and. (ip<=6)) &
          print*,'remove_file: Removed file <',trim(fname),'>'
!
    endsubroutine remove_file
!***********************************************************************
    function control_file_exists(fname,delete)
!
!  Does the given control file exist in either ./ or ./runtime/ ?
!  If DELETE is true, delete the file after checking for existence.
!
!  26-jul-09/wolf: coded
!
      logical :: control_file_exists
      character (len=*), intent(in) :: fname
      logical, optional, intent(in) :: delete
!
      logical :: ldelete
!
      if (present(delete)) then
        ldelete=delete
      else
        ldelete=.false.
      endif
!
      control_file_exists = parallel_file_exists(trim(fname), ldelete)
      if (.not. control_file_exists) &
          control_file_exists = parallel_file_exists(trim("runtime/"//fname), ldelete)
!
    endfunction control_file_exists
!***********************************************************************
    function read_line_from_file(fname)
!
!  Read the first line from a file; return empty string if file is empty
!
!  4-oct-02/wolf: coded
!
      use General, only : file_exists
!
      character (len=linelen) :: read_line_from_file
      character (len=*) :: fname
!
      integer :: unit=1
      integer :: ierr=0
!
      read_line_from_file=char(0)
      if (.not. lroot) return
!
      if (file_exists(fname)) then
        open(unit,FILE=fname,IOSTAT=ierr)
        read(unit,'(A)',IOSTAT=ierr) read_line_from_file
        close(unit)
        if (ierr /= 0) read_line_from_file=char(0)
      endif
!
    endfunction read_line_from_file
!***********************************************************************
    subroutine get_nseed(nseed)
!
!  Get length of state of random number generator. The current seed can
!  be represented by nseed (4-byte) integers.
!  Different compilers have different lengths:
!    NAG: 1, Compaq: 2, Intel: 47, SGI: 64, NEC: 256
!
      use Mpicomm, only: lroot
      use General, only: random_seed_wrapper
!
      integer :: nseed
!
      call random_seed_wrapper(SIZE=nseed)
!
!  Test whether mseed is large enough for this machine.
!
      if (nseed > mseed) then
        if (lroot) print*, 'This machine requires mseed >= ', nseed, &
            ', but you have only ', mseed
        call fatal_error('get_nseed','Need to increase mseed')
      endif
!
    endsubroutine get_nseed
!***********************************************************************
    subroutine get_where(mask, indices, status)
!
!  Get the indices where mask are .true.
!
!  Note: The pointer argument indices will be reassociated with a newly
!      allocated array.  It is the user's responsibility to pass in a
!      disassociated pointer and deallocate it after use.
!
!  10-feb-15/ccyang: coded.
!
      logical, dimension(:), intent(in) :: mask
      integer, dimension(:), pointer :: indices
      integer, intent(out), optional :: status
!
      integer :: i, n, stat
!
!  Allocate the index array.
!
      allocate(indices(count(mask)), stat=stat)
      error: if (stat /= 0) then
        soft: if (present(status)) then
          status = stat
          return
        endif soft
        call fatal_error('get_where', 'unable to allocate the index array. ')
      endif error
!
!  Scan through the elements of mask.
!
      n = 0
      scan: do i = 1, size(mask)
        true: if (mask(i)) then
          n = n + 1
          indices(n) = i
        endif true
      enddo scan
      if (present(status)) status = 0
!
      endsubroutine get_where
!***********************************************************************
    subroutine write_dx_general(file,x00,y00,z00)
!
!  Write .general file for data explorer (aka DX).
!
!  04-oct-02/wolf: coded
!  08-oct-02/tony: use safe_character_assign() to detect string overflows
!  12-sep-13/mcnallcp: make the endianness detection reflect the platform
!
      use General, only: safe_character_append, date_time_string
!
      real :: x00,y00,z00
      character (len=*) :: file
      character (len=datelen) :: date
      character (len=linelen) :: field='',struct='',type='',dep=''
!
!  This is True for big endian, False of little endian
      logical, parameter :: bigendian = ichar(transfer(1,'a')) == 0
!
      call date_time_string(date)
!
!  Accumulate a few lines.
!
      if (lhydro    ) then
        call safe_character_append(field,  'uu, '       )
        call safe_character_append(struct, '3-vector, ' )
        call safe_character_append(type,   'float, '    )
        call safe_character_append(dep,    'positions, ')
      endif
      if (ldensity  ) then
        call safe_character_append(field,  'lnrho, '    )
        call safe_character_append(struct, 'scalar, '   )
        call safe_character_append(type,   'float, '    )
        call safe_character_append(dep,    'positions, ')
      endif
      if (lentropy  ) then
        call safe_character_append(field,  'ss, '       )
        call safe_character_append(struct, 'scalar, '   )
        call safe_character_append(type,   'float, '    )
        call safe_character_append(dep,    'positions, ')
      endif
      if (lmagnetic ) then
        call safe_character_append(field,  'aa, '       )
        call safe_character_append(struct, '3-vector, ' )
        call safe_character_append(type,   'float, '    )
        call safe_character_append(dep,    'positions, ')
      endif
      if (lradiation) then
        call safe_character_append(field,  'e_rad, ff_rad, '       )
        call safe_character_append(struct, 'scalar, 3-vector, '    )
        call safe_character_append(type,   'float, float, '        )
        call safe_character_append(dep,    'positions, positions, ')
      endif
      if (lpscalar  ) then
        call safe_character_append(field,  'lncc, '     )
        call safe_character_append(struct, 'scalar, '   )
        call safe_character_append(type,   'float, '    )
        call safe_character_append(dep,    'positions, ')
      endif
!
!  Remove trailing comma.
!
      field  = field (1:len(trim(field ))-1)
      struct = struct(1:len(trim(struct))-1)
      type   = type  (1:len(trim(type  ))-1)
      dep    = dep   (1:len(trim(dep   ))-1)
!
!  Now write.
!
      open(1,FILE=file)
!
      write(1,'(A)'  ) '# Creator: The Pencil Code'
      write(1,'(A,A)') '# Date: ', trim(date)
      write(1,'(A,A)') 'file = ', trim(datadir)//'/proc0/var.dat'
      write(1,'(A,I4," x ",I4," x ",I4)') 'grid = ', mx, my, mz
      write(1,'(A,A)') '# NB: hardcoded assumption of ieee'
      if (bigendian) then
        write(1,'(A,A," ",A)') 'format = ', 'msb', 'ieee'
      else
        write(1,'(A,A," ",A)') 'format = ', 'lsb', 'ieee'
      endif
      write(1,'(A,A)') 'header = ', 'bytes 4'
      write(1,'(A,A)') 'interleaving = ', 'record'
      write(1,'(A,A)') 'majority = ', 'column'
      write(1,'(A,A)') 'field = ', trim(field)
      write(1,'(A,A)') 'structure = ', trim(struct)
      write(1,'(A,A)') 'type = ', trim(type)
      write(1,'(A,A)') 'dependency = ', trim(dep)
      write(1,'(A,A,6(", ",1PG12.4))') 'positions = ', &
           'regular, regular, regular', &
           x00, dx, y00, dy, z00, dz
      write(1,'(A)') ''
      write(1,'(A)') 'end'
!
      close(1)
!
    endsubroutine write_dx_general
!***********************************************************************
    subroutine write_zprof(fname,a)
!
!  Writes z-profile to a file.
!
!  10-jul-05/axel: coded
!
      real, dimension(:) :: a
      character (len=*) :: fname

      if (size(a)==mz) then
        call write_prof(fname,z,a,'z')
      else
        call write_prof(fname,z(n1:n2),a,'z')
      endif
 
    endsubroutine write_zprof
!***********************************************************************
    subroutine write_xprof(fname,a)
!
!  Writes x-profile to a file.
!
!  10-jul-05/axel: coded
!
      real, dimension(:) :: a
      character (len=*) :: fname

      if (size(a)==mx) then
        call write_prof(fname,x,a,'x')
      else
        call write_prof(fname,x(l1:l2),a,'x')
      endif
 
    endsubroutine write_xprof
!***********************************************************************
    subroutine write_prof(fname,coor,a,type,lsave_name)
!
!  Writes profile to a file.
!
!  10-jul-05/axel: coded
!
      use General, only: safe_character_assign, loptest
!
      real, dimension(:) :: coor, a
      character (len=*) :: fname
      character :: type
      logical, optional :: lsave_name
!
      integer, parameter :: unit=1
      character (len=fnlen) :: wfile
      integer :: i
!
!  If within a loop, do this only for the first step (indicated by lwrite_prof).
!
      if (lwrite_prof) then
!
!  Write zprofile file.
!
        call safe_character_assign(wfile, &
            trim(directory)//'/'//type//'prof_'//trim(fname)//'.dat')
        open(unit,file=wfile,position='append')
        do i=1,size(coor)
          write(unit,*) coor(i),a(i)
        enddo
        close(unit)
!
!  Add file name to list of f zprofile files if lsave_name *not* set or true.
!
        if (loptest(lsave_name,.true.)) then
          call safe_character_assign(wfile,trim(directory)//'/'//type//'prof_list.dat')
          open(unit,file=wfile,position='append')
          write(unit,*) fname
          close(unit)
        endif
      endif
!
    endsubroutine write_prof
!***********************************************************************
    subroutine read_zprof(fname,a)
!
!  Read z-profile from a file (taken from stratification, for example).
!
!  15-may-15/axel: coded
!
      use General, only: safe_character_assign
!
      real, dimension(:) :: a
      character (len=*) :: fname
!
      integer, parameter :: unit=1
      character (len=fnlen) :: wfile
      real, dimension(size(a)) :: zloc
!
!  Write zprofile file.
!
      call safe_character_assign(wfile, &
          trim(directory)//'/zprof_'//trim(fname)//'.dat')
      open(unit,file=wfile)
      do n=1,size(a)
        read(unit,*) zloc(n),a(n)
      enddo
      close(unit)
!
!  Should we check that zloc == z ?
!
    endsubroutine read_zprof
!***********************************************************************
    subroutine remove_zprof
!
!  Remove z-profile file.
!
!  10-jul-05/axel: coded
!
      use General, only: safe_character_assign
      use IO, only: lcollective_IO
!
      character (len=120) :: fname,wfile,listfile
      integer :: ierr, unit=2
!
      if (lcollective_IO .and. .not. lroot) return
!
!  Do this only for the first step.
!
      call safe_character_assign(listfile,trim(directory)//'/zprof_list.dat')
!
!  Read list of file and remove them one by one.
!
      open(unit,file=listfile,status='old',iostat=ierr)
      if (ierr /= 0) return
      do while ((it <= nt) .and. (ierr == 0))
        read(unit,*,iostat=ierr) fname
        if (ierr == 0) then
          call safe_character_assign(wfile, &
              trim(directory)//'/zprof_'//trim(fname)//'.dat')
          call remove_file(wfile)
        endif
      enddo
      close(unit)
!
!  Now delete this listfile altogether.
!
      call remove_file(listfile)
!
    endsubroutine remove_zprof
!***********************************************************************
    subroutine blob(ampl,f,i,radius,xblob,yblob,zblob)
!
!  Single blob.
!
!  27-jul-02/axel: coded
!
      integer :: i
      real, dimension (mx,my,mz,mfarray) :: f
      real, optional :: xblob,yblob,zblob
      real :: ampl,radius,x01=0.,y01=0.,z01=0.
!
!  Single  blob.
!
      if (present(xblob)) x01=xblob
      if (present(yblob)) y01=yblob
      if (present(zblob)) z01=zblob
      if (ampl==0) then
        if (lroot) print*,'ampl=0 in blob'
      else
        if (lroot.and.ip<14) print*,'blob: variable i,ampl=',i,ampl
        f(:,:,:,i)=f(:,:,:,i)+ampl*(&
           spread(spread(exp(-((x-x01)/radius)**2),2,my),3,mz)&
          *spread(spread(exp(-((y-y01)/radius)**2),1,mx),3,mz)&
          *spread(spread(exp(-((z-z01)/radius)**2),1,mx),2,my))
      endif
!
    endsubroutine blob
!***********************************************************************
    recursive function hypergeometric2F1(a,b,c,z,tol) result (hyp2F1)
!
!  DOCUMENT ME!!!!!!
!
      real, intent(in) :: a,b,c,z,tol
      real :: hyp2F1
      real :: fac
      integer :: n
!
      real :: aa,bb,cc
!
      aa=a; bb=b; cc=c
!
      fac=1
      hyp2F1=fac
      n=1
!
      if (z<=0.5) then
!
        do while (fac>tol)
          fac=fac*aa*bb*z/(cc*n)
          hyp2F1=hyp2F1+fac
          aa=aa+1
          bb=bb+1
          cc=cc+1
          n=n+1
        enddo
!
      else
!
        !!!!!!!! only valid for mu=-1 !!!!!!!!
        !hyp2F1=2*hypergeometric2F1(aa,bb,aa+bb-cc+1,1-z,tol)-sqrt(1-z)* &
               !2*hypergeometric2F1(cc-aa,cc-bb,cc-aa-bb+1,1-z,tol)
        hyp2F1=(gamma_function(cc)*gamma_function(cc-aa-bb))/ &
               (gamma_function(cc-aa)*gamma_function(cc-bb))* &
               hypergeometric2F1(aa,bb,aa+bb-cc+1,1-z,tol) &
              +(1-z)**(cc-aa-bb)* &
               (gamma_function(cc)*gamma_function(aa+bb-cc))/ &
               (gamma_function(aa)*gamma_function(bb))* &
               hypergeometric2F1(cc-aa,cc-bb,cc-aa-bb+1,1-z,tol)
!
      endif
!
    endfunction hypergeometric2F1
!***********************************************************************
    recursive function pi_function(x) result(pi_func)
!
!  Calculates the Pi-function using rational approximation.
!
!    Pi(x) = Gamma(x+1) = x!
!
!  Coefficients were determined using maple's minimax() function.
!
!  9-jun-04/tobi+wolf: coded
!
      real, intent(in) :: x
      real :: pi_func
      integer, parameter :: order=7
      real, dimension(order) :: coeff1,coeff2
      real :: enum,denom
      integer :: i
!
      coeff1=(/0.66761295020790986D00, &
               0.36946093910826145D00, &
               0.18669829780572704D00, &
               4.8801451277274492D-2, &
               1.36528684153155468D-2, &
               1.7488042503123817D-3, &
               3.6032044608268575D-4/)
!
      coeff2=(/0.66761295020791116D00, &
               0.754817592058897962D00, &
              -3.7915754844972276D-2, &
              -0.11379619871302534D00, &
               1.5035521280605477D-2, &
               3.1375176929984225D-3, &
              -5.5599617153443518D-4/)
!
      if (x>1) then
!
        pi_func=x*pi_function(x-1)
!
      elseif (x<0) then
!
        if (abs(x+1)<=epsilon(x)) then
          pi_func=pi_function(x+1)/epsilon(x)
        else
          pi_func=pi_function(x+1)/(x+1)
        endif
!
      else
!
        enum=coeff1(order)
        do i=order-1,1,-1
          enum=enum*x+coeff1(i)
        enddo
        denom=coeff2(order)
        do i=order-1,1,-1
          denom=denom*x+coeff2(i)
        enddo
        pi_func=enum/denom
!
      endif
!
    endfunction pi_function
!***********************************************************************
    function gamma_function(x)
!
!  Calculates the Gamma-function as
!
!    Gamma(x) = Pi(x-1)
!
!  9-jun-04/tobi+wolf: coded
!
      real, intent(in) :: x
      real :: gamma_function
!
      gamma_function=pi_function(x-1)
!
    endfunction gamma_function
!***********************************************************************
    subroutine tensor_diffusion_coef(gecr,ecr_ij,bij,bb,vKperp,vKpara,rhs,llog,gvKperp,gvKpara)
!
!  Calculates tensor diffusion with variable tensor (or constant tensor).
!  Calculates parts common to both variable and constant tensor first.
!
!  Note: ecr=lnecr in the below comment
!
!  Write diffusion tensor as K_ij = Kpara*ni*nj + (Kperp-Kpara)*del_ij.
!
!  vKperp*del2ecr + d_i(vKperp)d_i(ecr) + (vKpara-vKperp) d_i(n_i*n_j*d_j ecr)
!      + n_i*n_j*d_i(ecr)d_j(vKpara-vKperp)
!
!  = vKperp*del2ecr + gKperp.gecr + (vKpara-vKperp) (H.G + ni*nj*Gij)
!      + ni*nj*Gi*(vKpara_j - vKperp_j),
!  where H_i = (nj bij - 2 ni nj nk bk,j)/|b| and vKperp, vKpara are variable
!  diffusion coefficients.
!
!  Calculates (K.gecr).gecr
!  =  vKperp(gecr.gecr) + (vKpara-vKperp)*Gi(ni*nj*Gj)
!
!  Adds both parts into decr/dt.
!
!  10-oct-03/axel: adapted from pscalar
!  30-nov-03/snod: adapted from tensor_diff without variable diffusion
!  04-dec-03/snod: converted for evolution of lnecr (=ecr)
!   9-apr-04/axel: adapted for general purpose tensor diffusion
!  25-jun-05/bing:
!
      real, dimension (nx,3,3) :: ecr_ij,bij
      real, dimension (nx,3) :: gecr,bb,bunit,hhh,gvKperp1,gvKpara1,tmpv
      real, dimension (nx) :: abs_b,b1,del2ecr,gecr2,vKperp,vKpara
      real, dimension (nx) :: hhh2,quenchfactor,rhs,tmp,tmpi,tmpj,tmpk
      real :: limiter_tensordiff=3.
      integer :: i,j,k
      logical, optional :: llog
      real, optional, dimension (nx,3) :: gvKperp,gvKpara
!
      intent(in) :: bb,bij,gecr,ecr_ij
      intent(out) :: rhs
!
!  Calculate unit vector of bb.
!
!     call dot2_mn(bb,abs_b,PRECISE_SQRT=.true.)
      call dot2_mn(bb,abs_b,FAST_SQRT=.true.)
      b1=1./max(tini,abs_b)
      call multsv_mn(b1,bb,bunit)
!
!  Calculate first H_i.
!
      del2ecr=0.
      do i=1,3
        del2ecr=del2ecr+ecr_ij(:,i,i)
        hhh(:,i)=0.
        do j=1,3
          tmpj(:)=0.
          do k=1,3
            tmpj(:)=tmpj(:)-2.*bunit(:,k)*bij(:,k,j)
          enddo
          hhh(:,i)=hhh(:,i)+bunit(:,j)*(bij(:,i,j)+bunit(:,i)*tmpj(:))
        enddo
      enddo
      call multsv_mn(b1,hhh,tmpv)
!
!  Limit the length of H such that dxmin*H < 1, so we also multiply
!  by 1/sqrt(1.+dxmin^2*H^2).
!  And dot H with ecr gradient.
!
!     call dot2_mn(tmpv,hhh2,PRECISE_SQRT=.true.)
      call dot2_mn(tmpv,hhh2,FAST_SQRT=.true.)
      quenchfactor=1./max(1.,limiter_tensordiff*hhh2*dxmax)
      call multsv_mn(quenchfactor,tmpv,hhh)
      call dot_mn(hhh,gecr,tmp)
!
!  Dot Hessian matrix of ecr with bi*bj, and add into tmp.
!
      call multmv_mn(ecr_ij,bunit,hhh)
      call dot_mn(hhh,bunit,tmpj)
      tmp = tmp+tmpj
!
!  Calculate (Gi*ni)^2 needed for lnecr form; also add into tmp.
!
      gecr2=0.
      if (present(llog)) then
        call dot_mn(gecr,bunit,tmpi)
        tmp=tmp+tmpi**2
!
!  Calculate gecr2 - needed for lnecr form.
!
        call dot2_mn(gecr,gecr2)
      endif
!
!  If variable tensor, add extra terms and add result into decr/dt.
!
!  Set gvKpara, gvKperp.
!
     if (present(gvKpara)) then; gvKpara1=gvKpara; else; gvKpara1=0.; endif
     if (present(gvKperp)) then; gvKperp1=gvKperp; else; gvKperp1=0.; endif
!
!  Put d_i ecr d_i vKperp into tmpj.
!
      call dot_mn(gvKperp1,gecr,tmpj)
!
!  Nonuniform conductivities, add terms into tmpj.
!
      call dot(bunit,gvKpara1-gvKperp1,tmpi)
      call dot(bunit,gecr,tmpk)
      tmpj = tmpj+tmpi*tmpk
!
!  Calculate rhs.
!
      rhs=vKperp*(del2ecr+gecr2) + (vKpara-vKperp)*tmp + tmpj
!
    endsubroutine tensor_diffusion_coef
!***********************************************************************
    subroutine max_for_dt_nx_nx(f,maxf)
!
!  Like maxf = max(f,max), unless we have chosen to manipulate data
!  before taking the maximum value. Designed for calculation of time step,
!  where one may want to exclude certain regions, etc.
!
!  Would be nicer as an (assumed-size) array-valued function (as a plug-in
!  replacement for max), but this can be more than 2 times slower (NEC
!  SX-5, compared to about 15% slower with Intel F95) than a subroutine
!  call according to tests.
!
!  30-jan-04/wolf: coded
!
      real, dimension(nx) :: maxf,f
!
      intent(in)    :: f
      intent(inout) :: maxf
!
      maxf = max(f,maxf)
!
    endsubroutine max_for_dt_nx_nx
!***********************************************************************
    subroutine max_for_dt_1_nx(f,maxf)
!
!  Like max_for_dt_n_n, but with a different signature of argument shapes.
!
!  30-jan-04/wolf: coded
!
      real, dimension(nx) :: maxf
      real                :: f
!
      intent(in)    :: f
      intent(inout) :: maxf
!
      maxf = max(f,maxf)
!
    endsubroutine max_for_dt_1_nx
!***********************************************************************
    subroutine max_for_dt_1_1_1_nx(f1,f2,f3,maxf)
!
!  Like max_for_dt_n_n, but with a different signature of argument shapes.
!
!  30-jan-04/wolf: coded
!
      real, dimension(nx) :: maxf
      real                :: f1,f2,f3
!
      intent(in)    :: f1,f2,f3
      intent(inout) :: maxf
!
      maxf = max(f1,f2,f3,maxf)
!
    endsubroutine max_for_dt_1_1_1_nx
!***********************************************************************
    function pencil_multiply1(s,v)
!
!  The `*' operator may be extended through this function to allow
!  elementwise multiplication of a `pencil-scalar' with a `pencil-vector'.
!
!   6-Sep-05/tobi: coded
!
      real, dimension(nx), intent(in) :: s
      real, dimension(nx,3), intent(in) :: v
      real, dimension(nx,3) :: pencil_multiply1
!
      integer :: i
!
      do i=1,3; pencil_multiply1(:,i) = s(:) * v(:,i); enddo
!
    endfunction pencil_multiply1
!***********************************************************************
    function pencil_multiply2(v,s)
!
!  The `*' operator may be extended through this function to allow
!  elementwise multiplication of a `pencil-scalar' with a `pencil-vector'.
!
!   6-Sep-05/tobi: coded
!
      real, dimension(nx,3), intent(in) :: v
      real, dimension(nx), intent(in) :: s
      real, dimension(nx,3) :: pencil_multiply2
!
      integer :: i
!
      do i=1,3; pencil_multiply2(:,i) = v(:,i) * s(:); enddo
!
    endfunction pencil_multiply2
!***********************************************************************
    function pencil_add1(s,v)
!
!  The `+' operator may be extended through this function to allow
!  elementwise addition of a `pencil-scalar' to a `pencil-vector'.
!
!   6-Sep-05/tobi: coded
!
      real, dimension(nx), intent(in) :: s
      real, dimension(nx,3), intent(in) :: v
      real, dimension(nx,3) :: pencil_add1
!
      integer :: i
!
      do i=1,3; pencil_add1(:,i) = s(:) + v(:,i); enddo
!
    endfunction pencil_add1
!***********************************************************************
    function pencil_add2(v,s)
!
!  The `+' operator may be extended through this function to allow
!  elementwise addition of a `pencil-scalar' to a `pencil-vector'.
!
!   6-Sep-05/tobi: coded
!
      real, dimension(nx,3), intent(in) :: v
      real, dimension(nx), intent(in) :: s
      real, dimension(nx,3) :: pencil_add2
!
      integer :: i
!
      do i=1,3; pencil_add2(:,i) = v(:,i) + s(:); enddo
!
    endfunction pencil_add2
!***********************************************************************
    function pencil_divide1(s,v)
!
!  The `/' operator may be extended through this function to allow
!  elementwise division of a `pencil-scalar' by a `pencil-vector'.
!
!   6-Sep-05/tobi: coded
!
      real, dimension(nx), intent(in) :: s
      real, dimension(nx,3), intent(in) :: v
      real, dimension(nx,3) :: pencil_divide1
!
      integer :: i
!
      do i=1,3; pencil_divide1(:,i) = s(:) / v(:,i); enddo
!
    endfunction pencil_divide1
!***********************************************************************
    function pencil_divide2(v,s)
!
!  The `/' operator may be extended through this function to allow
!  elementwise division of a `pencil-vector' by a `pencil-scalar'.
!
!   6-Sep-05/tobi: coded
!
      real, dimension(nx,3), intent(in) :: v
      real, dimension(nx), intent(in) :: s
      real, dimension(nx,3) :: pencil_divide2
!
      integer :: i
!
      do i=1,3; pencil_divide2(:,i) = v(:,i) / s(:); enddo
!
    endfunction pencil_divide2
!***********************************************************************
    function pencil_subtract1(s,v)
!
!  The `-' operator may be extended through this function to allow
!  elementwise subtraction of a `pencil-vector' from a `pencil-scalar'.
!
!   6-Sep-05/tobi: coded
!
      real, dimension(nx), intent(in) :: s
      real, dimension(nx,3), intent(in) :: v
      real, dimension(nx,3) :: pencil_subtract1
!
      integer :: i
!
      do i=1,3; pencil_subtract1(:,i) = s(:) - v(:,i); enddo
!
    endfunction pencil_subtract1
!***********************************************************************
    function pencil_subtract2(v,s)
!
!  The `-' operator may be extended through this function to allow
!  elementwise subtraction of a `pencil-scalar' from a `pencil-vector'.
!
!   6-Sep-05/tobi: coded
!
      real, dimension(nx,3), intent(in) :: v
      real, dimension(nx), intent(in) :: s
      real, dimension(nx,3) :: pencil_subtract2
!
      integer :: i
!
      do i=1,3; pencil_subtract2(:,i) = v(:,i) - s(:); enddo
!
    endfunction pencil_subtract2
!***********************************************************************
    elemental real function one_minus_exp(x)
!
!  Returns 1 - exp(-x).
!
!  03-jan-15/ccyang: coded.
!
      real, intent(in) :: x
!
      if (x * x > epsilon(1.0)) then
        one_minus_exp = 1.0 - exp(-x)
      else
        one_minus_exp = x * (1.0 - 0.5 * x)
      endif
!
    endfunction one_minus_exp
!***********************************************************************
    elemental real function erfunc(x)
!
!  Error function from Numerical Recipes.
!
!  15-Jan-2007/dintrans: coded
!
      real, intent(in) :: x
!
      real :: dumerfc, t, z
!
      z = abs(x)
      t = 1.0 / ( 1.0 + 0.5 * z )
!
      dumerfc =  t * exp(-z * z - 1.26551223 + t *        &
          ( 1.00002368 + t * ( 0.37409196 + t *           &
          ( 0.09678418 + t * (-0.18628806 + t *           &
          ( 0.27886807 + t * (-1.13520398 + t *           &
          ( 1.48851587 + t * (-0.82215223 + t * 0.17087277 )))))))))
!
      if (x<0.0) dumerfc = 2.0 - dumerfc
      erfunc = 1.0 - dumerfc
!
    endfunction erfunc
!***********************************************************************
    subroutine power_law_mn(const,dist,power_law_index,output,xref)
!
!  General distance power law initial conditions.
!
!  24-feb-05/wlad: coded
!   4-jul-07/wlad: generalized for any power law case
!
      real, dimension(:) :: dist,output
      real :: const,power_law_index
      real, optional :: xref
!
      intent(in)  :: const,power_law_index
      intent(out) :: output
!
      if (present(xref)) dist=dist/xref
!
      if (rsmooth==0.) then
        output = const*dist**(-power_law_index)
      else
        output = const*(dist**2+rsmooth**2)**(-.5*power_law_index)
      endif
!
    endsubroutine power_law_mn
!***********************************************************************
    subroutine power_law_pt(const,dist,power_law_index,output,xref)
!
!  General distance power law initial conditions.
!
!  24-feb-05/wlad: coded
!   4-jul-07/wlad: generalized for any power law case
!
      real :: dist,output
      real :: const,power_law_index
      real, optional :: xref
!
      intent(in)  :: const,power_law_index
      intent(out) :: output
!
      if (present(xref)) dist=dist/xref
!
      if (rsmooth==0.) then
        output = const*dist**(-power_law_index)
      else
        output = const*(dist**2+rsmooth**2)**(-.5*power_law_index)
      endif
!
    endsubroutine power_law_pt
!***********************************************************************
    subroutine get_radial_distance(rrmn,rcylmn,e1_,e2_,e3_)
!
!  Calculate distance and its cylindrical projection for different
!  coordinate systems.
!
!  e1, e2, and e3 are the positions in the respective coordinate systems
!
!  15-mar-07/wlad : coded
!
      real, dimension(:),intent(out) :: rrmn,rcylmn
      real, dimension(size(rrmn,1)) :: xc
      real, intent(in), optional :: e1_,e2_,e3_
      real :: e1,e2,e3
      integer :: tmp
      logical :: lorigin
!
!  Check if we are dealing with distance from the origin.
!
      tmp=0 ; lorigin=.false.
      if (present(e1_)) then;e1=e1_;tmp=tmp+1;else;e1=0.;endif
      if (present(e2_)) then;e2=e2_;tmp=tmp+1;else;e2=0.;endif
      if (present(e3_)) then;e3=e3_;tmp=tmp+1;else;e3=0.;endif
      if (tmp==0) lorigin=.true.
!
!  Check if this array has size nx or mx.
!
      select case (size(rrmn))
      case (mx)
        xc=x
      case (nx)
        xc=x(l1:l2)
      case default
        print*,'get_radial_distance: '//&
             'the array has dimension=',size(rrmn),' is that correct?'
        call fatal_error('get_radial_distance','')
      endselect
!
!  Calculate the coordinate-free distance relative to the position (e1,e2,e3).
!
      if (lorigin) then
        if (coord_system=='cartesian') then
          rcylmn=sqrt(xc**2+y(m)**2)+tini
          rrmn  =sqrt(    rcylmn**2 +z(n)**2)
        elseif (coord_system=='cylindric') then
          rcylmn= xc                      +tini
          rrmn  =sqrt(  rcylmn**2+z(n)**2)
        elseif (coord_system=='spherical') then
          rcylmn=     xc*sinth(m)+tini
          rrmn  =     xc         +tini
        endif
      else
        if (coord_system=='cartesian') then
          rcylmn=sqrt((xc-e1)**2+(y(m)-e2)**2)+tini
          rrmn  =sqrt(       rcylmn**2+(z(n)-e3)**2)
        elseif (coord_system=='cylindric') then
          rcylmn=sqrt(xc**2+e1**2 - 2*xc*e1*cos(y(m)-e2))+tini
          rrmn  =sqrt(rcylmn**2+(z(n)-e3)**2)
        elseif (coord_system=='spherical') then
          rcylmn=sqrt((xc*sinth(m))**2 + (e1*sin(e2))**2 - &
               2*xc*e1*costh(m)*cos(e2))+tini
          rrmn  =sqrt(xc**2 + e1**2 - 2*xc*e1*&
               (costh(m)*cos(e2)+sinth(m)*sin(e2)*cos(z(n)-e3)))+tini
        endif
      endif
!
    endsubroutine get_radial_distance
!***********************************************************************
    function interp1(r,fr,nr,r0,ldescending)
!
!  20-dec-07/dintrans: coded
!  Note: if ldescending=T, then input arrays r and fr are in descending
!  order and we first reverse them.
!
      integer :: nr,istop,i,i1,i2
      real, dimension (nr) :: r,fr,r1,fr1
      real    :: r0,interp1
      logical, optional :: ldescending
!
      if (present(ldescending)) then
        if (ldescending) then
          r1=r(nr:1:-1)
          fr1=fr(nr:1:-1)
        else
          r1=r
          fr1=fr
        endif
      else
        r1=r
        fr1=fr
      endif
!
      if (r0 == r1(1)) then
        interp1=fr1(1)
        return
      elseif (r0 > r1(nr)) then
        interp1=fr1(nr)
        return
      else
        istop=0 ; i=1
        do while (istop /= 1)
          if (r1(i) >= r0) istop=1
          i=i+1
        enddo
        i1=i-2 ; i2=i-1
        interp1=(fr1(i1)*(r1(i2)-r0)+fr1(i2)*(r0-r1(i1)))/(r1(i2)-r1(i1))
      endif
!
    endfunction interp1
!***********************************************************************
    subroutine bspline_basis(k, x, b)
!
!  Computes the values of the non-zero B-spline basis functions
!  B_{i,k}(j+x) for i = j-k+1, j-k+2, ..., j.  The knot sequence {t_i)
!  is assumed to be infinite and be integers, i.e., t_i = i for all
!  integer i.
!
!  28-jul-15/ccyang: coded.
!
!  Input Arguments
!      k   Number of knot spans for each basis function, which has order
!          (k-1).
!      x   A number in [0,1).
!  Output Argument
!      b   An array of k elements, where b(i) = B_{j-k+i,k}(j+x).
!
      integer, intent(in) :: k
      real, intent(in) :: x
      real, dimension(k), intent(out) :: b
!
      integer :: i, j
!
!  Work up the order column by column.
!
      b = 0.0
      b(1) = 1.0
      order: do j = 2, k
        b(j) = x * b(j-1)
        do i = j - 1, 2, -1
          b(i) = (x - real(i-j)) * b(i-1) + (real(i) - x) * b(i)
        end do
        b(1) = (1.0 - x) * b(1)
        b(1:j) = b(1:j) / real(j-1)
      enddo order
!
    endsubroutine bspline_basis
!***********************************************************************
    subroutine bspline_interpolation(n, k, f, a, indx, shift)
!
!  Uses the B-spline interpolation to periodically shift a regular array
!  of data nodes.
!
!  31-jul-15/ccyang: coded.
!
!  Input/Output Argument
!      f   An array of node data; interpolated after shift on return.
!  Input Arguments
!      n   Number of nodes.
!      k   Number of knot spans for each basis function, which has order
!          (k-1).
!      a   Preconditioned by bspline_precondition() and then LU
!          decomposed by ludcmp().
!      indx
!          Index permutations returned by ludcmp().
!      shift
!          Shift in unit of array index.
!
      integer, intent(in) :: n, k
      real, dimension(n), intent(inout) :: f
      real, dimension(n,n), intent(in) :: a
      integer, dimension(n), intent(in) :: indx
      real, intent(in) :: shift
!
      real, dimension(:), allocatable, save :: bk
      integer :: k_old = -1
      real :: shift_old = 0.0
!
      real, dimension(n) :: b, c
      integer :: i, j
!
!  Solve for the coefficients for the B-spline basis functions.
!
      c = f
      call lubksb(a, indx, c)
!
!  Find the values of the basis functions at the interpolation points.
!
      j = ceiling(shift - 0.5)
      basis: if (k /= k_old .or. shift /= shift_old) then
        alloc: if (k /= k_old) then
          if (allocated(bk)) deallocate(bk)
          allocate(bk(k))
        endif alloc
        call bspline_basis(k, 0.5 - shift + real(j), bk)
        k_old = k
        shift_old = shift
      endif basis
      b = 0.0
      b(1:k) = bk
!
!  Make the interpolation.
!
      j = k + j
      do i = 1, n
        f(i) = sum(cshift(b, j - i) * c)
      enddo
!
    endsubroutine bspline_interpolation
!***********************************************************************
    subroutine bspline_precondition(n, k, a)
!
!  Sets up the linear system for the coefficients of the B-spline basis
!  functions, assuming periodic boundary conditions.  The linear system
!  reads
!
!      A x = f,
!
!  where
!          [ B_{1,k}(0.5)   B_{2,k}(0.5)   B_{3,k}(0.5)   ... B_{n,k}(0.5)   ]
!      A = [ B_{1,k}(1.5)   B_{2,k}(1.5)   B_{3,k}(1.5)   ... B_{n,k}(1.5)   ],
!          [ ...            ...            ...            ... ...            ]
!          [ B_{1,k}(n-0.5) B_{2,k}{n-0.5) B_{3,k}(n-0.5) ... B_{n,k}(n-0.5) ]
!
!      x = [ alpha_1, alpha_2, alpha_3, ..., alpha_n ]^T,
!
!  is the coefficients,
!
!      f = [ f_1, f_2, f_3, ..., f_n ],
!
!  is the node data.  The B-spline interpolation is then given by
!
!      f(x) = sum_i B_{i,k)(x)
!
!  with
!
!      f(i-0.5) = f_i, i = 1, 2, ..., n.
!
!  28-jul-15/ccyang: coded.
!
!  Input Arguments
!      n   Number of nodes.
!      k   Number of knot spans for each B-spline basis function, which
!          has order (k-1).
!  Output Argument
!      a   The square matrix for the linear system.
!
      integer, intent(in) :: n, k
      real, dimension(n,n), intent(out) :: a
!
      real, dimension(n) :: b
      integer :: i
!
!  Find the values of the B-spline basis functions.
!
      b = 0.0
      call bspline_basis(k, 0.5, b(1:k))
!
!  Cyclically assign the basis functions into the square matrix.
!
      do i = 1, n
        a(i,:) = cshift(b, k - i)
      enddo
!
    endsubroutine bspline_precondition
!***********************************************************************
    subroutine ludcmp(a,indx)
!
!  25-jun-09/rplasson: coded (adapted from numerical recipe)
!
!  Computes the LU decomposition of the matrix a.
!  The result is placed in the matrix a.
!  The row permutations are returned in indx.
!
      real, dimension(:,:), intent(INOUT) :: a
      integer, dimension(:), intent(OUT) :: indx
      real, dimension(size(a,1)) :: vv,swap
      integer :: j,n,imax
      integer, dimension(1) :: tmp
!
      n=size(a,1)
      if (n /= size(a,2)) call fatal_error('ludcmp','non square matrix')
      if (n /= size(indx)) call fatal_error('ludcmp','bad dimension for indx')
      vv=maxval(abs(a),dim=2)
      if (any(vv == 0.0)) call fatal_error('ludcmp','singular matrix')
      vv=1.0/vv
      do j=1,n
        tmp=maxloc(vv(j:n)*abs(a(j:n,j)))
        imax=(j-1)+tmp(1)
        if (j /= imax) then
          swap=a(imax,:)
          a(imax,:)=a(j,:)
          a(j,:)=swap
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if (a(j,j) == 0.0) a(j,j)=tiny(0.)
        a(j+1:n,j)=a(j+1:n,j)/a(j,j)
        a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-spread(a(j+1:n,j),dim=2,ncopies=(n-j)) * &
            spread(a(j,j+1:n),dim=1,ncopies=(n-j))
      enddo
!
    endsubroutine ludcmp
!***********************************************************************
    subroutine lubksb(a,indx,b)
!
!  25-jun-09/rplasson: coded (adapted from numerical recipe)
!
!  Solves the equation A.X=B .
!  'a' must contain the LU decomposition of matrix A obtained by ludcmp.
!  'indx' is the permutation vector obtained by ludcmp.
!  'b' contains B, and returns the solution vector X.
!
      real, dimension(:,:), intent(IN) :: a
      integer, dimension(:), intent(IN) :: indx
      real, dimension(:), intent(INOUT) :: b
      integer :: i,n,ii,ll
      real :: summ
!
      n=size(a,1)
      if (n /= size(a,2)) call fatal_error('lubksb','non square matrix')
      if (n /= size(indx)) call fatal_error('lubksb','bad dimension for indx')
      ii=0
      do i=1,n
        ll=indx(i)
        summ=b(ll)
        b(ll)=b(i)
        if (ii /= 0) then
          summ=summ-dot_product(a(i,ii:i-1),b(ii:i-1))
        else if (summ /= 0.0) then
          ii=i
        endif
        b(i)=summ
      enddo
      do i=n,1,-1
        b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i)
      enddo
!
    endsubroutine lubksb
!***********************************************************************
    subroutine gij_psi(psif,ee,g)
!
!  Calculate gradient of a scalar field multiplied by a constant vector),
!  return matrix.
!
!  31-jul-07/dhruba: adapted from gij
!
      use Deriv, only: der
!
      real, dimension (mx,my,mz) :: psif
      real, dimension (3) :: ee
      real, dimension (nx,3,3) :: g
      real, dimension (nx) :: tmp
      integer :: i,j
!
      intent(in) :: psif
      intent(out) :: g
!
      do i=1,3; do j=1,3
          call der(psif*ee(i),tmp,j)
          g(:,i,j) = tmp
      enddo; enddo
!
    endsubroutine gij_psi
!***********************************************************************
    subroutine gij_psi_etc(psif,ee,aa,aij,Bij,del2,graddiv)
!
!  Calculate B_i,j = eps_ikl A_l,jk and A_l,kk .
!
!  1-aug-07/dhruba : adapted from gij_etc
!
      use Deriv, only: der2,derij
!
      real, dimension (mx,my,mz), intent (in) :: psif
      real, dimension(3), intent(in) :: ee
      real, dimension (nx,3,3), intent (out) :: bij
      real, dimension (nx,3,3), intent (in), optional :: aij
      real, dimension (nx,3), intent (out), optional :: del2,graddiv
      real, dimension (nx,3), intent (in), optional :: aa
!
!  Locally used variables.
!
      real, dimension (nx,3,3,3) :: d2A
      real, dimension (nx) :: tmp
      integer :: i,j
!
!  Calculate all mixed and non-mixed second derivatives
!  of the vector potential (A_k,ij).
!
!  Do not calculate both d^2 A/(dx dy) and d^2 A/(d dx)
!  (This wasn't spotted by me but by a guy from SGI...)
!  Note: for non-cartesian coordinates there are different correction terms,
!  see below.
!
      do i=1,3
        do j=1,3
          call der2(psif*ee(i),tmp,j); d2A(:,j,j,i)=tmp
        enddo
!!DHRUBA
        call derij(psif*ee(i),tmp,2,3); d2A(:,2,3,i)=tmp; d2A(:,3,2,i)=tmp
        call derij(psif*ee(i),tmp,3,1); d2A(:,3,1,i)=tmp; d2A(:,1,3,i)=tmp
        call derij(psif*ee(i),tmp,1,2); d2A(:,1,2,i)=tmp; d2A(:,2,1,i)=tmp
      enddo
!
!  Corrections for spherical polars from swapping mixed derivatives:
!  Psi_{,theta^ r^} = Psi_{,r^ theta^} - Psi_{,\theta^}/r
!  Psi_{,phi^ r^} = Psi_{,r^ phi^} - Psi_{,\phi^}/r
!  Psi_{,phi^ theta^} = Psi_{,theta^ phi^} - Psi_{,\phi^}*r^{-1}*cot(theta)
!
      if (lspherical_coords) then
        do i=1,3
          d2A(:,2,1,i)=d2A(:,2,1,i)-aij(:,i,2)*r1_mn
          d2A(:,3,1,i)=d2A(:,3,1,i)-aij(:,i,3)*r1_mn
          d2A(:,3,2,i)=d2A(:,3,2,i)-aij(:,i,3)*r1_mn*cotth(m)
        enddo
      endif
!
!  For cylindrical, only
!  Psi_{,phi^ pom^} = Psi_{,pom^ phi^} - Psi_{,\phi^}/pom .
!
      if (lcylindrical_coords) then
         do i=1,3
            d2A(:,2,1,i)=d2A(:,2,1,i)-aij(:,i,2)*rcyl_mn1
         enddo
      endif
!
!  Calculate b_i,j = eps_ikl A_l,kj, as well as optionally,
!  del2_i = A_i,jj and graddiv_i = A_j,ji .
!
      bij(:,1,:)=d2A(:,2,:,3)-d2A(:,3,:,2)
      bij(:,2,:)=d2A(:,3,:,1)-d2A(:,1,:,3)
      bij(:,3,:)=d2A(:,1,:,2)-d2A(:,2,:,1)
!
!  Corrections for spherical coordinates.
!
      if (lspherical_coords) then
        bij(:,3,2)=bij(:,3,2)+aij(:,2,2)*r1_mn
        bij(:,2,3)=bij(:,2,3)-aij(:,3,3)*r1_mn
        bij(:,1,3)=bij(:,1,3)+aij(:,3,3)*r1_mn*cotth(m)
        bij(:,3,1)=bij(:,3,1)+aij(:,2,1)*r1_mn         -aa(:,2)*r2_mn
        bij(:,2,1)=bij(:,2,1)-aij(:,3,1)*r1_mn         +aa(:,3)*r2_mn
        bij(:,1,2)=bij(:,1,2)+aij(:,3,2)*r1_mn*cotth(m)-aa(:,3)*r2_mn*sin2th(m)
      endif
!
!  Corrections for cylindrical coordinates.
!
      if (lcylindrical_coords) then
        bij(:,3,2)=bij(:,3,2)+ aij(:,2,2)*r1_mn
        bij(:,3,1)=bij(:,3,1)+(aij(:,2,1)+aij(:,1,2))*rcyl_mn1-aa(:,2)*rcyl_mn2
      endif
!
!  Calculate del2 and graddiv, if requested.
!
      if (present(graddiv)) then
!--     graddiv(:,:)=d2A(:,:,1,1)+d2A(:,:,2,2)+d2A(:,:,3,3)
        graddiv(:,:)=d2A(:,1,:,1)+d2A(:,2,:,2)+d2A(:,3,:,3)
        if (lspherical_coords) then
          graddiv(:,1)=graddiv(:,1)+aij(:,1,1)*r1_mn*2+ &
             aij(:,2,1)*r1_mn*cotth(m)-aa(:,2)*r2_mn*cotth(m)-aa(:,1)*r2_mn*2
          graddiv(:,2)=graddiv(:,2)+aij(:,1,2)*r1_mn*2+ &
             aij(:,2,2)*r1_mn*cotth(m)-aa(:,2)*r2_mn*sin2th(m)
          graddiv(:,3)=graddiv(:,3)+aij(:,1,3)*r1_mn*2+ &
             aij(:,2,3)*r1_mn*cotth(m)
        endif
      endif
!
      if (present(del2)) then
        del2(:,:)=d2A(:,1,1,:)+d2A(:,2,2,:)+d2A(:,3,3,:)
        if (lspherical_coords.and.present(aij).and.present(aa)) then
          del2(:,1)= del2(:,1)+&
            r1_mn*(2.*(aij(:,1,1)-aij(:,2,2)-aij(:,3,3)&
            -r1_mn*aa(:,1)-cotth(m)*r1_mn*aa(:,2) ) &
            +cotth(m)*aij(:,1,2) )
          del2(:,2)=del2(:,2)+&
            r1_mn*(2.*(aij(:,2,1)-cotth(m)*aij(:,3,3)&
            +aij(:,1,2) )&
            +cotth(m)*aij(:,2,2)-r1_mn*sin2th(m)*aa(:,2) )
          del2(:,3)=del2(:,3)+&
            r1_mn*(2.*(aij(:,3,1)+aij(:,1,3)&
            +cotth(m)*aij(:,2,3) ) &
            +cotth(m)*aij(:,3,2)-r1_mn*sin2th(m)*aa(:,3) )
        else
        endif
        if (lcylindrical_coords) call fatal_error('gij_etc', &
            'use del2=graddiv-curlcurl for cylindrical coords')
      endif
!
    endsubroutine gij_psi_etc
!***********************************************************************
    logical function location_in_proc(xpos,ypos,zpos,lpos,mpos,npos)
!
! finds out if a points belongs to the location in the processor.
! If yes the also returns the nearest grid location of this point.
!
!  -dec-10/dhruba: coded
!  28-dec-10/MR: changed into function
!  13-jan-11/MR: made dependent on dimensionality,
!                irrelevant position indices set to 1
!
      real,    intent(in)  :: xpos,ypos,zpos
      integer, intent(out) :: lpos,mpos,npos
!
      logical :: linx,liny,linz
!
      lpos=1 ; mpos=1 ; npos=1
      linx=.true. ; liny=.true. ; linz=.true.
!
      if (dimensionality>0) then
        call xlocation(xpos,lpos,linx)
!
        if (dimensionality>1) then
          call ylocation(ypos,mpos,liny)
!
          if (dimensionality>2) &
            call zlocation(zpos,npos,linz)
        endif
      endif
!
      location_in_proc = linx.and.liny.and.linz
!
    endfunction location_in_proc
!***********************************************************************
    subroutine xlocation(xpos,ixpos,lproc)
!
!  If xpos lies within this processor, then lproc=T and xpos=x(ixpos).
!  Otherwise lproc=F and ixpos=1.
!
!  18-nov-06/axel: coded
!  14-oct-08/ccyang: use half-closed interval and include the top-most plane
!  03-dec-09/MR: moved here from module slices
!  16-dec-10/joern: adapted from zlocation
!  18-jan-11/axel+joern: choose between l and l+1 depending on sign of x
!
      real :: xpos
      integer :: ixpos,l
      logical :: lproc
!
!  Use ixpos for testing if position is found.
!
      ixpos=-1
!
!  Run through all x positions until we hit the right interval.
!  If the right interval is found, jump out of the loop.
!  To make the output symmetric about zero, we must swap upper
!  and lower boundaries depending on the sign of x.
!
      do l=l1,l2
        if (x(l)<=xpos.and.x(l+1)>xpos) then
          if (x(l)<0) then
            ixpos=l+1
          else
            ixpos=l
          endif
          lproc=.true.
          exit
        endif
      enddo
!
!  If nothing is found, we set lproc=.false. and and put ixpos=1
!
      if (ixpos==-1) then
        ixpos=1
        lproc=.false.
      endif
!
    endsubroutine xlocation
!***********************************************************************
    subroutine ylocation(ypos,iypos,lproc)
!
!  If ypos lies within this processor, then lproc=T and ypos=y(iypos).
!  Otherwise lproc=F and iypos=1.
!
!  18-nov-06/axel: coded
!  14-oct-08/ccyang: use half-closed interval and include the top-most plane
!  03-dec-09/MR: adapted from module slices
!  16-dec-10/joern: took it from zlocation
!  18-jan-11/axel+joern: choose between m and m+1 depending on sign of y
!
      real :: ypos
      integer :: iypos,m
      logical :: lproc
!
!  Use iypos for testing if position is found.
!
      iypos=-1
!
!  Run through all y positions until we hit the right interval.
!  If the right interval is found, jump out of the loop.
!  To make the output symmetric about zero, we must swap upper
!  and lower boundaries depending on the sign of y.
!
      do m=m1,m2
        if (y(m)<=ypos.and.y(m+1)>ypos) then
          if (y(m)<0) then
            iypos=m+1
          else
            iypos=m
          endif
          lproc=.true.
          exit
        endif
      enddo
!
!  If nothing is found, we set lproc=.false. and and put iypos=1
!
      if (iypos==-1) then
        iypos=1
        lproc=.false.
      endif
!
    endsubroutine ylocation
!***********************************************************************
    subroutine zlocation(zpos,izpos,lproc)
!
!  If zpos lies within this processor, then lproc=T and zpos=z(izpos).
!  Otherwise lproc=F and izpos=1.
!
!  18-nov-06/axel: coded
!  14-oct-08/ccyang: use half-closed interval and include the top-most plane
!  03-dec-09/MR: moved here from module slices
!  18-jan-11/axel+joern: choose between n and n+1 depending on sign of z
!
      real :: zpos
      integer :: izpos,n
      logical :: lproc
!
!  Use izpos for testing if position is found.
!
      izpos=-1
!
!  Run through all z positions until we hit the right interval.
!  If the right interval is found, jump out of the loop.
!  To make the output symmetric about zero, we must swap upper
!  and lower boundaries depending on the sign of z.
!
      do n=n1,n2
        if (z(n)<=zpos.and.z(n+1)>zpos) then
          if (z(n)<0) then
            izpos=n+1
          else
            izpos=n
          endif
!AB: old  izpos=n
          lproc=.true.
          exit
        endif
      enddo
!
!  If nothing is found, we set lproc=.false. and and put izpos=1
!
      if (izpos==-1) then
        izpos=1
        lproc=.false.
      endif
!
    endsubroutine zlocation
!***********************************************************************
    subroutine position(ind,ip,ngrid,ind_loc,flag)
!
!  Determines local position ind_loc with respect to processor ip corresponding to global position ind
!  if grid has global extent ngrid. flag is set if ind_loc lies within the local range.
!  On return, ind_loc is corrected for ghost zones, thus can be used to index the f array.
!  ind_loc and flag are not altered if ind <= 0.
!
!  21-apr-15/MR: coded
!
      integer, intent(in) :: ind,ip,ngrid
      integer, intent(out):: ind_loc
      logical, intent(out):: flag
     
      if (ind>0) then
        ind_loc=ind-ip*ngrid
        if (ind_loc>=1.and.ind_loc<=ngrid) then
          flag=.true.
          ind_loc=ind_loc+nghost
        else
          flag=.false.
        endif
      endif

    endsubroutine position
!***********************************************************************
    subroutine fourier_single_mode(arr,idims,k,idir,amps,l2nd)
!
!  No parallelization in x allowed here
!
!  08-dec-10/MR: coded
!
      use mpicomm, only: mpireduce_sum
!
      implicit none
!
      integer, dimension(2)             , intent(in)  :: idims
      real, dimension(idims(1),idims(2)), intent(in)  :: arr
      real                              , intent(in)  :: k
      integer                           , intent(in)  :: idir
      real   , dimension(2,*)           , intent(out) :: amps
      logical                , optional , intent(in)  :: l2nd
!
      integer :: n,i,idim
      real, dimension(:)  , allocatable :: cg,sg
      real, dimension(:,:), allocatable :: buffer
      real :: fac
      logical :: l2ndl
!
      if (present(l2nd)) then
        l2ndl=l2nd
      else
        l2ndl=.false.
      endif
!
      if (l2ndl) then
        idim=idims(1)
      else
        idim=idims(2)
      endif
!
      select case (idir)
      case (1)    ; n=nxgrid
      case (2)    ; n=ny
      case (3)    ; n=nz
      case default; n=nxgrid
      end select
!
      if (idims(1)+idims(2)-idim/=n) then
        amps(:,1:idim)=0.
        return
      endif
!
      allocate(cg(n),sg(n))
!
      select case (idir)
      case (1)    ; cg=cos(k*xgrid)   ;sg=sin(k*xgrid)   ; fac=dx
      case (2)    ; cg=cos(k*y(m1:m2));sg=sin(k*y(m1:m2)); fac=dy
      case (3)    ; cg=cos(k*z(n1:n2));sg=sin(k*z(n1:n2)); fac=dz
      case default; cg=cos(k*xgrid)   ;sg=sin(k*xgrid)   ; fac=dx
      end select
!
      do i=1,idim
        if (l2ndl) then
          amps(:,i) = fac*(/ sum(arr(i,:)*cg), sum(arr(i,:)*sg) /)
        else
          amps(:,i) = fac*(/ sum(arr(:,i)*cg), sum(arr(:,i)*sg) /)
        endif
      enddo
!
      if (ncpus>1) then
        allocate(buffer(2,idim))
        select case (idir)
        case (2)
          if (nprocy>1) then
            call mpireduce_sum(amps,buffer,(/2,idim/),idir=2)
            if (ipy==0) amps(:,1:idim)=buffer      !result is in root of y-beams
          endif
        case (3)
          if (nprocz>1) then
            call mpireduce_sum(amps,buffer,(/2,idim/),idir=3)
            if (ipz==0) amps(:,1:idim)=buffer      !result is in root of z-beams
          endif
        end select
        deallocate(buffer)
      endif
!
      deallocate(cg,sg)
!
    endsubroutine fourier_single_mode
!***********************************************************************
    subroutine register_report_aux(name, index, ind_aux1, ind_aux2, ind_aux3, communicated)
!
!  Registers aux variable named 'name' if not already registered
!  (i.e. if index==0). Variable is scalar if ind_aux1,ind_aux2,
!  ind_aux3 are missing, vector with number of components equal to
!  number of present ind_aux* parameters. Index of variable and its
!  components (if any) are returned in index,ind_aux1,ind_aux2,ind_aux3
!  If already registered: outputs indices in index.pro
!
!  13-jan-11/MR: coded
!  29-may-14/ccyang: add optional argument communicated
!
      use FArrayManager, only: farray_register_auxiliary
!
      implicit none
!
      integer,           intent(inout) :: index
      integer, optional, intent(inout) :: ind_aux1,ind_aux2,ind_aux3
      character (LEN=*), intent(in)    :: name
      logical, intent(in), optional :: communicated
!
      integer   :: vec
      character :: ch
      character (LEN=len(name)-2) :: tail
      character (LEN=32) :: equalto
!
      vec=-1
!
      if ( present(ind_aux1) ) then
        vec=1
        if ( present(ind_aux2) ) then
          if ( present(ind_aux3) ) then
            vec=3
          else
            vec=2
          endif
        endif
      endif
!
      if (index==0) then
!
        if (present(communicated)) then
          call farray_register_auxiliary(trim(name), index, vector=abs(vec), communicated=communicated)
        else
          call farray_register_auxiliary(trim(name), index, vector=abs(vec))
        endif
!
        if (vec>=1) then
          ind_aux1=index
          if (vec>=2) then
            ind_aux2=index+1
            if (vec==3) ind_aux3=index+2
          endif
        endif
!
      endif
!
      if (index/=0.and.lroot) then
!
        print*, 'register_report_aux: i'//trim(name)//' = ', index
        open(3,file=trim(datadir)//'/index.pro', POSITION='append')
        write(3,*) 'i'//trim(name)//'=',index
!
        if ( vec>=1 ) then
          tail='='
          if ( name(1:1)==name(2:2) ) then
            ch = name(1:1)
            if ( len_trim(name)>2 ) tail = trim(name(3:))//'='
          endif
!
          equalto='='
          write(3,*) 'i'//ch//'x'//trim(equalto),ind_aux1
          if ( vec>=2 ) then
            write(3,*) 'i'//ch//'y'//trim(equalto),ind_aux2
            if ( vec==3 ) write(3,*) 'i'//ch//'z'//trim(equalto),ind_aux3
          endif
        endif
!
        close(3)
!
      endif
!
    endsubroutine register_report_aux
!***********************************************************************
    subroutine unit_vector(bb,bb_hat)
!
!  Compute the unit vector for any given vector bb.
!  Tries to avoid division by zero.
!  Taken from http://nuclear.llnl.gov/CNP/apt/apt/aptvunb.html.
!
!  18-oct-11/bing: copied from bb_unitvec_shock in magnetic.f90
!
      real, dimension(nx,3) :: bb,bb_hat,bb2
      real, dimension(nx) :: a2,aerr2,bb_len
      integer :: j
      real :: tol
!
      intent(in) :: bb
      intent(out) :: bb_hat
!
!  Truncate small components to zero.
!
      tol=sqrt(tini)
!
      bb2 = bb**2
!
      aerr2 = tol * max(sum(bb2,2),1.)
!
      do j=1,3
        where (bb2(:,j) < aerr2)
          bb_hat(:,j) = 0.
        elsewhere
          bb_hat(:,j) = bb(:,j)
        endwhere
      enddo
!
!  Get unit vector.
!
      bb_len = sqrt(sum(bb_hat**2,2))
!
      do j=1,3; bb_hat(:,j) = bb_hat(:,j)/(bb_len+tini); enddo
!
!  Check if length is between 0. and 1.
!
      call dot2(bb_hat,a2)
!
      if (maxval(a2) > 1.+1e-6) &
          call fatal_error('unit_vector:','has not the length 1')
!
    endsubroutine unit_vector
!***********************************************************************
    subroutine doupwind(f,k,uu,ugradf,mask)
!
!  Calculates upwind correction, works incrementally on ugradf
!
!  26-mar-12/MR: outsourced from routines u_dot_grad_mat, u_dot_grad_scl, u_dot_grad_scl_alt
!   9-apr-12/MR: optional parameter plus added
!  12-apr-12/MR: optional parameter modified
!
      use Deriv, only: der6, deri_3d_inds
!
      real, dimension (mx,my,mz,mfarray), intent(IN)    :: f
      integer                                           :: k
      real, dimension (nx,3),             intent(IN)    :: uu
      real, dimension (nx),               intent(INOUT) :: ugradf
      integer,                            intent(IN), optional :: mask
!
      real, dimension (nx,3) :: del6f
      integer                :: ii,msk
      integer, dimension(nx) :: indxs
!
      msk=0
      if (present(mask)) then
        if ( mask>=1 .and. mask <=3 ) msk=mask
      endif
!
      do ii=1,3
!
        if (ii==msk) then
          del6f(:,ii) = 0.
        else
!
          if ( lequidist(ii) ) then
            call der6(f,k,del6f(1,ii),ii,UPWIND=.true.)
          else
            where( uu(:,ii)>=0 )
              indxs = 7
            elsewhere
              indxs = 8
            endwhere
            call deri_3d_inds(f(1,1,1,k),del6f(1,ii),indxs,ii,lnometric=.true.)
          endif
!
          del6f(:,ii) = abs(uu(:,ii))*del6f(:,ii)
!
        endif
      enddo
!
      if (lcylindrical_coords) &
        del6f(:,2) = rcyl_mn1*del6f(:,2)
!
      if (lspherical_coords) then
        del6f(:,2) = r1_mn*del6f(:,2)
        del6f(:,3) = r1_mn*sin1th(m)*del6f(:,3)
      endif
!
      if (msk>0) then
        ugradf = ugradf+sum(del6f,2)
      else
        ugradf = ugradf-sum(del6f,2)
      endif
!
    endsubroutine doupwind
!***********************************************************************
    subroutine global_mean(f,inda,mean,indep,lexp)
!
!  Calculate global mean for a (several) field(s) selected by the index inda
!  (the index range inda - indep) in f.
!
!  15-oct-12/MR: adapted from remove_mean
!  22-aug-13/MR: rewritten into subroutine because of internal compiler error of gfortran
!  25-aug-13/MR: removed allocatable attribute from mean to adhere to f95
!
      use Mpicomm, only: mpiallreduce_sum
      use General, only: ioptest, loptest
!
      real, dimension (mx,my,mz,*), intent(in)           :: f
      integer,                      intent(in)           :: inda
      real, dimension(inda:),       intent(out)          :: mean
      integer,                      intent(in), optional :: indep
      logical,                      intent(in), optional :: lexp
!
      real, allocatable, dimension(:) :: mean_tmp
      integer :: j, inde
      logical :: lexpl
!
      inde = ioptest(indep,inda)
!
      allocate(mean_tmp(inda:inde))
!
!  initialize mean
!
      mean = 0.0
!
!  Compute mean for each field.
!
      lexpl = loptest(lexp)
      do j=inda,inde
        if (lexpl) then
          mean(j) = mean(j) + sum(exp(f(l1:l2,m1:m2,n1:n2,j)))
        else
          mean(j) = mean(j) + sum(f(l1:l2,m1:m2,n1:n2,j))
        endif
      enddo
!
!  Compute total mean for all processors
!
      call mpiallreduce_sum(mean,mean_tmp,inde-inda+1)
      mean = mean_tmp/nwgrid
!
    endsubroutine global_mean
!***********************************************************************
    subroutine remove_mean(f,inda,indep)
!
!  Substract mean from a (several) field(s) selected by the index inda
!  (the index range inda - indep) in f.
!
!  08-may-12/MR: adapted from remove_mean_flow
!
      use Mpicomm, only: mpiallreduce_sum
      use General, only: ioptest
!
      real, dimension (mx,my,mz,*), intent (inout)        :: f
      integer,                      intent (in)           :: inda
      integer,                      intent (in), optional :: indep
!
      real, allocatable, dimension(:) :: mean, mean_tmp
      integer :: m,n,j, inde
!
      inde = ioptest(indep,inda)
      allocate( mean(inda:inde), mean_tmp(inda:inde) )
!
!  initialize mean
!
      mean = 0.0
!
!  Go through all pencils.
!
      do n = n1,n2
      do m = m1,m2
!
!  Compute mean for each field.
!
        do j=inda,inde
          mean(j) = mean(j) + sum(f(l1:l2,m,n,j))
        enddo
      enddo
      enddo
!
!  Compute total mean for all processors
!
      call mpiallreduce_sum(mean/nwgrid,mean_tmp,inde-inda+1)
      mean = mean_tmp
!
!  Go through all pencils and subtract out the mean
!  separately for each field.
!
      do n = n1,n2
      do m = m1,m2
        do j=inda,inde
          f(l1:l2,m,n,j) = f(l1:l2,m,n,j) - mean(j)
        enddo
      enddo
      enddo
!
      if (lroot.and.ip<6) print*,'remove_mean: mean=',mean
!
      deallocate( mean, mean_tmp )
!
    endsubroutine remove_mean
!***********************************************************************
    subroutine insert_carray( array, insert, index, leng )
!
! inserts string vector insert in string vector array at position index; array is assumed
! to be settled until position leng; leng is updated.
! Oobs: routine does not check whether there is enough free space in
!       array for the insertion.
!
! 15-feb-2013/MR: parameter leng_insert removed (now derived from insert)
! 21-aug-2013/MR: moved from Testflow
!
      character(LEN=*), dimension(*) :: array
      character(LEN=*), dimension(:) :: insert
      integer                        :: index, leng, leng_insert, i
!
      intent(in)    :: index, insert
      intent(inout) :: leng, array
!
!  insert only if position index is meaningful
!
      if ( index>0.and.index<=leng+1 ) then
!
        leng_insert = size(insert)
        do i=leng,index,-1
          array(i+leng_insert) = array(i)
        enddo
!
        array(index:index+leng_insert-1) = insert
!
        leng = leng+leng_insert

      endif
!
    endsubroutine insert_carray
!***********************************************************************
    subroutine insert_rarray( array, insert, index, leng )
!
! 21-aug-2013/MR: derived from insert_carray for vectors of reals.
!
      integer            :: index, leng, leng_insert, i
      real, dimension(*) :: array
      real, dimension(:) :: insert
!
      intent(in)    :: index, insert
      intent(inout) :: leng, array
!
!  insert only if position index is meaningful
!
      if ( index>0.and.index<=leng+1 ) then
!
        leng_insert = size(insert)
        do i=leng,index,-1
          array(i+leng_insert) = array(i)
        enddo
!
        array(index:index+leng_insert-1) = insert
!
        leng = leng+leng_insert

      endif
!
    endsubroutine insert_rarray
!***********************************************************************
    subroutine insert_carray_mult( array, insert, mult, index, leng )
!
! inserts string insert mult times in string vector array at position index; array is assumed
! to be settled until position leng; leng is updated.
! Oobs: routine does not check whether there is enough free space in
!       array for the insertion.
!
! 15-feb-2013/MR: derived from insert_carray
!
      character(LEN=*), dimension(*) :: array
      character(LEN=*)               :: insert
      integer                        :: index, leng, mult, i
!
      intent(in)    :: index, insert, mult
      intent(inout) :: leng, array
!
!  insert only if position index is meaningful
!
      if ( index>0.and.index<=leng+1 ) then
!
        do i=leng,index,-1
          array(i+mult) = array(i)
        enddo
!
        array(index:index+mult-1) = insert
!
        leng = leng+mult
!
      endif
    endsubroutine insert_carray_mult
!***********************************************************************
    real function find_max_fvec(f,iv)
!
!  Find the maximum of the modulus of a vector field.
!
!  19-aug-2011/ccyang: coded
!  12-sep-2013/MR: outsourced from hydro
!
      use Mpicomm, only: mpiallreduce_max, MPI_COMM_WORLD
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      integer, intent(in) :: iv
!
      real :: umax1
!
!  Find the maximum.
!
      umax1 = sqrt(maxval(  f(l1:l2,m1:m2,n1:n2,iv  )**2 &
                          + f(l1:l2,m1:m2,n1:n2,iv+1)**2 &
                          + f(l1:l2,m1:m2,n1:n2,iv+2)**2  ))
      call mpiallreduce_max(umax1, find_max_fvec, MPI_COMM_WORLD)
!
    endfunction find_max_fvec
!***********************************************************************
    real function find_rms_fvec(f, iv)
!
!  Find the root-mean-square of the modulus of a vector field.
!
!  25-aug-2015/ccyang: coded
!
      use Mpicomm, only: mpiallreduce_sum
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      integer, intent(in) :: iv
!
      real :: s1, s
!
!  Find the rms.
!
      s1 = sum(f(l1:l2,m1:m2,n1:n2,iv:iv+2)**2)
      call mpiallreduce_sum(s1, s)
      find_rms_fvec = sqrt(s / real(nwgrid))
!
    endfunction find_rms_fvec
!***********************************************************************
    function find_xyrms_fvec(f, iv) result(rms)
!
!  Find the rms magnitude of a vector field in each z.
!
!  29-jun-14/ccyang: coded
!
      use Mpicomm, only: mpiallreduce_sum
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      integer, intent(in) :: iv
      real, dimension(nz) :: rms
!
      real, dimension(nz) :: rms_loc
      integer :: k
!
      forall(k = n1:n2) rms_loc(k-nghost) = sum(f(l1:l2,m1:m2,k,iv:iv+2)**2)
      call mpiallreduce_sum(rms_loc, rms, nz)
      rms = sqrt(rms / nxygrid)
!
    endfunction find_xyrms_fvec
!***********************************************************************
    subroutine finalize_aver_3D(nproc,idir,arrm)
!
!  Finalize calculation of an average of a 3D array arrm, by MPI communication
!  in direction idir, nproc - number of processors in this direction,
!  direction can be a beam or a layer
!
!  12-sep-2013/MR: coded
!
      use Mpicomm, only: mpiallreduce_sum
      use Messages, only: fatal_error
      use Yinyang_mpi, only: reduce_zsum
!
      integer,                intent(IN)   :: nproc,idir
      real, dimension(:,:,:), intent(INOUT):: arrm

      real, dimension(:,:,:), allocatable :: temp
      integer, dimension(3) :: sz
!
!  communicate in direction idir
!
        if (nproc>1) then
!
          sz=(/size(arrm,1),size(arrm,2),size(arrm,3)/)
          allocate(temp(sz(1),sz(2),sz(3)))
          if (lyinyang) then
            if (idir/=3) & 
              call fatal_error('finalize_aver_3D', &
                               'Not implemented for other than phi direction on Yin-Yang grid')
            !call reduce_zsum(arrm,temp)
            !call mpibcast_z_yy(temp)
          else
            call mpiallreduce_sum(arrm,temp,sz,idir=idir)
            arrm=temp
          endif
!
        endif
!
    endsubroutine finalize_aver_3D
!***********************************************************************
    subroutine finalize_aver_1D(nproc,idir,arrm)
!
!  Finalize calculation of an average of a 1D array arrm, by MPI communication
!  in direction idir, nproc - number of processors in this direction,
!  direction can be a beam or a layer
!
!  12-sep-2013/MR: coded
!
      use Mpicomm, only: mpiallreduce_sum
!
      integer,            intent(IN)   :: nproc,idir
      real, dimension(:), intent(INOUT):: arrm

      real, dimension(:), allocatable :: temp
      integer :: sz
!
!  communicate over direction idir
!
        if (nproc>1) then
!
          sz = size(arrm)
          allocate(temp(sz))
          if (lyinyang.and.idir==3) then
            !call zaverages_xy
          else
            call mpiallreduce_sum(arrm,temp,sz,idir=idir)
          endif
          arrm=temp
!
        endif
!
    endsubroutine finalize_aver_1D
!***********************************************************************
    subroutine finalize_aver_2D(nproc,idir,arrm)
!
!  Finalize calculation of an average of a 2D array arrm, by MPI communication
!  in direction idir, nproc - number of processors in this direction,
!  direction can be a beam or a layer
!
!  12-sep-2013/MR: coded
!
      use Mpicomm, only: mpiallreduce_sum
!
      integer,              intent(IN)   :: nproc,idir
      real, dimension(:,:), intent(INOUT):: arrm

      real, dimension(:,:), allocatable :: temp
      integer, dimension(2) :: sz
!
!  communicate over direction idir
!
        if (nproc>1) then
!
          sz=(/size(arrm,1),size(arrm,2)/)
          allocate(temp(sz(1),sz(2)))
          call mpiallreduce_sum(arrm,temp,sz,idir=idir)
          arrm=temp
!
        endif
!
    endsubroutine finalize_aver_2D
!***********************************************************************
    subroutine finalize_aver_4D(nproc,idir,arrm)
!
!  Finalize calculation of an average of a 4D array arrm, by MPI communication
!  in direction idir, nproc - number of processors in this direction,
!  direction can be a beam or a layer
!
!  12-sep-2013/MR: coded
!
      use Mpicomm, only: mpiallreduce_sum, mpireduce_sum   !!, mpicomm
!
!      include 'mpif.h'
!
      integer,                  intent(IN)   :: nproc,idir
      real, dimension(:,:,:,:), intent(INOUT):: arrm

      real, dimension(:,:,:,:), allocatable :: temp

      integer, dimension(4) :: sz
!
!  communicate over direction idir
!
        if (nproc>1) then
!
          sz=(/size(arrm,1),size(arrm,2),size(arrm,3),size(arrm,4)/)
          allocate(temp(sz(1),sz(2),sz(3),sz(4)))
          call mpiallreduce_sum(arrm,temp,sz,idir=idir)
!
          !!call MPI_BCAST(temp, product(sz), MPI_REAL, iprocx+nprocx*ipy, &
          !!mpicomm(idir),mpierr)
          arrm=temp
!
        endif
!
    endsubroutine finalize_aver_4D
!***********************************************************************
    subroutine fseek_pos(unit, rec_len, num_rec, reference)
!
!  Seeks to a given position in an opened file relative to a reference point.
!  If reference=0, this is relative to the beginning of the file,
!  if reference=1, this is relative to the current position in the file,
!  and if reference=2, this is relative to the end of the file.
!  'rec_len' and 'num_rec' are referring to a record length and a given number
!  of records that one likes to seek, boths must be representable in integer.
!  If 'num_rec' is negative, seeking is done backwards.
!
!  20-Feb-2012/Bourdin.KIS: coded
!
      use General, only: itoa

      integer, intent(in) :: unit
      integer(kind=8) :: rec_len, num_rec
      integer, intent(in) :: reference
!
      integer :: i, num, len
!
      if (num_rec < 0) then
        num_rec = -num_rec
        rec_len = -rec_len
      endif
!
      ! all numbers must be representable as integer(kind=4)
      len = rec_len
      num = num_rec
      if (len /= rec_len) call fatal_error ('fseek_pos on unit '//trim (itoa (unit)), &
          "rec_len is not representable as integer(kind=4).", .true.)
      if (num /= num_rec) call fatal_error ('fseek_pos on unit '//trim (itoa (unit)), &
          "num_rec is not representable as integer(kind=4).", .true.)
!
! WORKAROUND:
! Even though the ifort manual states that ifort would be able to fseek
! with an 64-bit integer argument, this is NOT working!
! Therefore, we have to iterate the fseek with a 32-bit integer to be save.
! Note: gfortran would be able to seek with a 64-bit integer value, though.
! (20-Feb-2012, Bourdin.KIS)
!
      !!!call fseek (unit, rec_len, reference)
      if (num >= 2) then
        do i = 2, num
          !!!call fseek (unit, rec_len, 1)
        enddo
      endif
!
    endsubroutine fseek_pos
!***********************************************************************
    function parallel_count_lines(file,ignore_comments)
!
!  Determines in parallel the number of lines in a file.
!
!  Returns:
!  * Integer containing the number of lines in a given file
!  * -1 on error
!
!  23-mar-10/Bourdin.KIS: implemented
!  26-aug-13/MR: optional parameter ignore_comments added
!  28-May-2015/Bourdin.KIS: reworked
!
      use General, only: count_lines
      use Mpicomm, only: mpibcast_int, MPI_COMM_WORLD

      integer :: parallel_count_lines
      character (len=*), intent(in) :: file
      logical, optional, intent(in) :: ignore_comments
!
      if (lroot) parallel_count_lines = count_lines(file,ignore_comments)
      call mpibcast_int(parallel_count_lines,comm=MPI_COMM_WORLD)
!
    endfunction parallel_count_lines
!***********************************************************************
    function parallel_file_exists(file, delete)
!
!  Determines in parallel if a given file exists.
!  If delete is true, deletes the file.
!
!  Returns:
!  * Integer containing the number of lines in a given file
!  * -1 on error
!
!  23-mar-10/Bourdin.KIS: implemented
!
      use General, only: file_exists, loptest
      use Mpicomm, only: mpibcast_logical, MPI_COMM_WORLD

      character (len=*) :: file
      logical :: parallel_file_exists
      logical, optional :: delete
!
      ! Let the root node do the dirty work
      if (lroot) parallel_file_exists = file_exists(file,loptest(delete))
!
      call mpibcast_logical(parallel_file_exists,comm=MPI_COMM_WORLD)
!
    endfunction parallel_file_exists
!***********************************************************************
    subroutine read_namelist(reader,name,lactive)
!
!  Encapsulates reading of pars + error handling.
!
!  31-oct-13/MR: coded
!  16-dec-13/MR: handling of optional ierr corrected
!  18-dec-13/MR: changed handling of ierr to avoid compiler trouble
!  19-dec-13/MR: changed ierr into logical lierr to avoid compiler trouble
!  11-jul-14/MR: end-of-file caught to avoid program crash when a namelist is missing
!  18-aug-15/PABourdin: reworked to simplify code and display all errors at once
!  19-aug-15/PABourdin: renamed from 'read_pars' to 'read_namelist'
!
      use File_io, only: parallel_rewind, find_namelist
      use General, only: loptest, itoa
      use Messages, only: warning
!
      interface
        subroutine reader(iostat)
          integer, intent(out) :: iostat
        endsubroutine reader
      endinterface
!
      character(len=*), intent(in) :: name
      logical, optional, intent(in) :: lactive
!
      integer :: ierr
      logical :: found
      character(len=5) :: type, suffix
!
      if (.not. loptest (lactive, .true.)) return
!
      if (lstart .or. lparam_nml) then
        type = 'init'
      else
        type = 'run'
      endif
      if (name /= '') type = '_'//type
      suffix = '_pars'
      if (name == 'initial_condition_pars') then
        type = ''
        suffix = ''
      endif
!
      !if (.not. find_namelist (trim(name)//trim(type)//trim(suffix))) then
      call find_namelist (trim(name)//trim(type)//trim(suffix),found)

      if (.not. found) then
        if (.not. lparam_nml) lnamelist_error = .true.
        return
      endif
!
      ierr = 0 ! G95 complains 'ierr' is used but not set, even though 'reader' has intent(out).
      call reader(ierr)
!
      if (ierr /= 0) then
        lnamelist_error = .true.
        if (lroot) then
          if (ierr == -1) then
            call warning ('read_namelist', 'namelist "'//trim(name)//trim(type)//trim(suffix)//'" is missing!')
          else
            call warning ('read_namelist', 'namelist "'//trim(name)//trim(type)//trim(suffix)// &
                          '" has an error ('//trim(itoa(ierr))//')!')
          endif
        endif
      endif
!
      call parallel_rewind
!
    endsubroutine read_namelist
!***********************************************************************
    subroutine calc_diffusive_flux(diffs,c_char,islope_limiter,h_slope_limited,flux)
!
!  Calculates diffusive flux from variable differences diffs acc. to Eqs. (6)-(10) in Rempel (2014).
!  c_char is vector of characteristic velocities.
!
!  23-sep-15/MR,joern,fred,petri: coded
!
      use General, only: notanumber

      real, dimension(:),intent(in) :: diffs,c_char
      real,              intent(in) :: h_slope_limited
      character(LEN=*),  intent(in) :: islope_limiter
      real, dimension(:),intent(out):: flux

      real, dimension(size(diffs)-1) :: slope
      real, dimension(size(diffs)-2) :: phi
      integer :: len

      len=size(diffs)

      call slope_limiter(diffs(2:),diffs(:len-1),slope,islope_limiter)
      flux = diffs(2:len-1) - 0.5*(slope(2:) + slope(1:len-2))        ! = u_r - u_l

      call diff_flux(h_slope_limited, diffs(2:len-1), flux, phi)
      flux = -0.5*c_char*phi*flux
if (notanumber(c_char)) then
   print*, 'CALC_DIFFUSIVE_FLUX: c_char=', len
   stop
endif 

    endsubroutine calc_diffusive_flux
!***********************************************************************
    subroutine slope_limiter(diff_right,diff_left,limited,type)
!
!  Returns limited slope in parameter limited, see Rempel (2014). 
!  Presently only limiter minmod is implemented. 
!
!  25-sep-15/MR,joern: coded
!  27-jan-16/MR: converted into non-elemental subroutine
!
      real, dimension(:), intent(IN) :: diff_left, diff_right
      real, dimension(:), intent(OUT):: limited
      character(LEN=*),   intent(IN) :: type
      integer :: len, ii

      len=size(diff_right)

      select case (type)
      case ('minmod-mine')
        where( sign(1.,diff_left)*sign(1.,diff_right)>0)
          limited = sign(min( 2.*abs(diff_left), 2.*abs(diff_right), &
                             0.5*abs(diff_right+diff_left) ), diff_left )
        elsewhere
          limited = 0.
        endwhere
      case ('minmod') 
         do ii=1,len
           if (diff_left(ii)*diff_right(ii)>0) then
             if (diff_left(ii)>0) limited(ii)=min(2.*diff_left(ii),2.*diff_right(ii), &
                                                  0.5*(diff_right(ii)+diff_left(ii)))
             if (diff_left(ii)<0) limited(ii)=max(2.*diff_left(ii),2.*diff_right(ii), &
                                                  0.5*(diff_right(ii)+diff_left(ii)))
           else
             limited(ii)=0
           endif
         enddo
      case ('superbee')
        limited=0.
        call fatal_error('slope_limiter','limiter not implemented')
      case ('')
        limited=0.
        call fatal_error('slope_limiter','limiter not implemented')
      case default
        limited=0.
        call fatal_error('slope_limiter','limiter not implemented')
      end select
    
    endsubroutine slope_limiter
!***********************************************************************
     subroutine diff_flux(h, diff_right, diff_lr, phi)
!
!  Calculates diffusive flux for one coordinate direction from u_i+1-u_i and u_r-u_l 
!  and returns it in parameter phi, see Rempel (2014). 
!
!  25-sep-15/MR,joern: coded
!  27-jan-16/MR: converted into non-elemental subroutine, because of malcompilation by gcc version 4.6.3
!
      real, intent(IN) :: h
      real, dimension(:), intent(IN)  :: diff_lr, diff_right
      real, dimension(:), intent(OUT) :: phi

      where (diff_right*diff_lr>0.)
        phi=max(0.,1.+h*(diff_lr/diff_right-1.))
      elsewhere
        phi=0.
      endwhere

    endsubroutine diff_flux
!***********************************************************************
    subroutine periodic_fold_back(dd,Boxsize)
      real,dimension(3),intent(in) :: Boxsize
      real,dimension(3),intent(inout) :: dd
      integer :: p,q,idim
      do idim=1,3
        q = floor(dd(idim)/(Boxsize(idim)/2))
        p = 0
        if (q.eq.1)  p = -1
        if (q.eq.-2) p = 1
        dd(idim) = dd(idim) + Boxsize(idim)*p
      enddo
    endsubroutine periodic_fold_back
!***********************************************************************
    subroutine calc_all_diff_fluxes(f,k,islope_limiter,h_slope_limited)
!
!  Calculates all <=3 components of the diffusive flux according to dimensionality. 
!
!  8-oct-15/MR: carved out from viscosity
!
      use General, only: notanumber
!
      real, dimension (:,:,:,:), intent(INOUT):: f
      integer,                   intent(IN)   :: k
      character(LEN=*),          intent(IN)   :: islope_limiter
      real,                      intent(IN)   :: h_slope_limited
             
      integer :: iff, nn, mm, ll
      real, dimension(mx-1) :: tmpx
      real, dimension(my-1) :: tmpy
      real, dimension(mz-1) :: tmpz

      if (.not.lcartesian_coords) &
        call fatal_error('calc_all_diff_fluxes', &
                         'Slope-limited diffusion not implemented for curvilinear coordinates.')

      f(:,:,:,iFF_diff1:iFF_diff2)=0.

      iff=iFF_diff

      if (nxgrid>1) then
        do nn=n1,n2; do mm=m1,m2
          tmpx = f(2:,mm,nn,k)-f(:mx-1,mm,nn,k)
if (notanumber(tmpx)) print*, 'TMPX:k,mm,nn=', k,mm,nn
!if (j==1) print*, 'TMPX:', tmpx
          call calc_diffusive_flux(tmpx,f(2:mx-2,mm,nn,iFF_char_c),islope_limiter,h_slope_limited,f(2:mx-2,mm,nn,iff))
if (notanumber(f(2:mx-2,mm,nn,iff))) print*, 'DIFFX:k,mm,nn=', k,mm,nn
        enddo; enddo
        iff=iff+1
      endif

      if (nygrid>1) then
        do nn=n1,n2; do ll=l1,l2
          tmpy = f(ll,2:,nn,k)-f(ll,:my-1,nn,k)
if (notanumber(tmpy)) print*, 'TMPY:k,mm,nn=', k,mm,nn
!if (j==2) print*, 'UY, iproc=:', iproc    !, '1-5,my-5:my-1'
!if (j==2) print'(16(e13.6,1x))', f(ll,:,nn,iuu+j-1)
          call calc_diffusive_flux(tmpy,f(ll,2:my-2,nn,iFF_char_c),islope_limiter,h_slope_limited,f(ll,2:my-2,nn,iff))
if (notanumber(f(ll,2:my-2,nn,iff))) print*, 'DIFFY:k,ll,mm=', k,ll,mm
        enddo; enddo
        iff=iff+1
      endif

      if (nzgrid>1) then
        do mm=m1,m2; do ll=l1,l2
if (notanumber(f(ll,mm,n1:n2,k))) print*, 'FARRAY(n1:n2):k,ll,mm=', k,ll,mm
if (notanumber(f(ll,mm,:,k))) print*, 'FARRAY(:):k,ll,mm=', k,ll,mm
          tmpz = f(ll,mm,2:,k)-f(ll,mm,:mz-1,k)
if (notanumber(tmpz)) print*, 'TMPZ:k,ll,mm=', k,ll,mm
          call calc_diffusive_flux(tmpz,f(ll,mm,2:mz-2,iFF_char_c),islope_limiter,h_slope_limited,f(ll,mm,2:mz-2,iff))
if (notanumber(f(ll,mm,2:mz-2,iff))) print*, 'DIFFZ:k,ll,mm=', k,ll,mm
        enddo; enddo
      endif
!
    endsubroutine calc_all_diff_fluxes
!***********************************************************************
endmodule Sub