mod_viscosity.t

Go to the documentation of this file.

!>   The module add viscous source terms and check time step
!>
!>   Viscous forces in the momentum equations:
!>   d m_i/dt +=  - div (vc_mu * PI)
!>   !! Viscous work in the energy equation:
!>   !! de/dt    += - div (v . vc_mu * PI)
!>   where the PI stress tensor is
!>   PI_i,j = - (dv_j/dx_i + dv_i/dx_j) + Kronecker delta_i,j*(2/3)*Sum_k dv_k/dx_k
!>   where vc_mu is the dynamic viscosity coefficient (g cm^-1 s^-1).
module mod_viscosity
  use mod_comm_lib, only: mpistop
  implicit none
 
  !> Viscosity coefficient
  double precision, public :: vc_mu = 1.d0
 
  !> Index of the density (in the w array)
  integer, private, parameter              :: rho_ = 1
 
  !> Indices of the momentum density
  integer, allocatable, private, protected :: mom(:)
 
  !> Index of the energy density (-1 if not present)
  integer, private, protected              :: e_
  !> Indices of the velocity for the form of better vectorization
  integer, public, protected              :: v1_,v2_,v3_
 
  !> fourth order
  logical :: vc_4th_order = .false.
 
  !> source split or not
  logical :: vc_split= .false.
 
  !> whether to compute the viscous terms as
  !> fluxes (ie in the div on the LHS), or not (by default)
  logical :: viscindiv= .false.
 
  procedure(sub_add_source), pointer :: viscosity_add_source => null()
  ! Public methods
  public :: viscosity_add_source
  public :: visc_get_flux_prim
 
