! $Id$
!
!  This modules implements viscous heating and diffusion terms
!  here for shock viscosity
!    nu_total = nu + nu_shock*dx^2*smooth(max5(-(div u))))
!
!  NOTE: this works and has been tested for periodic boundaries.
!  With the current version, if your shock fronts reach a non-periodic
!  boundary, unexpected things may happen, so you should monitor the
!  behavior on the boundaries in this case.
!
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lshock = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 1
! COMMUNICATED AUXILIARIES 1
!
! PENCILS PROVIDED shock; gshock(3); shock_perp; gshock_perp(3)
!
!***************************************************************
module Shock
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
!
  implicit none
!
  include 'shock.h'
!
  integer :: ishock_max=1
  logical :: lshock_linear = .false.
  logical :: lgaussian_smooth=.false.
  logical :: lforce_periodic_shockviscosity=.false.
  logical :: lfix_Re_mesh=.false.
  real    :: div_threshold=0.0
  real    :: shock_linear = 0.01
  real    :: shock_div_pow = 1.
  logical :: lrewrite_shock_boundary=.false., lconvergence_only=.true.
!
  namelist /shock_run_pars/ &
      ishock_max, lgaussian_smooth, lforce_periodic_shockviscosity, &
      div_threshold, lrewrite_shock_boundary, lfix_Re_mesh, lshock_linear, shock_linear, &
      shock_div_pow, lconvergence_only
!
!  Diagnostic variables for print.in
! (needs to be consistent with reset list below)
!
  integer :: idiag_shockm=0        ! DIAG_DOC:
  integer :: idiag_shockmin=0      ! DIAG_DOC:
  integer :: idiag_shockmax=0      ! DIAG_DOC:
  integer :: idiag_gshockmax = 0   ! DIAG_DOC: $\max\left|\nabla\nu_{shock}\right|$
!
! xy averaged diagnostics given in xyaver.in written every it1d timestep
!
  integer :: idiag_shockmz=0       ! XYAVG_DOC:
!
! xz averaged diagnostics given in xzaver.in
!
  integer :: idiag_shockmy=0       ! XZAVG_DOC:
!
! yz averaged diagnostics given in yzaver.in
!
  integer :: idiag_shockmx=0       ! YZAVG_DOC:
!
  real, dimension (-nghost:nghost,-nghost:nghost,-nghost:nghost) :: smooth_factor
!
  interface shock_divu_perp
    module procedure shock_divu_perp_pencil
  endinterface
!
  contains
!***********************************************************************
    subroutine register_shock()
!
!  19-nov-02/tony: coded
!  24-jan-05/tony: modified from visc_shock.f90
!
      use FArrayManager
      use Messages, only: svn_id
!
      call farray_register_auxiliary('shock',ishock,communicated=.true.)
!
!  Identify version number.
!
      if (lroot) call svn_id( &
           "$Id$")
!
!  Writing files for use with IDL
!
      if (naux+naux_com <  maux+maux_com) aux_var(aux_count)=',shock $'
      if (naux+naux_com  == maux+maux_com) aux_var(aux_count)=',shock'
      aux_count=aux_count+1
      if (lroot) write(15,*) 'shock = fltarr(mx,my,mz)*one'
!
    endsubroutine register_shock
!***********************************************************************
    subroutine initialize_shock(f)
!
!  20-nov-02/tony: coded
!
      use Messages, only: fatal_error
      use Sub, only: register_report_aux, smoothing_kernel
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      real, dimension (-3:3) :: weights
      integer :: i,j,k
!
!  Initialize shock profile to zero
!
      f(:,:,:,ishock)=0.0
      if (ldivu_perp) then
        call register_report_aux('shock_perp',ishock_perp,communicated=.true.)
        if (.not. lmagnetic) call fatal_error('initialize_shock', &
                    'ldivu_perp requires magnetic field module for B_perp')
        f(:,:,:,ishock_perp)=0.0
      endif
!
!  Calculate the smoothing factors
!
      call smoothing_kernel(smooth_factor,lgaussian_smooth)
