mod_geometry.t

Go to the documentation of this file.

!> Module with geometry-related routines (e.g., divergence, curl)
module mod_geometry
  use mod_comm_lib, only: mpistop
  implicit none
  public
 
  integer :: coordinate=-1
  integer, parameter :: cartesian          = 0
  integer, parameter :: cartesian_stretched= 1
  integer, parameter :: cylindrical        = 2
  integer, parameter :: spherical          = 3
  integer, parameter :: cartesian_expansion= 4
 
  integer :: type_curl=0
  integer, parameter :: central=1
  integer, parameter :: gaussbased=2
  integer, parameter :: stokesbased=3
 
contains
 
  !> Set the coordinate system to be used
  subroutine set_coordinate_system(geom)
    use mod_global_parameters
    character(len=*), intent(in) :: geom !< Name of the coordinate system
 
    ! Store the geometry name
    geometry_name = geom
 
    select case (geom)
    case ("Cartesian","Cartesian_1D","Cartesian_2D","Cartesian_3D")
      ndir = ndim
      coordinate=cartesian
    case ("Cartesian_1D_expansion")
      if (ndim /= 1) call mpistop("Geometry Cartesian_1D_expansion but ndim /= 1")
      ndir = ndim
      coordinate=cartesian_expansion
    case ("Cartesian_1.5D")
      if (ndim /= 1) call mpistop("Geometry Cartesian_1.5D but ndim /= 1")
      coordinate=cartesian
      ndir = 2
    case ("Cartesian_1.75D")
      if (ndim /= 1) call mpistop("Geometry Cartesian_1.75D but ndim /= 1")
      coordinate=cartesian
      ndir = 3
    case ("Cartesian_2.5D")
      if (ndim /= 2) call mpistop("Geometry Cartesian_2.5D but ndim /= 2")
      coordinate=cartesian
      ndir = 3
    case ("cylindrical","cylindrical_2D","cylindrical_3D")
      ndir = ndim
      r_   = 1
      z_   = 2
      if(ndir==3) phi_ = 3
      coordinate=cylindrical
    case ("cylindrical_2.5D")
      if (ndim /= 2) call mpistop("Geometry cylindrical_2.5D but ndim /= 2")
      ndir = 3
      r_   = 1
      z_   = 2
      phi_ = 3
      coordinate=cylindrical
    case ("polar","polar_2D","polar_3D")
      ndir = ndim
      r_   = 1
      phi_ = 2
      if(ndir==3) z_ = 3
      coordinate=cylindrical
    case ("polar_1.5D")
       if (ndim /= 1) call mpistop("Geometry polar_1.5D but ndim /= 1")
       ndir = 2
       r_   = 1
       phi_ = 2
       coordinate=cylindrical
    case ("polar_2.5D")
      if (ndim /= 2) call mpistop("Geometry polar_2.5D but ndim /= 2")
      ndir = 3
      r_   = 1
      phi_ = 2
      z_   = 3
      coordinate=cylindrical
    case ("spherical","spherical_2D","spherical_3D")
      ndir = ndim
      r_   = 1
      if(ndir==3) phi_ = 3
      z_   = -1
      coordinate=spherical
    case ("spherical_2.5D")
      if (ndim /= 2) &
           call mpistop("Geometry spherical_2.5D requires ndim == 2")
      ndir = 3
      r_   = 1
      phi_ = 3
      z_   = -1
      coordinate=spherical
    case default
      call mpistop("Unknown geometry specified")
    end select
  end subroutine set_coordinate_system
 
  subroutine set_pole
    use mod_global_parameters
 
    select case (coordinate)
    case (spherical) {^ifthreed
      ! For spherical grid, check whether phi-direction is periodic
      if(periodb(ndim)) then
        if(phi_/=3) call mpistop("phi_ should be 3 in 3D spherical coord!")
        if(mod(ng3(1),2)/=0) &
             call mpistop("Number of meshes in phi-direction should be even!")
        if(abs(xprobmin2)<smalldouble) then
          if(mype==0) write(unitterm,*) &
               "Will apply pi-periodic conditions at northpole!"
          poleb(1,2)=.true.
        else
          if(mype==0) write(unitterm,*) "There is no northpole!"
        end if
        if(abs(xprobmax2-dpi)<smalldouble) then
          if(mype==0) write(unitterm,*) &
               "Will apply pi-periodic conditions at southpole!"
          poleb(2,2)=.true.
        else
          if(mype==0) write(unitterm,*) "There is no southpole!"
        end if
      end if}
    case (cylindrical)
      {
      if (^d == phi_ .and. periodb(^d)) then
        if(mod(ng^d(1),2)/=0) then
          call mpistop("Number of meshes in phi-direction should be even!")
        end if
 
        if(abs(xprobmin1)<smalldouble) then
          if (mype==0) then
            write(unitterm,*) "Will apply pi-periodic conditions at r=0"
          end if
          poleb(1,1)=.true.
        else
          if (mype==0) then
            write(unitterm,*) "There is no cylindrical axis!"
          end if
        end if
      end if\}
    end select
 
  end subroutine set_pole
 
  !> Deallocate geometry-related variables
  subroutine putgridgeo(igrid)
    use mod_global_parameters
    integer, intent(in) :: igrid
 
    deallocate(ps(igrid)%surfaceC,ps(igrid)%surface,ps(igrid)%dvolume,ps(igrid)%dx,&
         psc(igrid)%dx,ps(igrid)%ds,psc(igrid)%dvolume)
 
  end subroutine putgridgeo
 
  !> calculate area of surfaces of cells
  subroutine get_surface_area(s,ixG^L)
    use mod_global_parameters
    use mod_usr_methods, only: usr_set_surface
 
    type(state) :: s
    integer, intent(in) :: ixG^L
 
    double precision :: x(ixG^S,ndim), xext(ixGmin1-1:ixGmax1,1), drs(ixG^S), drs_ext(ixGmin1-1:ixGmax1), dx2(ixG^S), dx3(ixG^S)
    double precision :: exp_factor_ext(ixGmin1-1:ixGmax1),del_exp_factor_ext(ixGmin1-1:ixGmax1),exp_factor_primitive_ext(ixGmin1-1:ixGmax1)
    double precision :: exp_factor(ixGmin1:ixGmax1),del_exp_factor(ixGmin1:ixGmax1),exp_factor_primitive(ixGmin1:ixGmax1)
 
    select case (coordinate)
 
    case (cartesian_expansion)
      drs(ixg^s)=s%dx(ixg^s,1)
      x(ixg^s,1)=s%x(ixg^s,1)
      {^ifoned
      if(associated(usr_set_surface))then
           call usr_set_surface(ixgmin1,ixgmax1,x,drs,exp_factor,del_exp_factor,exp_factor_primitive)
           if (any(exp_factor <= zero)) call mpistop("The area must always be strictly positive!")
      endif
      s%surface(ixg^s,1)=exp_factor(ixg^s)
      xext(ixgmin1-1,1)=x(1,1)-half*drs(1)
      xext(ixg^s,1)=x(ixg^s,1)+half*drs(ixg^s)
      drs_ext(ixgmin1-1)=drs(1)
      drs_ext(ixg^s)=drs(ixg^s)
      if(associated(usr_set_surface)) call usr_set_surface(ixgmin1-1,ixgmax1,xext,drs_ext,exp_factor_ext,del_exp_factor_ext,exp_factor_primitive_ext)
      s%surfaceC(ixgmin1-1:ixgmax1,1)=exp_factor_ext(ixgmin1-1:ixgmax1)
      }
 
    case (cartesian,cartesian_stretched)
      drs(ixg^s)=s%dx(ixg^s,1)
      {^nooned
      dx2(ixg^s)=s%dx(ixg^s,2)}
      {^ifthreed
      dx3(ixg^s)=s%dx(ixg^s,3)}
 
      {^ifoned
      s%surfaceC(ixg^s,1)=1.d0
      s%surface(ixg^s,1) =1.d0
      }
      {^iftwod
      s%surfaceC(ixg^s,1)=dx2(ixg^s)
      s%surfaceC(ixg^s,2)=drs(ixg^s)
      s%surface(ixg^s,1) =dx2(ixg^s)
      s%surface(ixg^s,2)=drs(ixg^s)
      }
      {^ifthreed
      s%surfaceC(ixg^s,1)= dx2(ixg^s)*dx3(ixg^s)
      s%surfaceC(ixg^s,2)= drs(ixg^s)*dx3(ixg^s)
      s%surfaceC(ixg^s,3)= drs(ixg^s)*dx2(ixg^s)
      s%surface(ixg^s,1)=s%surfaceC(ixg^s,1)
      s%surface(ixg^s,2)=s%surfaceC(ixg^s,2)
      s%surface(ixg^s,3)=s%surfaceC(ixg^s,3)
      }
      {s%surfaceC(0^d%ixG^s,^d)=s%surfaceC(1^d%ixG^s,^d);\}
    case (spherical)
      x(ixg^s,1)=s%x(ixg^s,1)
      {^nooned
      x(ixg^s,2)=s%x(ixg^s,2)
      }
      drs(ixg^s)=s%dx(ixg^s,1)
      {^nooned
      dx2(ixg^s)=s%dx(ixg^s,2)
      }
      {^ifthreed
      dx3(ixg^s)=s%dx(ixg^s,3)
      }
 
      s%surfaceC(ixg^s,1)=(x(ixg^s,1)+half*drs(ixg^s))**2 {^nooned &
           *two*dabs(dsin(x(ixg^s,2)))*dsin(half*dx2(ixg^s))}{^ifthreed*dx3(ixg^s)}
      ! change negative theta to positive to preserve positive area near pole
 
      {^nooned
      s%surfaceC(ixg^s,2)=x(ixg^s,1)*drs(ixg^s)&
           *dabs(dsin(x(ixg^s,2)+half*dx2(ixg^s)))}{^ifthreed*dx3(ixg^s)}
 
      {^ifthreed
      s%surfaceC(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)
      }
 
      {^ifoned
      s%surfaceC(0,1)=(x(1,1)-half*drs(1))**2
      }
      {^iftwod
      s%surfaceC(0,ixgmin2:ixgmax2,1)=(x(1,ixgmin2:ixgmax2,1)-half*drs(1,&
         ixgmin2:ixgmax2))**2*two*dabs(dsin(x(1,ixgmin2:ixgmax2,2)))*dsin(half*dx2(1,&
         ixgmin2:ixgmax2))
      s%surfaceC(ixgmin1:ixgmax1,0,2)=x(ixgmin1:ixgmax1,1,&
         1)*drs(ixgmin1:ixgmax1,1)*dabs(dsin(x(ixgmin1:ixgmax1,1,&
         2)-half*dx2(ixgmin1:ixgmax1,1)))
      }
      {^ifthreed
      s%surfaceC(0,ixgmin2:ixgmax2,ixgmin3:ixgmax3,1)=(x(1,ixgmin2:ixgmax2,&
         ixgmin3:ixgmax3,1)-half*drs(1,ixgmin2:ixgmax2,&
         ixgmin3:ixgmax3))**2*two*dabs(dsin(x(1,ixgmin2:ixgmax2,ixgmin3:ixgmax3,&
         2)))*dsin(half*dx2(1,ixgmin2:ixgmax2,ixgmin3:ixgmax3))*dx3(1,&
         ixgmin2:ixgmax2,ixgmin3:ixgmax3)
      s%surfaceC(ixgmin1:ixgmax1,0,ixgmin3:ixgmax3,2)=x(ixgmin1:ixgmax1,1,&
         ixgmin3:ixgmax3,1)*drs(ixgmin1:ixgmax1,1,&
         ixgmin3:ixgmax3)*dabs(dsin(x(ixgmin1:ixgmax1,1,ixgmin3:ixgmax3,&
         2)-half*dx2(ixgmin1:ixgmax1,1,ixgmin3:ixgmax3)))*dx3(ixgmin1:ixgmax1,1,&
         ixgmin3:ixgmax3)
      s%surfaceC(ixgmin1:ixgmax1,ixgmin2:ixgmax2,0,3)=&
         s%surfaceC(ixgmin1:ixgmax1,ixgmin2:ixgmax2,1,3)
      }
 
      s%surface(ixg^s,1)=x(ixg^s,1)**2 {^nooned &
           *two*dabs(dsin(x(ixg^s,2)))*dsin(half*dx2(ixg^s))}{^ifthreed*dx3(ixg^s)}
      {^nooned
      s%surface(ixg^s,2)=x(ixg^s,1)*drs(ixg^s)&
           *dabs(dsin(x(ixg^s,2)))}{^ifthreed*dx3(ixg^s)}
 
      {^ifthreed
      s%surface(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)}
 
    case (cylindrical)
      x(ixg^s,1)=s%x(ixg^s,1)
      drs(ixg^s)=s%dx(ixg^s,1)
      {^nooned
      dx2(ixg^s)=s%dx(ixg^s,2)}
      {^ifthreed
      dx3(ixg^s)=s%dx(ixg^s,3)}
 
      s%surfaceC(ixg^s,1)=dabs(x(ixg^s,1)+half*drs(ixg^s)){^de&*dx^de(ixg^s) }
      ! change negative r to positive to preserve positive area near pole
      {^nooned
      if (z_==2) s%surfaceC(ixg^s,2)=dabs(x(ixg^s,1))*drs(ixg^s){^ifthreed*dx3(ixg^s)}
      if (phi_==2) s%surfaceC(ixg^s,2)=drs(ixg^s){^ifthreed*dx3(ixg^s)}
      }
      {^ifthreed
      if (z_==3) s%surfaceC(ixg^s,3)=dabs(x(ixg^s,1))*drs(ixg^s)*dx2(ixg^s)
      if (phi_==3) s%surfaceC(ixg^s,3)=drs(ixg^s)*dx2(ixg^s)
      }
      {^ifoned
      s%surfaceC(0,1)=dabs(x(1,1)-half*drs(1))
      }
      {^iftwod
      s%surfaceC(0,ixgmin2:ixgmax2,1)=dabs(x(1,ixgmin2:ixgmax2,1)-half*drs(1,&
         ixgmin2:ixgmax2))*dx2(1,ixgmin2:ixgmax2)
      }
      {^ifthreed
      s%surfaceC(0,ixgmin2:ixgmax2,ixgmin3:ixgmax3,1)=dabs(x(1,ixgmin2:ixgmax2,&
         ixgmin3:ixgmax3,1)-half*drs(1,ixgmin2:ixgmax2,ixgmin3:ixgmax3))*dx2(1,&
         ixgmin2:ixgmax2,ixgmin3:ixgmax3)*dx3(1,ixgmin2:ixgmax2,&
         ixgmin3:ixgmax3)
      }
      {s%surfaceC(0^de%ixG^s,^de)=s%surfaceC(1^de%ixG^s,^de);\}
 
      s%surface(ixg^s,1)=dabs(x(ixg^s,1)){^de&*dx^de(ixg^s) }
      {^nooned
      if (z_==2) s%surface(ixg^s,2)=dabs(x(ixg^s,1))*drs(ixg^s){^ifthreed*dx3(ixg^s)}
      if (phi_==2) s%surface(ixg^s,2)=drs(ixg^s){^ifthreed*dx3(ixg^s)}}
      {^ifthreed
      if (z_==3) s%surface(ixg^s,3)=dabs(x(ixg^s,1))*drs(ixg^s)*dx2(ixg^s)
      if (phi_==3) s%surface(ixg^s,3)=drs(ixg^s)*dx2(ixg^s)}
 
    case default
      call mpistop("Sorry, coordinate unknown")
    end select
  end subroutine get_surface_area
 
  !**************************************************************************
  ! Purpose: Computes the gradient of a scalar field q(ixI^S) at cell 
  !          centers in the idir direction and outputs it in gradq(ixO^S), 
  !          which is also defined at cell **centers**. The gradient is 
  !          computed using a 2nd-order central difference scheme, 
  !          utilizing cell center values on both sides. 
  !
  !          For uniform Cartesian coordinates, an optional input 
  !          parameter nth_in allows increasing the order of the 
  !          central difference scheme to 2*nth_in.
  !**************************************************************************
  subroutine gradient(q,ixI^L,ixO^L,idir,gradq,nth_in)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L, idir
    integer, intent(in), optional   :: nth_in
    double precision, intent(in)    :: q(ixI^S)
    double precision, intent(inout) :: gradq(ixI^S)
    integer                         :: jxO^L, hxO^L, nth
 
    if(present(nth_in)) then
      nth = nth_in
    else
      nth = 1
    endif
    if(nth .gt. nghostcells) then 
      call mpistop("gradient stencil too wide")
    endif
 
    hxo^l=ixo^l-kr(idir,^d);
    jxo^l=ixo^l+kr(idir,^d);
    select case(coordinate)
    case(cartesian)
      select case(nth)
      case(1)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))*half
      case(2)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))*2.d0/3.d0
        hxo^l=ixo^l-2*kr(idir,^d);
        jxo^l=ixo^l+2*kr(idir,^d);
        gradq(ixo^s)=gradq(ixo^s)-(q(jxo^s)-q(hxo^s))/12.d0
      case(3)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))*3.d0/4.d0
        hxo^l=ixo^l-2*kr(idir,^d);
        jxo^l=ixo^l+2*kr(idir,^d);
        gradq(ixo^s)=gradq(ixo^s)-(q(jxo^s)-q(hxo^s))*3.d0/20.d0
        hxo^l=ixo^l-3*kr(idir,^d);
        jxo^l=ixo^l+3*kr(idir,^d);
        gradq(ixo^s)=gradq(ixo^s)+(q(jxo^s)-q(hxo^s))/60.d0
      case default
        call mpistop("unknown stencil gradient")
      end select
      gradq(ixo^s)=gradq(ixo^s)/dxlevel(idir)
    case(cartesian_stretched,cartesian_expansion)
      gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(block%x(jxo^s,idir)-block%x(hxo^s,idir))
    case(spherical)
      select case(idir)
      case(1)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((block%x(jxo^s,1)-block%x(hxo^s,1)))
        {^nooned
      case(2)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
        }
        {^ifthreed
      case(3)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1)*dsin(block%x(ixo^s,2)))
        }
      end select
    case(cylindrical)
      if(idir==phi_) then
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((block%x(jxo^s,phi_)-block%x(hxo^s,phi_))*block%x(ixo^s,r_))
      else
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(block%x(jxo^s,idir)-block%x(hxo^s,idir))
      end if
    case default
      call mpistop('Unknown geometry')
    end select
  end subroutine gradient
 
  !**************************************************************************
  ! Purpose: Originally named as gradientS.
  !          Computes the gradient of a scalar field q(ixI^S) at cell 
  !          centers in the idir direction and outputs it in gradq(ixO^S), 
  !          which is also defined at cell **centers**. Unlike the standard 
  !          gradient computation, this routine first reconstructs cell 
  !          face values using a slope **Limiter** and then computes the 
  !          gradient from these face values.
  !**************************************************************************
  subroutine gradientl(q,ixI^L,ixO^L,idir,gradq)
    use mod_global_parameters
    use mod_limiter
    use mod_ppm
    integer, intent(in)                :: ixI^L, ixO^L, idir
    double precision, intent(in)       :: q(ixI^S)
    double precision, intent(inout)    :: gradq(ixI^S)
    double precision ,dimension(ixI^S) :: qC,qL,qR,dqC,ldq,rdq
    double precision :: x(ixI^S,1:ndim)
    double precision :: invdx
    integer          :: hxO^L,ixC^L,jxC^L,gxC^L,hxC^L
 
    x(ixi^s,1:ndim)=block%x(ixi^s,1:ndim)
    invdx=1.d0/dxlevel(idir)
    hxo^l=ixo^l-kr(idir,^d);
    ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
    jxc^l=ixc^l+kr(idir,^d);
    gxcmin^d=ixcmin^d-kr(idir,^d);gxcmax^d=jxcmax^d;
    hxc^l=gxc^l+kr(idir,^d);
    ! set the gradient limiter here
    qr(gxc^s) = q(hxc^s)
    ql(gxc^s) = q(gxc^s)
    if(type_gradient_limiter(block%level)/=limiter_ppm) then
      dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
      call dwlimiter2(dqc,ixi^l,gxc^l,idir,type_gradient_limiter(block%level),ldw=ldq,rdw=rdq)
      ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
      qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
    else
      call ppmlimitervar(ixi^l,ixo^l,idir,q,q,ql,qr)
    endif
 
    select case(coordinate)
    case(cartesian)
      gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))*invdx
    case(cartesian_stretched,cartesian_expansion)
      gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/block%dx(ixo^s,idir)
    case(spherical)
      gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/block%dx(ixo^s,idir)
      select case(idir)
      case(2)
        gradq(ixo^s)=gradq(ixo^s)/x(ixo^s,1)
        {^ifthreed
      case(3)
        gradq(ixo^s)=gradq(ixo^s)/(x(ixo^s,1)*dsin(x(ixo^s,2)))
        }
      end select
    case(cylindrical)
      gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/block%dx(ixo^s,idir)
      if(idir==phi_) gradq(ixo^s)=gradq(ixo^s)/x(ixo^s,1)
    end select
  end subroutine gradientl
 
  !**************************************************************************
  ! Purpose: Merged from gradientx, gradientq and gradientC. 
  !          Computes the gradient of a scalar field q(ixI^S) at cell 
  !          centers in the idir direction and outputs it in gradq(ixO^S), 
  !          which is defined at cell **Faces**. The gradient is computed 
  !          using a 2nd-order central difference scheme, utilizing 
  !          the cell center values on both sides of a given cell face.
  !
  !          For uniform Cartesian coordinates, an optional input 
  !          parameter nth_in allows increasing the order of the 
  !          central difference scheme to 2*nth_in.
  !
  !          By default, the index is shifted one position to the left. 
  !          This structure also allows the routine to be used for 
  !          one-sided difference calculations of the gradient at cell 
  !          **centers**. In such cases, the optional logical parameter pm_in
  !          controls whether the index shifts to the left (default) 
  !          or to the right.
  !**************************************************************************
  subroutine gradientf(q,x,ixI^L,ixO^L,idir,gradq,nth_in,pm_in)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L, idir
    double precision, intent(in)    :: q(ixI^S),x(ixI^S,1:ndim)
    double precision, intent(inout) :: gradq(ixI^S)
    integer, intent(in), optional   :: nth_in
    logical, intent(in), optional   :: pm_in
    logical                         :: pm
    integer                         :: nth,jxO^L,hxO^L
 
    if(present(nth_in)) then
      nth = nth_in
    else
      nth = 1
    endif
    if(nth .gt. nghostcells) then
      call mpistop("gradient stencil too wide")
    endif
    if(present(pm_in)) then
      pm = pm_in
    else
      pm = .true.
    endif
 
    if(pm) then
      hxo^l=ixo^l;
      jxo^l=ixo^l+kr(idir,^d);
    else
      hxo^l=ixo^l-kr(idir,^d);
      jxo^l=ixo^l;
    endif
    select case(coordinate)
    case(cartesian)
      select case(nth)
      case(1)
        gradq(ixo^s)=q(jxo^s)-q(hxo^s)
      case(2)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))*27.d0/24.d0
        if(pm) then
          hxo^l=ixo^l-kr(idir,^d);
          jxo^l=ixo^l+2*kr(idir,^d);
        else
          hxo^l=ixo^l-2*kr(idir,^d);
          jxo^l=ixo^l+kr(idir,^d);
        endif
        gradq(ixo^s)=gradq(ixo^s)-(q(jxo^s)-q(hxo^s))/24.d0
      case(3)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))*225.d0/192.d0
        if(pm) then
          hxo^l=ixo^l-kr(idir,^d);
          jxo^l=ixo^l+2*kr(idir,^d);
        else
          hxo^l=ixo^l-2*kr(idir,^d);
          jxo^l=ixo^l+kr(idir,^d);
        endif
        gradq(ixo^s)=gradq(ixo^s)-(q(jxo^s)-q(hxo^s))*25.d0/384.d0
        if(pm) then
          hxo^l=ixo^l-2*kr(idir,^d);
          jxo^l=ixo^l+3*kr(idir,^d);
        else
          hxo^l=ixo^l-3*kr(idir,^d);
          jxo^l=ixo^l+2*kr(idir,^d);
        endif
        gradq(ixo^s)=gradq(ixo^s)+(q(jxo^s)-q(hxo^s))*3.0/640.d0
      case default
        call mpistop("unknown stencil gradient")
      end select
      gradq(ixo^s)=gradq(ixo^s)/dxlevel(idir)
    case(cartesian_stretched)
      gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
    case(spherical)
      select case(idir)
      case(1)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,1)-x(hxo^s,1))
        {^nooned
      case(2)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,2)-x(hxo^s,2))*x(ixo^s,1))
        }
        {^ifthreed
      case(3)
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,3)-x(hxo^s,3))*x(ixo^s,1)*dsin(x(ixo^s,2)))
        }
      end select
    case(cylindrical)
      if(idir==phi_) then
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,phi_)-x(hxo^s,phi_))*x(ixo^s,r_))
      else
        gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
      end if
    case default
      call mpistop('Unknown geometry')
    end select
  end subroutine gradientf
 
  !**************************************************************************
  ! Purpose: Computes the divergence of a vector field qvec(ixI^S,1:ndim) 
  !          at cell centers and stores the result in q(ixO^S), which is 
  !          also defined at cell centers. The divergence is computed using 
  !          a 2nd-order central difference scheme, utilizing the vector 
  !          components at adjacent cell centers.
  !
  !          For uniform Cartesian coordinates, an optional input parameter 
  !          nth_in allows increasing the order of the central difference 
  !          scheme to 2*nth_in, improving the accuracy of the divergence 
  !          computation.
  !**************************************************************************
  subroutine divvector(qvec,ixI^L,ixO^L,divq,nth_in)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L,ixO^L
    double precision, intent(in)    :: qvec(ixI^S,1:ndir)
    double precision, intent(inout) :: divq(ixI^S)
    integer, intent(in), optional   :: nth_in
    double precision                :: qC(ixI^S), invdx(1:ndim)
    integer                         :: jxO^L, hxO^L, ixC^L, jxC^L
    integer                         :: idims, gxO^L, kxO^L
    integer                         :: lxO^L, fxO^L, nth
 
    if(present(nth_in)) then
      nth = nth_in
    else
      nth = 1
    endif
    if(nth .gt. nghostcells) then
      call mpistop("divvector stencil too wide")
    endif
 
    invdx=1.d0/dxlevel
    divq(ixo^s)=0.0d0
    if (slab_uniform) then
      do idims=1,ndim
        select case(nth)
        case(1)
          jxo^l=ixo^l+kr(idims,^d);
          hxo^l=ixo^l-kr(idims,^d);
          divq(ixo^s)=divq(ixo^s)+half*(qvec(jxo^s,idims) &
               - qvec(hxo^s,idims))*invdx(idims)
        case(2)
          kxo^l=ixo^l+2*kr(idims,^d);
          jxo^l=ixo^l+kr(idims,^d);
          hxo^l=ixo^l-kr(idims,^d);
          gxo^l=ixo^l-2*kr(idims,^d);
          divq(ixo^s)=divq(ixo^s)+&
                    (-qvec(kxo^s,idims)+8.d0*qvec(jxo^s,idims)-8.d0*&
                      qvec(hxo^s,idims)+qvec(gxo^s,idims))/(12.d0*dxlevel(idims))
        case(3)
          lxo^l=ixo^l+3*kr(idims,^d);
          kxo^l=ixo^l+2*kr(idims,^d);
          jxo^l=ixo^l+kr(idims,^d);
          hxo^l=ixo^l-kr(idims,^d);
          gxo^l=ixo^l-2*kr(idims,^d);
          fxo^l=ixo^l-3*kr(idims,^d);
          divq(ixo^s)=divq(ixo^s)+&
                    (-qvec(fxo^s,idims)+9.d0*qvec(gxo^s,idims)-45.d0*qvec(hxo^s,idims)&
                     +qvec(lxo^s,idims)-9.d0*qvec(kxo^s,idims)+45.d0*qvec(jxo^s,idims))&
                     /(60.d0*dxlevel(idims))
        case default
          call mpistop("unknown stencil in divvector")
        end select
      enddo
    else
      do idims=1,ndim
        hxo^l=ixo^l-kr(idims,^d);
        ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
        jxc^l=ixc^l+kr(idims,^d);
        if(stretched_dim(idims) .and. stretch_uncentered) then
          ! linear interpolation at cell interface along stretched dimension
          qc(ixc^s)=block%surfaceC(ixc^s,idims)*(qvec(ixc^s,idims)+0.5d0*block%dx(ixc^s,idims)*&
               (qvec(jxc^s,idims)-qvec(ixc^s,idims))/(block%x(jxc^s,idims)-block%x(ixc^s,idims)))
        else
          qc(ixc^s)=block%surfaceC(ixc^s,idims)*half*(qvec(ixc^s,idims)+qvec(jxc^s,idims))
        end if
        divq(ixo^s)=divq(ixo^s)+qc(ixo^s)-qc(hxo^s)
      end do
      divq(ixo^s)=divq(ixo^s)/block%dvolume(ixo^s)
    end if
  end subroutine divvector
 
  !> Calculate divergence of a vector qvec within ixL
  !> using limited extrapolation to cell edges
  subroutine divvectors(qvec,ixI^L,ixO^L,divq)
    use mod_global_parameters
    use mod_limiter
    use mod_ppm
    integer, intent(in)                :: ixI^L,ixO^L
    double precision, intent(in)       :: qvec(ixI^S,1:ndir)
    double precision, intent(inout)    :: divq(ixI^S)
    double precision, dimension(ixI^S) :: qL,qR,dqC,ldq,rdq
 
    double precision :: invdx(1:ndim)
    integer          :: hxO^L,ixC^L,jxC^L,idims,ix^L,gxC^L,hxC^L
 
    ix^l=ixo^l^ladd2;
 
    if (iximin^d>ixmin^d.or.iximax^d<ixmax^d|.or.) &
         call mpistop("Error in divvectorS: Non-conforming input limits")
 
    invdx=1.d0/dxlevel
    divq(ixo^s)=zero
    do idims=1,ndim
      hxo^l=ixo^l-kr(idims,^d);
      ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
      jxc^l=ixc^l+kr(idims,^d);
      gxcmin^d=ixcmin^d-kr(idims,^d);gxcmax^d=jxcmax^d;
      hxc^l=gxc^l+kr(idims,^d);
      qr(gxc^s) = qvec(hxc^s,idims)
      ql(gxc^s) = qvec(gxc^s,idims)
      if(type_gradient_limiter(block%level)/=limiter_ppm) then
        dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
        call dwlimiter2(dqc,ixi^l,gxc^l,idims,type_gradient_limiter(block%level),ldw=ldq,rdw=rdq)
        ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
        qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
      else
        dqc(ixi^s)=qvec(ixi^s,idims)
        call ppmlimitervar(ixi^l,ixo^l,idims,dqc,dqc,ql,qr)
      endif
 
      if (slab_uniform) then
        divq(ixo^s)=divq(ixo^s)+half*(qr(ixo^s)-ql(hxo^s))*invdx(idims)
      else
        qr(ixc^s)=block%surfaceC(ixc^s,idims)*qr(ixc^s)
        ql(ixc^s)=block%surfaceC(ixc^s,idims)*ql(ixc^s)
        divq(ixo^s)=divq(ixo^s)+qr(ixo^s)-ql(hxo^s)
      end if
    end do
    if(.not.slab_uniform) divq(ixo^s)=divq(ixo^s)/block%dvolume(ixo^s)
  end subroutine divvectors
 
  !**************************************************************************
  ! Purpose: Computes the Laplacian of a scalar field q(ixI^S) at cell 
  !          centers and outputs it in laplq(ixO^S), 
  !          which is also defined at cell **centers**. 
  !
  !          For uniform Cartesian coordinates, an optional input 
  !          parameter nth_in allows increasing the order of the 
  !          central difference scheme to 2*nth_in.
  !**************************************************************************
  subroutine laplacian(q,ixI^L,ixO^L,laplq,nth_in)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    integer, intent(in), optional   :: nth_in
    double precision, intent(in)    :: q(ixI^S)
    double precision, intent(inout) :: laplq(ixI^S)
    integer                         :: lxO^L, jxO^L, hxO^L, kxO^L, nth, idim
 
    if(present(nth_in)) then
      nth = nth_in
    else
      nth = 1
    endif
    if(nth .gt. nghostcells) then 
      call mpistop("laplacian stencil too wide")
    endif
 
    select case(coordinate)
    case(cartesian)
      select case(nth)
      case(1)
        laplq(ixo^s)=zero
        do idim=1,ndim
           jxo^l=ixo^l+kr(idim,^d);
           hxo^l=ixo^l-kr(idim,^d);
           laplq(ixo^s)=laplq(ixo^s)+&
                 (q(jxo^s)-2.0d0*q(ixo^s)+q(hxo^s))/dxlevel(idim)**2
         end do
      case(2)
        laplq(ixo^s)=zero
        do idim=1,ndim
           lxo^l=ixo^l+2*kr(idim,^d);
           jxo^l=ixo^l+kr(idim,^d);
           hxo^l=ixo^l-kr(idim,^d);
           kxo^l=ixo^l-2*kr(idim,^d);
           laplq(ixo^s)=laplq(ixo^s)+&
                 (-q(lxo^s)+16.0d0*q(jxo^s)-30.0d0*q(ixo^s)+16.0d0*q(hxo^s)-q(kxo^s)) &
                 /(12.0d0 * dxlevel(idim)**2)
        end do
      case default
        call mpistop("unknown stencil laplacian")
      end select
    case(cartesian_stretched,cartesian_expansion)
        call mpistop("No laplacian available")
    case(spherical)
        call mpistop("No laplacian available")
    case(cylindrical)
        call mpistop("No laplacian available")
    case default
      call mpistop('Unknown geometry')
    end select
  end subroutine laplacian
 
  !**************************************************************************
  ! Purpose: Computes the Laplacian of a vector field qvec(ixI^S,1:ndir) at cell 
  !          centers and outputs it in lapl_qvec(ixO^S,1:ndir), 
  !          which is also defined at cell **centers**. 
  !**************************************************************************
  subroutine laplacian_of_vector(qvec,ixI^L,ixO^L,lapl_qvec)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qvec(ixI^S,1:ndir)
    double precision, intent(inout) :: lapl_qvec(ixI^S,1:ndir)
    integer :: idir
    double precision :: tmp(ixI^S)
 
    do idir=1,ndir
       call laplacian(qvec(ixi^s,idir),ixi^l,ixo^l,tmp)
       lapl_qvec(ixo^s,idir)=tmp(ixo^s)
    enddo
 
  end subroutine laplacian_of_vector
 
 
  !> Calculate curl of a vector qvec within ixL
  !> Options to
  !>        employ standard second order CD evaluations
  !>        use Gauss theorem for non-Cartesian grids
  !>        use Stokes theorem for non-Cartesian grids
  subroutine curlvector(qvec,ixI^L,ixO^L,curlvec,idirmin,idirmin0,ndir0,fourthorder)
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L,ixO^L
    integer, intent(in)             :: ndir0, idirmin0
    integer, intent(inout)          :: idirmin
    double precision, intent(in)    :: qvec(ixI^S,1:ndir0)
    double precision, intent(inout) :: curlvec(ixI^S,idirmin0:3)
    logical, intent(in), optional   :: fourthorder !< Default: false
 
    double precision :: invdx(1:ndim)
    double precision :: tmp(ixI^S),tmp2(ixI^S),xC(ixI^S),surface(ixI^S)
    integer          :: ixA^L,ixC^L,jxC^L,idir,jdir,kdir,hxO^L,jxO^L,kxO^L,gxO^L
    logical          :: use_4th_order
 
    ! Calculate curl within ixL: CurlV_i=eps_ijk*d_j V_k
    ! Curl can have components (idirmin:3)
    ! Determine exact value of idirmin while doing the loop.
 
    use_4th_order = .false.
    if (present(fourthorder)) use_4th_order = fourthorder
 
    if (use_4th_order) then
      if (.not. slab_uniform) &
           call mpistop("curlvector: 4th order only supported for slab geometry")
      ! Fourth order, stencil width is two
      ixa^l=ixo^l^ladd2;
    else
      ! Second order, stencil width is one
      ixa^l=ixo^l^ladd1;
    end if
 
    if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
         call mpistop("Error in curlvector: Non-conforming input limits")
 
    idirmin=4
    curlvec(ixo^s,idirmin0:3)=zero
 
    ! all non-Cartesian cases
    select case(coordinate)
      case(cartesian) ! Cartesian grids
        invdx=1.d0/dxlevel
        do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
          if(lvc(idir,jdir,kdir)/=0)then
            if (.not. use_4th_order) then
              ! Use second order scheme
              jxo^l=ixo^l+kr(jdir,^d);
              hxo^l=ixo^l-kr(jdir,^d);
              tmp(ixo^s)=half*(qvec(jxo^s,kdir) &
                   - qvec(hxo^s,kdir))*invdx(jdir)
            else
              ! Use fourth order scheme
              kxo^l=ixo^l+2*kr(jdir,^d);
              jxo^l=ixo^l+kr(jdir,^d);
              hxo^l=ixo^l-kr(jdir,^d);
              gxo^l=ixo^l-2*kr(jdir,^d);
              tmp(ixo^s)=(-qvec(kxo^s,kdir) + 8.0d0 * qvec(jxo^s,kdir) - 8.0d0 * &
                   qvec(hxo^s,kdir) + qvec(gxo^s,kdir))/(12.0d0 * dxlevel(jdir))
            end if
            if(lvc(idir,jdir,kdir)==1)then
              curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp(ixo^s)
            else
              curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp(ixo^s)
            endif
            if(idir<idirmin)idirmin=idir
          endif
        enddo; enddo; enddo;
      case(cartesian_stretched) ! stretched Cartesian grids
        do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
          if(lvc(idir,jdir,kdir)/=0)then
            select case(type_curl)
              case(central)
                tmp(ixa^s)=qvec(ixa^s,kdir)
                hxo^l=ixo^l-kr(jdir,^d);
                jxo^l=ixo^l+kr(jdir,^d);
                ! second order centered differencing
                tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,jdir)-block%x(hxo^s,jdir))
              case(gaussbased)
                hxo^l=ixo^l-kr(jdir,^d);
                ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                jxc^l=ixc^l+kr(jdir,^d);
                tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
              case(stokesbased)
                hxo^l=ixo^l-kr(jdir,^d);
                ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                jxc^l=ixc^l+kr(jdir,^d);
                if(kdir<=ndim)then
                  tmp(ixc^s)=block%ds(ixo^s,kdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                else
                  tmp(ixc^s)=(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                endif
                if(idir<=ndim)then
                  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%surface(ixo^s,idir)
                else ! essentially for 2.5D case, idir=3 and jdir,kdir<=2
                  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/(block%ds(ixo^s,jdir)*block%ds(ixo^s,kdir))
                endif
              case default
                call mpistop('no such curl evaluator')
            end select
            if(lvc(idir,jdir,kdir)==1)then
              curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
            else
              curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
            endif
            if(idir<idirmin)idirmin=idir
          endif
        enddo; enddo; enddo;
      case(spherical) ! possibly stretched spherical grids
        select case(type_curl)
          case(central) ! ok for any dimensionality
            do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
                tmp(ixa^s)=qvec(ixa^s,kdir)
                hxo^l=ixo^l-kr(jdir,^d);
                jxo^l=ixo^l+kr(jdir,^d);
                select case(jdir)
                case(1)
                tmp(ixa^s)=tmp(ixa^s)*block%x(ixa^s,1)
                tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,1)-block%x(hxo^s,1))*block%x(ixo^s,1))
                {^nooned    case(2)
                if(idir==1) tmp(ixa^s)=tmp(ixa^s)*dsin(block%x(ixa^s,2))
                tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
                if(idir==1) tmp2(ixo^s)=tmp2(ixo^s)/dsin(block%x(ixo^s,2))
                }
                {^ifthreed  case(3)
                tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1)*dsin(block%x(ixo^s,2)))
                }
                end select
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
                else
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
                endif
                if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo; enddo;
          case(gaussbased)
            do idir=idirmin0,3;
            do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
                hxo^l=ixo^l-kr(jdir,^d);
                ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                jxc^l=ixc^l+kr(jdir,^d);
                tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
                else
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
                 endif
                 if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo;
            ! geometric terms
            if(idir==2.and.phi_>0) curlvec(ixo^s,2)=curlvec(ixo^s,2)+qvec(ixo^s,phi_)/block%x(ixo^s,r_)
            {^nooned
            if(idir==phi_) curlvec(ixo^s,phi_)=curlvec(ixo^s,phi_)-qvec(ixo^s,2)/block%x(ixo^s,r_) &
                 +qvec(ixo^s,r_)*dcos(block%x(ixo^s,2))/(block%x(ixo^s,r_)*dsin(block%x(ixo^s,2)))
            }
            enddo;
          case(stokesbased)
            !if(ndim<3) call mpistop("Stokesbased for 3D spherical only")
            do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
             if(lvc(idir,jdir,kdir)/=0)then
              select case(idir)
              case(1)
                if(jdir<kdir) then
                  ! idir=1,jdir=2,kdir=3
                  !! integral along 3rd dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(3) at cell interface along 2nd dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  ! 2nd coordinate at cell interface along 2nd dimension
                  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
                  curlvec(ixo^s,idir)=(dsin(xc(ixo^s))*tmp(ixo^s)-&
                       dsin(xc(hxo^s))*tmp(hxo^s))*block%dx(ixo^s,kdir)
                  !! integral along 2nd dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(2) at cell interface along 3rd dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir))&
                       /block%surface(ixo^s,idir)*block%x(ixo^s,idir)
                end if
              case(2)
                if(jdir<kdir) then
                  ! idir=2,jdir=1,kdir=3
                  !! integral along 1st dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(1) at cell interface along 3rd dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,1)
                  !! integral along 3rd dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(3) at cell interface along 1st dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  ! 1st coordinate at cell interface along 1st dimension
                  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(hxo^s)*tmp(hxo^s)-xc(ixo^s)*tmp(ixo^s))*&
                       dsin(block%x(ixo^s,idir))*block%dx(ixo^s,kdir))/block%surface(ixo^s,idir)
                end if
              case(3)
                if(jdir<kdir) then
                  ! idir=3,jdir=1,kdir=2
                  !! integral along 1st dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(1) at cell interface along 2nd dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
                  !! integral along 2nd dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(2) at cell interface along 1st dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  ! 1st coordinate at cell interface along 1st dimension
                  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
                  if(ndim==3) then
                    surface(ixo^s)=block%surface(ixo^s,idir)
                  else
                    surface(ixo^s)=block%x(ixo^s,jdir)*block%dx(ixo^s,kdir)*block%dx(ixo^s,jdir)
                  end if
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*block%dx(ixo^s,kdir))&
                       /surface(ixo^s)
                end if
              end select
              if(idir<idirmin)idirmin=idir
             endif
            enddo; enddo; enddo;
          case default
            call mpistop('no such curl evaluator')
        end select
      case(cylindrical) ! possibly stretched cylindrical grids
        select case(type_curl)
          case(central)  ! works for any dimensionality, polar/cylindrical
            do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
                tmp(ixa^s)=qvec(ixa^s,kdir)
                hxo^l=ixo^l-kr(jdir,^d);
                jxo^l=ixo^l+kr(jdir,^d);
                if(z_==3.or.z_==-1) then
                  ! Case Polar_2D, Polar_2.5D or Polar_3D, i.e. R,phi,Z
                  select case(jdir)
                  case(1)
                  if(idir==z_) tmp(ixa^s)=tmp(ixa^s)*block%x(ixa^s,1) ! R V_phi
                  ! computes d(R V_phi)/dR or d V_Z/dR
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
                  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
                  {^nooned      case(2)
                  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
                  }
                  {^ifthreed    case(3)
                  ! handles d V_phi/dZ or d V_R/dZ
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,3)-block%x(hxo^s,3))
                  }
                  end select
                end if
                if(phi_==3.or.phi_==-1) then
                  ! Case Cylindrical_2D, Cylindrical_2.5D or Cylindrical_3D, i.e. R,Z,phi
                  select case(jdir)
                  case(1)
                  if(idir==z_) tmp(ixa^s)=tmp(ixa^s)*block%x(ixa^s,1) ! R V_phi
                  ! computes d(R V_phi)/dR or d V_Z/dR
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
                  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
                  {^nooned      case(2)
                  ! handles d V_phi/dZ or d V_R/dZ
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,2)-block%x(hxo^s,2))
                  }
                  {^ifthreed    case(3)
                  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1))
                  }
                  end select
                end if
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
                else
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
                 endif
                 if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo; enddo;
          case(gaussbased) ! works for any dimensionality, polar/cylindrical
            if(ndim<2) call mpistop("Gaussbased for 2D, 2.5D or 3D polar or cylindrical only")
            do idir=idirmin0,3;
            do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
                hxo^l=ixo^l-kr(jdir,^d);
                ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                jxc^l=ixc^l+kr(jdir,^d);
                tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
                else
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
                 endif
                 if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo;
            ! geometric term from d e_R/d phi= e_phi for unit vectors e_R, e_phi
            !       but minus sign appears due to R,Z,phi ordering (?)
            ! note that in cylindrical 2D (RZ), phi_ is -1
            ! note that in polar 2D     (Rphi), z_ is -1
            if((idir==phi_.or.(phi_==-1.and.idir==3)).and.z_>0) then
              ! cylindrical
              if(  z_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-qvec(ixo^s,z_)/block%x(ixo^s,r_)
              ! polar
              if(phi_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+qvec(ixo^s,z_)/block%x(ixo^s,r_)
            endif
            enddo;
          case(stokesbased)
            if(ndim<3) call mpistop("Stokesbased for 3D cylindrical only")
            do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
               if(idir==r_) then
                if(jdir==phi_) then
                  ! idir=r,jdir=phi,kdir=z
                  !! integral along z dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(z) at cell interface along phi dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,kdir)
                  !! integral along phi dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(phi) at cell interface along z dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%x(ixo^s,idir)*block%dx(ixo^s,jdir))&
                       /block%surface(ixo^s,idir)
                end if
               else if(idir==phi_) then
                if(jdir<kdir) then
                  ! idir=phi,jdir=r,kdir=z
                  !! integral along r dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(r) at cell interface along z dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,jdir)
                  !! integral along z dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(z) at cell interface along r dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,kdir))&
                       /block%surface(ixo^s,idir)
                end if
               else ! idir==z_
                if(jdir<kdir) then
                  ! idir=z,jdir=r,kdir=phi
                  !! integral along r dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(r) at cell interface along phi dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
                  !! integral along phi dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(phi) at cell interface along r dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  ! r coordinate at cell interface along r dimension
                  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*block%dx(ixo^s,kdir))&
                       /block%surface(ixo^s,idir)
                end if
               end if
               if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo; enddo;
          case default
            call mpistop('no such curl evaluator')
        end select
        case default
          call mpistop('not possible to calculate curl')
    end select
 
  end subroutine curlvector
 
  !> Calculate idim transverse components of curl of a vector qvec within ixL
  !> Options to
  !>        employ standard second order CD evaluations
  !>        use Gauss theorem for non-Cartesian grids
  !>        use Stokes theorem for non-Cartesian grids
  subroutine curlvector_trans(qvec,qvecc,ixI^L,ixO^L,curlvec,idim,idirmin,idirmin0,ndir0)
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L,ixO^L
    integer, intent(in)             :: idim, ndir0, idirmin0
    integer, intent(inout)          :: idirmin
    double precision, intent(in)    :: qvec(ixI^S,1:ndir0),qvecc(ixI^S,1:ndir0)
    double precision, intent(inout) :: curlvec(ixI^S,idirmin0:3)
 
    double precision :: invdx(1:ndim)
    double precision :: tmp(ixI^S),tmp2(ixI^S),xC(ixI^S),surface(ixI^S)
    integer          :: ixA^L,ixC^L,jxC^L,idir,jdir,kdir,hxO^L,jxO^L
 
    ! Calculate curl within ixL: CurlV_i=eps_ijk*d_j V_k
    ! Curl can have components (idirmin:3)
    ! Determine exact value of idirmin while doing the loop.
 
    idirmin=4
    curlvec(ixo^s,idirmin0:3)=zero
    ! Second order, stencil width is one
    ixa^l=ixo^l^ladd1;
 
    ! all non-Cartesian cases
    select case(coordinate)
      case(cartesian) ! Cartesian grids
        invdx=1.d0/dxlevel
        do idir=idirmin0,3
          do jdir=1,ndim; do kdir=1,ndir0
            if(lvc(idir,jdir,kdir)/=0)then
              if(jdir/=idim) then
                tmp(ixi^s)=qvec(ixi^s,kdir)
                hxo^l=ixo^l-kr(jdir,^d);
                jxo^l=ixo^l+kr(jdir,^d);
              else
                ! use two cell-center values for gradient at interface of the two cells
                ! because left and right neighbor interface values is unavailable at block boundary faces
                tmp(ixi^s)=qvecc(ixi^s,kdir)
                hxo^l=ixo^l;
                jxo^l=ixo^l+kr(jdir,^d);
              end if
              ! second order centered differencing
              tmp(ixo^s)=half*(tmp(jxo^s)-tmp(hxo^s))*invdx(jdir)
              if(lvc(idir,jdir,kdir)==1)then
                curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp(ixo^s)
              else
                curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp(ixo^s)
              endif
              if(idir<idirmin)idirmin=idir
            endif
          enddo; enddo;
        end do
      case(cartesian_stretched) ! stretched Cartesian grids
        do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
          if(lvc(idir,jdir,kdir)/=0)then
            select case(type_curl)
              case(central)
                tmp(ixi^s)=qvec(ixi^s,kdir)
                hxo^l=ixo^l-kr(jdir,^d);
                jxo^l=ixo^l+kr(jdir,^d);
                ! second order centered differencing
                tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,jdir)-block%x(hxo^s,jdir))
              case(gaussbased)
                hxo^l=ixo^l-kr(jdir,^d);
                ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                jxc^l=ixc^l+kr(jdir,^d);
                tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
              case(stokesbased)
                hxo^l=ixo^l-kr(jdir,^d);
                ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                jxc^l=ixc^l+kr(jdir,^d);
                if(kdir<=ndim)then
                  tmp(ixc^s)=block%ds(ixo^s,kdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                else
                  tmp(ixc^s)=(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                endif
                if(idir<=ndim)then
                  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%surface(ixo^s,idir)
                else ! essentially for 2.5D case, idir=3 and jdir,kdir<=2
                  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/(block%ds(ixo^s,jdir)*block%ds(ixo^s,kdir))
                endif
              case default
                call mpistop('no such curl evaluator')
            end select
            if(lvc(idir,jdir,kdir)==1)then
              curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
            else
              curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
            endif
            if(idir<idirmin)idirmin=idir
          endif
        enddo; enddo; enddo;
      case(spherical) ! possibly stretched spherical grids
        select case(type_curl)
          case(central) ! ok for any dimensionality
            do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
                tmp(ixi^s)=qvec(ixi^s,kdir)
                hxo^l=ixo^l-kr(jdir,^d);
                jxo^l=ixo^l+kr(jdir,^d);
                select case(jdir)
                case(1)
                tmp(ixa^s)=tmp(ixa^s)*block%x(ixa^s,1)
                tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,1)-block%x(hxo^s,1))*block%x(ixo^s,1))
                {^nooned    case(2)
                if(idir==1) tmp(ixa^s)=tmp(ixa^s)*dsin(block%x(ixa^s,2))
                tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
                if(idir==1) tmp2(ixo^s)=tmp2(ixo^s)/dsin(block%x(ixo^s,2))
                }
                {^ifthreed  case(3)
                tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1)*dsin(block%x(ixo^s,2)))
                }
                end select
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
                else
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
                endif
                if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo; enddo;
          case(gaussbased)
            do idir=idirmin0,3;
            do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
                hxo^l=ixo^l-kr(jdir,^d);
                ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                jxc^l=ixc^l+kr(jdir,^d);
                tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
                else
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
                 endif
                 if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo;
            ! geometric terms
            if(idir==2.and.phi_>0) curlvec(ixo^s,2)=curlvec(ixo^s,2)+qvec(ixo^s,phi_)/block%x(ixo^s,r_)
            {^nooned
            if(idir==phi_) curlvec(ixo^s,phi_)=curlvec(ixo^s,phi_)-qvec(ixo^s,2)/block%x(ixo^s,r_) &
                 +qvec(ixo^s,r_)*dcos(block%x(ixo^s,2))/(block%x(ixo^s,r_)*dsin(block%x(ixo^s,2)))
            }
            enddo;
          case(stokesbased)
            !if(ndim<3) call mpistop("Stokesbased for 3D spherical only")
            do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
             if(lvc(idir,jdir,kdir)/=0)then
              select case(idir)
              case(1)
                if(jdir<kdir) then
                  ! idir=1,jdir=2,kdir=3
                  !! integral along 3rd dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(3) at cell interface along 2nd dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  ! 2nd coordinate at cell interface along 2nd dimension
                  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
                  curlvec(ixo^s,idir)=(dsin(xc(ixo^s))*tmp(ixo^s)-&
                       dsin(xc(hxo^s))*tmp(hxo^s))*block%dx(ixo^s,kdir)
                  !! integral along 2nd dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(2) at cell interface along 3rd dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir))&
                       /block%surface(ixo^s,idir)*block%x(ixo^s,idir)
                end if
              case(2)
                if(jdir<kdir) then
                  ! idir=2,jdir=1,kdir=3
                  !! integral along 1st dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(1) at cell interface along 3rd dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,1)
                  !! integral along 3rd dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(3) at cell interface along 1st dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  ! 1st coordinate at cell interface along 1st dimension
                  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(hxo^s)*tmp(hxo^s)-xc(ixo^s)*tmp(ixo^s))*&
                       dsin(block%x(ixo^s,idir))*block%dx(ixo^s,kdir))/block%surface(ixo^s,idir)
                end if
              case(3)
                if(jdir<kdir) then
                  ! idir=3,jdir=1,kdir=2
                  !! integral along 1st dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(1) at cell interface along 2nd dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
                  !! integral along 2nd dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(2) at cell interface along 1st dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  ! 1st coordinate at cell interface along 1st dimension
                  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
                  if(ndim==3) then
                    surface(ixo^s)=block%surface(ixo^s,idir)
                  else
                    surface(ixo^s)=block%x(ixo^s,jdir)*block%dx(ixo^s,kdir)*block%dx(ixo^s,jdir)
                  end if
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*block%dx(ixo^s,kdir))&
                       /surface(ixo^s)
                end if
              end select
              if(idir<idirmin)idirmin=idir
             endif
            enddo; enddo; enddo;
          case default
            call mpistop('no such curl evaluator')
        end select
      case(cylindrical) ! possibly stretched cylindrical grids
        select case(type_curl)
          case(central)  ! works for any dimensionality, polar/cylindrical
            do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
                tmp(ixi^s)=qvec(ixi^s,kdir)
                hxo^l=ixo^l-kr(jdir,^d);
                jxo^l=ixo^l+kr(jdir,^d);
                if(z_==3.or.z_==-1) then
                  ! Case Polar_2D, Polar_2.5D or Polar_3D, i.e. R,phi,Z
                  select case(jdir)
                  case(1)
                  if(idir==z_) tmp(ixa^s)=tmp(ixa^s)*block%x(ixa^s,1) ! R V_phi
                  ! computes d(R V_phi)/dR or d V_Z/dR
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
                  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
                  {^nooned      case(2)
                  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
                  }
                  {^ifthreed    case(3)
                  ! handles d V_phi/dZ or d V_R/dZ
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,3)-block%x(hxo^s,3))
                  }
                  end select
                end if
                if(phi_==3.or.phi_==-1) then
                  ! Case Cylindrical_2D, Cylindrical_2.5D or Cylindrical_3D, i.e. R,Z,phi
                  select case(jdir)
                  case(1)
                  if(idir==z_) tmp(ixa^s)=tmp(ixa^s)*block%x(ixa^s,1) ! R V_phi
                  ! computes d(R V_phi)/dR or d V_Z/dR
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
                  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
                  {^nooned      case(2)
                  ! handles d V_phi/dZ or d V_R/dZ
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,2)-block%x(hxo^s,2))
                  }
                  {^ifthreed    case(3)
                  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
                  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1))
                  }
                  end select
                end if
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
                else
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
                 endif
                 if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo; enddo;
          case(gaussbased) ! works for any dimensionality, polar/cylindrical
            if(ndim<2) call mpistop("Gaussbased for 2D, 2.5D or 3D polar or cylindrical only")
            do idir=idirmin0,3;
            do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
                hxo^l=ixo^l-kr(jdir,^d);
                ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                jxc^l=ixc^l+kr(jdir,^d);
                tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                       (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
                tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
                else
                  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
                 endif
                 if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo;
            ! geometric term from d e_R/d phi= e_phi for unit vectors e_R, e_phi
            !       but minus sign appears due to R,Z,phi ordering (?)
            ! note that in cylindrical 2D (RZ), phi_ is -1
            ! note that in polar 2D     (Rphi), z_ is -1
            if((idir==phi_.or.(phi_==-1.and.idir==3)).and.z_>0) then
              ! cylindrical
              if(  z_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-qvec(ixo^s,z_)/block%x(ixo^s,r_)
              ! polar
              if(phi_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+qvec(ixo^s,z_)/block%x(ixo^s,r_)
            endif
            enddo;
          case(stokesbased)
            if(ndim<3) call mpistop("Stokesbased for 3D cylindrical only")
            do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
              if(lvc(idir,jdir,kdir)/=0)then
               if(idir==r_) then
                if(jdir==phi_) then
                  ! idir=r,jdir=phi,kdir=z
                  !! integral along z dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(z) at cell interface along phi dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,kdir)
                  !! integral along phi dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(phi) at cell interface along z dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%x(ixo^s,idir)*block%dx(ixo^s,jdir))&
                       /block%surface(ixo^s,idir)
                end if
               else if(idir==phi_) then
                if(jdir<kdir) then
                  ! idir=phi,jdir=r,kdir=z
                  !! integral along r dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(r) at cell interface along z dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,jdir)
                  !! integral along z dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(z) at cell interface along r dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,kdir))&
                       /block%surface(ixo^s,idir)
                end if
               else ! idir==z_
                if(jdir<kdir) then
                  ! idir=z,jdir=r,kdir=phi
                  !! integral along r dimension
                  hxo^l=ixo^l-kr(kdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(kdir,^d);
                  ! qvec(r) at cell interface along phi dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
                         (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
                  end if
                  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
                  !! integral along phi dimension
                  hxo^l=ixo^l-kr(jdir,^d);
                  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
                  jxc^l=ixc^l+kr(jdir,^d);
                  ! qvec(phi) at cell interface along r dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
                         (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
                  else
                    tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
                  end if
                  ! r coordinate at cell interface along r dimension
                  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
                  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*block%dx(ixo^s,kdir))&
                       /block%surface(ixo^s,idir)
                end if
               end if
               if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo; enddo;
          case default
            call mpistop('no such curl evaluator')
        end select
        case default
          call mpistop('not possible to calculate curl')
    end select
 
  end subroutine curlvector_trans
 
end module mod_geometry