contains
  !> Read this module"s parameters from a file
  subroutine vc_params_read(files)
    use mod_global_parameters, only: unitpar
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /vc_list/ vc_mu, vc_4th_order, vc_split, viscindiv
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, vc_list, end=111)
111    close(unitpar)
    end do
 
 
  end subroutine vc_params_read
 
  !> Initialize the module
  subroutine viscosity_init(phys_wider_stencil)
    use mod_global_parameters
    use mod_geometry
    integer, intent(inout) :: phys_wider_stencil
    integer :: nwx,idir
 
    call vc_params_read(par_files)
 
    if(vc_split) any_source_split=.true.
 
    ! Determine flux variables
    nwx = 1                  ! rho (density)
 
    allocate(mom(ndir))
    do idir = 1, ndir
       nwx    = nwx + 1
       mom(idir) = nwx       ! momentum density
    end do
    v^c_=mom(^c);
 
    nwx = nwx + 1
    e_     = nwx          ! energy density
 
    if (viscindiv) then
      ! to compute the derivatives from left and right upwinded values
      phys_wider_stencil = 1
    end if
    select case (coordinate)
      case (cartesian,cartesian_stretched,cartesian_expansion)
        viscosity_add_source => viscosity_add_source_cartesian
      case (cylindrical)
        viscosity_add_source => viscosity_add_source_cylinder
      case (spherical)
        viscosity_add_source => viscosity_add_source_sphere
    end select
 
  end subroutine viscosity_init
 
  subroutine viscosity_add_source_cartesian(qdt,ixI^L,ixO^L,wCT,wp,w,x,&
       energy,qsourcesplit,active)
  ! Add viscosity source in isotropic Newtonian fluids to w within ixO
  ! neglecting bulk viscosity
  ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wp(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
 
    double precision:: lambda(ixI^S,ndir,ndir),tmp(ixI^S),vlambda(ixI^S,ndir),qdtmu,divv23
    integer:: ix^L,ix^D
 
    if (viscindiv) return
 
    if(qsourcesplit .eqv. vc_split) then
      active = .true.
      ! standard case, textbook viscosity
      ! Calculating viscosity sources
      if(.not.vc_4th_order) then
        ! involves second derivatives, two extra layers
        ix^l=ixo^l^ladd2;
        if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
          call mpistop("error for viscous source addition, 2 layers needed")
        ix^l=ixo^l^ladd1;
      else
        ! involves second derivatives, four extra layers
        ix^l=ixo^l^ladd4;
        if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
          call mpistop("error for viscous source addition"//&
          "requested fourth order gradients: 4 layers needed")
        ix^l=ixo^l^ladd2;
      end if
      qdtmu=qdt*vc_mu
      {^ifthreed
      {do ix^db=ixmin^db,ixmax^db\}
        ! idim=1, idir=1
        lambda(ix^d,1,1)=(wp(ix1+1,ix2,ix3,v1_)-wp(ix1-1,ix2,ix3,v1_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))*two*qdtmu
        ! idim=1, idir=2
        lambda(ix^d,1,2)=(wp(ix1+1,ix2,ix3,v2_)-wp(ix1-1,ix2,ix3,v2_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
        ! idim=1, idir=3
        lambda(ix^d,1,3)=(wp(ix1+1,ix2,ix3,v3_)-wp(ix1-1,ix2,ix3,v3_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
        ! idim=2, idir=1
        lambda(ix^d,2,1)=(wp(ix1,ix2+1,ix3,v1_)-wp(ix1,ix2-1,ix3,v1_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))
        ! idim=2, idir=2
        lambda(ix^d,2,2)=(wp(ix1,ix2+1,ix3,v2_)-wp(ix1,ix2-1,ix3,v2_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*two*qdtmu
        ! idim=2, idir=3
        lambda(ix^d,2,3)=(wp(ix1,ix2+1,ix3,v3_)-wp(ix1,ix2-1,ix3,v3_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))
        ! idim=3, idir=1
        lambda(ix^d,3,1)=(wp(ix1,ix2,ix3+1,v1_)-wp(ix1,ix2,ix3-1,v1_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))
        ! idim=3, idir=2
        lambda(ix^d,3,2)=(wp(ix1,ix2,ix3+1,v2_)-wp(ix1,ix2,ix3-1,v2_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))
        ! idim=3, idir=3
        lambda(ix^d,3,3)=(wp(ix1,ix2,ix3+1,v3_)-wp(ix1,ix2,ix3-1,v3_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*two*qdtmu
        ! dv_i/d_j + dv_j/d_i
        lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*qdtmu
        lambda(ix^d,2,1)=lambda(ix^d,1,2)
        lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*qdtmu
        lambda(ix^d,3,1)=lambda(ix^d,1,3)
        lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*qdtmu
        lambda(ix^d,3,2)=lambda(ix^d,2,3)
        divv23=third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
        lambda(ix^d,1,1)=lambda(ix^d,1,1)-divv23
        lambda(ix^d,2,2)=lambda(ix^d,2,2)-divv23
        lambda(ix^d,3,3)=lambda(ix^d,3,3)-divv23
        if(energy) then
          vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
          vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
          vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
        end if
      {end do\}
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      {do ix^db=ixomin^db,ixomax^db\}
        w(ix^d,mom(1))=(lambda(ix1+1,ix2,ix3,1,1)-lambda(ix1-1,ix2,ix3,1,1))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
                       (lambda(ix1,ix2+1,ix3,2,1)-lambda(ix1,ix2-1,ix3,2,1))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+&
                       (lambda(ix1,ix2,ix3+1,3,1)-lambda(ix1,ix2,ix3-1,3,1))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+w(ix^d,mom(1))
        w(ix^d,mom(2))=(lambda(ix1+1,ix2,ix3,1,2)-lambda(ix1-1,ix2,ix3,1,2))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
                       (lambda(ix1,ix2+1,ix3,2,2)-lambda(ix1,ix2-1,ix3,2,2))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+&
                       (lambda(ix1,ix2,ix3+1,3,2)-lambda(ix1,ix2,ix3-1,3,2))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+w(ix^d,mom(2))
        w(ix^d,mom(3))=(lambda(ix1+1,ix2,ix3,1,3)-lambda(ix1-1,ix2,ix3,1,3))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
                       (lambda(ix1,ix2+1,ix3,2,3)-lambda(ix1,ix2-1,ix3,2,3))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+&
                       (lambda(ix1,ix2,ix3+1,3,3)-lambda(ix1,ix2,ix3-1,3,3))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+w(ix^d,mom(3))
      {end do\}
      }
      {^iftwod
      {do ix^db=ixmin^db,ixmax^db\}
        ! idim=1, idir=1
        lambda(ix^d,1,1)=(wp(ix1+1,ix2,v1_)-wp(ix1-1,ix2,v1_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))*two*qdtmu
        ! idim=1, idir=2
        lambda(ix^d,1,2)=(wp(ix1+1,ix2,v2_)-wp(ix1-1,ix2,v2_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
        ! idim=2, idir=1
        lambda(ix^d,2,1)=(wp(ix1,ix2+1,v1_)-wp(ix1,ix2-1,v1_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))
        ! idim=2, idir=2
        lambda(ix^d,2,2)=(wp(ix1,ix2+1,v2_)-wp(ix1,ix2-1,v2_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*two*qdtmu
        ! dv_i/d_j + dv_j/d_i
        lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*qdtmu
        lambda(ix^d,2,1)=lambda(ix^d,1,2)
        divv23=third*(lambda(ix^d,1,1)+lambda(ix^d,2,2))
        lambda(ix^d,1,1)=lambda(ix^d,1,1)-divv23
        lambda(ix^d,2,2)=lambda(ix^d,2,2)-divv23
        if(ndir==3) then
          ! idim=1, idir=3
          lambda(ix1,ix2,1,3)=(wp(ix1+1,ix2,v3_)-wp(ix1-1,ix2,v3_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))*qdtmu
          ! idim=2, idir=3
          lambda(ix1,ix2,2,3)=(wp(ix1,ix2+1,v3_)-wp(ix1,ix2-1,v3_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*qdtmu
          lambda(ix1,ix2,3,1)=lambda(ix1,ix2,1,3)
          lambda(ix1,ix2,3,2)=lambda(ix1,ix2,2,3)
          lambda(ix1,ix2,3,3)=-divv23
          if(energy) then
            vlambda(ix1,ix2,1)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,1)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,1)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,1)
            vlambda(ix1,ix2,2)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,2)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,2)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,2)
            vlambda(ix1,ix2,3)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,3)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,3)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,3)
          end if
        else if(energy) then
          vlambda(ix1,ix2,1)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,1)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,1)
          vlambda(ix1,ix2,2)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,2)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,2)
        end if
      {end do\}
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      {do ix^db=ixomin^db,ixomax^db\}
        w(ix^d,mom(1))=(lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                       (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+w(ix^d,mom(1))
        w(ix^d,mom(2))=(lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                       (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+w(ix^d,mom(2))
        if(ndir==3) then
          w(ix1,ix2,mom(3))=(lambda(ix1+1,ix2,1,3)-lambda(ix1-1,ix2,1,3))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                         (lambda(ix1,ix2+1,2,3)-lambda(ix1,ix2-1,2,3))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+w(ix1,ix2,mom(3))
        end if
      {end do\}
      }
      {^ifoned
      do ix1=ixmin1,ixmax1
        ! idim=1, idir=1
        lambda(ix1,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))*two*qdtmu
        ! dv_i/d_j + dv_j/d_i
        divv23=third*lambda(ix1,1,1)
        lambda(ix1,1,1)=lambda(ix1,1,1)-divv23
        if(ndir==2) then
          ! idim=1, idir=2
          lambda(ix1,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))*qdtmu
          ! dv_i/d_j + dv_j/d_i
          lambda(ix1,2,1)=lambda(ix1,1,2)
          lambda(ix1,2,2)=-divv23
          if(energy) then
            vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)+wp(ix1,v2_)*lambda(ix1,2,1)
            vlambda(ix1,2)=wp(ix1,v1_)*lambda(ix1,1,2)+wp(ix1,v2_)*lambda(ix1,2,2)
          end if
        else if(ndir==3) then
          ! idim=1, idir=2
          lambda(ix1,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))*qdtmu
          ! dv_i/d_j + dv_j/d_i
          lambda(ix1,2,1)=lambda(ix1,1,2)
          lambda(ix1,2,2)=-divv23
          ! idim=1, idir=3
          lambda(ix1,1,3)=(wp(ix1+1,v3_)-wp(ix1-1,v3_))/(x(ix1+1,1)-x(ix1-1,1))*qdtmu
          lambda(ix1,3,1)=lambda(ix1,1,3)
          lambda(ix1,3,3)=-divv23
          if(energy) then
            vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)+wp(ix1,v2_)*lambda(ix1,2,1)+wp(ix1,v3_)*lambda(ix1,3,1)
            vlambda(ix1,2)=wp(ix1,v1_)*lambda(ix1,1,2)+wp(ix1,v2_)*lambda(ix1,2,2)+wp(ix1,v3_)*lambda(ix1,3,2)
            vlambda(ix1,3)=wp(ix1,v1_)*lambda(ix1,1,3)+wp(ix1,v2_)*lambda(ix1,2,3)+wp(ix1,v3_)*lambda(ix1,3,3)
          end if
        else if(energy) then
          vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)
        end if
      end do
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      do ix1=ixomin1,ixomax1
        if(ndir==1) then
          w(ix1,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+w(ix1,mom(1))
        else if(ndir==2) then
          w(ix1,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+w(ix1,mom(1))
          w(ix1,mom(2))=(lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+w(ix1,mom(2))
        else
          w(ix1,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+w(ix1,mom(1))
          w(ix1,mom(2))=(lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+w(ix1,mom(2))
          w(ix1,mom(3))=(lambda(ix1+1,1,3)-lambda(ix1-1,1,3))/(x(ix1+1,1)-x(ix1-1,1))+w(ix1,mom(3))
        end if
      end do
      }
      if(energy) then
        ! de/dt= +div(v.dot.[mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v *kr])
        ! thus e=e+d_i v_j tensor_ji
        call divvector(vlambda,ixi^l,ixo^l,tmp)
        w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)
      end if
    end if
 
  end subroutine viscosity_add_source_cartesian
 
  subroutine viscosity_add_source_sphere(qdt,ixI^L,ixO^L,wCT,wp,w,x,&
       energy,qsourcesplit,active)
  ! Add viscosity source in isotropic Newtonian fluids to w within ixO
  ! neglecting bulk viscosity
  ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wp(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
 
    double precision :: lambda(ixI^S,ndir,ndir),vlambda(ixI^S,ndir),invr(ixI^S),tctan(ixI^S),invrsin(ixI^S)
    double precision :: qdtmu,tsin,tcos,divv23
    integer:: ix^L,ix^D
 
    if (viscindiv) return
 
    if(qsourcesplit .eqv. vc_split) then
      active = .true.
      ! standard case, textbook viscosity
      ! Calculating viscosity sources
      if(.not.vc_4th_order) then
        ! involves second derivatives, two extra layers
        ix^l=ixo^l^ladd2;
        if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
          call mpistop("error for viscous source addition, 2 layers needed")
        ix^l=ixo^l^ladd1;
      else
        ! involves second derivatives, four extra layers
        ix^l=ixo^l^ladd4;
        if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
          call mpistop("error for viscous source addition"//&
          "requested fourth order gradients: 4 layers needed")
        ix^l=ixo^l^ladd2;
      end if
 
      ! construct lambda tensor: lambda_ij = gradv_ij + gradv_ji
      ! initialize
      qdtmu=qdt*vc_mu
      {^ifthreed
      {do ix^db=ixmin^db,ixmax^db\}
        invr(ix^d)=1.d0/x(ix^d,1)
        tcos=dcos(x(ix^d,2))
        ! idim=1, idir=1
        lambda(ix^d,1,1)=(wp(ix1+1,ix2,ix3,v1_)-wp(ix1-1,ix2,ix3,v1_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))*two*qdtmu
        ! idim=1, idir=2
        lambda(ix^d,1,2)=(wp(ix1+1,ix2,ix3,v2_)-wp(ix1-1,ix2,ix3,v2_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
        ! idim=1, idir=3
        lambda(ix^d,1,3)=(wp(ix1+1,ix2,ix3,v3_)-wp(ix1-1,ix2,ix3,v3_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
        ! idim=2, idir=1
        lambda(ix^d,2,1)=((wp(ix1,ix2+1,ix3,v1_)-wp(ix1,ix2-1,ix3,v1_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))-wp(ix^d,v2_))*invr(ix^d)
        ! idim=2, idir=2
        lambda(ix^d,2,2)=((wp(ix1,ix2+1,ix3,v2_)-wp(ix1,ix2-1,ix3,v2_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+wp(ix^d,v1_))*invr(ix^d)*&
                         two*qdtmu
        ! idim=2, idir=3
        lambda(ix^d,2,3)=((wp(ix1,ix2+1,ix3,v3_)-wp(ix1,ix2-1,ix3,v3_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+wp(ix^d,v2_)*tcos)*invr(ix^d)
        tsin=dsin(x(ix^d,2))
        invrsin(ix^d)=invr(ix^d)/tsin
        tctan(ix^d)=tcos/tsin
        ! idim=3, idir=1
        lambda(ix^d,3,1)=((wp(ix1,ix2,ix3+1,v1_)-wp(ix1,ix2,ix3-1,v1_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))-wp(ix^d,v3_)*tsin)*invrsin(ix^d)
        ! idim=3, idir=2
        lambda(ix^d,3,2)=((wp(ix1,ix2,ix3+1,v2_)-wp(ix1,ix2,ix3-1,v2_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))-wp(ix^d,v3_)*tcos)*invrsin(ix^d)
        ! idim=3, idir=3
        lambda(ix^d,3,3)=((wp(ix1,ix2,ix3+1,v3_)-wp(ix1,ix2,ix3-1,v3_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+wp(ix^d,v1_)*tsin+&
            wp(ix^d,v2_)*tcos)*invrsin(ix^d)*two*qdtmu
        ! dv_i/d_j + dv_j/d_i
        lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*qdtmu
        lambda(ix^d,2,1)=lambda(ix^d,1,2)
        lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*qdtmu
        lambda(ix^d,3,1)=lambda(ix^d,1,3)
        lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*qdtmu
        lambda(ix^d,3,2)=lambda(ix^d,2,3)
        divv23=third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
        lambda(ix^d,1,1)=lambda(ix^d,1,1)-divv23
        lambda(ix^d,2,2)=lambda(ix^d,2,2)-divv23
        lambda(ix^d,3,3)=lambda(ix^d,3,3)-divv23
        if(energy) then
          vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
          vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
          vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
        end if
      {end do\}
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      {do ix^db=ixomin^db,ixomax^db\}
        w(ix^d,mom(1))=(lambda(ix1+1,ix2,ix3,1,1)-lambda(ix1-1,ix2,ix3,1,1))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
                       (lambda(ix1,ix2+1,ix3,2,1)-lambda(ix1,ix2-1,ix3,2,1))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
                       (lambda(ix1,ix2,ix3+1,3,1)-lambda(ix1,ix2,ix3-1,3,1))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invrsin(ix^d)+&
                    two*lambda(ix^d,1,1)*invr(ix^d)+lambda(ix^d,2,1)*invr(ix^d)*tctan(ix^d)+w(ix^d,mom(1))
        w(ix^d,mom(2))=(lambda(ix1+1,ix2,ix3,1,2)-lambda(ix1-1,ix2,ix3,1,2))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
                       (lambda(ix1,ix2+1,ix3,2,2)-lambda(ix1,ix2-1,ix3,2,2))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
                       (lambda(ix1,ix2,ix3+1,3,2)-lambda(ix1,ix2,ix3-1,3,2))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invrsin(ix^d)+&
                    two*lambda(ix^d,1,2)*invr(ix^d)+lambda(ix^d,2,2)*invr(ix^d)*tctan(ix^d)-lambda(ix^d,3,3)*tctan(ix^d)*invr(ix^d)**2+&
                    w(ix^d,mom(2))
        w(ix^d,mom(3))=(lambda(ix1+1,ix2,ix3,1,3)-lambda(ix1-1,ix2,ix3,1,3))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
                       (lambda(ix1,ix2+1,ix3,2,3)-lambda(ix1,ix2-1,ix3,2,3))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
                       (lambda(ix1,ix2,ix3+1,3,3)-lambda(ix1,ix2,ix3-1,3,3))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invrsin(ix^d)+&
                    two*lambda(ix^d,1,3)*invr(ix^d)+lambda(ix^d,2,3)*(invr(ix^d)+invr(ix^d)**2)*tctan(ix^d)+w(ix^d,mom(3))
      {end do\}
      }
      {^iftwod
      {do ix^db=ixmin^db,ixmax^db\}
        invr(ix1,ix2)=1.d0/x(ix1,ix2,1)
        ! idim=1, idir=1
        lambda(ix1,ix2,1,1)=(wp(ix1+1,ix2,v1_)-wp(ix1-1,ix2,v1_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))*two*qdtmu
        ! idim=1, idir=2
        lambda(ix1,ix2,1,2)=(wp(ix1+1,ix2,v2_)-wp(ix1-1,ix2,v2_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
        ! idim=2, idir=1
        lambda(ix1,ix2,2,1)=((wp(ix1,ix2+1,v1_)-wp(ix1,ix2-1,v1_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))-wp(ix1,ix2,v2_))*invr(ix1,ix2)
        ! idim=2, idir=2
        lambda(ix1,ix2,2,2)=((wp(ix1,ix2+1,v2_)-wp(ix1,ix2-1,v2_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+wp(ix1,ix2,v1_))*invr(ix1,ix2)*&
                         two*qdtmu
        ! dv_i/d_j + dv_j/d_i
        lambda(ix1,ix2,1,2)=(lambda(ix1,ix2,1,2)+lambda(ix1,ix2,2,1))*qdtmu
        lambda(ix1,ix2,2,1)=lambda(ix1,ix2,1,2)
        divv23=third*(lambda(ix1,ix2,1,1)+lambda(ix1,ix2,2,2))
        tcos=dcos(x(ix1,ix2,2))
        tsin=dsin(x(ix1,ix2,2))
        tctan(ix1,ix2)=tcos/tsin
        if(ndir==2) then
          lambda(ix1,ix2,1,1)=lambda(ix1,ix2,1,1)-divv23
          lambda(ix1,ix2,2,2)=lambda(ix1,ix2,2,2)-divv23
          if(energy) then
            vlambda(ix1,ix2,1)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,1)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,1)
            vlambda(ix1,ix2,2)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,2)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,2)
          end if
        else
          invrsin(ix1,ix2)=invr(ix1,ix2)/tsin
          ! idim=1, idir=3
          lambda(ix1,ix2,1,3)=(wp(ix1+1,ix2,v3_)-wp(ix1-1,ix2,v3_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
          ! idim=2, idir=3
          lambda(ix1,ix2,2,3)=((wp(ix1,ix2+1,v3_)-wp(ix1,ix2-1,v3_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+wp(ix1,ix2,v2_)*tcos)*invr(ix1,ix2)
          ! idim=3, idir=1
          lambda(ix1,ix2,3,1)=-wp(ix1,ix2,v3_)*tsin*invrsin(ix1,ix2)
          ! idim=3, idir=2
          lambda(ix1,ix2,3,2)=-wp(ix1,ix2,v3_)*tcos*invrsin(ix1,ix2)
          ! idim=3, idir=3
          lambda(ix1,ix2,3,3)=(wp(ix1,ix2,v1_)*tsin+wp(ix1,ix2,v2_)*tcos)*invrsin(ix1,ix2)*two*qdtmu
          lambda(ix1,ix2,1,3)=(lambda(ix1,ix2,1,3)+lambda(ix1,ix2,3,1))*qdtmu
          lambda(ix1,ix2,3,1)=lambda(ix1,ix2,1,3)
          lambda(ix1,ix2,2,3)=(lambda(ix1,ix2,2,3)+lambda(ix1,ix2,3,2))*qdtmu
          lambda(ix1,ix2,3,2)=lambda(ix1,ix2,2,3)
          lambda(ix1,ix2,1,1)=lambda(ix1,ix2,1,1)-divv23
          lambda(ix1,ix2,2,2)=lambda(ix1,ix2,2,2)-divv23
          lambda(ix1,ix2,3,3)=-divv23
          if(energy) then
            vlambda(ix1,ix2,1)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,1)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,1)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,1)
            vlambda(ix1,ix2,2)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,2)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,2)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,2)
            vlambda(ix1,ix2,3)=wp(ix1,ix2,v1_)*lambda(ix1,ix2,1,3)+wp(ix1,ix2,v2_)*lambda(ix1,ix2,2,3)+wp(ix1,ix2,v3_)*lambda(ix1,ix2,3,3)
          end if
        end if
      {end do\}
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      {do ix^db=ixomin^db,ixomax^db\}
        if(ndir==2) then
          w(ix1,ix2,mom(1))=(lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                            (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
                         two*lambda(ix1,ix2,1,1)*invr(ix1,ix2)+lambda(ix1,ix2,2,1)*invr(ix1,ix2)*tctan(ix1,ix2)+w(ix1,ix2,mom(1))
          w(ix1,ix2,mom(2))=(lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                         (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
                      two*lambda(ix1,ix2,1,2)*invr(ix1,ix2)+lambda(ix1,ix2,2,2)*invr(ix1,ix2)*tctan(ix1,ix2)+w(ix1,ix2,mom(2))
        else
          w(ix1,ix2,mom(1))=(lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                         (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
                      two*lambda(ix1,ix2,1,1)*invr(ix1,ix2)+lambda(ix1,ix2,2,1)*invr(ix1,ix2)*tctan(ix1,ix2)+w(ix1,ix2,mom(1))
          w(ix1,ix2,mom(2))=(lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                         (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
                      two*lambda(ix1,ix2,1,2)*invr(ix1,ix2)+lambda(ix1,ix2,2,2)*invr(ix1,ix2)*tctan(ix1,ix2)-&
                      lambda(ix1,ix2,3,3)*tctan(ix1,ix2)*invr(ix1,ix2)**2+w(ix1,ix2,mom(2))
          w(ix1,ix2,mom(3))=(lambda(ix1+1,ix2,1,3)-lambda(ix1-1,ix2,1,3))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                         (lambda(ix1,ix2+1,2,3)-lambda(ix1,ix2-1,2,3))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix1,ix2)+&
                      two*lambda(ix1,ix2,1,3)*invr(ix1,ix2)+lambda(ix1,ix2,2,3)*(invr(ix1,ix2)+invr(ix1,ix2)**2)*tctan(ix1,ix2)+w(ix1,ix2,mom(3))
        end if
      {end do\}
      }
      {^ifoned
      do ix1=ixmin1,ixmax1
        invr(ix1)=1.d0/x(ix1,1)
        ! idim=1, idir=1
        lambda(ix1,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))*two*qdtmu
        ! dv_i/d_j + dv_j/d_i
        divv23=third*lambda(ix1,1,1)
        lambda(ix1,1,1)=lambda(ix1,1,1)-divv23
        if(ndir==1) then
          if(energy) then
            vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)
          end if
        else if(ndir==2) then
          ! idim=1, idir=2
          lambda(ix1,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))
          lambda(ix1,1,2)=lambda(ix1,1,2)*qdtmu
          lambda(ix1,2,1)=lambda(ix1,1,2)
          lambda(ix1,2,2)=-divv23
          if(energy) then
            vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)+wp(ix1,v2_)*lambda(ix1,2,1)
            vlambda(ix1,2)=wp(ix1,v1_)*lambda(ix1,1,2)+wp(ix1,v2_)*lambda(ix1,2,2)
          end if
        else
          tcos=0.d0
          tsin=1.d0
          tctan(ix1)=tcos/tsin
          invrsin(ix1)=invr(ix1)/tsin
          ! idim=1, idir=3
          lambda(ix1,1,3)=(wp(ix1+1,v3_)-wp(ix1-1,v3_))/(x(ix1+1,1)-x(ix1-1,1))
          ! idim=3, idir=1
          lambda(ix1,3,1)=-wp(ix1,v3_)*tsin*invrsin(ix1)
          ! idim=3, idir=2
          lambda(ix1,3,2)=-wp(ix1,v3_)*tcos*invrsin(ix1)
          ! idim=3, idir=3
          lambda(ix1,3,3)=(wp(ix1,v1_)*tsin+wp(ix1,v2_)*tcos)*invrsin(ix1)*two*qdtmu
          lambda(ix1,1,3)=(lambda(ix1,1,3)+lambda(ix1,3,1))*qdtmu
          lambda(ix1,3,1)=lambda(ix1,1,3)
          lambda(ix1,3,2)=lambda(ix1,3,2)*qdtmu
          lambda(ix1,2,3)=lambda(ix1,3,2)
          lambda(ix1,3,3)=lambda(ix1,3,3)-divv23
          if(energy) then
            vlambda(ix1,1)=wp(ix1,v1_)*lambda(ix1,1,1)+wp(ix1,v2_)*lambda(ix1,2,1)+wp(ix1,v3_)*lambda(ix1,3,1)
            vlambda(ix1,2)=wp(ix1,v1_)*lambda(ix1,1,2)+wp(ix1,v2_)*lambda(ix1,2,2)+wp(ix1,v3_)*lambda(ix1,3,2)
            vlambda(ix1,3)=wp(ix1,v1_)*lambda(ix1,1,3)+wp(ix1,v2_)*lambda(ix1,2,3)+wp(ix1,v3_)*lambda(ix1,3,3)
          end if
        end if
      end do
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      do ix1=ixomin1,ixomax1
        if(ndir==1) then
          w(ix1,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
                      two*lambda(ix1,1,1)*invr(ix1)+w(ix1,mom(1))
        else if(ndir==2) then
          w(ix1,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
                      two*lambda(ix1,1,1)*invr(ix1)+lambda(ix1,2,1)*invr(ix1)*tctan(ix1)+w(ix1,mom(1))
          w(ix1,mom(2))=(lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
                      two*lambda(ix1,1,2)*invr(ix1)+lambda(ix1,2,2)*invr(ix1)*tctan(ix1)+w(ix1,mom(2))
        else
          w(ix1,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
                      two*lambda(ix1,1,1)*invr(ix1)+lambda(ix1,2,1)*invr(ix1)*tctan(ix1)+w(ix1,mom(1))
          w(ix1,mom(2))=(lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
                      two*lambda(ix1,1,2)*invr(ix1)+lambda(ix1,2,2)*invr(ix1)*tctan(ix1)-lambda(ix1,3,3)*tctan(ix1)*invr(ix1)**2+&
                      w(ix1,mom(2))
          w(ix1,mom(3))=(lambda(ix1+1,1,3)-lambda(ix1-1,1,3))/(x(ix1+1,1)-x(ix1-1,1))+&
                      two*lambda(ix1,1,3)*invr(ix1)+lambda(ix1,2,3)*(invr(ix1)+invr(ix1)**2)*tctan(ix1)+w(ix1,mom(3))
        end if
      end do
      }
      if(energy) then
        ! de/dt= +div(v.dot.[mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v *kr])
        ! thus e=e+d_i v_j tensor_ji
        call divvector(vlambda,ixi^l,ixo^l,invr)
        w(ixo^s,e_)=w(ixo^s,e_)+invr(ixo^s)
      end if
 
    end if
 
  end subroutine viscosity_add_source_sphere
 
  subroutine viscosity_add_source_cylinder(qdt,ixI^L,ixO^L,wCT,wp,w,x,&
       energy,qsourcesplit,active)
  ! Add viscosity source in isotropic Newtonian fluids to w within ixO
  ! neglecting bulk viscosity
  ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wp(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
 
    double precision:: lambda(ixI^S,ndir,ndir),vlambda(ixI^S,ndir),invr(ixI^S)
    double precision :: qdtmu,divv23
    integer:: ix^L,ix^D
 
    if (viscindiv) return
 
    if(qsourcesplit .eqv. vc_split) then
      active = .true.
      ! standard case, textbook viscosity
      ! Calculating viscosity sources
      if(.not.vc_4th_order) then
        ! involves second derivatives, two extra layers
        ix^l=ixo^l^ladd2;
        if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
          call mpistop("error for viscous source addition, 2 layers needed")
        ix^l=ixo^l^ladd1;
      else
        ! involves second derivatives, four extra layers
        ix^l=ixo^l^ladd4;
        if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
          call mpistop("error for viscous source addition"//&
          "requested fourth order gradients: 4 layers needed")
        ix^l=ixo^l^ladd2;
      end if
 
      ! construct lambda tensor: lambda_ij = gradv_ij + gradv_ji
      ! initialize
      qdtmu=qdt*vc_mu
      {^ifthreed
      {do ix^db=ixmin^db,ixmax^db\}
        invr(ix^d)=1.d0/x(ix^d,1)
        ! idim=1, idir=1
        lambda(ix^d,1,1)=(wp(ix1+1,ix2,ix3,v1_)-wp(ix1-1,ix2,ix3,v1_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))*two*qdtmu
        ! idim=1, idir=2
        lambda(ix^d,1,2)=(wp(ix1+1,ix2,ix3,v2_)-wp(ix1-1,ix2,ix3,v2_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
        ! idim=1, idir=3
        lambda(ix^d,1,3)=(wp(ix1+1,ix2,ix3,v3_)-wp(ix1-1,ix2,ix3,v3_))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))
        ! idim=2, idir=1
        lambda(ix^d,2,1)=((wp(ix1,ix2+1,ix3,v1_)-wp(ix1,ix2-1,ix3,v1_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))-wp(ix^d,v2_))*invr(ix^d)
        ! idim=2, idir=2
        lambda(ix^d,2,2)=((wp(ix1,ix2+1,ix3,v2_)-wp(ix1,ix2-1,ix3,v2_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))+wp(ix^d,v1_))*invr(ix^d)*&
                         two*qdtmu
        ! idim=2, idir=3
        lambda(ix^d,2,3)=(wp(ix1,ix2+1,ix3,v3_)-wp(ix1,ix2-1,ix3,v3_))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)
        ! idim=3, idir=1
        lambda(ix^d,3,1)=(wp(ix1,ix2,ix3+1,v1_)-wp(ix1,ix2,ix3-1,v1_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invr(ix^d)
        ! idim=3, idir=2
        lambda(ix^d,3,2)=(wp(ix1,ix2,ix3+1,v2_)-wp(ix1,ix2,ix3-1,v2_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invr(ix^d)
        ! idim=3, idir=3
        lambda(ix^d,3,3)=(wp(ix1,ix2,ix3+1,v3_)-wp(ix1,ix2,ix3-1,v3_))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))*invr(ix^d)*two*qdtmu
        ! dv_i/d_j + dv_j/d_i
        lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*qdtmu
        lambda(ix^d,2,1)=lambda(ix^d,1,2)
        lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*qdtmu
        lambda(ix^d,3,1)=lambda(ix^d,1,3)
        lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*qdtmu
        lambda(ix^d,3,2)=lambda(ix^d,2,3)
        divv23=third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
        lambda(ix^d,1,1)=lambda(ix^d,1,1)-divv23
        lambda(ix^d,2,2)=lambda(ix^d,2,2)-divv23
        lambda(ix^d,3,3)=lambda(ix^d,3,3)-divv23
        if(energy) then
          vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
          vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
          vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
        end if
      {end do\}
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      {do ix^db=ixomin^db,ixomax^db\}
        w(ix^d,mom(1))=(lambda(ix1+1,ix2,ix3,1,1)-lambda(ix1-1,ix2,ix3,1,1))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
                       (lambda(ix1,ix2+1,ix3,2,1)-lambda(ix1,ix2-1,ix3,2,1))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
                       (lambda(ix1,ix2,ix3+1,3,1)-lambda(ix1,ix2,ix3-1,3,1))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+&
                       (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d)+w(ix^d,mom(1))
        w(ix^d,mom(2))=(lambda(ix1+1,ix2,ix3,1,2)-lambda(ix1-1,ix2,ix3,1,2))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
                       (lambda(ix1,ix2+1,ix3,2,2)-lambda(ix1,ix2-1,ix3,2,2))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
                       (lambda(ix1,ix2,ix3+1,3,2)-lambda(ix1,ix2,ix3-1,3,2))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+&
                       two*lambda(ix^d,1,2)*invr(ix^d)+w(ix^d,mom(2))
        w(ix^d,mom(3))=(lambda(ix1+1,ix2,ix3,1,3)-lambda(ix1-1,ix2,ix3,1,3))/(x(ix1+1,ix2,ix3,1)-x(ix1-1,ix2,ix3,1))+&
                       (lambda(ix1,ix2+1,ix3,2,3)-lambda(ix1,ix2-1,ix3,2,3))/(x(ix1,ix2+1,ix3,2)-x(ix1,ix2-1,ix3,2))*invr(ix^d)+&
                       (lambda(ix1,ix2,ix3+1,3,3)-lambda(ix1,ix2,ix3-1,3,3))/(x(ix1,ix2,ix3+1,3)-x(ix1,ix2,ix3-1,3))+&
                        lambda(ix^d,1,3)*invr(ix^d)+w(ix^d,mom(3))
      {end do\}
      }
      {^iftwod
      {do ix^db=ixmin^db,ixmax^db\}
        invr(ix^d)=1.d0/x(ix^d,1)
        ! idim=1, idir=1
        lambda(ix^d,1,1)=(wp(ix1+1,ix2,v1_)-wp(ix1-1,ix2,v1_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))*two*qdtmu
        ! idim=1, idir=2
        lambda(ix^d,1,2)=(wp(ix1+1,ix2,v2_)-wp(ix1-1,ix2,v2_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
        if(phi_==2) then
          ! idim=2, idir=1
          lambda(ix^d,2,1)=((wp(ix1,ix2+1,v1_)-wp(ix1,ix2-1,v1_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))-wp(ix^d,v2_))*invr(ix^d)
          ! idim=2, idir=2
          lambda(ix^d,2,2)=((wp(ix1,ix2+1,v2_)-wp(ix1,ix2-1,v2_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+wp(ix^d,v1_))*invr(ix^d)*&
                           two*qdtmu
          if(ndir==3) then
            ! idim=2, idir=3
            lambda(ix^d,2,3)=(wp(ix1,ix2+1,v3_)-wp(ix1,ix2-1,v3_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)
            ! idim=3, idir=1
            lambda(ix^d,3,1)=0.d0
            ! idim=3, idir=2
            lambda(ix^d,3,2)=0.d0
            ! idim=3, idir=3
            lambda(ix^d,3,3)=0.d0
          end if
        else
          ! idim=2, idir=1
          lambda(ix^d,2,1)=(wp(ix1,ix2+1,v1_)-wp(ix1,ix2-1,v1_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)
          ! idim=2, idir=2
          lambda(ix^d,2,2)=(wp(ix1,ix2+1,v2_)-wp(ix1,ix2-1,v2_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)*two*qdtmu
          if(ndir==3) then
            ! idim=2, idir=3
            lambda(ix^d,2,3)=(wp(ix1,ix2+1,v3_)-wp(ix1,ix2-1,v3_))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)
            ! idim=3, idir=1
            lambda(ix^d,3,1)=-wp(ix^d,v3_)*invr(ix^d)
            ! idim=3, idir=2
            lambda(ix^d,3,2)=0.d0
            ! idim=3, idir=3
            lambda(ix^d,3,3)=wp(ix^d,v1_)*invr(ix^d)
          end if
        end if
        ! dv_i/d_j + dv_j/d_i
        lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*qdtmu
        lambda(ix^d,2,1)=lambda(ix^d,1,2)
        divv23=third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
        lambda(ix^d,1,1)=lambda(ix^d,1,1)-divv23
        lambda(ix^d,2,2)=lambda(ix^d,2,2)-divv23
        if(ndir==3) then
          ! idim=1, idir=3
          lambda(ix^d,1,3)=(wp(ix1+1,ix2,v3_)-wp(ix1-1,ix2,v3_))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))
          lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*qdtmu
          lambda(ix^d,3,1)=lambda(ix^d,1,3)
          lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*qdtmu
          lambda(ix^d,3,2)=lambda(ix^d,2,3)
          lambda(ix^d,3,3)=lambda(ix^d,3,3)-divv23
          if(energy) then
            vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
            vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
            vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
          end if
        else if(energy) then
          vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)
          vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)
        end if
      {end do\}
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      {do ix^db=ixomin^db,ixomax^db\}
        if(ndir==2) then
          if(phi_==2) then
            w(ix^d,mom(1))=(lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
                           (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d)+w(ix^d,mom(1))
            w(ix^d,mom(2))=(lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
                        two*lambda(ix^d,1,2)*invr(ix^d)+w(ix^d,mom(2))
          else
            w(ix^d,mom(1))=(lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
                            lambda(ix^d,1,1)*invr(ix^d)+w(ix^d,mom(1))
            w(ix^d,mom(2))=(lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
                            lambda(ix^d,1,2)*invr(ix^d)+w(ix^d,mom(2))
          end if
        else
          if(phi_==2) then
            w(ix^d,mom(1))=(lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
                           (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d)+w(ix^d,mom(1))
            w(ix^d,mom(2))=(lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
                           two*lambda(ix^d,1,2)*invr(ix^d)+w(ix^d,mom(2))
            w(ix^d,mom(3))=(lambda(ix1+1,ix2,1,3)-lambda(ix1-1,ix2,1,3))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,3)-lambda(ix1,ix2-1,2,3))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))*invr(ix^d)+&
                            lambda(ix^d,1,3)*invr(ix^d)+w(ix^d,mom(3))
          else
            w(ix^d,mom(1))=(lambda(ix1+1,ix2,1,1)-lambda(ix1-1,ix2,1,1))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,1)-lambda(ix1,ix2-1,2,1))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
                           (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d)+w(ix^d,mom(1))
            w(ix^d,mom(2))=(lambda(ix1+1,ix2,1,3)-lambda(ix1-1,ix2,1,3))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,3)-lambda(ix1,ix2-1,2,3))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
                           two*lambda(ix^d,1,3)*invr(ix^d)+w(ix^d,mom(2))
            w(ix^d,mom(3))=(lambda(ix1+1,ix2,1,2)-lambda(ix1-1,ix2,1,2))/(x(ix1+1,ix2,1)-x(ix1-1,ix2,1))+&
                           (lambda(ix1,ix2+1,2,2)-lambda(ix1,ix2-1,2,2))/(x(ix1,ix2+1,2)-x(ix1,ix2-1,2))+&
                            lambda(ix^d,1,2)*invr(ix^d)+w(ix^d,mom(3))
          end if
        end if
      {end do\}
      }
      {^ifoned
      {do ix^db=ixmin^db,ixmax^db\}
        invr(ix^d)=1.d0/x(ix^d,1)
        if(ndir==1) then
          ! idim=1, idir=1
          lambda(ix^d,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))*two*qdtmu
          divv23=third*lambda(ix^d,1,1)
          lambda(ix^d,1,1)=lambda(ix^d,1,1)-divv23
          if(energy) then
            vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)
          end if
        else if(ndir==2) then
          ! idim=1, idir=1
          lambda(ix^d,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))*two*qdtmu
          ! idim=1, idir=2
          lambda(ix^d,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))
          if(phi_==2) then
            ! idim=2, idir=1
            lambda(ix^d,2,1)=-wp(ix^d,v2_)*invr(ix^d)
            ! idim=2, idir=2
            lambda(ix^d,2,2)=+wp(ix^d,v1_)*invr(ix^d)*two*qdtmu
          else
            ! idim=2, idir=1
            lambda(ix^d,2,1)=0.d0
            ! idim=2, idir=2
            lambda(ix^d,2,2)=0.d0
          end if
          ! dv_i/d_j + dv_j/d_i
          lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*qdtmu
          lambda(ix^d,2,1)=lambda(ix^d,1,2)
          divv23=third*(lambda(ix^d,1,1)+lambda(ix^d,2,2))
          lambda(ix^d,1,1)=lambda(ix^d,1,1)-divv23
          lambda(ix^d,2,2)=lambda(ix^d,2,2)-divv23
          if(energy) then
            vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)
            vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)
          end if
        else
          ! idim=1, idir=1
          lambda(ix^d,1,1)=(wp(ix1+1,v1_)-wp(ix1-1,v1_))/(x(ix1+1,1)-x(ix1-1,1))*two*qdtmu
          ! idim=1, idir=2
          lambda(ix^d,1,2)=(wp(ix1+1,v2_)-wp(ix1-1,v2_))/(x(ix1+1,1)-x(ix1-1,1))
          ! idim=1, idir=3
          lambda(ix^d,1,3)=(wp(ix1+1,v3_)-wp(ix1-1,v3_))/(x(ix1+1,1)-x(ix1-1,1))
          if(phi_==2) then
            ! idim=2, idir=1
            lambda(ix^d,2,1)=-wp(ix^d,v2_)*invr(ix^d)
            ! idim=2, idir=2
            lambda(ix^d,2,2)=+wp(ix^d,v1_)*invr(ix^d)*two*qdtmu
            ! idim=2, idir=3
            lambda(ix^d,2,3)=0.d0
            ! idim=3, idir=1
            lambda(ix^d,3,1)=0.d0
            ! idim=3, idir=2
            lambda(ix^d,3,2)=0.d0
            ! idim=3, idir=3
            lambda(ix^d,3,3)=0.d0
          else
            ! idim=2, idir=1
            lambda(ix^d,2,1)=0.d0
            ! idim=2, idir=2
            lambda(ix^d,2,2)=0.d0
            ! idim=2, idir=3
            lambda(ix^d,2,3)=0.d0
            ! idim=3, idir=1
            lambda(ix^d,3,1)=-wp(ix^d,v3_)*invr(ix^d)
            ! idim=3, idir=2
            lambda(ix^d,3,2)=0.d0
            ! idim=3, idir=3
            lambda(ix^d,3,3)=wp(ix^d,v1_)*invr(ix^d)*two*qdtmu
          end if
          ! dv_i/d_j + dv_j/d_i
          lambda(ix^d,1,2)=(lambda(ix^d,1,2)+lambda(ix^d,2,1))*qdtmu
          lambda(ix^d,2,1)=lambda(ix^d,1,2)
          lambda(ix^d,1,3)=(lambda(ix^d,1,3)+lambda(ix^d,3,1))*qdtmu
          lambda(ix^d,3,1)=lambda(ix^d,1,3)
          lambda(ix^d,2,3)=(lambda(ix^d,2,3)+lambda(ix^d,3,2))*qdtmu
          lambda(ix^d,3,2)=lambda(ix^d,2,3)
          divv23=third*(lambda(ix^d,1,1)+lambda(ix^d,2,2)+lambda(ix^d,3,3))
          lambda(ix^d,1,1)=lambda(ix^d,1,1)-divv23
          lambda(ix^d,2,2)=lambda(ix^d,2,2)-divv23
          lambda(ix^d,3,3)=lambda(ix^d,3,3)-divv23
          if(energy) then
            vlambda(ix^d,1)=wp(ix^d,v1_)*lambda(ix^d,1,1)+wp(ix^d,v2_)*lambda(ix^d,2,1)+wp(ix^d,v3_)*lambda(ix^d,3,1)
            vlambda(ix^d,2)=wp(ix^d,v1_)*lambda(ix^d,1,2)+wp(ix^d,v2_)*lambda(ix^d,2,2)+wp(ix^d,v3_)*lambda(ix^d,3,2)
            vlambda(ix^d,3)=wp(ix^d,v1_)*lambda(ix^d,1,3)+wp(ix^d,v2_)*lambda(ix^d,2,3)+wp(ix^d,v3_)*lambda(ix^d,3,3)
          end if
        end if
      {end do\}
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      {do ix^db=ixomin^db,ixomax^db\}
        if(ndir==1) then
          w(ix^d,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+w(ix^d,mom(1))
        else if(ndir==2) then
          if(phi_==2) then
            w(ix^d,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
                           (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d)+w(ix^d,mom(1))
            w(ix^d,mom(2))=(lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
                           two*lambda(ix^d,1,2)*invr(ix^d)+w(ix^d,mom(2))
          else
            w(ix^d,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
                           (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d)+w(ix^d,mom(1))
          end if
        else
          if(phi_==2) then
            w(ix^d,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
                           (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d)+w(ix^d,mom(1))
            w(ix^d,mom(2))=(lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
                           two*lambda(ix^d,1,2)*invr(ix^d)+w(ix^d,mom(2))
            w(ix^d,mom(3))=(lambda(ix1+1,1,3)-lambda(ix1-1,1,3))/(x(ix1+1,1)-x(ix1-1,1))+&
                            lambda(ix^d,1,3)*invr(ix^d)+w(ix^d,mom(3))
          else
            w(ix^d,mom(1))=(lambda(ix1+1,1,1)-lambda(ix1-1,1,1))/(x(ix1+1,1)-x(ix1-1,1))+&
                           (lambda(ix^d,1,1)-lambda(ix^d,2,2))*invr(ix^d)+w(ix^d,mom(1))
            w(ix^d,mom(2))=(lambda(ix1+1,1,3)-lambda(ix1-1,1,3))/(x(ix1+1,1)-x(ix1-1,1))+&
                           two*lambda(ix^d,1,3)*invr(ix^d)+w(ix^d,mom(2))
            w(ix^d,mom(3))=(lambda(ix1+1,1,2)-lambda(ix1-1,1,2))/(x(ix1+1,1)-x(ix1-1,1))+&
                            lambda(ix^d,1,2)*invr(ix^d)+w(ix^d,mom(3))
          end if
        end if
      {end do\}
      }
      if(energy) then
        ! de/dt= +div(v.dot.[mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v *kr])
        ! thus e=e+d_i v_j tensor_ji
        call divvector(vlambda,ixi^l,ixo^l,invr)
        w(ixo^s,e_)=w(ixo^s,e_)+invr(ixo^s)
      end if
 
    end if
 
  end subroutine viscosity_add_source_cylinder
 
  subroutine viscosity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    ! Check diffusion time limit for dt < dtdiffpar * dx**2 / (mu/rho)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: dx^D, x(ixI^S,1:ndim)
    double precision, intent(in) :: w(ixI^S,1:nw)
    double precision, intent(inout) :: dtnew
 
    double precision :: tmp(ixI^S)
    double precision:: dtdiff_visc, dxinv2(1:ndim), max_mu
    integer:: idim
 
    ! Calculate the kinematic viscosity tmp=mu/rho
    ! here vc_mu must be non-zero!!!
 
    tmp(ixo^s)=vc_mu/w(ixo^s,rho_)
 
    if(slab_uniform)then
      ^d&dxinv2(^d)=one/dx^d**2;
      do idim=1,ndim
         dtdiff_visc=dtdiffpar/maxval(tmp(ixo^s)*dxinv2(idim))
         ! limit the time step
         dtnew=min(dtnew,dtdiff_visc)
      enddo
    else
      do idim=1,ndim
         max_mu=maxval(tmp(ixo^s)/block%ds(ixo^s,idim)**2)
         dtdiff_visc=dtdiffpar/max_mu
         ! limit the time step
         dtnew=min(dtnew,dtdiff_visc)
      enddo
    endif
 
  end subroutine viscosity_get_dt
 
  ! viscInDiv
  ! Get the viscous stress tensor terms in the idim direction
  ! Beware : a priori, won't work for ndir /= ndim
  ! Rq : we work with primitive w variables here
  ! Rq : ixO^L is already extended by 1 unit in the direction we work on
 
  subroutine visc_get_flux_prim(w, x, ixI^L, ixO^L, idim, f, energy)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)             :: ixi^l, ixo^l, idim
    double precision, intent(in)    :: w(ixi^s, 1:nw), x(ixi^s, 1:^nd)
    double precision, intent(inout) :: f(ixi^s, nwflux)
    logical, intent(in) :: energy
    integer                         :: idir, i
    double precision :: v(ixi^s,1:ndir)
 
    double precision                :: divergence(ixi^s)
 
    double precision:: lambda(ixi^s,ndir) !, tmp(ixI^S) !gradV(ixI^S,ndir,ndir)
 
    if (.not. viscindiv) return
 
    do i=1,ndir
     v(ixi^s,i)=w(ixi^s,i+1)
    enddo
    call divvector(v,ixi^l,ixo^l,divergence)
 
    call get_crossgrad(ixi^l,ixo^l,x,w,idim,lambda)
    lambda(ixo^s,idim) = lambda(ixo^s,idim) - (2.d0/3.d0) * divergence(ixo^s)
 
    ! Compute the idim-th row of the viscous stress tensor
    do idir = 1, ndir
      f(ixo^s, mom(idir)) = f(ixo^s, mom(idir)) - vc_mu * lambda(ixo^s,idir)
      if (energy) f(ixo^s, e_) = f(ixo^s, e_) - vc_mu * lambda(ixo^s,idir) * v(ixi^s,idir)
    enddo
 
  end subroutine visc_get_flux_prim
 
  ! Compute the cross term ( d_i v_j + d_j v_i in Cartesian BUT NOT IN
  ! CYLINDRICAL AND SPHERICAL )
  subroutine get_crossgrad(ixI^L,ixO^L,x,w,idim,cross)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: w(ixI^S, 1:nw), x(ixI^S, 1:ndim)
    double precision, intent(out)   :: cross(ixI^S,ndir)
 
    double precision :: tmp(ixI^S), v(ixI^S)
    integer :: idir
 
    if (ndir/=ndim) call mpistop("This formula are probably wrong for ndim/=ndir")
    ! Beware also, we work w/ the angle as the 3rd component in cylindrical
    ! and the colatitude as the 2nd one in spherical
    cross(ixi^s,:)=zero
    tmp(ixi^s)=zero
    select case(coordinate)
    case (cartesian,cartesian_stretched)
      call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
    case (cylindrical)
      if (idim==1) then
        ! for rr and rz
        call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
        ! then we overwrite rth w/ the correct expression
        {^nooned
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
        cross(ixi^s,2)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
        }
      elseif (idim==2) then
        ! thr (idem as above)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        {^nooned
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,1)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
        ! thth
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,2)=two*tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
        !thz
        v(ixi^s)=w(ixi^s,mom(3)) ! v_z
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        }
        cross(ixi^s,3)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))  ! v_th
        {^ifthreed
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_z
        }
        cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)
        {^ifthreed
      elseif (idim==3) then
        ! for zz and rz
        call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
        ! then we overwrite zth w/ the correct expression
        !thz
        v(ixi^s)=w(ixi^s,mom(3)) ! v_z
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,2)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))  ! v_th
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_z
        cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)
        }
      endif
    case (spherical)
      if (idim==1) then
        ! rr (normal, simply 2 * dr vr)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,1)=two*tmp(ixi^s)
        !rth
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        {^nooned
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
        cross(ixi^s,2)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
        }
        {^ifthreed
        !rph
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
        cross(ixi^s,3)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*x(ixi^s,1)
        }
      elseif (idim==2) then
        ! thr
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        {^nooned
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
        cross(ixi^s,1)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
        ! thth
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,2)=two*tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
        }
        {^ifthreed
        !thph
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
        cross(ixi^s,3)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*dsin(x(ixi^s,2))
        }
        {^ifthreed
      elseif (idim==3) then
        !phr
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
        cross(ixi^s,1)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
        !phth
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
        cross(ixi^s,2)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*dsin(x(ixi^s,2))
        !phph
        v(ixi^s)=w(ixi^s,mom(3)) ! v_ph
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
        cross(ixi^s,3)=two*tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        cross(ixi^s,3)=cross(ixi^s,3)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        cross(ixi^s,3)=cross(ixi^s,3)+two*v(ixi^s)/(x(ixi^s,1)*dtan(x(ixi^s,2))) ! + 2 vth/(rtan(th))
        }
      endif
    case default
      call mpistop("Unknown geometry specified")
    end select
 
  end subroutine get_crossgrad
 
  !> yields d_i v_j + d_j v_i for a given i, OK in Cartesian and for some
  !> tensor terms in cylindrical (rr & rz) and in spherical (rr)
  subroutine cart_cross_grad(ixI^L,ixO^L,x,w,idim,cross)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: w(ixI^S, 1:nw), x(ixI^S, 1:^ND)
    double precision, intent(out)   :: cross(ixI^S,ndir)
 
    double precision :: tmp(ixI^S), v(ixI^S)
    integer :: idir
 
    v(ixi^s)=w(ixi^s,mom(idim))
    do idir=1,ndir
      call gradient(v,ixi^l,ixo^l,idir,tmp)
      cross(ixo^s,idir)=tmp(ixo^s)
    enddo
    do idir=1,ndir
      v(ixi^s)=w(ixi^s,mom(idir))
      call gradient(v,ixi^l,ixo^l,idim,tmp)
      cross(ixo^s,idir)=cross(ixo^s,idir)+tmp(ixo^s)
    enddo
 
  end subroutine cart_cross_grad
 
  subroutine visc_add_source_geom(qdt, ixI^L, ixO^L, wCT, w, x)
    use mod_global_parameters
    use mod_geometry
    ! w and wCT conservative variables here
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt, x(ixI^S, 1:ndim)
    double precision, intent(inout) :: wCT(ixI^S, 1:nw), w(ixI^S, 1:nw)
    ! to change and to set as a parameter in the parfile once the possibility to
    ! solve the equations in an angular momentum conserving form has been
    ! implemented (change tvdlf.t eg)
    double precision :: vv(ixI^S), divergence(ixI^S)
    double precision :: tmp(ixI^S),tmp1(ixI^S)
    integer          :: i
 
    if (.not. viscindiv) return
 
    select case (coordinate)
    case (cylindrical)
      ! thth tensor term - - -
        ! 1st the cross grad term
{^nooned
      call gradient(wct(ixi^s,mom(2)),ixi^l,ixo^l,2,tmp1) ! d_th
      tmp(ixo^s)=two*(tmp1(ixo^s)+wct(ixo^s,mom(1))/x(ixo^s,1)) ! 2 ( d_th v_th / r + vr/r )
        ! 2nd the divergence
      call divvector(wct(ixi^s,mom(1:ndir)),ixi^l,ixo^l,divergence)
      tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
      ! s[mr]=-thth/radius
      w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
      ! rth tensor term - - -
      call gradient(wct(ixi^s,mom(1)),ixi^l,ixo^l,2,tmp) ! d_th
      vv(ixi^s)=wct(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
      call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
      tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
      ! s[mphi]=+rth/radius
      w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
}
    case (spherical)
      ! thth tensor term - - -
      ! 1st the cross grad term
{^nooned
      call gradient(wct(ixi^s,mom(2)),ixi^l,ixo^l,2,tmp1) ! d_th
      tmp(ixo^s)=two*(tmp1(ixo^s)+wct(ixo^s,mom(1))/x(ixo^s,1)) ! 2 ( 1/r * d_th v_th + vr/r )
      ! 2nd the divergence
      call divvector(wct(ixi^s,mom(1:ndir)),ixi^l,ixo^l,divergence)
      tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
      ! s[mr]=-thth/radius
      w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
}
      ! phiphi tensor term - - -
      ! 1st the cross grad term
{^ifthreed
      call gradient(wct(ixi^s,mom(3)),ixi^l,ixo^l,3,tmp1) ! d_phi
      tmp(ixo^s)=two*(tmp1(ixo^s)+wct(ixo^s,mom(1))/x(ixo^s,1)+wct(ixo^s,mom(2))/(x(ixo^s,1)*dtan(x(ixo^s,2)))) ! 2 ( 1/rsinth * d_ph v_ph + vr/r + vth/rtanth )
}
      ! 2nd the divergence
      tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
      ! s[mr]=-phiphi/radius
      w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
      ! s[mth]=-cotanth*phiphi/radius
{^nooned
      w(ixo^s,mom(2))=w(ixo^s,mom(2))-qdt*vc_mu*tmp(ixo^s)/(x(ixo^s,1)*dtan(x(ixo^s,2)))
}
      ! rth tensor term - - -
      call gradient(wct(ixi^s,mom(1)),ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
      vv(ixi^s)=wct(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
      call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
      tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
      ! s[mth]=+rth/radius
      w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
      ! rphi tensor term - - -
{^ifthreed
      call gradient(wct(ixi^s,mom(1)),ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
}
      vv(ixi^s)=wct(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
      call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
      tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
      ! s[mphi]=+rphi/radius
      w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
      ! phith tensor term - - -
{^ifthreed
      call gradient(wct(ixi^s,mom(2)),ixi^l,ixo^l,3,tmp) ! d_phi
}
{^nooned
      vv(ixi^s)=wct(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
      call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
      tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*dsin(x(ixo^s,2))
      ! s[mphi]=+cotanth*phith/radius
      w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*vc_mu*tmp(ixo^s)/(x(ixo^s,1)*dtan(x(ixo^s,2)))
}
    end select
 
  end subroutine visc_add_source_geom
 
  subroutine sub_add_source(qdt,ixI^L,ixO^L,wCT,wp,w,x,&
       energy,qsourcesplit,active)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wp(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
  end subroutine sub_add_source
 
end module mod_viscosity