!
!      if (lgaussian_smooth) then
!        weights = (/1.,9.,45.,70.,45.,9.,1./)
!      else
!        weights = (/1.,6.,15.,20.,15.,6.,1./)
!      endif
!!
!      if (nxgrid > 1) then
!        do i = -3,3
!          smooth_factor(i,:,:) = smooth_factor(i,:,:)*weights(i)
!        enddo
!      else
!        smooth_factor(-3:-1,:,:) = 0.
!        smooth_factor(+1:+3,:,:) = 0.
!      endif
!!
!      if (nygrid > 1) then
!        do j = -3,3
!          smooth_factor(:,j,:) = smooth_factor(:,j,:)*weights(j)
!        enddo
!      else
!        smooth_factor(:,-3:-1,:) = 0.
!        smooth_factor(:,+1:+3,:) = 0.
!      endif
!!
!      if (nzgrid > 1) then
!        do k = -3,3
!          smooth_factor(:,:,k) = smooth_factor(:,:,k)*weights(k)
!        enddo
!      else
!        smooth_factor(:,:,-3:-1) = 0.
!        smooth_factor(:,:,+1:+3) = 0.
!      endif
!!
!      smooth_factor = smooth_factor / sum(smooth_factor)
!
!  Check that smooth order is within bounds
!
      if (lroot) then
        if (ishock_max < 1.or.ishock_max > nghost) &
          call fatal_error('initialize_shock', 'ishock_max needs to be between 1 and nghost.')
      endif
!
!  Die if periodic boundary condition for shock viscosity, but not for
!  velocity field. It can lead to subtle numerical errors near non-
!  periodic boundaries if the shock viscosity is assumed periodic.
!
      if (lrun .and. .not. lforce_periodic_shockviscosity) then
        if ((bcx(ishock) == 'p') .and. .not. all(bcx(iux:iuz) == 'p')) then
          if (lroot) then
            print*, 'initialize_shock: shock viscosity has bcx=''p'', but the velocity field is not'
            print*, '                  periodic! (you must set a proper boundary condition for the'
            print*, '                  shock viscosity)'
            print*, 'initialize_shock: bcx=', bcx
            print*, 'initialize_shock: to suppress this error,'
            print*, '                  set lforce_periodic_shockviscosity=T in &shock_run_pars'
          endif
          call fatal_error('initialize_shock','')
        endif
        if (bcy(ishock)=='p' .and. .not. all(bcy(iux:iuz)=='p')) then
          if (lroot) then
            print*, 'initialize_shock: shock viscosity has bcy=''p'', but the velocity field is not'
            print*, '                  periodic! (you must set a proper boundary condition for the'
            print*, '                  shock viscosity)'
            print*, 'initialize_shock: bcy=', bcy
            print*, 'initialize_shock: to suppress this error,'
            print*, '                  set lforce_periodic_shockviscosity=T in &shock_run_pars'
          endif
          call fatal_error('initialize_shock','')
        endif
        if (bcz(ishock)=='p' .and. .not. all(bcz(iux:iuz)=='p')) then
          if (lroot) then
            print*, 'initialize_shock: shock viscosity has bcz=''p'', but the velocity field is not'
            print*, '                  periodic! (you must set a proper boundary condition for the'
            print*, '                  shock viscosity)'
            print*, 'initialize_shock: bcz=', bcz
            print*, 'initialize_shock: to suppress this error,'
            print*, '                  set lforce_periodic_shockviscosity=T in &shock_run_pars'
          endif
          call fatal_error('initialize_shock','')
        endif
      endif
!
    endsubroutine initialize_shock
!***********************************************************************
    subroutine read_shock_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=shock_run_pars, IOSTAT=iostat)
!
    endsubroutine read_shock_run_pars
!***********************************************************************
    subroutine write_shock_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=shock_run_pars)
!
    endsubroutine write_shock_run_pars
!***********************************************************************
    subroutine rprint_shock(lreset,lwrite)
!
!  Writes ishock to index.pro file
!
!  24-nov-03/tony: adapted from rprint_ionization
!
      use Diagnostics
      use FArrayManager, only: farray_index_append
!
      logical :: lreset
      logical, optional :: lwrite
!
      integer :: iname, inamex, inamey, inamez
!
!  Reset everything in case of reset.
!  (this needs to be consistent with what is defined above!)
!
      if (lreset) then
        idiag_shockm=0; idiag_shockmin=0; idiag_shockmax=0
        idiag_gshockmax = 0
        idiag_shockmx=0; idiag_shockmy=0; idiag_shockmz=0
      endif
