curlvector Subroutine

  subroutine curlvector(qvec,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
     ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,curlvec,idirmin,idirmin0,&
     ndir0,fourthorder)
    use mod_global_parameters

    integer, intent(in)             :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
       ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
    integer, intent(in)             :: ndir0, idirmin0
    integer, intent(inout)          :: idirmin
    double precision, intent(in)    :: qvec(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3,1:ndir0)
    double precision, intent(inout) :: curlvec(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3,idirmin0:3)
    logical, intent(in), optional   :: fourthorder
    !< Default: false

    integer          :: ixAmin1,ixAmin2,ixAmin3,ixAmax1,ixAmax2,ixAmax3,&
       ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,jxCmin1,jxCmin2,jxCmin3,&
       jxCmax1,jxCmax2,jxCmax3,idir,jdir,kdir,hxOmin1,hxOmin2,hxOmin3,hxOmax1,&
       hxOmax2,hxOmax3,jxOmin1,jxOmin2,jxOmin3,jxOmax1,jxOmax2,jxOmax3,kxOmin1,&
       kxOmin2,kxOmin3,kxOmax1,kxOmax2,kxOmax3,gxOmin1,gxOmin2,gxOmin3,gxOmax1,&
       gxOmax2,gxOmax3
    double precision :: invdx(1:ndim)
    double precision :: tmp(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
       tmp2(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
       xC(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
       surface(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3)
    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.

    ! The Stokesbased discretisation divides by block%surface, the cell-centre
    ! face areas.  Those are no longer built: no other routine reads them, and
    ! carrying a block-sized array per direction for 1:max_blocks purely for
    ! this branch is not worth the memory.  Fail here rather than divide by an
    ! unassociated pointer.
    if (type_curl == Stokesbased) call mpistop(&
       "curlvector: typecurl='Stokesbased' needs block%surface, which is no &
       &longer built - see geo_t in mod_physicaldata")

    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
      ixAmin1=ixOmin1-2;ixAmin2=ixOmin2-2;ixAmin3=ixOmin3-2;ixAmax1=ixOmax1+2
      ixAmax2=ixOmax2+2;ixAmax3=ixOmax3+2;
    else
      ! Second order, stencil width is one
      ixAmin1=ixOmin1-1;ixAmin2=ixOmin2-1;ixAmin3=ixOmin3-1;ixAmax1=ixOmax1+1
      ixAmax2=ixOmax2+1;ixAmax3=ixOmax3+1;
    end if

    if (ixImin1>ixAmin1.or.ixImax1<ixAmax1.or.ixImin2>ixAmin2.or.&
       ixImax2<ixAmax2.or.ixImin3>ixAmin3.or.ixImax3<ixAmax3) call &
       mpistop("Error in curlvector: Non-conforming input limits")

    idirmin=4
    curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,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
              jxOmin1=ixOmin1+kr(jdir,1);jxOmin2=ixOmin2+kr(jdir,2)
              jxOmin3=ixOmin3+kr(jdir,3);jxOmax1=ixOmax1+kr(jdir,1)
              jxOmax2=ixOmax2+kr(jdir,2);jxOmax3=ixOmax3+kr(jdir,3);
              hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
              hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
              hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
              tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3)=half*(qvec(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                 jxOmin3:jxOmax3,kdir) - qvec(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                 hxOmin3:hxOmax3,kdir))*invdx(jdir)
            else
              ! Use fourth order scheme
              kxOmin1=ixOmin1+2*kr(jdir,1);kxOmin2=ixOmin2+2*kr(jdir,2)
              kxOmin3=ixOmin3+2*kr(jdir,3);kxOmax1=ixOmax1+2*kr(jdir,1)
              kxOmax2=ixOmax2+2*kr(jdir,2);kxOmax3=ixOmax3+2*kr(jdir,3);
              jxOmin1=ixOmin1+kr(jdir,1);jxOmin2=ixOmin2+kr(jdir,2)
              jxOmin3=ixOmin3+kr(jdir,3);jxOmax1=ixOmax1+kr(jdir,1)
              jxOmax2=ixOmax2+kr(jdir,2);jxOmax3=ixOmax3+kr(jdir,3);
              hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
              hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
              hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
              gxOmin1=ixOmin1-2*kr(jdir,1);gxOmin2=ixOmin2-2*kr(jdir,2)
              gxOmin3=ixOmin3-2*kr(jdir,3);gxOmax1=ixOmax1-2*kr(jdir,1)
              gxOmax2=ixOmax2-2*kr(jdir,2);gxOmax3=ixOmax3-2*kr(jdir,3);
              tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3)=(-qvec(kxOmin1:kxOmax1,kxOmin2:kxOmax2,&
                 kxOmin3:kxOmax3,kdir) + 8.0d0 * qvec(jxOmin1:jxOmax1,&
                 jxOmin2:jxOmax2,jxOmin3:jxOmax3,&
                 kdir) - 8.0d0 * qvec(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                 hxOmin3:hxOmax3,kdir) + qvec(gxOmin1:gxOmax1,gxOmin2:gxOmax2,&
                 gxOmin3:gxOmax3,kdir))/(12.0d0 * dxlevel(jdir))
            end if
            if(lvc(idir,jdir,kdir)==1)then
              curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                 idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                 idir)+tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)
            else
              curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                 idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                 idir)-tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)
            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(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3)=qvec(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3,kdir)
                hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                jxOmin1=ixOmin1+kr(jdir,1);jxOmin2=ixOmin2+kr(jdir,2)
                jxOmin3=ixOmin3+kr(jdir,3);jxOmax1=ixOmax1+kr(jdir,1)
                jxOmax2=ixOmax2+kr(jdir,2);jxOmax3=ixOmax3+kr(jdir,3);
                ! second order centered differencing
                tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                   jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3))/(block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                   jxOmin3:jxOmax3,jdir)-block%x(hxOmin1:hxOmax1,&
                   hxOmin2:hxOmax2,hxOmin3:hxOmax3,jdir))
              case(Gaussbased)
                hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3)=block%surfaceC(ixCmin1:ixCmax1,&
                   ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir)*(qvec(ixCmin1:ixCmax1,&
                   ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                   kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3,jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                   jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3,kdir))/(block%x(jxCmin1:jxCmax1,&
                   jxCmin2:jxCmax2,jxCmin3:jxCmax3,&
                   jdir)-block%x(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3,jdir)))
                tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)=(tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3))/block%dvolume(ixOmin1:ixOmax1,&
                   ixOmin2:ixOmax2,ixOmin3:ixOmax3)
              case(Stokesbased)
                hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                if(kdir<=ndim)then
                  tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3)=block%ds(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,kdir)*(qvec(ixCmin1:ixCmax1,&
                     ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                     kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3,jdir)*(qvec(jxCmin1:jxCmax1,&
                     jxCmin2:jxCmax2,jxCmin3:jxCmax3,&
                     kdir)-qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3,kdir))/(block%x(jxCmin1:jxCmax1,&
                     jxCmin2:jxCmax2,jxCmin3:jxCmax3,&
                     jdir)-block%x(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3,jdir)))
                else
                  tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3)=(qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3,kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                     ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                     jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                     jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,&
                     ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                     kdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                     jxCmin3:jxCmax3,jdir)-block%x(ixCmin1:ixCmax1,&
                     ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir)))
                endif
                if(idir<=ndim)then
                  tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=(tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))/block%surface(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,idir)
                else ! essentially for 2.5D case, idir=3 and jdir,kdir<=2
                  tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=(tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))/(block%ds(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     jdir)*block%ds(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,kdir))
                endif
              case default
                call mpistop('no such curl evaluator')
            end select
            if(lvc(idir,jdir,kdir)==1)then
              curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                 idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                 idir)+tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)
            else
              curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                 idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                 idir)-tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)
            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(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3)=qvec(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3,kdir)
                hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                jxOmin1=ixOmin1+kr(jdir,1);jxOmin2=ixOmin2+kr(jdir,2)
                jxOmin3=ixOmin3+kr(jdir,3);jxOmax1=ixOmax1+kr(jdir,1)
                jxOmax2=ixOmax2+kr(jdir,2);jxOmax3=ixOmax3+kr(jdir,3);
                select case(jdir)
                case(1)
                tmp(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3)=tmp(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3)*block%x(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3,1)
                tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                   jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3))/((block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                   jxOmin3:jxOmax3,1)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3,1))*block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3,1))
                    case(2)
                if(idir==1) tmp(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3)=tmp(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3)*dsin(block%x(ixAmin1:ixAmax1,&
                   ixAmin2:ixAmax2,ixAmin3:ixAmax3,2))
                tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                   jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3))/((block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                   jxOmin3:jxOmax3,2)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3,2))*block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3,1))
                if(idir==1) tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)=tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)/dsin(block%x(ixOmin1:ixOmax1,&
                   ixOmin2:ixOmax2,ixOmin3:ixOmax3,2))
               
                  case(3)
                tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                   jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3))/((block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                   jxOmin3:jxOmax3,3)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3,3))*block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3,1)*dsin(block%x(ixOmin1:ixOmax1,&
                   ixOmin2:ixOmax2,ixOmin3:ixOmax3,2)))
               
                end select
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+tmp2(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                else
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)-tmp2(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                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
                hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3)=block%surfaceC(ixCmin1:ixCmax1,&
                   ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir)*(qvec(ixCmin1:ixCmax1,&
                   ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                   kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3,jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                   jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3,kdir))/(block%x(jxCmin1:jxCmax1,&
                   jxCmin2:jxCmax2,jxCmin3:jxCmax3,&
                   jdir)-block%x(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3,jdir)))
                tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)=(tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3))/block%dvolume(ixOmin1:ixOmax1,&
                   ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+tmp2(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                else
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)-tmp2(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                 endif
                 if(idir<idirmin)idirmin=idir
              endif
            enddo; enddo;
            ! geometric terms
            if(idir==2.and.phi_>0) curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,2)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,2)+qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,phi_)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,r_)
            
            if(idir==phi_) curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,phi_)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,phi_)-qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,2)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,r_) +qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,r_)*dcos(block%x(ixOmin1:ixOmax1,&
               ixOmin2:ixOmax2,ixOmin3:ixOmax3,2))/(block%x(ixOmin1:ixOmax1,&
               ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
               r_)*dsin(block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
               ixOmin3:ixOmax3,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
                  hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                  hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                  hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                  jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                  jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                  ! qvec(3) at cell interface along 2nd dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir))
                  end if
                  ! 2nd coordinate at cell interface along 2nd dimension
                  xC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3)=block%x(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                     ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir)
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(dsin(xC(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3))*tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)-dsin(xC(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))*tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))*block%dx(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,kdir)
                  !! integral along 2nd dimension
                  hxOmin1=ixOmin1-kr(kdir,1);hxOmin2=ixOmin2-kr(kdir,2)
                  hxOmin3=ixOmin3-kr(kdir,3);hxOmax1=ixOmax1-kr(kdir,1)
                  hxOmax2=ixOmax2-kr(kdir,2);hxOmax3=ixOmax3-kr(kdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(kdir,1);jxCmin2=ixCmin2+kr(kdir,2)
                  jxCmin3=ixCmin3+kr(kdir,3);jxCmax1=ixCmax1+kr(kdir,1)
                  jxCmax2=ixCmax2+kr(kdir,2);jxCmax3=ixCmax3+kr(kdir,3);
                  ! qvec(2) at cell interface along 3rd dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,kdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir))
                  end if
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+(tmp(hxOmin1:hxOmax1,&
                     hxOmin2:hxOmax2,hxOmin3:hxOmax3)-tmp(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3))*block%dx(&
                     ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     jdir))/block%surface(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)*block%x(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,idir)
                end if
              case(2)
                if(jdir<kdir) then
                  ! idir=2,jdir=1,kdir=3
                  !! integral along 1st dimension
                  hxOmin1=ixOmin1-kr(kdir,1);hxOmin2=ixOmin2-kr(kdir,2)
                  hxOmin3=ixOmin3-kr(kdir,3);hxOmax1=ixOmax1-kr(kdir,1)
                  hxOmax2=ixOmax2-kr(kdir,2);hxOmax3=ixOmax3-kr(kdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(kdir,1);jxCmin2=ixCmin2+kr(kdir,2)
                  jxCmin3=ixCmin3+kr(kdir,3);jxCmax1=ixCmax1+kr(kdir,1)
                  jxCmax2=ixCmax2+kr(kdir,2);jxCmax3=ixCmax3+kr(kdir,3);
                  ! qvec(1) at cell interface along 3rd dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,kdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir))
                  end if
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))*block%dx(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,1)
                  !! integral along 3rd dimension
                  hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                  hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                  hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                  jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                  jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                  ! qvec(3) at cell interface along 1st dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir))
                  end if
                  ! 1st coordinate at cell interface along 1st dimension
                  xC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3)=block%x(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                     ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir)
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+(xC(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3)*tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3)-xC(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)*tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3))*dsin(block%x(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir))*block%dx(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,kdir))/block%surface(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,idir)
                end if
              case(3)
                if(jdir<kdir) then
                  ! idir=3,jdir=1,kdir=2
                  !! integral along 1st dimension
                  hxOmin1=ixOmin1-kr(kdir,1);hxOmin2=ixOmin2-kr(kdir,2)
                  hxOmin3=ixOmin3-kr(kdir,3);hxOmax1=ixOmax1-kr(kdir,1)
                  hxOmax2=ixOmax2-kr(kdir,2);hxOmax3=ixOmax3-kr(kdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(kdir,1);jxCmin2=ixCmin2+kr(kdir,2)
                  jxCmin3=ixCmin3+kr(kdir,3);jxCmax1=ixCmax1+kr(kdir,1)
                  jxCmax2=ixCmax2+kr(kdir,2);jxCmax3=ixCmax3+kr(kdir,3);
                  ! qvec(1) at cell interface along 2nd dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,kdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir))
                  end if
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3)-tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3))*block%dx(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,jdir)
                  !! integral along 2nd dimension
                  hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                  hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                  hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                  jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                  jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                  ! qvec(2) at cell interface along 1st dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir))
                  end if
                  ! 1st coordinate at cell interface along 1st dimension
                  xC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3)=block%x(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                     ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir)
                  if(ndim==3) then
                    surface(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                       ixOmin3:ixOmax3)=block%surface(ixOmin1:ixOmax1,&
                       ixOmin2:ixOmax2,ixOmin3:ixOmax3,idir)
                  else
                    surface(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                       ixOmin3:ixOmax3)=block%x(ixOmin1:ixOmax1,&
                       ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                       jdir)*block%dx(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                       ixOmin3:ixOmax3,kdir)*block%dx(ixOmin1:ixOmax1,&
                       ixOmin2:ixOmax2,ixOmin3:ixOmax3,jdir)
                  end if
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+(xC(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)*tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)-xC(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3)*tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))*block%dx(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     kdir))/surface(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)
                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(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3)=qvec(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                   ixAmin3:ixAmax3,kdir)
                hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                jxOmin1=ixOmin1+kr(jdir,1);jxOmin2=ixOmin2+kr(jdir,2)
                jxOmin3=ixOmin3+kr(jdir,3);jxOmax1=ixOmax1+kr(jdir,1)
                jxOmax2=ixOmax2+kr(jdir,2);jxOmax3=ixOmax3+kr(jdir,3);
                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(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                     ixAmin3:ixAmax3)=tmp(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                     ixAmin3:ixAmax3)*block%x(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                     ixAmin3:ixAmax3,1) !R V_phi
                  ! computes d(R V_phi)/dR or d V_Z/dR
                  tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                     jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))/(block%x(jxOmin1:jxOmax1,&
                     jxOmin2:jxOmax2,jxOmin3:jxOmax3,&
                     1)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3,1))
                  if(idir==z_) tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,1) !(1/R)*d(R V_phi)/dR
                        case(2)
                  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
                  tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                     jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))/((block%x(jxOmin1:jxOmax1,&
                     jxOmin2:jxOmax2,jxOmin3:jxOmax3,&
                     2)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3,2))*block%x(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,1))
                 
                      case(3)
                  ! handles d V_phi/dZ or d V_R/dZ
                  tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                     jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))/(block%x(jxOmin1:jxOmax1,&
                     jxOmin2:jxOmax2,jxOmin3:jxOmax3,&
                     3)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3,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(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                     ixAmin3:ixAmax3)=tmp(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                     ixAmin3:ixAmax3)*block%x(ixAmin1:ixAmax1,ixAmin2:ixAmax2,&
                     ixAmin3:ixAmax3,1) !R V_phi
                  ! computes d(R V_phi)/dR or d V_Z/dR
                  tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                     jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))/(block%x(jxOmin1:jxOmax1,&
                     jxOmin2:jxOmax2,jxOmin3:jxOmax3,&
                     1)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3,1))
                  if(idir==z_) tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,1) !(1/R)*d(R V_phi)/dR
                        case(2)
                  ! handles d V_phi/dZ or d V_R/dZ
                  tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                     jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))/(block%x(jxOmin1:jxOmax1,&
                     jxOmin2:jxOmax2,jxOmin3:jxOmax3,&
                     2)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3,2))
                 
                      case(3)
                  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
                  tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)=(tmp(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
                     jxOmin3:jxOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))/((block%x(jxOmin1:jxOmax1,&
                     jxOmin2:jxOmax2,jxOmin3:jxOmax3,&
                     3)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3,3))*block%x(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,1))
                 
                  end select
                end if
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+tmp2(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                else
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)-tmp2(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                 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
                hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3)=block%surfaceC(ixCmin1:ixCmax1,&
                   ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir)*(qvec(ixCmin1:ixCmax1,&
                   ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                   kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3,jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                   jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3,kdir))/(block%x(jxCmin1:jxCmax1,&
                   jxCmin2:jxCmax2,jxCmin3:jxCmax3,&
                   jdir)-block%x(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                   ixCmin3:ixCmax3,jdir)))
                tmp2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)=(tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                   hxOmin3:hxOmax3))/block%dvolume(ixOmin1:ixOmax1,&
                   ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                if(lvc(idir,jdir,kdir)==1)then
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+tmp2(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                else
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)-tmp2(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3)
                 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(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3,idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3,idir)-qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3,z_)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3,r_)
              ! polar
              if(phi_==2) curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3,idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3,idir)+qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3,z_)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                 ixOmin3:ixOmax3,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
                  hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                  hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                  hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                  jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                  jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                  ! qvec(z) at cell interface along phi dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir))
                  end if
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))*block%dx(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,kdir)
                  !! integral along phi dimension
                  hxOmin1=ixOmin1-kr(kdir,1);hxOmin2=ixOmin2-kr(kdir,2)
                  hxOmin3=ixOmin3-kr(kdir,3);hxOmax1=ixOmax1-kr(kdir,1)
                  hxOmax2=ixOmax2-kr(kdir,2);hxOmax3=ixOmax3-kr(kdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(kdir,1);jxCmin2=ixCmin2+kr(kdir,2)
                  jxCmin3=ixCmin3+kr(kdir,3);jxCmax1=ixCmax1+kr(kdir,1)
                  jxCmax2=ixCmax2+kr(kdir,2);jxCmax3=ixCmax3+kr(kdir,3);
                  ! qvec(phi) at cell interface along z dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,kdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir))
                  end if
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+(tmp(hxOmin1:hxOmax1,&
                     hxOmin2:hxOmax2,hxOmin3:hxOmax3)-tmp(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3))*block%x(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)*block%dx(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,jdir))/block%surface(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,idir)
                end if
               else if(idir==phi_) then
                if(jdir<kdir) then
                  ! idir=phi,jdir=r,kdir=z
                  !! integral along r dimension
                  hxOmin1=ixOmin1-kr(kdir,1);hxOmin2=ixOmin2-kr(kdir,2)
                  hxOmin3=ixOmin3-kr(kdir,3);hxOmax1=ixOmax1-kr(kdir,1)
                  hxOmax2=ixOmax2-kr(kdir,2);hxOmax3=ixOmax3-kr(kdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(kdir,1);jxCmin2=ixCmin2+kr(kdir,2)
                  jxCmin3=ixCmin3+kr(kdir,3);jxCmax1=ixCmax1+kr(kdir,1)
                  jxCmax2=ixCmax2+kr(kdir,2);jxCmax3=ixCmax3+kr(kdir,3);
                  ! qvec(r) at cell interface along z dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,kdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir))
                  end if
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))*block%dx(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,jdir)
                  !! integral along z dimension
                  hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                  hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                  hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                  jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                  jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                  ! qvec(z) at cell interface along r dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir))
                  end if
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+(tmp(hxOmin1:hxOmax1,&
                     hxOmin2:hxOmax2,hxOmin3:hxOmax3)-tmp(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3))*block%dx(&
                     ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     kdir))/block%surface(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)
                end if
               else ! idir==z_
                if(jdir<kdir) then
                  ! idir=z,jdir=r,kdir=phi
                  !! integral along r dimension
                  hxOmin1=ixOmin1-kr(kdir,1);hxOmin2=ixOmin2-kr(kdir,2)
                  hxOmin3=ixOmin3-kr(kdir,3);hxOmax1=ixOmax1-kr(kdir,1)
                  hxOmax2=ixOmax2-kr(kdir,2);hxOmax3=ixOmax3-kr(kdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(kdir,1);jxCmin2=ixCmin2+kr(kdir,2)
                  jxCmin3=ixCmin3+kr(kdir,3);jxCmax1=ixCmax1+kr(kdir,1)
                  jxCmax2=ixCmax2+kr(kdir,2);jxCmax3=ixCmax3+kr(kdir,3);
                  ! qvec(r) at cell interface along phi dimension
                  if(stretched_dim(kdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,kdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir))
                  end if
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3)-tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3))*block%dx(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,jdir)
                  !! integral along phi dimension
                  hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2)
                  hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1)
                  hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3);
                  ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3
                  ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3;
                  jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2)
                  jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1)
                  jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3);
                  ! qvec(phi) at cell interface along r dimension
                  if(stretched_dim(jdir) .and. stretch_uncentered) then
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3,kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,jdir)-block%x(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir))
                  else
                    tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                       ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,&
                       ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
                       kdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
                       jxCmin3:jxCmax3,kdir))
                  end if
                  ! r coordinate at cell interface along r dimension
                  xC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3)=block%x(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
                     ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,&
                     ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir)
                  curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     idir)=(curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,idir)+(xC(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)*tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3)-xC(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3)*tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
                     hxOmin3:hxOmax3))*block%dx(ixOmin1:ixOmax1,&
                     ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                     kdir))/block%surface(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                     ixOmin3:ixOmax3,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