!
!  iname runs through all possible names that may be listed in print.in
!
      if (lroot.and.ip<14) print*,'rprint_shock: run through parse list'
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),&
            'shockm',idiag_shockm)
        call parse_name(iname,cname(iname),cform(iname),&
            'shockmin',idiag_shockmin)
        call parse_name(iname,cname(iname),cform(iname),&
            'shockmax',idiag_shockmax)
        call parse_name(iname, cname(iname), cform(iname), 'gshockmax', idiag_gshockmax)
      enddo
!
!  Check for those quantities for which we want yz-averages.
!
      do inamex=1,nnamex
        call parse_name(inamex,cnamex(inamex),cformx(inamex),'shockmx',idiag_shockmx)
      enddo
!
!  Check for those quantities for which we want xz-averages.
!
      do inamey=1,nnamey
        call parse_name(inamey,cnamey(inamey),cformy(inamey),'shockmy',idiag_shockmy)
      enddo
!
!  Check for those quantities for which we want xy-averages.
!
      do inamez=1,nnamez
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'shockmz',idiag_shockmz)
      enddo
!
!  check for those quantities for which we want video slices
!
      if (lwrite_slices) then 
        where(cnamev=='shock') cformv='DEFINED'
      endif
!
!  Write column where which shock variable is stored.
!
      if (present(lwrite)) then
        if (lwrite) then
          call farray_index_append('ishock',ishock)
          call farray_index_append('ishock_perp',ishock_perp)
        endif
      endif
!
    endsubroutine rprint_shock
!***********************************************************************
    subroutine get_slices_shock(f,slices)
!
!  Write slices for animation of shock variable.
!
!  26-jul-06/tony: coded
!
      use Slices_methods, only: assign_slices_scal

      real, dimension (mx,my,mz,mfarray) :: f
      type (slice_data) :: slices
!
!  Loop over slices
!
      select case (trim(slices%name))
!
!  Shock profile
!
        case ('shock'); call assign_slices_scal(slices,f,ishock)
!
      endselect
!
    endsubroutine get_slices_shock
!***********************************************************************
    subroutine pencil_criteria_shock()
!
!  All pencils that the Viscosity module depends on are specified here.
!
!  20-11-04/anders: coded
!
      if (idiag_shockm/=0 .or. idiag_shockmin/=0 .or. idiag_shockmax/=0 .or. &
          idiag_shockmx/=0 .or. idiag_shockmy/=0 .or. idiag_shockmz/=0) then
        lpenc_diagnos(i_shock)=.true.
      endif
!
      if (idiag_gshockmax /= 0) lpenc_diagnos(i_gshock) = .true.
!
    endsubroutine pencil_criteria_shock
!***********************************************************************
    subroutine pencil_interdep_shock(lpencil_in)
!
!  Interdependency among pencils from the Viscosity module is specified here.
!
!  20-11-04/anders: coded
!
      logical, dimension (npencils) :: lpencil_in
!
      call keep_compiler_quiet(lpencil_in)
!
    endsubroutine pencil_interdep_shock
!***********************************************************************
    subroutine calc_pencils_shock(f,p)
!
!  Calculate Viscosity pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  20-11-04/anders: coded
!
      use Diagnostics
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      intent(in) :: f
      intent(inout) :: p
! shock
      if (lpencil(i_shock)) p%shock=f(l1:l2,m,n,ishock)
! gshock
      if (lpencil(i_gshock)) call grad(f,ishock,p%gshock)
      if (ldivu_perp) then
! shock_perp
        if (lpencil(i_shock_perp)) p%shock_perp=f(l1:l2,m,n,ishock_perp)
! gshock_perp
        if (lpencil(i_gshock_perp)) call grad(f,ishock_perp,p%gshock_perp)
      endif
!
      if (ldiagnos) then
        if (idiag_shockm/=0)   call sum_mn_name( p%shock,idiag_shockm)
        if (idiag_shockmin/=0) call max_mn_name(-p%shock,idiag_shockmin,lneg=.true.)
        if (idiag_shockmax/=0) call max_mn_name( p%shock,idiag_shockmax)
        if (idiag_gshockmax /= 0) call max_mn_name(sqrt(sum(p%gshock**2, dim=2)), idiag_gshockmax)
      endif
!
      if (l1davgfirst) then
        call yzsum_mn_name_x(p%shock,idiag_shockmx)
        call xzsum_mn_name_y(p%shock,idiag_shockmy)
        call xysum_mn_name_z(p%shock,idiag_shockmz)
      endif
!
    endsubroutine calc_pencils_shock
!***********************************************************************
    subroutine calc_shock_profile(f)
!
!  Calculate divu based shock profile to be used in viscosity and
!  diffusion type terms.
!
!  23-nov-02/tony: coded
!  17-dec-08/ccyang: add divergence threshold
!
      use Boundcond, only: boundconds_x, boundconds_y, boundconds_z
      use Mpicomm, only: initiate_isendrcv_bdry, finalize_isendrcv_bdry, &
                         mpiallreduce_max, MPI_COMM_WORLD
      use Magnetic, only: bb_unitvec_shock
      use Sub, only: div
!
      real, dimension (mx,my,mz,mfarray), intent (inout) :: f
!
      real, dimension (mx,my,mz) :: tmp
      real, dimension (nx) :: penc, penc1, penc_perp
      real, dimension (mx,3) :: bb_hat
      integer :: imn
      integer :: i,j,k
      integer :: ni,nj,nk
      real :: shock_max, a=0., a1
!
!  Initialize shock to impossibly large number to force code crash in case of
!  inconsistencies in boundary conditions.
!
      f(:,:,:,ishock)=impossible
!
!  Compute divergence of velocity field.
!
      call boundconds_x(f,iux,iux)
!
      call initiate_isendrcv_bdry(f,iux,iuz)
!
      if (lshock_linear) a1 = shock_linear / dxmax
!
      do imn=1,ny*nz
!
        n = nn(imn)
        m = mm(imn)
!
        if (necessary(imn)) then
          call finalize_isendrcv_bdry(f,iux,iuz)
          call boundconds_y(f,iuy,iuy)
          call boundconds_z(f,iuz,iuz)
        endif
!
! The following will calculate div u for any coordinate system.
!
        call div(f,iuu,penc)
        if (lconvergence_only) then
          f(l1:l2,m,n,ishock) = max(0.0,-penc)
        else
          f(l1:l2,m,n,ishock) = abs(penc)
        endif
        if (shock_div_pow /= 1.) f(l1:l2,m,n,ishock)=f(l1:l2,m,n,ishock)**shock_div_pow
!
!  Add the linear term if requested
!
        if (lshock_linear) f(l1:l2,m,n,ishock) = f(l1:l2,m,n,ishock) + a1 * wave_speed(f)
!
      enddo
!
!  Cut off small divergence if requested.
!
      if (div_threshold > 0.0) &
        where(abs(f(:,:,:,ishock)) < div_threshold) f(:,:,:,ishock) = 0.0
!
!  Take maximum over a number of grid cells
!
      ni = merge(ishock_max,0,nxgrid > 1)
      nj = merge(ishock_max,0,nygrid > 1)
      nk = merge(ishock_max,0,nzgrid > 1)
!
!  Because of a bug in the shearing boundary conditions we must first manually
!  set the y boundary conditions on the shock profile.
!
      if (lshear) then
        call boundconds_y(f,ishock,ishock)
        call initiate_isendrcv_bdry(f,ishock,ishock)
        call finalize_isendrcv_bdry(f,ishock,ishock)
      endif
!
      call boundconds_x(f,ishock,ishock)
!
      call initiate_isendrcv_bdry(f,ishock,ishock)
!
      tmp = 0.0
!
      do imn=1,ny*nz
!
        n = nn(imn)
        m = mm(imn)
!
        if (necessary(imn)) then
          call finalize_isendrcv_bdry(f,ishock,ishock)
          call boundconds_y(f,ishock,ishock)
          call boundconds_z(f,ishock,ishock)
        endif
!
        penc = 0.0
!
        do k=-nk,nk
        do j=-nj,nj
        do i=-ni,ni
          penc = max(penc,f(l1+i:l2+i,m+j,n+k,ishock))
        enddo
        enddo
        enddo
!
        tmp(l1:l2,m,n) = penc
!
      enddo
!
      f(:,:,:,ishock) = tmp
!
!  Smooth with a Gaussian profile
!
      ni = merge(3,0,nxgrid > 1)
      nj = merge(3,0,nygrid > 1)
      nk = merge(3,0,nzgrid > 1)
!
!  Because of a bug in the shearing boundary conditions we must first manually
!  set the y boundary conditions on the shock profile.
!
      if (lshear) then
        call boundconds_y(f,ishock,ishock)
        call initiate_isendrcv_bdry(f,ishock,ishock)
        call finalize_isendrcv_bdry(f,ishock,ishock)
      endif
!
      call boundconds_x(f,ishock,ishock)
      call initiate_isendrcv_bdry(f,ishock,ishock)
!
      tmp = 0.0
!
      do imn=1,ny*nz
!
        n = nn(imn)
        m = mm(imn)
!
        if (necessary(imn)) then
          call finalize_isendrcv_bdry(f,ishock,ishock)
          call boundconds_y(f,ishock,ishock)
          call boundconds_z(f,ishock,ishock)
        endif
!
        penc = 0.0
!
        do k=-nk,nk
        do j=-nj,nj
        do i=-ni,ni
          penc = penc + smooth_factor(i,j,k)*f(l1+i:l2+i,m+j,n+k,ishock)
        enddo
        enddo
        enddo
!
        tmp(l1:l2,m,n) = penc
!
      enddo
!
      fix_Re: if (lfix_Re_mesh) then
!
!  Scale given a fixed mesh Reynolds number.
!
        if (headtt) print *, 'Shock: fix mesh Reynolds number'
        call mpiallreduce_max(maxval(tmp(l1:l2,m1:m2,n1:n2)), shock_max, comm=MPI_COMM_WORLD)
        shock: if (shock_max > 0.) then
          a1 = 0.
          scan_z: do n = n1, n2
            scan_y: do m = m1, m2
              penc = 0.
              penc1 = 0.
              xdir: if (nxgrid > 1) then
                penc = penc + abs(f(l1:l2,m,n,iux) * dx_1(l1:l2))
                penc1 = penc1 + dx_1(l1:l2)**2
              endif xdir
              ydir: if (nygrid > 1) then
                penc = penc + abs(f(l1:l2,m,n,iuy) * dy_1(m))
                penc1 = penc1 + dy_1(m)**2
              endif ydir
              zdir: if (nzgrid > 1) then
                penc = penc + abs(f(l1:l2,m,n,iuz) * dz_1(n))
                penc1 = penc1 + dz_1(n)**2
              endif zdir
              a1 = max(a1, maxval(penc / penc1))
            enddo scan_y
          enddo scan_z
          call mpiallreduce_max(a1, a, comm=MPI_COMM_WORLD)
          a = a / (re_mesh * pi * shock_max)
        else shock
          a = dxmin**2
        endif shock
        f(:,:,:,ishock) = a * tmp
      else fix_Re
!
!  Scale by dxmin**2.
!
!  The shearing boundary conditions have a bug that can cause nonsensical
!  numbers in the corners (bug #61). We can overwrite this bug by defining the
!  shock viscosity in the ghost zones too. THIS WOULD NOT BE NECESSARY
!  IF THE BUG IN THE SHEARING BOUNDARY CONDITIONS WOULD BE FIXED.
!
        if (.not.lrewrite_shock_boundary) then
          f(l1:l2,m1:m2,n1:n2,ishock) = tmp(l1:l2,m1:m2,n1:n2)*dxmin**2
        else
          f(:,:,:,ishock) = tmp*dxmin**2
        endif
      endif fix_Re
!
      if (ldivu_perp) then
        call boundconds_x(f,ishock_perp,ishock_perp)
        call initiate_isendrcv_bdry(f,ishock_perp,ishock_perp)
!
        do imn=1,ny*nz
!
          n = nn(imn)
          m = mm(imn)
!
          if (necessary(imn)) then
            call finalize_isendrcv_bdry(f,ishock_perp,ishock_perp)
            call boundconds_y(f,ishock_perp,ishock_perp)
            call boundconds_z(f,ishock_perp,ishock_perp)
          endif
!
          call div(f,iuu,penc)
          call bb_unitvec_shock(f,bb_hat)
          call shock_divu_perp(f,bb_hat,penc,penc_perp)
          f(l1:l2,m,n,ishock_perp)=max(0.,-penc_perp)
!
        enddo
        tmp = 0.0
        do imn=1,ny*nz
!
          n = nn(imn)
          m = mm(imn)
!
          if (necessary(imn)) then
            call finalize_isendrcv_bdry(f,ishock_perp,ishock_perp)
            call boundconds_y(f,ishock_perp,ishock_perp)
            call boundconds_z(f,ishock_perp,ishock_perp)
          endif
!
          penc = 0.0
!
          do k=-nk,nk
          do j=-nj,nj
          do i=-ni,ni
            penc = penc + smooth_factor(i,j,k)*f(l1+i:l2+i,m+j,n+k,ishock_perp)
          enddo
          enddo
          enddo
!
          tmp(l1:l2,m,n) = penc
!
        enddo
        if (.not.lrewrite_shock_boundary) then
          f(l1:l2,m1:m2,n1:n2,ishock_perp) = tmp(l1:l2,m1:m2,n1:n2) * dxmin**2
        else
          f(:,:,:,ishock_perp) = tmp * dxmin**2
        endif
!
      endif
!
    endsubroutine calc_shock_profile
!***********************************************************************
    subroutine shock_divu_perp_pencil(f,bb_hat,divu,divu_perp)
!
!  Calculate `perpendicular divergence' of u.
!  nabla_perp.uu = nabla.uu - (1/b2)*bb.(bb.nabla)*uu
!
!  16-aug-06/tobi: coded
!  07-jun-18/fred: revised to include higher order gradient u
!
      real, dimension (mx,my,mz,mfarray), intent (in) :: f
      real, dimension (mx,3), intent (in) :: bb_hat
      real, dimension (nx), intent (in) :: divu
      real, dimension (nx), intent (out) :: divu_perp
!
      real, dimension (nx) :: fac
!
      divu_perp = divu
!
      if (nxgrid/=1) then
         fac=bb_hat(l1:l2,1)/(60*dx)
         divu_perp(:) = divu_perp(:)                          &
           - fac*(bb_hat(l1:l2,1)*( f(l1+3:l2+3,m  ,n  ,iux)   &
                                   + 45.0*(f(l1+1:l2+1,m  ,n  ,iux)   &
                                         - f(l1-1:l2-1,m  ,n  ,iux))  &
                                   -  9.0*(f(l1+2:l2+2,m  ,n  ,iux)   &
                                         - f(l1-2:l2-2,m  ,n  ,iux))  &
                                         - f(l1-3:l2-3,m  ,n  ,iux) ) &
                       + bb_hat(l1:l2,2)*( f(l1+3:l2+3,m  ,n  ,iuy)   &
                                   + 45.0*(f(l1+1:l2+1,m  ,n  ,iuy)   &
                                         - f(l1-1:l2-1,m  ,n  ,iuy))  &
                                   -  9.0*(f(l1+2:l2+2,m  ,n  ,iuy)   &
                                         - f(l1-2:l2-2,m  ,n  ,iuy))  &
                                         - f(l1-3:l2-3,m  ,n  ,iuy) ) &
                       + bb_hat(l1:l2,3)*( f(l1+3:l2+3,m  ,n  ,iuz)   &
                                   + 45.0*(f(l1+1:l2+1,m  ,n  ,iuz)   &
                                         - f(l1-1:l2-1,m  ,n  ,iuz))  &
                                   -  9.0*(f(l1+2:l2+2,m  ,n  ,iuz)   &
                                         - f(l1-2:l2-2,m  ,n  ,iuz))  &
                                         - f(l1-3:l2-3,m  ,n  ,iuz) ) )
      endif
!
      if (nygrid/=1) then
         fac=bb_hat(l1:l2,2)/(60*dy)
         divu_perp(:) = divu_perp(:)                          &
           - fac*(bb_hat(l1:l2,1)*( f(  l1:l2  ,m+3,n  ,iux)   &
                                   + 45.0*(f(  l1:l2  ,m+1,n  ,iux)   &
                                         - f(  l1:l2  ,m-1,n  ,iux))  &
                                   -  9.0*(f(  l1:l2  ,m+2,n  ,iux)   &
                                         - f(  l1:l2  ,m-2,n  ,iux))  &
                                         - f(  l1:l2  ,m-3,n  ,iux) ) &
                       + bb_hat(l1:l2,2)*( f(  l1:l2  ,m+3,n  ,iuy)   &
                                   + 45.0*(f(  l1:l2  ,m+1,n  ,iuy)   &
                                         - f(  l1:l2  ,m-1,n  ,iuy))  &
                                   -  9.0*(f(  l1:l2  ,m+2,n  ,iuy)   &
                                         - f(  l1:l2  ,m-2,n  ,iuy))  &
                                         - f(  l1:l2  ,m-3,n  ,iuy) ) &
                       + bb_hat(l1:l2,3)*( f(  l1:l2  ,m+3,n  ,iuz)   &
                                   + 45.0*(f(  l1:l2  ,m+1,n  ,iuz)   &
                                         - f(  l1:l2  ,m-1,n  ,iuz))  &
                                   -  9.0*(f(  l1:l2  ,m+2,n  ,iuz)   &
                                         - f(  l1:l2  ,m-2,n  ,iuz))  &
                                         - f(  l1:l2  ,m-3,n  ,iuz) ) )
      endif
!
      if (nzgrid/=1) then
         fac=bb_hat(l1:l2,3)/(60*dz)
         divu_perp(:) = divu_perp(:)                          &
           - fac*(bb_hat(l1:l2,1)*( f(  l1:l2  ,m  ,n+3,iux)   &
                                   + 45.0*(f(  l1:l2  ,m,  n+1,iux)   &
                                         - f(  l1:l2  ,m,  n-1,iux))  &
                                   -  9.0*(f(  l1:l2  ,m,  n+2,iux)   &
                                         - f(  l1:l2  ,m,  n-2,iux))  &
                                         - f(  l1:l2  ,m  ,n-3,iux) ) &
                       + bb_hat(l1:l2,2)*( f(  l1:l2  ,m  ,n+3,iuy)   &
                                   + 45.0*(f(  l1:l2  ,m,  n+1,iuy)   &
                                         - f(  l1:l2  ,m,  n-1,iuy))  &
                                   -  9.0*(f(  l1:l2  ,m,  n+2,iuy)   &
                                         - f(  l1:l2  ,m,  n-2,iuy))  &
                                         - f(  l1:l2  ,m  ,n-3,iuy) ) &
                       + bb_hat(l1:l2,3)*( f(  l1:l2  ,m  ,n+3,iuz)   &
                                   + 45.0*(f(  l1:l2  ,m,  n+1,iuz)   &
                                         - f(  l1:l2  ,m,  n-1,iuz))  &
                                   -  9.0*(f(  l1:l2  ,m,  n+2,iuz)   &
                                         - f(  l1:l2  ,m,  n-2,iuz))  &
                                         - f(  l1:l2  ,m  ,n-3,iuz) ) )
      endif
!
    endsubroutine shock_divu_perp_pencil
!***********************************************************************
    subroutine calc_shock_profile_simple(f)
!
!  Calculate divu based shock profile to be used in viscosity and
!  diffusion type terms.
!
!  12-apr-05/tony: coded
!
      real, dimension (mx,my,mz,mfarray), intent (in) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine calc_shock_profile_simple
!***********************************************************************
    function wave_speed(f) result(speed)
!
!  Calculate the wave speeds along one pencil.
!
!  23-nov-13/ccyang: coded.
!
      use EquationOfState, only: eoscalc, rho0
      use Magnetic, only: get_bext
      use Sub, only: gij, curl_mn
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      real, dimension(nx) :: speed
!
      real, dimension(nx,3,3) :: aij
      real, dimension(nx,3) :: bb
      real, dimension(nx) :: b2
      real, dimension(3) :: B_ext
!
!  Get the sound speed.
!
      call eoscalc(f, nx, cs2=speed)
!
!  Add the Alfven speed.
!
      alfven: if (lmagnetic) then
!
        bfield: if (lbfield) then
          bb = f(l1:l2,m,n,ibx:ibz)
        else bfield
          call gij(f, iaa, aij, 1)
          call curl_mn(aij, bb, f(:,m,n,iax:iaz))
        endif bfield
        call get_bext(B_ext)
        b2 = sum((bb + spread(B_ext,1,nx))**2, dim=2)
!
        density: if (ldensity) then
          if (ldensity_nolog) then
            speed = speed + mu01 * b2 / f(l1:l2,m,n,irho)
          else
            speed = speed + mu01 * b2 / exp(f(l1:l2,m,n,ilnrho))
          endif
        else density
          speed = speed + mu01 / rho0 * b2
        endif density
!
      endif alfven
!
      speed = sqrt(speed)
!
    endfunction wave_speed
!***********************************************************************
endmodule Shock