mod_convert_files.t

Go to the documentation of this file.

module mod_convert_files
  use mod_comm_lib, only: mpistop
 
  implicit none
  public
 
contains
 
  subroutine generate_plotfile
    use mod_usr_methods, only: usr_special_convert
    use mod_global_parameters
    use mod_ghostcells_update
    use mod_physics, only: phys_te_images
    use mod_convert, only: convert_all
    use mod_thermal_emission
    use mod_magnetic_helicity, only: mh_run_task
    use mod_magnetic_topology, only: mt_run_topology_task
 
    character(len=std_len) :: convert_type_elem
    integer :: i
 
    select case(convert_type)
      case('tecplot','tecplotCC','tecline')
       call tecplot(unitconvert)
      case('tecplotmpi','tecplotCCmpi','teclinempi')
       call tecplot_mpi(unitconvert)
      case('vtu','vtuCC')
       call unstructuredvtk(unitconvert)
      case('vtumpi','vtuCCmpi')
       call unstructuredvtk_mpi(unitconvert)
      case('vtuB','vtuBCC','vtuBmpi','vtuBCCmpi')
       call unstructuredvtkb(unitconvert)
      case('vtuB64','vtuBCC64','vtuBmpi64','vtuBCCmpi64')
       call unstructuredvtkb64(unitconvert)
      {^iftwod
      case('vtuB23','vtuBCC23')
       call unstructuredvtkb23(unitconvert)
      case('vtuBsym23','vtuBCCsym23')
       call unstructuredvtkbsym23(unitconvert)
      }
      case('pvtumpi','pvtuCCmpi')
       call punstructuredvtk_mpi(unitconvert)
      case('pvtuBmpi','pvtuBCCmpi')
       call punstructuredvtkb_mpi(unitconvert)
      case('vtimpi','vtiCCmpi')
       call imagedatavtk_mpi(unitconvert)
      case('onegrid','onegridmpi')
       call onegrid(unitconvert)
      case('oneblock','oneblockB')
       call oneblock(unitconvert)
      case('EIvtiCCmpi','ESvtiCCmpi','SIvtiCCmpi','WIvtiCCmpi','EIvtuCCmpi','ESvtuCCmpi','SIvtuCCmpi','WIvtuCCmpi')
        ! output synthetic euv emission
        if (ndim==3 .and. associated(phys_te_images)) then
          call phys_te_images
        endif
      case('dat_generic_mpi')
        call convert_all()
      case('magnetic_helicity')
        call mh_run_task()
      case('magnetic_topology')
        call mt_run_topology_task()
      case('user','usermpi')
         if (.not. associated(usr_special_convert)) then
            call mpistop("usr_special_convert not defined")
         else
            call usr_special_convert(unitconvert)
         end if
      case default
       call mpistop("Error in generate_plotfile: Unknown convert_type")
    end select
 
  end subroutine generate_plotfile
 
  subroutine oneblock(qunit)
    ! this is for turning an AMR run into a single block
    ! the data will be all on selected level level_io
    ! this version should work for any dimension
    ! only writes w_write selected 1:nw variables, also nwauxio
    ! may use saveprim to switch to primitives
    ! this version can not work on multiple CPUs
    ! does not renormalize variables
    ! header info differs from onegrid below 
    ! ASCII or binary output
 
    use mod_forest, only: morton_start, morton_stop, sfc_to_igrid, igrid_to_node
    use mod_global_parameters
    use mod_usr_methods, only: usr_aux_output
    use mod_physics
    use mod_calculate_xw
    use mod_amr_solution_node, only: alloc_state_output
    integer, intent(in) :: qunit
 
    double precision:: normconv(0:nw+nwauxio)
    integer             :: Morton_no,igrid,ix^D,ig^D,level
    integer, pointer    :: ig_to_igrid(:^D&,:)
    integer             :: filenr,ncells^D,ncellx^D,jg^D,jig^D
    integer           :: iw,iiw,writenw,iwrite(1:nw+nwauxio),iigrid,idim
    logical             :: fileopen,writeblk(max_blocks)
    logical :: patchw(ixG^T)
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
    character(len=80)   :: filename
 
    if(level_io<1)then
     call mpistop('please specify level_io>0 for usage with oneblock')
    end if
 
    if(autoconvert)then
     call mpistop('Set autoconvert=F and convert oneblock data manually')
    end if
 
    if(npe>1)then
     if(mype==0) print *,'ONEBLOCK as yet to be parallelized'
     call mpistop('npe>1, oneblock')
    end if
 
    ! only variables selected by w_write will be written out
    normconv(0:nw+nwauxio)=one
    normconv(0) = length_convert_factor
    normconv(1:nw) = w_convert_factor
    writenw=count(w_write(1:nw))+nwauxio
    iiw=0
    do iw =1,nw
     if (.not.w_write(iw))cycle
     iiw=iiw+1
     iwrite(iiw)=iw
    end do
    if(nwauxio>0)then
      do iw =nw+1,nw+nwauxio
       iiw=iiw+1
       iwrite(iiw)=iw
      end do
    end if
 
    allocate(ig_to_igrid(ng^d(level_io),0:npe-1))
    ig_to_igrid=-1
    writeblk=.false.
    do morton_no=morton_start(mype),morton_stop(mype)
      igrid=sfc_to_igrid(morton_no)
      level=node(plevel_,igrid)
      ig^d=igrid_to_node(igrid,mype)%node%ig^d;
      ig_to_igrid(ig^d,mype)=igrid
      if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
            *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
           <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
        writeblk(igrid)=.true.
      end if
    end do
 
    call getheadernames(wnamei,xandwnamei,outfilehead)
    ncells^d=0;
    ncellx^d=ixmhi^d-ixmlo^d+1\
    {do ig^d=1,ng^d(level_io)\}
      igrid=ig_to_igrid(ig^d,mype)
      if(writeblk(igrid)) go to 20
    {end do\}
    20 continue
    jg^d=ig^d;
    {
    jig^dd=jg^dd;
    do ig^d=1,ng^d(level_io)
      jig^d=ig^d
      igrid=ig_to_igrid(jig^dd,mype)
      if(writeblk(igrid)) ncells^d=ncells^d+ncellx^d
    end do
    \}
 
    do iigrid=1,igridstail; igrid=igrids(iigrid)
      if(.not.writeblk(igrid)) cycle
      block=>ps(igrid)
      call alloc_state_output(igrid,ps1(igrid),ixg^ll)
      ps1(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
      if(nwauxio > 0) then
        if (.not. associated(usr_aux_output)) then
          call mpistop("usr_aux_output not defined")
        else
          call usr_aux_output(ixg^ll,ixm^ll^ladd1,ps1(igrid)%w,ps(igrid)%x,normconv)
        end if
      end if
      if(saveprim) then
        call phys_to_primitive(ixg^ll,ixg^ll^lsub2,ps1(igrid)%w,ps(igrid)%x)
        if(b0field) ps1(igrid)%w(ixg^t,iw_mag(:))=ps1(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
      else
        if(b0field) then
          ! add background magnetic field B0 to B
          if(phys_total_energy) then
            ps1(igrid)%w(ixg^t,iw_e)=ps1(igrid)%w(ixg^t,iw_e)+0.5d0*sum(ps(igrid)%B0(ixg^t,:,0)**2,dim=ndim+1) &
                + sum(ps1(igrid)%w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=ndim+1)
          end if
          ps1(igrid)%w(ixg^t,iw_mag(:))=ps1(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
        end if
      end if
    end do
 
    master_cpu_open : if (mype == 0) then
     inquire(qunit,opened=fileopen)
     if (.not.fileopen) then
       ! generate filename
        filenr=snapshotini
        if (autoconvert) filenr=snapshotnext
       write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
       select case(convert_type)
        case("oneblock")
         open(qunit,file=filename,status='unknown')
         write(qunit,*) trim(outfilehead)
         write(qunit,*) ncells^d
         write(qunit,*) real(global_time*time_convert_factor)
        case("oneblockB")
         open(qunit,file=filename,form='unformatted',status='unknown')
         write(qunit) outfilehead
         write(qunit) ncells^d
         write(qunit) real(global_time*time_convert_factor)
       end select
     end if
    end if master_cpu_open
 
    {^ifthreed
    do ig3=1,ng3(level_io)
     do ix3=ixmlo3,ixmhi3}
 
       {^nooned
       do ig2=1,ng2(level_io)
         do ix2=ixmlo2,ixmhi2}
 
           do ig1=1,ng1(level_io)
             igrid=ig_to_igrid(ig^d,mype)
             if(.not.writeblk(igrid)) cycle
             do ix1=ixmlo1,ixmhi1
               master_write : if(mype==0) then
                 select case(convert_type)
                   case("oneblock")
                     write(qunit,fmt="(100(e14.6))") &
                      ps(igrid)%x(ix^d,1:ndim)*normconv(0),&
                      (ps1(igrid)%w(ix^d,iwrite(iw))*normconv(iwrite(iw)),iw=1,writenw)
                   case("oneblockB")
                     write(qunit) real(ps(igrid)%x(ix^d,1:ndim)*normconv(0)),&
                      (real(ps1(igrid)%w(ix^d,iwrite(iw))*normconv(iwrite(iw))),iw=1,writenw)
                 end select
               end if master_write
             end do
           end do
        {^nooned
         end do
       end do}
     {^ifthreed
     end do
    end do}
 
    close(qunit)
 
  end subroutine oneblock
 
  subroutine onegrid(qunit)
    ! this is for turning an AMR run into a single grid
    ! this version should work for any dimension, can be in parallel
    ! in 1D, should behave much like oneblock, except for header info
    
    ! only writes all 1:nw variables, no nwauxio
    ! may use saveprim to switch to primitives
    ! this version can work on multiple CPUs
    ! does not renormalize variables
    ! ASCII output
 
    use mod_forest, only: morton_start, morton_stop, sfc_to_igrid
    use mod_global_parameters
    use mod_physics
    use mod_calculate_xw
    integer, intent(in) :: qunit
 
    double precision  :: w_recv(ixG^T,1:nw),x_recv(ixG^T,1:ndim)
    integer             :: itag,Morton_no,igrid,ix^D,iw
    integer             :: filenr
    !.. MPI variables ..
    integer           :: igrid_recv,ipe
    integer, allocatable :: intstatus(:,:)
    logical             :: fileopen
    character(len=80)   :: filename
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
 
 
    if(nwauxio>0)then
      if(mype==0) print *,'ONEGRID to be used without nwauxio'
      call mpistop('nwauxio>0, onegrid')
    end if
 
    if(saveprim)then
      if(mype==0.and.nwaux>0) print *,'warning: ONEGRID used with saveprim, check auxiliaries'
    end if
 
    master_cpu_open : if (mype == 0) then
      call getheadernames(wnamei,xandwnamei,outfilehead)
      write(outfilehead,'(a)') "#"//" "//trim(outfilehead)
      inquire(qunit,opened=fileopen)
      if (.not.fileopen) then
        ! generate filename
         filenr=snapshotini
         if (autoconvert) filenr=snapshotnext
        write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
        open(qunit,file=filename,status='unknown')
      end if
      write(qunit,"(a)")outfilehead
      write(qunit,"(i7)") ( {^d&(ixmhi^d-ixmlo^d+1)*} )*(morton_stop(npe-1)-morton_start(0)+1)
    end if master_cpu_open
 
    do morton_no=morton_start(mype),morton_stop(mype)
      igrid=sfc_to_igrid(morton_no)
      if(saveprim) call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
      if(mype/=0)then
        itag=morton_no
        call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
        call mpi_send(ps(igrid)%x,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
        itag=igrid
        call mpi_send(ps(igrid)%w,1,type_block_io, 0,itag,icomm,ierrmpi)
      else
       {do ix^db=ixmlo^db,ixmhi^db\}
          do iw=1,nw
            if( dabs(ps(igrid)%w(ix^d,iw)) < 1.0d-32 ) ps(igrid)%w(ix^d,iw) = zero
          end do
          write(qunit,fmt="(100(e14.6))") ps(igrid)%x(ix^d,1:ndim),ps(igrid)%w(ix^d,1:nw)
       {end do\}
      end if
    end do
 
    if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
 
    manycpu : if (npe>1) then
      if (mype==0) then
       loop_cpu : do ipe =1, npe-1
        loop_morton : do morton_no=morton_start(ipe),morton_stop(ipe)
              itag=morton_no
              call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              call mpi_recv(x_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              itag=igrid_recv
              call mpi_recv(w_recv,1,type_block_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              {do ix^db=ixmlo^db,ixmhi^db\}
                 do iw=1,nw
                   if( dabs(ps(igrid)%w(ix^d,iw)) < smalldouble ) ps(igrid)%w(ix^d,iw) = zero
                 end do
                 write(qunit,fmt="(100(e14.6))") x_recv(ix^d,1:ndim),w_recv(ix^d,1:nw)
              {end do\}
        end do loop_morton
       end do loop_cpu
      end if
    end if manycpu
 
    if (npe>1) then
      call mpi_barrier(icomm,ierrmpi)
      if(mype==0)deallocate(intstatus)
    end if
 
    if(mype==0) close(qunit)
  end subroutine onegrid 
 
  subroutine tecplot(qunit)
 
    ! output for tecplot (ASCII format)
    ! not parallel, uses calc_grid to compute nwauxio variables
    ! allows renormalizing using convert factors
 
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) :: qunit
 
    double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC_TMP
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC_TMP
    double precision, dimension(0:nw+nwauxio)                   :: normconv
    integer::               igrid,iigrid,level,igonlevel,iw,idim,ix^D
    integer::               NumGridsOnLevel(1:nlevelshi)
    integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
    integer ::              nodes, elems
    integer  :: filenr
    logical :: fileopen,first
    character(len=80) :: filename
    !!! possible length conflict
    character(len=1024) :: tecplothead
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
 
    if(npe>1)then
      if(mype==0) print *,'tecplot not parallel, use tecplotmpi'
      call mpistop('npe>1, tecplot')
    end if
 
    if(nw/=count(w_write(1:nw)))then
      if(mype==0) print *,'tecplot does not use w_write=F'
      call mpistop('w_write, tecplot')
    end if
 
    if(nocartesian)then
      if(mype==0) print *,'tecplot with nocartesian'
    end if
 
    inquire(qunit,opened=fileopen)
    if(.not.fileopen) then
      ! generate filename
      filenr=snapshotini
      if (autoconvert) filenr=snapshotnext
      write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
      open(qunit,file=filename,status='unknown')
    end if
 
    call getheadernames(wnamei,xandwnamei,outfilehead)
 
    write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
    write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
 
    numgridsonlevel(1:nlevelshi)=0
    do level=levmin,levmax
      numgridsonlevel(level)=0
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        if (node(plevel_,igrid)/=level) cycle
        numgridsonlevel(level)=numgridsonlevel(level)+1
      end do
    end do
 
    nx^d=ixmhi^d-ixmlo^d+1;
    nxc^d=nx^d+1;
 
    {^ifoned
    if(convert_type=='tecline') then
      nodes=0
      elems=0
      do level=levmin,levmax
        nodes=nodes + numgridsonlevel(level)*{nxc^d*}
        elems=elems + numgridsonlevel(level)*{nx^d*}
      end do
 
      write(qunit,"(a,i7,a,1pe12.5,a)") &
            'ZONE T="all levels", I=',elems, &
            ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT' 
 
      igonlevel=0
      do iigrid=1,igridstail; igrid=igrids(iigrid);
         block=>ps(igrid)
         call calc_x(igrid,xc,xcc)
         call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
        {do ix^db=ixccmin^db,ixccmax^db\}
           x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
           w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
           write(qunit,fmt="(100(e24.16))") x_tec, w_tec
        {end do\}
       end do
       close(qunit)
    else
    }
    do level=levmin,levmax
      nodesonlevel=numgridsonlevel(level)*{nxc^d*}
      elemsonlevel=numgridsonlevel(level)*{nx^d*}
      ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
      ! with the AMR grid LEVEL. Other options would be
      !    let each grid define a zone: inefficient for TECPLOT internal workings
      !       hence not implemented
      !    let entire octree define 1 zone: no difference in interpolation 
      !       properties across TECPLOT zones detected as yet, hence not done
      select case(convert_type)
        case('tecplot')
          ! in this option, we store the corner coordinates, as well as the corner
          ! values of all variables (obtained by averaging). This allows POINT packaging, 
          ! and thus we can save full grid info by using one call to calc_grid
          write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
               'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
               ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
            {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
          do iigrid=1,igridstail; igrid=igrids(iigrid);
            if (node(plevel_,igrid)/=level) cycle
            block=>ps(igrid)
            call calc_x(igrid,xc,xcc)
            call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                           ixc^l,ixcc^l,.true.)
            {do ix^db=ixcmin^db,ixcmax^db\}
               x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
               w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
               write(qunit,fmt="(100(e14.6))") x_tec, w_tec
            {end do\}
          end do
        case('tecplotCC')
          ! in this option, we store the corner coordinates, and the cell center
          ! values of all variables. Due to this mix of corner/cell center, we must 
          ! use BLOCK packaging, and thus we have enormous overhead by using 
          ! calc_grid repeatedly to merely fill values of cell corner coordinates 
          ! and cell center values per dimension, per variable
          if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
          if(nw+nwauxio==1)then
            ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
            ! and just set [ndim+1]
            write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
               'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
               ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
               ndim+1,']=CELLCENTERED), ZONETYPE=', &
            {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
          else
            if(ndim+nw+nwauxio<10) then
             ! difference only in length of integer format specification for ndim+nw+nwauxio
             write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
                'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
                ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
                ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
             {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
            else
             write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
                'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
                ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
                ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
             {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
            end if
          end if
          do idim=1,ndim
            first=(idim==1) 
            do iigrid=1,igridstail; igrid=igrids(iigrid);
              if (node(plevel_,igrid)/=level) cycle
              block=>ps(igrid)
              call calc_x(igrid,xc,xcc)
              call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                             ixc^l,ixcc^l,first)
              write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
            end do
          end do
          do iw=1,nw+nwauxio
            do iigrid=1,igridstail; igrid=igrids(iigrid);
               if (node(plevel_,igrid)/=level) cycle
               block=>ps(igrid)
               call calc_x(igrid,xc,xcc)
               call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                              ixc^l,ixcc^l,.true.)
               write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
            enddo
          enddo
        case default
          call mpistop('no such tecplot type')
      end select
      igonlevel=0
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        if (node(plevel_,igrid)/=level) cycle
        block=>ps(igrid)
        igonlevel=igonlevel+1
        call save_conntec(qunit,igrid,igonlevel)
      end do
    end do
    {^ifoned endif}
 
    close(qunit)
 
  end subroutine tecplot
 
  subroutine save_conntec(qunit,igrid,igonlevel)
 
    ! this saves the basic line, quad and brick connectivity,
    ! as used by TECPLOT file outputs for unstructured grid
    use mod_global_parameters
 
    integer, intent(in) :: qunit, igrid, igonlevel
 
    integer :: nx^D, nxC^D, ix^D
 
    nx^d=ixmhi^d-ixmlo^d+1;
    nxc^d=nx^d+1;
 
    ! connectivity list
    {do ix^db=1,nx^db\}
       {^ifthreed
       ! basic brick connectivity
       write(qunit,'(8(i7,1x))') &
          nodenumbertec3d(ix1,  ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
          nodenumbertec3d(ix1+1,ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
          nodenumbertec3d(ix1+1,ix2  ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
          nodenumbertec3d(ix1  ,ix2  ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
          nodenumbertec3d(ix1  ,ix2-1,ix3  ,nxc1,nxc2,nxc3,igonlevel,igrid),&
          nodenumbertec3d(ix1+1,ix2-1,ix3  ,nxc1,nxc2,nxc3,igonlevel,igrid),&
          nodenumbertec3d(ix1+1,ix2  ,ix3  ,nxc1,nxc2,nxc3,igonlevel,igrid),&
          nodenumbertec3d(ix1  ,ix2  ,ix3  ,nxc1,nxc2,nxc3,igonlevel,igrid)
       }
       {^iftwod
       ! basic quadrilateral connectivity
       write(qunit,'(4(i7,1x))') &
          nodenumbertec2d(ix1,  ix2-1,nxc1,nxc2,igonlevel,igrid),&
          nodenumbertec2d(ix1+1,ix2-1,nxc1,nxc2,igonlevel,igrid),&
          nodenumbertec2d(ix1+1,ix2  ,nxc1,nxc2,igonlevel,igrid),&
          nodenumbertec2d(ix1  ,ix2  ,nxc1,nxc2,igonlevel,igrid)
       }
       {^ifoned
       ! basic line connectivity
       write(qunit,'(2(i7,1x))') nodenumbertec1d(ix1,nxc1,igonlevel,igrid),&
                                 nodenumbertec1d(ix1+1,nxc1,igonlevel,igrid)
       }
    {end do\}
 
  end subroutine save_conntec
 
  integer function nodenumbertec1d(i1,nx1,ig,igrid)
    use mod_comm_lib, only: mpistop
    integer, intent(in):: i1,nx1,ig,igrid
 
    nodenumbertec1d=i1+(ig-1)*nx1
    if(nodenumbertec1d>9999999)call mpistop("too large nodenumber")
 
  end function nodenumbertec1d
 
  integer function nodenumbertec2d(i1,i2,nx1,nx2,ig,igrid)
 
    integer, intent(in):: i1,i2,nx1,nx2,ig,igrid
 
    nodenumbertec2d=i1+i2*nx1+(ig-1)*nx1*nx2
    if(nodenumbertec2d>9999999)call mpistop("too large nodenumber")
 
  end function nodenumbertec2d
 
  integer function nodenumbertec3d(i1,i2,i3,nx1,nx2,nx3,ig,igrid)
 
    integer, intent(in):: i1,i2,i3,nx1,nx2,nx3,ig,igrid
 
    nodenumbertec3d=i1+i2*nx1+i3*nx1*nx2+(ig-1)*nx1*nx2*nx3
    if(nodenumbertec3d>9999999)call mpistop("too large nodenumber")
 
  end function nodenumbertec3d
 
  subroutine unstructuredvtk(qunit)
 
    ! output for vtu format to paraview
    ! not parallel, uses calc_grid to compute nwauxio variables
    ! allows renormalizing using convert factors
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) ::    qunit
 
    double precision ::  x_VTK(1:3)
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC_TMP
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC_TMP
    double precision, dimension(0:nw+nwauxio)                   :: normconv
    integer::               igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,iw
    integer::               NumGridsOnLevel(1:nlevelshi)
    integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
    integer          :: filenr
    character(len=80)::  filename
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
    logical :: fileopen
 
    if(npe>1)then
      if(mype==0) print *,'unstructuredvtk not parallel, use vtumpi'
      call mpistop('npe>1, unstructuredvtk')
    end if
 
    inquire(qunit,opened=fileopen)
    if(.not.fileopen)then
      ! generate filename 
      filenr=snapshotini
      if (autoconvert) filenr=snapshotnext
      write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
      ! Open the file for the header part
      open(qunit,file=filename,status='unknown')
    end if
 
    call getheadernames(wnamei,xandwnamei,outfilehead)
 
    ! generate xml header
    write(qunit,'(a)')'<?xml version="1.0"?>'
    write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
    write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
    write(qunit,'(a)')'<UnstructuredGrid>'
    write(qunit,'(a)')'<FieldData>'
    write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                       'NumberOfTuples="1" format="ascii">'
    write(qunit,*) dble(dble(global_time)*time_convert_factor)
    write(qunit,'(a)')'</DataArray>'
    write(qunit,'(a)')'</FieldData>'
 
    ! number of cells, number of corner points, per grid.
    nx^d=ixmhi^d-ixmlo^d+1;
    nxc^d=nx^d+1;
    nc={nx^d*}
    np={nxc^d*}
 
    ! Note: using the w_write, writelevel, writespshift
    ! we can clip parts of the grid away, select variables, levels etc.
    do level=levmin,levmax
      if (writelevel(level)) then
        do iigrid=1,igridstail; igrid=igrids(iigrid);
          if (node(plevel_,igrid)/=level) cycle
          block=>ps(igrid)
          ! only output a grid when fully within clipped region selected
          ! by writespshift array
          if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
               *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
              <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
            call calc_x(igrid,xc,xcc)
            call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                          ixc^l,ixcc^l,.true.)
            select case(convert_type)
            case('vtu')
              ! we write out every grid as one VTK PIECE
              write(qunit,'(a,i7,a,i7,a)') &
                 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
              write(qunit,'(a)')'<PointData>'
              do iw=1,nw+nwauxio
                if(iw<=nw) then 
                  if(.not.w_write(iw)) cycle
                end if
                write(qunit,'(a,a,a)')&
                '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
                write(qunit,'(200(1pe14.6))') {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
                write(qunit,'(a)')'</DataArray>'
              end do
              write(qunit,'(a)')'</PointData>'
              write(qunit,'(a)')'<Points>'
              write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
              ! write cell corner coordinates in a backward dimensional loop, always 3D output
              {do ix^db=ixcmin^db,ixcmax^db \}
                 x_vtk(1:3)=zero;
                 x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
                 write(qunit,'(3(1pe14.6))') x_vtk
              {end do \}
              write(qunit,'(a)')'</DataArray>'
              write(qunit,'(a)')'</Points>'
            case('vtuCC')
              ! we write out every grid as one VTK PIECE
              write(qunit,'(a,i7,a,i7,a)') &
                 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
              write(qunit,'(a)')'<CellData>'
              do iw=1,nw+nwauxio
                if(iw<=nw) then 
                  if(.not.w_write(iw)) cycle
                end if
                write(qunit,'(a,a,a)')&
                '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
                write(qunit,'(200(1pe14.6))') {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
                write(qunit,'(a)')'</DataArray>'
              end do
              write(qunit,'(a)')'</CellData>'
              write(qunit,'(a)')'<Points>'
              write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
              ! write cell corner coordinates in a backward dimensional loop, always 3D output
              {do ix^db=ixcmin^db,ixcmax^db \}
                 x_vtk(1:3)=zero;
                 x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
                 write(qunit,'(3(1pe14.6))') x_vtk
              {end do \}
              write(qunit,'(a)')'</DataArray>'
              write(qunit,'(a)')'</Points>'
            end select
    
            write(qunit,'(a)')'<Cells>'
            ! connectivity part
            write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
            call save_connvtk(qunit,igrid)
            write(qunit,'(a)')'</DataArray>'
    
            ! offsets data array
            write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
            do icel=1,nc
              write(qunit,'(i7)') icel*(2**^nd)
            end do
            write(qunit,'(a)')'</DataArray>'
    
            ! VTK cell type data array
            write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
            vtk_type=vtk_cell_type()
            do icel=1,nc
              write(qunit,'(i2)') vtk_type
            end do
            write(qunit,'(a)')'</DataArray>'
    
            write(qunit,'(a)')'</Cells>'
    
            write(qunit,'(a)')'</Piece>'
          end if
        end do
      end if
    end do
 
    write(qunit,'(a)')'</UnstructuredGrid>'
    write(qunit,'(a)')'</VTKFile>'
    close(qunit)
 
  end subroutine unstructuredvtk
 
  subroutine unstructuredvtkb(qunit)
 
    ! output for vtu format to paraview, binary version output
    ! not parallel, uses calc_grid to compute nwauxio variables
    ! allows renormalizing using convert factors
    use mod_forest, only: morton_start, morton_stop, sfc_to_igrid
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) ::    qunit
 
    double precision ::  x_VTK(1:3)
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC_TMP
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio):: wC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)  :: wCC_TMP
    double precision :: normconv(0:nw+nwauxio)
    integer, allocatable :: intstatus(:,:)
    integer*8 :: offset
    integer :: itag,ipe,igrid,level,icel,ixC^L,ixCC^L,Morton_no,Morton_length
    integer :: nx^D,nxC^D,nc,np,VTK_type,ix^D,filenr
    integer::  k,iw
    integer::  length,lengthcc,length_coords,length_conn,length_offsets
    character::  buf
    character(len=80)::  filename
    character(len=19)::  offset_char
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
    logical ::   fileopen,cell_corner=.false.
    logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
 
    normconv=one
    morton_length=morton_stop(npe-1)-morton_start(0)+1
    allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
    allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
    morton_aim=.false.
    morton_aim_p=.false.
    do morton_no=morton_start(mype),morton_stop(mype)
      igrid=sfc_to_igrid(morton_no)
      level=node(plevel_,igrid)
      ! we can clip parts of the grid away, select variables, levels etc.
      if(writelevel(level)) then
        ! only output a grid when fully within clipped region selected
        ! by writespshift array
        if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
              *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
             <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
          morton_aim_p(morton_no)=.true.
        end if
      end if
    end do
    call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
                             icomm,ierrmpi)
    select case(convert_type)
     case('vtuB','vtuBmpi')
       cell_corner=.true.
     case('vtuBCC','vtuBCCmpi')
       cell_corner=.false.
    end select
    if (mype /= 0) then
      do morton_no=morton_start(mype),morton_stop(mype)
        if(.not. morton_aim(morton_no)) cycle
        igrid=sfc_to_igrid(morton_no)
        call calc_x(igrid,xc,xcc)
        call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                          ixc^l,ixcc^l,.true.)
        itag=morton_no
        call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
        if(cell_corner) then
          call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
        else
          call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
        end if
      end do
 
    else
      ! mype==0
      offset=0
      inquire(qunit,opened=fileopen)
      if(.not.fileopen)then
        ! generate filename 
        filenr=snapshotini
        if (autoconvert) filenr=snapshotnext
        write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
        ! Open the file for the header part
        open(qunit,file=filename,status='replace')
      end if
      call getheadernames(wnamei,xandwnamei,outfilehead)
      ! generate xml header
      write(qunit,'(a)')'<?xml version="1.0"?>'
      write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
      write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
      write(qunit,'(a)')'<UnstructuredGrid>'
      write(qunit,'(a)')'<FieldData>'
      write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                         'NumberOfTuples="1" format="ascii">'
      write(qunit,*) real(global_time*time_convert_factor)
      write(qunit,'(a)')'</DataArray>'
      write(qunit,'(a)')'</FieldData>'
 
      ! number of cells, number of corner points, per grid.
      nx^d=ixmhi^d-ixmlo^d+1;
      nxc^d=nx^d+1;
      nc={nx^d*}
      np={nxc^d*}
      length=np*size_real
      lengthcc=nc*size_real
      length_coords=3*length
      length_conn=2**^nd*size_int*nc
      length_offsets=nc*size_int
 
      ! Note: using the w_write, writelevel, writespshift
      do morton_no=morton_start(0),morton_stop(0)
        if(.not. morton_aim(morton_no)) cycle
        if(cell_corner) then
          ! we write out every grid as one VTK PIECE
          write(qunit,'(a,i7,a,i7,a)') &
             '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
          write(qunit,'(a)')'<PointData>'
          do iw=1,nw+nwauxio
            if(iw<=nw) then 
              if(.not.w_write(iw)) cycle
            endif
            write(offset_char,'(i19)') offset
            write(qunit,'(a,a,a,a,a)')&
                '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
                '" format="appended" offset="',trim(adjustl(offset_char)),'">'
            write(qunit,'(a)')'</DataArray>'
            offset=offset+length+size_int
          end do
          write(qunit,'(a)')'</PointData>'
          write(qunit,'(a)')'<Points>'
          write(offset_char,'(i19)') offset
          write(qunit,'(a,a,a)') &
          '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
          ! write cell corner coordinates in a backward dimensional loop, always 3D output
          offset=offset+length_coords+size_int
          write(qunit,'(a)')'</Points>'
        else
          ! we write out every grid as one VTK PIECE
          write(qunit,'(a,i7,a,i7,a)') &
             '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
          write(qunit,'(a)')'<CellData>'
          do iw=1,nw+nwauxio
            if(iw<=nw) then 
               if(.not.w_write(iw)) cycle
            end if
            write(offset_char,'(i19)') offset
            write(qunit,'(a,a,a,a,a)')&
                '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
                '" format="appended" offset="',trim(adjustl(offset_char)),'">'
            write(qunit,'(a)')'</DataArray>'
            offset=offset+lengthcc+size_int
          end do
          write(qunit,'(a)')'</CellData>'
          write(qunit,'(a)')'<Points>'
          write(offset_char,'(i19)') offset
          write(qunit,'(a,a,a)') &
          '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
          ! write cell corner coordinates in a backward dimensional loop, always 3D output
          offset=offset+length_coords+size_int
          write(qunit,'(a)')'</Points>'
        end if
        write(qunit,'(a)')'<Cells>'
        ! connectivity part
        write(offset_char,'(i19)') offset
        write(qunit,'(a,a,a)')&
          '<DataArray type="Int32" Name="connectivity" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
        offset=offset+length_conn+size_int    
        ! offsets data array
        write(offset_char,'(i19)') offset
        write(qunit,'(a,a,a)') &
          '<DataArray type="Int32" Name="offsets" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
        offset=offset+length_offsets+size_int    
        ! VTK cell type data array
        write(offset_char,'(i19)') offset
        write(qunit,'(a,a,a)') &
          '<DataArray type="Int32" Name="types" format="appended" offset="',trim(adjustl(offset_char)),'"/>' 
        offset=offset+size_int+nc*size_int
        write(qunit,'(a)')'</Cells>'
        write(qunit,'(a)')'</Piece>'
      end do
      ! write metadata communicated from other processors
      if(npe>1)then
        do ipe=1, npe-1
          do morton_no=morton_start(ipe),morton_stop(ipe)
            if(.not. morton_aim(morton_no)) cycle
            if(cell_corner) then
              ! we write out every grid as one VTK PIECE
              write(qunit,'(a,i7,a,i7,a)') &
                 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
              write(qunit,'(a)')'<PointData>'
              do iw=1,nw+nwauxio
                if(iw<=nw) then 
                  if(.not.w_write(iw)) cycle
                end if
                write(offset_char,'(i19)') offset
                write(qunit,'(a,a,a,a,a)')&
                    '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
                    '" format="appended" offset="',trim(adjustl(offset_char)),'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+length+size_int
              end do
              write(qunit,'(a)')'</PointData>'
              write(qunit,'(a)')'<Points>'
              write(offset_char,'(i19)') offset
              write(qunit,'(a,a,a)') &
              '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
              ! write cell corner coordinates in a backward dimensional loop, always 3D output
              offset=offset+length_coords+size_int
              write(qunit,'(a)')'</Points>'
            else
              ! we write out every grid as one VTK PIECE
              write(qunit,'(a,i7,a,i7,a)') &
                 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
              write(qunit,'(a)')'<CellData>'
              do iw=1,nw+nwauxio
                if(iw<=nw) then 
                  if(.not.w_write(iw)) cycle
                end if
                write(offset_char,'(i19)') offset
                write(qunit,'(a,a,a,a,a)')&
                    '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
                    '" format="appended" offset="',trim(adjustl(offset_char)),'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+lengthcc+size_int
              end do
              write(qunit,'(a)')'</CellData>'
              write(qunit,'(a)')'<Points>'
              write(offset_char,'(i19)') offset
              write(qunit,'(a,a,a)') &
              '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
              ! write cell corner coordinates in a backward dimensional loop, always 3D output
              offset=offset+length_coords+size_int
              write(qunit,'(a)')'</Points>'
            end if
            write(qunit,'(a)')'<Cells>'
            ! connectivity part
            write(offset_char,'(i19)') offset
            write(qunit,'(a,a,a)')&
              '<DataArray type="Int32" Name="connectivity" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
            offset=offset+length_conn+size_int    
            ! offsets data array
            write(offset_char,'(i19)') offset
            write(qunit,'(a,a,a)') &
              '<DataArray type="Int32" Name="offsets" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
            offset=offset+length_offsets+size_int    
            ! VTK cell type data array
            write(offset_char,'(i19)') offset
            write(qunit,'(a,a,a)') &
              '<DataArray type="Int32" Name="types" format="appended" offset="',trim(adjustl(offset_char)),'"/>' 
            offset=offset+size_int+nc*size_int
            write(qunit,'(a)')'</Cells>'
            write(qunit,'(a)')'</Piece>'
          end do
        end do
      end if
 
      write(qunit,'(a)')'</UnstructuredGrid>'
      write(qunit,'(a)')'<AppendedData encoding="raw">'
      close(qunit)
      open(qunit,file=filename,access='stream',form='unformatted',position='append')
      buf='_'
      write(qunit) trim(buf)
 
      do morton_no=morton_start(0),morton_stop(0)
        if(.not. morton_aim(morton_no)) cycle
        igrid=sfc_to_igrid(morton_no)
        call calc_x(igrid,xc,xcc)
        call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                       ixc^l,ixcc^l,.true.)
        do iw=1,nw+nwauxio
          if(iw<=nw) then 
            if(.not.w_write(iw)) cycle
          end if
          if(cell_corner) then
            write(qunit) length
            write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
          else
            write(qunit) lengthcc
            write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
          end if
        end do
 
        write(qunit) length_coords
        {do ix^db=ixcmin^db,ixcmax^db \}
          x_vtk(1:3)=zero;
          x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
          do k=1,3
            write(qunit) real(x_vtk(k))
          end do
        {end do \}
 
        write(qunit) length_conn
        call write_connvtk_binary(qunit)
 
        write(qunit) length_offsets
        do icel=1,nc
          write(qunit) icel*(2**^nd)
        end do
 
        vtk_type=vtk_cell_type()
        write(qunit) size_int*nc
        do icel=1,nc
          write(qunit) vtk_type
        end do
      end do
      allocate(intstatus(mpi_status_size,1))
      if(npe>1)then
        ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d;
        ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
        do ipe=1, npe-1
          do morton_no=morton_start(ipe),morton_stop(ipe)
            if(.not. morton_aim(morton_no)) cycle
            itag=morton_no
            call mpi_recv(xc_tmp,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            if(cell_corner) then
              call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            else
              call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            end if
            do iw=1,nw+nwauxio
              if(iw<=nw) then 
                if(.not.w_write(iw)) cycle
              end if
              if(cell_corner) then
                write(qunit) length
                write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
              else
                write(qunit) lengthcc
                write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
              end if
            end do
            write(qunit) length_coords
            {do ix^db=ixcmin^db,ixcmax^db \}
              x_vtk(1:3)=zero;
              x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
              do k=1,3
               write(qunit) real(x_vtk(k))
              end do
            {end do \}
            write(qunit) length_conn
            call write_connvtk_binary(qunit)
            write(qunit) length_offsets
            do icel=1,nc
              write(qunit) icel*(2**^nd)
            end do
            vtk_type=vtk_cell_type()
            write(qunit) size_int*nc
            do icel=1,nc
              write(qunit) vtk_type
            end do
          end do
        end do
      end if
      close(qunit)
      open(qunit,file=filename,status='unknown',form='formatted',position='append')
      write(qunit,'(a)')'</AppendedData>'
      write(qunit,'(a)')'</VTKFile>'
      close(qunit)
      deallocate(intstatus)
    end if
 
    deallocate(morton_aim,morton_aim_p)
    if (npe>1) then
      call mpi_barrier(icomm,ierrmpi)
    end if
 
  end subroutine unstructuredvtkb
 
  subroutine unstructuredvtkb64(qunit)
    ! output for vtu format to paraview, binary version output
    ! not parallel, uses calc_grid to compute nwauxio variables
    ! allows renormalizing using convert factors
    use mod_forest, only: morton_start, morton_stop, sfc_to_igrid
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) ::    qunit
 
    double precision ::  x_VTK(1:3)
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC_TMP
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio):: wC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)  :: wCC_TMP
    double precision :: normconv(0:nw+nwauxio)
    integer, allocatable :: intstatus(:,:)
    integer*8 :: offset
    integer :: itag,ipe,igrid,level,icel,ixC^L,ixCC^L,Morton_no,Morton_length
    integer :: nx^D,nxC^D,nc,np,VTK_type,ix^D,filenr
    integer::  k,iw
    integer::  length,lengthcc,length_coords,length_conn,length_offsets
    character::  buf
    character(len=80)::  filename
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
    logical ::   fileopen,cell_corner=.false.
    logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
 
    normconv=one
    morton_length=morton_stop(npe-1)-morton_start(0)+1
    allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
    allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
    morton_aim=.false.
    morton_aim_p=.false.
    do morton_no=morton_start(mype),morton_stop(mype)
      igrid=sfc_to_igrid(morton_no)
      level=node(plevel_,igrid)
      ! we can clip parts of the grid away, select variables, levels etc.
      if(writelevel(level)) then
        ! only output a grid when fully within clipped region selected
        ! by writespshift array
        if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
              *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
             <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
          morton_aim_p(morton_no)=.true.
        end if
      end if
    end do
    call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
                             icomm,ierrmpi)
    select case(convert_type)
     case('vtuB64','vtuBmpi64')
       cell_corner=.true.
     case('vtuBCC64','vtuBCCmpi64')
       cell_corner=.false.
    end select
    if (mype /= 0) then
      do morton_no=morton_start(mype),morton_stop(mype)
        if(.not. morton_aim(morton_no)) cycle
        igrid=sfc_to_igrid(morton_no)
        call calc_x(igrid,xc,xcc)
        call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                          ixc^l,ixcc^l,.true.)
        itag=morton_no
        call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
        if(cell_corner) then
          call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
        else
          call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
        end if
      end do
    else
      ! mype==0
      offset=0
      inquire(qunit,opened=fileopen)
      if(.not.fileopen)then
        ! generate filename 
        filenr=snapshotini
        if (autoconvert) filenr=snapshotnext
        write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
        ! Open the file for the header part
        open(qunit,file=filename,status='replace')
      end if
      call getheadernames(wnamei,xandwnamei,outfilehead)
      ! generate xml header
      write(qunit,'(a)')'<?xml version="1.0"?>'
      write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
      write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
      write(qunit,'(a)')'<UnstructuredGrid>'
      write(qunit,'(a)')'<FieldData>'
      write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                         'NumberOfTuples="1" format="ascii">'
      write(qunit,*) real(global_time*time_convert_factor)
      write(qunit,'(a)')'</DataArray>'
      write(qunit,'(a)')'</FieldData>'
      ! number of cells, number of corner points, per grid.
      nx^d=ixmhi^d-ixmlo^d+1;
      nxc^d=nx^d+1;
      nc={nx^d*}
      np={nxc^d*}
      length=np*size_double
      lengthcc=nc*size_double
      length_coords=3*length
      length_conn=2**^nd*size_int*nc
      length_offsets=nc*size_int
      ! Note: using the w_write, writelevel, writespshift
      do morton_no=morton_start(0),morton_stop(0)
        if(.not. morton_aim(morton_no)) cycle
        if(cell_corner) then
          ! we write out every grid as one VTK PIECE
          write(qunit,'(a,i7,a,i7,a)') &
             '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
          write(qunit,'(a)')'<PointData>'
          do iw=1,nw+nwauxio
            if(iw<=nw) then 
              if(.not.w_write(iw)) cycle
            end if
            write(qunit,'(a,a,a,i16,a)')&
                '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
                '" format="appended" offset="',offset,'">'
            write(qunit,'(a)')'</DataArray>'
            offset=offset+length+size_int
          end do
          write(qunit,'(a)')'</PointData>'
          write(qunit,'(a)')'<Points>'
          write(qunit,'(a,i16,a)') &
          '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
          ! write cell corner coordinates in a backward dimensional loop, always 3D output
          offset=offset+length_coords+size_int
          write(qunit,'(a)')'</Points>'
        else
          ! we write out every grid as one VTK PIECE
          write(qunit,'(a,i7,a,i7,a)') &
             '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
          write(qunit,'(a)')'<CellData>'
          do iw=1,nw+nwauxio
            if(iw<=nw) then 
               if(.not.w_write(iw)) cycle
            end if
            write(qunit,'(a,a,a,i16,a)')&
                '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
                '" format="appended" offset="',offset,'">'
            write(qunit,'(a)')'</DataArray>'
            offset=offset+lengthcc+size_int
          end do
          write(qunit,'(a)')'</CellData>'
          write(qunit,'(a)')'<Points>'
          write(qunit,'(a,i16,a)') &
          '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
          ! write cell corner coordinates in a backward dimensional loop, always 3D output
          offset=offset+length_coords+size_int
          write(qunit,'(a)')'</Points>'
        end if
        write(qunit,'(a)')'<Cells>'
        ! connectivity part
        write(qunit,'(a,i16,a)')&
          '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
        offset=offset+length_conn+size_int    
        ! offsets data array
        write(qunit,'(a,i16,a)') &
          '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
        offset=offset+length_offsets+size_int    
        ! VTK cell type data array
        write(qunit,'(a,i16,a)') &
          '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>' 
        offset=offset+size_int+nc*size_int
        write(qunit,'(a)')'</Cells>'
        write(qunit,'(a)')'</Piece>'
      end do
      ! write metadata communicated from other processors
      if(npe>1)then
        do ipe=1, npe-1
          do morton_no=morton_start(ipe),morton_stop(ipe)
            if(.not. morton_aim(morton_no)) cycle
            if(cell_corner) then
              ! we write out every grid as one VTK PIECE
              write(qunit,'(a,i7,a,i7,a)') &
                 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
              write(qunit,'(a)')'<PointData>'
              do iw=1,nw+nwauxio
                if(iw<=nw) then 
                  if(.not.w_write(iw)) cycle
                end if
                write(qunit,'(a,a,a,i16,a)')&
                    '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
                    '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+length+size_int
              end do
              write(qunit,'(a)')'</PointData>'
              write(qunit,'(a)')'<Points>'
              write(qunit,'(a,i16,a)') &
              '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
              ! write cell corner coordinates in a backward dimensional loop, always 3D output
              offset=offset+length_coords+size_int
              write(qunit,'(a)')'</Points>'
            else
              ! we write out every grid as one VTK PIECE
              write(qunit,'(a,i7,a,i7,a)') &
                 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
              write(qunit,'(a)')'<CellData>'
              do iw=1,nw+nwauxio
                if(iw<=nw) then 
                  if(.not.w_write(iw)) cycle
                end if
                write(qunit,'(a,a,a,i16,a)')&
                    '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
                    '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+lengthcc+size_int
              end do
              write(qunit,'(a)')'</CellData>'
              write(qunit,'(a)')'<Points>'
              write(qunit,'(a,i16,a)') &
              '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
              ! write cell corner coordinates in a backward dimensional loop, always 3D output
              offset=offset+length_coords+size_int
              write(qunit,'(a)')'</Points>'
            end if
            write(qunit,'(a)')'<Cells>'
            ! connectivity part
            write(qunit,'(a,i16,a)')&
              '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
            offset=offset+length_conn+size_int    
            ! offsets data array
            write(qunit,'(a,i16,a)') &
              '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
            offset=offset+length_offsets+size_int    
            ! VTK cell type data array
            write(qunit,'(a,i16,a)') &
              '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>' 
            offset=offset+size_int+nc*size_int
            write(qunit,'(a)')'</Cells>'
            write(qunit,'(a)')'</Piece>'
          end do
        end do
      end if
      write(qunit,'(a)')'</UnstructuredGrid>'
      write(qunit,'(a)')'<AppendedData encoding="raw">'
      close(qunit)
      open(qunit,file=filename,access='stream',form='unformatted',position='append')
      buf='_'
      write(qunit) trim(buf)
      do morton_no=morton_start(0),morton_stop(0)
        if(.not. morton_aim(morton_no)) cycle
        igrid=sfc_to_igrid(morton_no)
        call calc_x(igrid,xc,xcc)
        call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                       ixc^l,ixcc^l,.true.)
        do iw=1,nw+nwauxio
          if(iw<=nw) then 
            if(.not.w_write(iw)) cycle
          end if
          if(cell_corner) then
            write(qunit) length
            write(qunit) {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
          else
            write(qunit) lengthcc
            write(qunit) {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
          end if
        end do
        write(qunit) length_coords
        {do ix^db=ixcmin^db,ixcmax^db \}
          x_vtk(1:3)=zero;
          x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
          do k=1,3
            write(qunit) x_vtk(k)
          end do
        {end do \}
        write(qunit) length_conn
        call write_connvtk_binary(qunit)
        write(qunit) length_offsets
        do icel=1,nc
          write(qunit) icel*(2**^nd)
        end do
        vtk_type=vtk_cell_type()
        write(qunit) size_int*nc
        do icel=1,nc
          write(qunit) vtk_type
        end do
      end do
      allocate(intstatus(mpi_status_size,1))
      if(npe>1)then
        ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d;
        ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
        do ipe=1, npe-1
          do morton_no=morton_start(ipe),morton_stop(ipe)
            if(.not. morton_aim(morton_no)) cycle
            itag=morton_no
            call mpi_recv(xc_tmp,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            if(cell_corner) then
              call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            else
              call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            end if
            do iw=1,nw+nwauxio
              if(iw<=nw) then 
                if(.not.w_write(iw)) cycle
              end if
              if(cell_corner) then
                write(qunit) length
                write(qunit) {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
              else
                write(qunit) lengthcc
                write(qunit) {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
              end if
            end do
            write(qunit) length_coords
            {do ix^db=ixcmin^db,ixcmax^db \}
              x_vtk(1:3)=zero;
              x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
              do k=1,3
               write(qunit) x_vtk(k)
              end do
            {end do \}
            write(qunit) length_conn
            call write_connvtk_binary(qunit)
            write(qunit) length_offsets
            do icel=1,nc
              write(qunit) icel*(2**^nd)
            end do
            vtk_type=vtk_cell_type()
            write(qunit) size_int*nc
            do icel=1,nc
              write(qunit) vtk_type
            end do
          end do
        end do
      end if
      close(qunit)
      open(qunit,file=filename,status='unknown',form='formatted',position='append')
      write(qunit,'(a)')'</AppendedData>'
      write(qunit,'(a)')'</VTKFile>'
      close(qunit)
      deallocate(intstatus)
    end if
    deallocate(morton_aim,morton_aim_p)
    if (npe>1) then
      call mpi_barrier(icomm,ierrmpi)
    end if
 
  end subroutine unstructuredvtkb64
 
  subroutine save_connvtk(qunit,igrid)
    ! this saves the basic line, pixel and voxel connectivity,
    ! as used by VTK file outputs for unstructured grid
    use mod_global_parameters
    use mod_geometry, only: cartesian, cartesian_expansion, cartesian_stretched, coordinate
 
    integer, intent(in) :: qunit, igrid
 
    integer :: nx^D, nxC^D, ix^D
    logical :: transformed_coordinates
 
    nx^d=ixmhi^d-ixmlo^d+1;
    nxc^d=nx^d+1;
    transformed_coordinates=vtk_coordinates_transformed()
    {do ix^db=1,nx^db\}
            {^ifoned write(qunit,'(2(i7,1x))')ix1-1,ix1 \}
            {^iftwod
            if(transformed_coordinates) then
              write(qunit,'(4(i7,1x))')(ix2-1)*nxc1+ix1-1, &
                   (ix2-1)*nxc1+ix1,ix2*nxc1+ix1,ix2*nxc1+ix1-1
            else
              write(qunit,'(4(i7,1x))')(ix2-1)*nxc1+ix1-1, &
                   (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
            end if
            \}
            {^ifthreed
            if(transformed_coordinates) then
              write(qunit,'(8(i7,1x))')&
                       (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
                       (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
                       (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1,&
                       (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1-1,&
                           ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
                           ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
                           ix3*nxc2*nxc1+    ix2*nxc1+ix1,&
                           ix3*nxc2*nxc1+    ix2*nxc1+ix1-1
            else
              write(qunit,'(8(i7,1x))')&
                       (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
                       (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
                       (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1-1,&
                       (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1,&
                           ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
                           ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
                           ix3*nxc2*nxc1+    ix2*nxc1+ix1-1,&
                           ix3*nxc2*nxc1+    ix2*nxc1+ix1
                          end if
            \}
    {end do\}
 
  end subroutine save_connvtk
 
  integer function vtk_cell_type()
    logical :: transformed_coordinates
 
    transformed_coordinates=vtk_coordinates_transformed()
    {^ifoned vtk_cell_type=3 \} !VTK_LINE
    {^iftwod
    if(transformed_coordinates) then
      vtk_cell_type=9 !VTK_QUAD
    else
      vtk_cell_type=8 !VTK_PIXEL
    end if
    \}
    {^ifthreed
    if(transformed_coordinates) then
      vtk_cell_type=12 !VTK_HEXAHEDRON
    else
      vtk_cell_type=11 !VTK_VOXEL
    end if
    \}
  end function vtk_cell_type
 
  logical function vtk_coordinates_transformed()
    use mod_global_parameters, only: ndim, nocartesian, phi_
    use mod_geometry, only: cartesian, cartesian_expansion, cartesian_stretched, cylindrical, &
                            coordinate, spherical
 
    vtk_coordinates_transformed=.false.
    if (nocartesian) return
 
    select case(coordinate)
    case(cylindrical)
      vtk_coordinates_transformed=(ndim==3) .or. (ndim==2 .and. phi_==2)
    case(spherical)
      vtk_coordinates_transformed=ndim>1
    case(cartesian, cartesian_stretched, cartesian_expansion)
      vtk_coordinates_transformed=.false.
    end select
  end function vtk_coordinates_transformed
 
  subroutine write_connvtk_binary(qunit)
    use mod_global_parameters
    use mod_geometry, only: cartesian, cartesian_expansion, cartesian_stretched, coordinate
 
    integer, intent(in) :: qunit
    integer :: nx^D, nxC^D, ix^D
    logical :: transformed_coordinates
 
    nx^d=ixmhi^d-ixmlo^d+1;
    nxc^d=nx^d+1;
    transformed_coordinates=vtk_coordinates_transformed()
    {do ix^db=1,nx^db\}
      {^ifoned write(qunit)ix1-1,ix1 \}
      {^iftwod
      if(transformed_coordinates) then
        write(qunit)(ix2-1)*nxc1+ix1-1, &
                    (ix2-1)*nxc1+ix1,ix2*nxc1+ix1,ix2*nxc1+ix1-1
      else
        write(qunit)(ix2-1)*nxc1+ix1-1, &
                    (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
      end if
      \}
      {^ifthreed
      if(transformed_coordinates) then
        write(qunit)&
        (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
        (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
        (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1,&
        (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1-1,&
         ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
         ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
         ix3*nxc2*nxc1+    ix2*nxc1+ix1,&
         ix3*nxc2*nxc1+    ix2*nxc1+ix1-1
      else
        write(qunit)&
        (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
        (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
        (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1-1,&
        (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1,&
         ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
         ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
         ix3*nxc2*nxc1+    ix2*nxc1+ix1-1,&
         ix3*nxc2*nxc1+    ix2*nxc1+ix1
      end if
      \}
    {end do\}
  end subroutine write_connvtk_binary
 
  subroutine imagedatavtk_mpi(qunit)
    ! output for vti format to paraview, non-binary version output
    ! parallel, uses calc_grid to compute nwauxio variables
    ! allows renormalizing using convert factors
    ! allows skipping of w_write selected variables
    ! implementation such that length of ASCII output is identical when 
    ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
    use mod_forest, only: morton_start, morton_stop, tree_node_ptr, igrid_to_node, sfc_to_igrid
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) ::    qunit
 
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC_TMP,xCC_TMP_recv
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC_TMP,wC_TMP_recv
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC_TMP,wCC_TMP_recv
    double precision, dimension(0:nw+nwauxio)                   :: normconv
    double precision :: origin(1:3), spacing(1:3)
    integer :: igrid,iigrid,level,ixC^L,ixCC^L
    integer :: NumGridsOnLevel(1:nlevelshi)
    integer :: nx^D
    integer :: filenr
    integer :: itag,ipe,Morton_no,Morton_length
    integer :: ixrvC^L, ixrvCC^L, siz_ind, ind_send(5*^ND), ind_recv(5*^ND)
    integer :: wholeExtent(1:6), ig^D
    integer, allocatable :: intstatus(:,:)
    logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
    logical :: fileopen
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
    character(len=80)::  filename
    type(tree_node_ptr) :: tree
 
    if(levmin/=levmax) call mpistop('ImageData can only be used when levmin=levmax')
    normconv(0) = length_convert_factor
    normconv(1:nw) = w_convert_factor
    siz_ind=5*^nd
    morton_length=morton_stop(npe-1)-morton_start(0)+1
    allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
    allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
    morton_aim=.false.
    morton_aim_p=.false.
    do morton_no=morton_start(mype),morton_stop(mype)
      igrid=sfc_to_igrid(morton_no)
      level=node(plevel_,igrid)
      ! we can clip parts of the grid away, select variables, levels etc.
      if(writelevel(level)) then
       ! only output a grid when fully within clipped region selected
       ! by writespshift array
       if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
             *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
            <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
         morton_aim_p(morton_no)=.true.
       end if
      end if
    end do
    call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
                             icomm,ierrmpi)
    if(mype /= 0) then
      do morton_no=morton_start(mype),morton_stop(mype)
        if(.not. morton_aim(morton_no)) cycle
        igrid=sfc_to_igrid(morton_no)
        call calc_x(igrid,xc,xcc)
        call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                          ixc^l,ixcc^l,.true.)
        tree%node => igrid_to_node(igrid, mype)%node
        {^d& ig^d = tree%node%ig^d; }
        itag=morton_no
        ind_send=(/ ixc^l,ixcc^l, ig^d /)
        call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
        call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
        call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
      end do
    else
      inquire(qunit,opened=fileopen)
      if(.not.fileopen)then
        ! generate filename 
        filenr=snapshotini
        if (autoconvert) filenr=snapshotnext
        write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vti"
        ! Open the file for the header part
        open(qunit,file=filename,status='unknown',form='formatted')
      end if
      call getheadernames(wnamei,xandwnamei,outfilehead)
      ! number of cells per grid.
      nx^d=ixmhi^d-ixmlo^d+1;
      origin      = 0
      {^d& origin(^d) = xprobmin^d*normconv(0); }
      spacing     = zero
      {^d&spacing(^d) = dxlevel(^d)*normconv(0); }
      wholeextent = 0
      ! if we use writespshift, the whole extent has to be calculated:
      {^d&wholeextent(^d*2-1) = nx^d * ceiling(((xprobmax^d-xprobmin^d)*writespshift(^d,1)) &
           /(nx^d*dxlevel(^d))) \}
      {^d&wholeextent(^d*2)   = nx^d * floor(((xprobmax^d-xprobmin^d)*(1.0d0-writespshift(^d,2))) &
           /(nx^d*dxlevel(^d))) \}
      
      ! generate xml header
      write(qunit,'(a)')'<?xml version="1.0"?>'
      write(qunit,'(a)',advance='no') '<VTKFile type="ImageData"'
      write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
      write(qunit,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')'  <ImageData Origin="',&
           origin,'" WholeExtent="',wholeextent,'" Spacing="',spacing,'">'
      write(qunit,'(a)')'<FieldData>'
      write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                         'NumberOfTuples="1" format="ascii">'
      write(qunit,*) real(global_time*time_convert_factor)
      write(qunit,'(a)')'</DataArray>'
      write(qunit,'(a)')'</FieldData>'
 
      ! write the data from proc 0
      do morton_no=morton_start(0),morton_stop(0)
        if(.not. morton_aim(morton_no)) cycle
        igrid=sfc_to_igrid(morton_no)
        tree%node => igrid_to_node(igrid, 0)%node
        {^d& ig^d = tree%node%ig^d; }
        call calc_x(igrid,xc,xcc)
        call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
             ixc^l,ixcc^l,.true.)
        call write_vti(qunit,ixg^ll,ixc^l,ixcc^l,ig^d,&
             nx^d,normconv,wnamei,wc_tmp,wcc_tmp)   
      end do
 
      if(npe>1)then
        allocate(intstatus(mpi_status_size,1))
        do ipe=1, npe-1
          do morton_no=morton_start(ipe),morton_stop(ipe)
            if(.not. morton_aim(morton_no)) cycle
            itag=morton_no
            call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
            ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
            ig^d=ind_recv(4*^nd+^d);
            call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            call write_vti(qunit,ixg^ll,ixrvc^l,ixrvcc^l,ig^d,&
                 nx^d,normconv,wnamei,wc_tmp,wcc_tmp)   
          end do
        end do
      end if
      write(qunit,'(a)')'</ImageData>'
      write(qunit,'(a)')'</VTKFile>'
      close(qunit)
      if(npe>1) deallocate(intstatus)
    end if
 
    deallocate(morton_aim,morton_aim_p)
    if (npe>1) then
      call mpi_barrier(icomm,ierrmpi)
    endif
 
  end subroutine imagedatavtk_mpi
 
  subroutine punstructuredvtk_mpi(qunit)
    ! Write one pvtu and vtu files for each processor
    ! Otherwise like unstructuredvtk_mpi
    use mod_forest, only: morton_start, morton_stop, sfc_to_igrid
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) ::    qunit
 
    double precision, dimension(0:nw+nwauxio)                   :: normconv
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim)         :: xC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)           :: xCC_TMP
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim)         :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)           :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC_TMP
    integer             :: nx^D,nxC^D,nc,np, igrid,ixC^L,ixCC^L,level,Morton_no
    integer             :: filenr
    logical             :: fileopen,conv_grid
    character(len=name_len)   :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
    character(len=80)   :: pfilename
 
    ! Write pvtu-file:
    if (mype==0) then
       call write_pvtu(qunit)
    endif
    ! Now write the Source files:
    inquire(qunit,opened=fileopen)
    if(.not.fileopen)then
      ! generate filename 
      filenr=snapshotini
      if (autoconvert) filenr=snapshotnext
      ! Open the file for the header part
      write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
      open(qunit,file=pfilename,status='unknown',form='formatted')
    end if
    ! generate xml header
    write(qunit,'(a)')'<?xml version="1.0"?>'
    write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
    write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
    write(qunit,'(a)')'  <UnstructuredGrid>'
    write(qunit,'(a)')'<FieldData>'
    write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                       'NumberOfTuples="1" format="ascii">'
    write(qunit,*) real(global_time*time_convert_factor)
    write(qunit,'(a)')'</DataArray>'
    write(qunit,'(a)')'</FieldData>'
 
    call getheadernames(wnamei,xandwnamei,outfilehead)
 
    ! number of cells, number of corner points, per grid.
    nx^d=ixmhi^d-ixmlo^d+1;
    nxc^d=nx^d+1;
    nc={nx^d*}
    np={nxc^d*}
 
    ! Note: using the w_write, writelevel, writespshift
    ! we can clip parts of the grid away, select variables, levels etc.
    do level=levmin,levmax
      if (.not.writelevel(level)) cycle
      do morton_no=morton_start(mype),morton_stop(mype)
        igrid=sfc_to_igrid(morton_no)
        if (node(plevel_,igrid)/=level) cycle
        ! only output a grid when fully within clipped region selected
        ! by writespshift array
        conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
              *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
             <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
        if (.not.conv_grid) cycle
        call calc_x(igrid,xc,xcc)
        call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
             ixc^l,ixcc^l,.true.)
        call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
             normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
      end do ! Morton_no loop
    end do ! level loop
 
    write(qunit,'(a)')'  </UnstructuredGrid>'
    write(qunit,'(a)')'</VTKFile>'
    close(qunit)
 
    if (npe>1) then
      call mpi_barrier(icomm,ierrmpi)
    end if
 
  end subroutine punstructuredvtk_mpi
 
  subroutine unstructuredvtk_mpi(qunit)
    ! output for vtu format to paraview, non-binary version output
    ! parallel, uses calc_grid to compute nwauxio variables
    ! allows renormalizing using convert factors
    ! allows skipping of w_write selected variables
    ! implementation such that length of ASCII output is identical when 
    ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
    use mod_forest, only: morton_start, morton_stop, sfc_to_igrid
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) ::    qunit
 
    double precision ::  x_VTK(1:3)
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC_TMP,xCC_TMP_recv
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC_TMP,wC_TMP_recv
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC_TMP,wCC_TMP_recv
    double precision, dimension(0:nw+nwauxio)                   :: normconv
    integer::               igrid,iigrid,level,ixC^L,ixCC^L
    integer::               NumGridsOnLevel(1:nlevelshi)
    integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,ix^D
    integer ::           filenr
    integer :: itag,ipe,Morton_no,siz_ind
    integer :: ind_send(4*^ND),ind_recv(4*^ND)
    integer :: levmin_recv,levmax_recv,level_recv,igrid_recv,ixrvC^L,ixrvCC^L
    integer, allocatable :: intstatus(:,:)
    logical :: fileopen,conv_grid,cond_grid_recv
    character(len=80)::  filename
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
 
    if(mype==0) then
      inquire(qunit,opened=fileopen)
      if(.not.fileopen)then
        ! generate filename 
        filenr=snapshotini
        if (autoconvert) filenr=snapshotnext
        write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
        ! Open the file for the header part
        open(qunit,file=filename,status='unknown',form='formatted')
      end if
      ! generate xml header
      write(qunit,'(a)')'<?xml version="1.0"?>'
      write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
      write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
      write(qunit,'(a)')'<UnstructuredGrid>'
      write(qunit,'(a)')'<FieldData>'
      write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                         'NumberOfTuples="1" format="ascii">'
      write(qunit,*) real(global_time*time_convert_factor)
      write(qunit,'(a)')'</DataArray>'
      write(qunit,'(a)')'</FieldData>'
    end if
 
    call getheadernames(wnamei,xandwnamei,outfilehead)
    ! number of cells, number of corner points, per grid.
    nx^d=ixmhi^d-ixmlo^d+1;
    nxc^d=nx^d+1;
    nc={nx^d*}
    np={nxc^d*}
    ! all slave processors send their minmal/maximal levels
    if  (mype/=0) then
     if (morton_stop(mype)==0) call mpistop("nultag")
     itag=1000*morton_stop(mype)
     !print *,'ype,itag for levmin=',mype,itag,levmin
     call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
     itag=2000*morton_stop(mype)
     !print *,'mype,itag for levmax=',mype,itag,levmax
     call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
    end if
    ! Note: using the w_write, writelevel, writespshift
    ! we can clip parts of the grid away, select variables, levels etc.
    do level=levmin,levmax
      if (.not.writelevel(level)) cycle
      do morton_no=morton_start(mype),morton_stop(mype)
        igrid=sfc_to_igrid(morton_no)
        if (mype/=0)then
          itag=morton_no
          call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
          itag=igrid
          call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
        end if
        if (node(plevel_,igrid)/=level) cycle
        ! only output a grid when fully within clipped region selected
        ! by writespshift array
        conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
              *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
             <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
        if (mype/=0)then
          call mpi_send(conv_grid,1,mpi_logical,0,itag,icomm,ierrmpi)
        end if
        if (.not.conv_grid) cycle
        call calc_x(igrid,xc,xcc)
        call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                         ixc^l,ixcc^l,.true.)
        if(mype/=0) then
          itag=morton_no
          ind_send=(/ ixc^l,ixcc^l /)
          siz_ind=4*^nd
          call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
          call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
          call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
          call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
          itag=igrid
          call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
          call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
        else
          call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
                             normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
        end if
      end do ! Morton_no loop
    end do ! level loop
 
    if(mype==0) then
      allocate(intstatus(mpi_status_size,1))
      if(npe>1)then
        do ipe=1,npe-1
          itag=1000*morton_stop(ipe)
          call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
          !!print *,'mype RECEIVES,itag for levmin=',mype,itag,levmin_recv
          itag=2000*morton_stop(ipe)
          call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
          !!print *,'mype RECEIVES itag for levmax=',mype,itag,levmax_recv
          do level=levmin_recv,levmax_recv
            if (.not.writelevel(level)) cycle
            do  morton_no=morton_start(ipe),morton_stop(ipe)
              itag=morton_no
              call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              itag=igrid_recv
              call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              if (level_recv/=level) cycle
              call mpi_recv(cond_grid_recv,1,mpi_logical, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              if(.not.cond_grid_recv)cycle
              itag=morton_no
              siz_ind=4*^nd
              call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
              ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
              call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
                           ,icomm,intstatus(:,1),ierrmpi)
              call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              itag=igrid_recv
              call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
              call write_vtk(qunit,ixg^ll,ixrvc^l,ixrvcc^l,igrid_recv,&
                             nc,np,nx^d,nxc^d,normconv,wnamei,&
                             xc_tmp_recv,xcc_tmp_recv,wc_tmp_recv,wcc_tmp_recv)
            end do ! Morton_no loop
          end do ! level loop
        end do ! processor loop
      end if ! multiple processors
      write(qunit,'(a)')'</UnstructuredGrid>'
      write(qunit,'(a)')'</VTKFile>'
      close(qunit)
    end if
    if (npe>1) then
      call mpi_barrier(icomm,ierrmpi)
      if(mype==0)deallocate(intstatus)
    end if
 
  end subroutine unstructuredvtk_mpi
 
  subroutine write_vtk(qunit,ixI^L,ixC^L,ixCC^L,igrid,nc,np,nx^D,nxC^D,&
                       normconv,wnamei,xC,xCC,wC,wCC)
    use mod_global_parameters
 
    integer, intent(in) :: qunit
    integer, intent(in) :: ixI^L,ixC^L,ixCC^L
    integer, intent(in) :: igrid,nc,np,nx^D,nxC^D
    double precision, intent(in) :: normconv(0:nw+nwauxio) 
    character(len=name_len), intent(in)::  wnamei(1:nw+nwauxio)
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC
 
    double precision ::  x_VTK(1:3)
    integer :: iw,ix^D,icel,VTK_type
 
    select case(convert_type)
      case('vtumpi','pvtumpi')
           ! we write out every grid as one VTK PIECE
        write(qunit,'(a,i7,a,i7,a)') &
              '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
        write(qunit,'(a)')'<PointData>'
        do iw=1,nw+nwauxio
          if(iw<=nw) then 
            if(.not.w_write(iw)) cycle
          end if
          write(qunit,'(a,a,a)')&
           '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
          write(qunit,'(200(1pe14.6))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
          write(qunit,'(a)')'</DataArray>'
        end do
        write(qunit,'(a)')'</PointData>'
        write(qunit,'(a)')'<Points>'
        write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
           ! write cell corner coordinates in a backward dimensional loop, always 3D output
        {do ix^db=ixcmin^db,ixcmax^db \}
              x_vtk(1:3)=zero;
              x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
              write(qunit,'(3(1pe14.6))') x_vtk
        {end do \}
        write(qunit,'(a)')'</DataArray>'
        write(qunit,'(a)')'</Points>'
 
      case('vtuCCmpi','pvtuCCmpi')
        ! we write out every grid as one VTK PIECE
        write(qunit,'(a,i7,a,i7,a)') &
           '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
        write(qunit,'(a)')'<CellData>'
        do iw=1,nw+nwauxio
          if(iw<=nw) then 
            if(.not.w_write(iw)) cycle
          end if
          write(qunit,'(a,a,a)')&
           '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
          write(qunit,'(200(1pe14.6))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
          write(qunit,'(a)')'</DataArray>'
        end do
        write(qunit,'(a)')'</CellData>'
        write(qunit,'(a)')'<Points>'
        write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
        ! write cell corner coordinates in a backward dimensional loop, always 3D output
        {do ix^db=ixcmin^db,ixcmax^db \}
           x_vtk(1:3)=zero;
           x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
           write(qunit,'(3(1pe14.6))') x_vtk
        {end do \}
        write(qunit,'(a)')'</DataArray>'
        write(qunit,'(a)')'</Points>'
    end select
 
    write(qunit,'(a)')'<Cells>'
    ! connectivity part
    write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
    call save_connvtk(qunit,igrid)
    write(qunit,'(a)')'</DataArray>'
    ! offsets data array
    write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
    do icel=1,nc
      write(qunit,'(i7)') icel*(2**^nd)
    end do
    write(qunit,'(a)')'</DataArray>'
    ! VTK cell type data array
    write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
    vtk_type=vtk_cell_type()
    do icel=1,nc
      write(qunit,'(i2)') vtk_type
    end do
    write(qunit,'(a)')'</DataArray>'
    write(qunit,'(a)')'</Cells>'
    write(qunit,'(a)')'</Piece>'
 
  end subroutine write_vtk
 
  subroutine write_vti(qunit,ixI^L,ixC^L,ixCC^L,ig^D,nx^D,&
                       normconv,wnamei,wC,wCC)
    use mod_global_parameters
 
    integer, intent(in) :: qunit
    integer, intent(in) :: ixI^L,ixC^L,ixCC^L
    integer, intent(in) :: ig^D,nx^D
    double precision, intent(in) :: normconv(0:nw+nwauxio) 
    character(len=name_len), intent(in)::  wnamei(1:nw+nwauxio)
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC
 
    integer :: iw,ix^D
    integer :: extent(1:6)
 
    extent = 0
    {^d& extent(^d*2-1) = (ig^d-1) * nx^d; }
    {^d& extent(^d*2)   = (ig^d)   * nx^d; }
 
    select case(convert_type)
      case('vtimpi','pvtimpi')
        ! we write out every grid as one VTK PIECE
        write(qunit,'(a,6(i10),a)') &
              '<Piece Extent="',extent,'">'
        write(qunit,'(a)')'<PointData>'
        do iw=1,nw+nwauxio
          if(iw<=nw) then 
             if(.not.w_write(iw)) cycle
          end if
          write(qunit,'(a,a,a)')&
           '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
          write(qunit,'(200(1pe20.12))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
          write(qunit,'(a)')'</DataArray>'
        end do
        write(qunit,'(a)')'</PointData>'
      case('vtiCCmpi','pvtiCCmpi')
        ! we write out every grid as one VTK PIECE
        write(qunit,'(a,6(i10),a)') &
              '<Piece Extent="',extent,'">'
        write(qunit,'(a)')'<CellData>'
        do iw=1,nw+nwauxio
          if(iw<=nw) then 
            if(.not.w_write(iw)) cycle
          end if
          write(qunit,'(a,a,a)')&
           '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
          write(qunit,'(200(1pe20.12))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
          write(qunit,'(a)')'</DataArray>'
        end do
        write(qunit,'(a)')'</CellData>'
    end select
 
    write(qunit,'(a)')'</Piece>'
 
  end subroutine write_vti
 
  subroutine write_pvtu(qunit)
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) :: qunit
 
    integer             :: filenr,iw,ipe,iscalars
    logical             :: fileopen
    character(len=name_len)   :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio),outtype
    character(len=1024) :: outfilehead
    character(len=80)   :: filename,pfilename
 
    select case(convert_type)
    case('pvtumpi','pvtuBmpi')
       outtype="PPointData"
    case('pvtuCCmpi','pvtuBCCmpi')
       outtype="PCellData"
    end select
    inquire(qunit,opened=fileopen)
    if(.not.fileopen)then
       ! generate filename 
       filenr=snapshotini
       if (autoconvert) filenr=snapshotnext
       write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".pvtu"
       ! Open the file
       open(qunit,file=filename,status='unknown',form='formatted')
    end if
 
    call getheadernames(wnamei,xandwnamei,outfilehead)
    ! Get the default selection:
    iscalars=1
    do iw=nw,1, -1
       if (w_write(iw)) iscalars=iw
    end do
    ! generate xml header
    write(qunit,'(a)')'<?xml version="1.0"?>'
    write(qunit,'(a)',advance='no') '<VTKFile type="PUnstructuredGrid"'
    write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
    write(qunit,'(a)')'  <PUnstructuredGrid GhostLevel="0">'
    ! Either celldata or pointdata:
    write(qunit,'(a,a,a,a,a)')&
         '    <',trim(outtype),' Scalars="',trim(wnamei(iscalars))//'">'
    do iw=1,nw
      if(.not.w_write(iw))cycle
      write(qunit,'(a,a,a)')&
           '      <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
    end do
    do iw=nw+1,nw+nwauxio
      write(qunit,'(a,a,a)')&
           '      <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
    end do
    write(qunit,'(a,a,a)')'    </',trim(outtype),'>'
    write(qunit,'(a)')'    <PPoints>'
    write(qunit,'(a)')'      <PDataArray type="Float32" NumberOfComponents="3"/>'
    write(qunit,'(a)')'    </PPoints>'
 
    do ipe=0,npe-1
      write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename(&
           index(base_filename, '/', back = .true.)+1:&
           len(base_filename))),filenr,"p",&
           ipe,".vtu"
      write(qunit,'(a,a,a)')'    <Piece Source="',trim(pfilename),'"/>'
    end do
    write(qunit,'(a)')'  </PUnstructuredGrid>'
    write(qunit,'(a)')'</VTKFile>'
    close(qunit)
 
  end subroutine write_pvtu
 
  subroutine tecplot_mpi(qunit)
    ! output for tecplot (ASCII format)
    ! parallel, uses calc_grid to compute nwauxio variables
    ! allows renormalizing using convert factors
    ! the current implementation is such that tecplotmpi and tecplotCCmpi will 
    ! create different length output ASCII files when used on 1 versus multiple CPUs
    ! in fact, on 1 CPU, there will be as many zones as there are levels
    ! on multiple CPUs, there will be a number of zones up to the number of
    ! levels times the number of CPUs (can be less, when some level not on a CPU)
    use mod_forest, only: morton_start, morton_stop, sfc_to_igrid
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) :: qunit
 
    double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC_TMP,xCC_TMP_recv
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC_TMP,wC_TMP_recv
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC_TMP,wCC_TMP_recv
    double precision, dimension(0:nw+nwauxio)                   :: normconv
    integer::               igrid,iigrid,level,igonlevel,iw,idim,ix^D
    integer::               NumGridsOnLevel(1:nlevelshi)
    integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
    integer :: nodesonlevelmype,elemsonlevelmype
    integer ::              nodes, elems
    integer, allocatable :: intstatus(:,:)
    integer :: itag,Morton_no,ipe,levmin_recv,levmax_recv,igrid_recv,level_recv
    integer :: ixrvC^L,ixrvCC^L
    integer :: ind_send(2*^ND),ind_recv(2*^ND),siz_ind,igonlevel_recv
    integer :: NumGridsOnLevel_mype(1:nlevelshi,0:npe-1)
    integer ::           filenr
    logical :: fileopen,first
    character(len=80) :: filename
    character(len=1024) :: tecplothead
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
 
    if(nw/=count(w_write(1:nw)))then
     if(mype==0) print *,'tecplot_mpi does not use w_write=F'
     call mpistop('w_write, tecplot')
    end if
 
    if(nocartesian)then
      if(mype==0) print *,'tecplot_mpi with nocartesian'
    end if
 
    master_cpu_open : if (mype == 0) then
      inquire(qunit,opened=fileopen)
      if (.not.fileopen) then
        ! generate filename
         filenr=snapshotini
         if (autoconvert) filenr=snapshotnext
        write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
        open(qunit,file=filename,status='unknown')
      end if
      call getheadernames(wnamei,xandwnamei,outfilehead)
      write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
      write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
    end if  master_cpu_open
 
    ! determine overall number of grids per level, and the same info per CPU
    numgridsonlevel(1:nlevelshi)=0
    do level=levmin,levmax
      numgridsonlevel(level)=0
      do morton_no=morton_start(mype),morton_stop(mype)
        igrid = sfc_to_igrid(morton_no)
        if (node(plevel_,igrid)/=level) cycle
        numgridsonlevel(level)=numgridsonlevel(level)+1
      end do
      numgridsonlevel_mype(level,0:npe-1)=0
      numgridsonlevel_mype(level,mype) = numgridsonlevel(level)
      call mpi_allreduce(mpi_in_place,numgridsonlevel_mype(level,0:npe-1),npe,mpi_integer,&
                    mpi_max,icomm,ierrmpi)
      call mpi_allreduce(mpi_in_place,numgridsonlevel(level),1,mpi_integer,mpi_sum, &
                      icomm,ierrmpi)
    end do
 
    nx^d=ixmhi^d-ixmlo^d+1;
    nxc^d=nx^d+1;
 
    if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
 
    {^ifoned
    if(convert_type=='teclinempi') then
      nodes=0
      elems=0
      do level=levmin,levmax
        nodes=nodes + numgridsonlevel(level)*{nxc^d*}
        elems=elems + numgridsonlevel(level)*{nx^d*}
      end do
 
      if (mype==0) write(qunit,"(a,i7,a,1pe12.5,a)") &
            'ZONE T="all levels", I=',elems, &
            ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT' 
    
      igonlevel=0
      do morton_no=morton_start(mype),morton_stop(mype)
        igrid = sfc_to_igrid(morton_no)
        call calc_x(igrid,xc,xcc)
        call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
        if (mype==0) then
         {do ix^db=ixccmin^db,ixccmax^db\}
              x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
              w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
             write(qunit,fmt="(100(e14.6))") x_tec, w_tec
         {end do\}
        else if (mype/=0) then
          itag=morton_no
          call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
          call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision,0,itag,icomm,ierrmpi)
          call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
          call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
        end if
      end do
      if(mype==0) then
        do ipe=1,npe-1
          do  morton_no=morton_start(ipe),morton_stop(ipe)
            itag=morton_no
            call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,&
                          itag,icomm,intstatus(:,1),ierrmpi)
            call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,&
                               icomm,intstatus(:,1),ierrmpi)
            call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,&
                               icomm,intstatus(:,1),ierrmpi)
           {do ix^db=ixccmin^db,ixccmax^db\}
              x_tec(1:ndim)=xcc_tmp_recv(ix^d,1:ndim)*normconv(0)
              w_tec(1:nw+nwauxio)=wcc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
              write(qunit,fmt="(100(e14.6))") x_tec, w_tec
           {end do\}
          end do 
        end do
        close(qunit)
      end if
    else
    }
    if(mype/=0) then
      itag=1000*morton_stop(mype)
      call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
      itag=2000*morton_stop(mype)
      call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
    end if
 
    do level=levmin,levmax
       nodesonlevelmype=numgridsonlevel_mype(level,mype)*{nxc^d*}
       elemsonlevelmype=numgridsonlevel_mype(level,mype)*{nx^d*}
       nodesonlevel=numgridsonlevel(level)*{nxc^d*}
       elemsonlevel=numgridsonlevel(level)*{nx^d*}
       ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
       ! with the AMR grid LEVEL. Other options would be
       !    let each grid define a zone: inefficient for TECPLOT internal workings
       !       hence not implemented
       !    let entire octree define 1 zone: no difference in interpolation 
       !       properties across TECPLOT zones detected as yet, hence not done
       select case(convert_type)
         case('tecplotmpi')
           ! in this option, we store the corner coordinates, as well as the corner
           ! values of all variables (obtained by averaging). This allows POINT packaging, 
           ! and thus we can save full grid info by using one call to calc_grid
           if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
            write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
                 'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
                 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
              {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
          do morton_no=morton_start(mype),morton_stop(mype)
             igrid = sfc_to_igrid(morton_no)
             if (mype/=0)then
               itag=morton_no
               call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
               itag=igrid
               call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
             end if
             if (node(plevel_,igrid)/=level) cycle
             call calc_x(igrid,xc,xcc)
             call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                            ixc^l,ixcc^l,.true.)
             if (mype/=0) then
                itag=morton_no
                ind_send=(/ ixc^l /)
                siz_ind=2*^nd
                call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
                call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
    
                call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
                call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
             else  
               {do ix^db=ixcmin^db,ixcmax^db\}
                  x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
                  w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
                  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
               {end do\}
             end if
           end do
         case('tecplotCCmpi')
           ! in this option, we store the corner coordinates, and the cell center
           ! values of all variables. Due to this mix of corner/cell center, we must 
           ! use BLOCK packaging, and thus we have enormous overhead by using 
           ! calc_grid repeatedly to merely fill values of cell corner coordinates 
           ! and cell center values per dimension, per variable
           if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
           if(nw+nwauxio==1)then
             ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
             ! and just set [ndim+1]
             if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
              write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
                'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
                ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
                ndim+1,']=CELLCENTERED), ZONETYPE=', &
             {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
           else
             if(ndim+nw+nwauxio<10) then
               ! difference only in length of integer format specification for ndim+nw+nwauxio
               if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
                write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
                  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
                  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
                  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
               {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
             else
               if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
                write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
                  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
                  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
                  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
               {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
             end if
           end if
 
           do idim=1,ndim
             first=(idim==1)
             do morton_no=morton_start(mype),morton_stop(mype)
               igrid = sfc_to_igrid(morton_no)
               if (mype/=0)then
                 itag=morton_no*idim
                 call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
                 itag=igrid*idim
                 call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
               end if
               if (node(plevel_,igrid)/=level) cycle
               call calc_x(igrid,xc,xcc)
               call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                                ixc^l,ixcc^l,first)
               if (mype/=0)then
                 ind_send=(/ ixc^l /)
                 siz_ind=2*^nd
                 itag=igrid*idim
                 call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
                 call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
                 call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
               else
                 write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
               end if
             end do
           end do
           do iw=1,nw+nwauxio
             do morton_no=morton_start(mype),morton_stop(mype)
               igrid = sfc_to_igrid(morton_no)
               if(mype/=0)then
                 itag=morton_no*(ndim+iw)
                 call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
                 itag=igrid*(ndim+iw)
                 call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
               end if
               if (node(plevel_,igrid)/=level) cycle
               call calc_x(igrid,xc,xcc)
               call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                                 ixc^l,ixcc^l,.true.)
               if(mype/=0)then
                 ind_send=(/ ixcc^l /)
                 siz_ind=2*^nd
                 itag=igrid*(ndim+iw)
                 call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
                 call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
                 call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
               else
                 write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
               end if
             end do
           end do
         case default
           call mpistop('no such tecplot type')
       end select
 
       igonlevel=0
       do morton_no=morton_start(mype),morton_stop(mype)
         igrid = sfc_to_igrid(morton_no)
         if(mype/=0)then
           itag=morton_no
           call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
           itag=igrid
           call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
         end if
         if(node(plevel_,igrid)/=level) cycle
         igonlevel=igonlevel+1
         if(mype/=0)then
           itag=igrid
           call mpi_send(igonlevel,1,mpi_integer, 0,itag,icomm,ierrmpi)
         end if
         if(mype==0)then
           call save_conntec(qunit,igrid,igonlevel)
         end if
       end do
    end do
 
    if(mype==0 .and.npe>1) then
      do ipe=1,npe-1
        itag=1000*morton_stop(ipe)
        call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
        itag=2000*morton_stop(ipe)
        call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
        do level=levmin_recv,levmax_recv
          nodesonlevelmype=numgridsonlevel_mype(level,ipe)*{nxc^d*}
          elemsonlevelmype=numgridsonlevel_mype(level,ipe)*{nx^d*}
          nodesonlevel=numgridsonlevel(level)*{nxc^d*}
          elemsonlevel=numgridsonlevel(level)*{nx^d*}
          select case(convert_type)
           case('tecplotmpi')
              ! in this option, we store the corner coordinates, as well as the corner
              ! values of all variables (obtained by averaging). This allows POINT packaging, 
              ! and thus we can save full grid info by using one call to calc_grid
              if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
              write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
                   'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
                   ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
                {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
              do  morton_no=morton_start(ipe),morton_stop(ipe)
               itag=morton_no
               call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
               itag=igrid_recv
               call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
               if (level_recv/=level) cycle
               itag=morton_no
               siz_ind=2*^nd
               call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,&
                              icomm,intstatus(:,1),ierrmpi)
               ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
               call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
                        ,icomm,intstatus(:,1),ierrmpi)
               call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,&
                              icomm,intstatus(:,1),ierrmpi)
               call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,&
                              icomm,intstatus(:,1),ierrmpi)
               {do ix^db=ixrvcmin^db,ixrvcmax^db\}
                    x_tec(1:ndim)=xc_tmp_recv(ix^d,1:ndim)*normconv(0)
                    w_tec(1:nw+nwauxio)=wc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
                    write(qunit,fmt="(100(e14.6))") x_tec, w_tec
               {end do\}
              end do
           case('tecplotCCmpi')
             ! in this option, we store the corner coordinates, and the cell center
             ! values of all variables. Due to this mix of corner/cell center, we must 
             ! use BLOCK packaging, and thus we have enormous overhead by using 
             ! calc_grid repeatedly to merely fill values of cell corner coordinates 
             ! and cell center values per dimension, per variable
             if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
             if(nw+nwauxio==1)then
               ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
               ! and just set [ndim+1]
               if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
               write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
                  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
                  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
                  ndim+1,']=CELLCENTERED), ZONETYPE=', &
               {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
             else
               if(ndim+nw+nwauxio<10) then
                 ! difference only in length of integer format specification for ndim+nw+nwauxio
                 if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
                 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
                    'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
                    ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
                    ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
                 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
               else
                 if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
                 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
                    'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
                    ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
                    ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
                 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
               end if
             end if
 
             do idim=1,ndim
               do  morton_no=morton_start(ipe),morton_stop(ipe)
                 itag=morton_no*idim
                 call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
                 itag=igrid_recv*idim
                 call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
                 if (level_recv/=level) cycle
                 siz_ind=2*^nd
                 itag=igrid_recv*idim
                 call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
                 ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);     
                 call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
                        ,icomm,intstatus(:,1),ierrmpi)
                 call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
                 write(qunit,fmt="(100(e14.6))") xc_tmp_recv(ixrvc^s,idim)*normconv(0)
               end do
             end do
             do iw=1,nw+nwauxio
               do morton_no=morton_start(ipe),morton_stop(ipe)
                 itag=morton_no*(ndim+iw)
                 call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
                 itag=igrid_recv*(ndim+iw)
                 call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
                 if (level_recv/=level) cycle
                 siz_ind=2*^nd
                 itag=igrid_recv*(ndim+iw)
                 call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
                 ixrvccmin^d=ind_recv(^d);ixrvccmax^d=ind_recv(^nd+^d);
                 call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
                        ,icomm,intstatus(:,1),ierrmpi)
                 call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
                 write(qunit,fmt="(100(e14.6))") wcc_tmp_recv(ixrvcc^s,iw)*normconv(iw)
               end do
             end do
           case default
             call mpistop('no such tecplot type')
          end select
 
          do morton_no=morton_start(ipe),morton_stop(ipe)
            itag=morton_no
            call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            itag=igrid_recv
            call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            if (level_recv/=level) cycle
            itag=igrid_recv
            call mpi_recv(igonlevel_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
            call save_conntec(qunit,igrid_recv,igonlevel_recv)
          end do ! morton loop
        end do ! level loop
      end do ! ipe loop
    end if ! mype=0 if
    {^ifoned endif}
 
    if (npe>1) then
      call mpi_barrier(icomm,ierrmpi)
      if(mype==0)deallocate(intstatus)
    end if
 
  end subroutine tecplot_mpi
 
  subroutine punstructuredvtkb_mpi(qunit)
    ! Write one pvtu and vtu files for each processor
    ! Otherwise like unstructuredvtk_mpi
    ! output for vtu format to paraview, binary version output
    ! uses calc_grid to compute nwauxio variables
    ! allows renormalizing using convert factors
    use mod_forest, only: morton_start, morton_stop, sfc_to_igrid
    use mod_global_parameters
    use mod_calculate_xw
 
    integer, intent(in) ::    qunit
 
    double precision ::  x_VTK(1:3)
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC_TMP
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC_TMP
    double precision :: normconv(0:nw+nwauxio)
    integer*8 :: offset
    integer :: igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,Morton_no
    integer ::               NumGridsOnLevel(1:nlevelshi)
    integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
    integer::  recsep,k,iw,filenr
    integer::  length,lengthcc,offset_points,offset_cells, &
               length_coords,length_conn,length_offsets
    logical ::   fileopen
    character::  buf
    character(len=6)::  bufform
    character(len=80) :: pfilename
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
 
    ! Write pvtu-file:
    if (mype==0) then
      call write_pvtu(qunit)
    end if
    ! Now write the Source files:
    inquire(qunit,opened=fileopen)
    if(.not.fileopen)then
      ! generate filename 
      filenr=snapshotnext-1
      if (autoconvert) filenr=snapshotnext
      ! Open the file for the header part
      write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
      open(qunit,file=pfilename,status='unknown',form='formatted')
    end if
    ! generate xml header
    write(qunit,'(a)')'<?xml version="1.0"?>'
    write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
    write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
    write(qunit,'(a)')'  <UnstructuredGrid>'
    write(qunit,'(a)')'<FieldData>'
    write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                       'NumberOfTuples="1" format="ascii">'
    write(qunit,*) real(global_time*time_convert_factor)
    write(qunit,'(a)')'</DataArray>'
    write(qunit,'(a)')'</FieldData>'
    offset=0
    recsep=4
 
    call getheadernames(wnamei,xandwnamei,outfilehead)
 
    ! number of cells, number of corner points, per grid.
    nx^d=ixmhi^d-ixmlo^d+1;
    nxc^d=nx^d+1;
    nc={nx^d*}
    np={nxc^d*}
 
    length=np*size_real
    lengthcc=nc*size_real
 
    length_coords=3*length
    length_conn=2**^nd*size_int*nc
    length_offsets=nc*size_int
 
    ! Note: using the w_write, writelevel, writespshift
    ! we can clip parts of the grid away, select variables, levels etc.
    do level=levmin,levmax
     if (writelevel(level)) then
       do morton_no=morton_start(mype),morton_stop(mype)
        igrid=sfc_to_igrid(morton_no)
        if (node(plevel_,igrid)/=level) cycle
        ! only output a grid when fully within clipped region selected
        ! by writespshift array
        if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
              *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
             <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
          select case(convert_type)
           case('pvtuBmpi')
             ! we write out every grid as one VTK PIECE
             write(qunit,'(a,i7,a,i7,a)') &
                '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
             write(qunit,'(a)')'<PointData>'
             do iw=1,nw
                if(.not.w_write(iw))cycle
    
                write(qunit,'(a,a,a,i16,a)')&
                    '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
                    '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+length+size_int
             enddo
             do iw=nw+1,nw+nwauxio
    
                write(qunit,'(a,a,a,i16,a)')&
                    '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
                    '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+length+size_int
             enddo
             write(qunit,'(a)')'</PointData>'
    
             write(qunit,'(a)')'<Points>'
             write(qunit,'(a,i16,a)') &
         '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
             ! write cell corner coordinates in a backward dimensional loop, always 3D output
             offset=offset+length_coords+size_int
             write(qunit,'(a)')'</Points>'
           case('pvtuBCCmpi')
             ! we write out every grid as one VTK PIECE
             write(qunit,'(a,i7,a,i7,a)') &
                '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
             write(qunit,'(a)')'<CellData>'
             do iw=1,nw
                if(.not.w_write(iw))cycle
    
                write(qunit,'(a,a,a,i16,a)')&
                    '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
                    '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+lengthcc+size_int
             enddo
             do iw=nw+1,nw+nwauxio
    
                write(qunit,'(a,a,a,i16,a)')&
                    '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
                    '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+lengthcc+size_int
             enddo
             write(qunit,'(a)')'</CellData>'
             write(qunit,'(a)')'<Points>'
             write(qunit,'(a,i16,a)') &
         '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
             ! write cell corner coordinates in a backward dimensional loop, always 3D output
             offset=offset+length_coords+size_int
             write(qunit,'(a)')'</Points>'
          end select
          write(qunit,'(a)')'<Cells>'
          ! connectivity part
          write(qunit,'(a,i16,a)')&
            '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
          offset=offset+length_conn+size_int    
          ! offsets data array
          write(qunit,'(a,i16,a)') &
            '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
          offset=offset+length_offsets+size_int    
          ! VTK cell type data array
          write(qunit,'(a,i16,a)') &
            '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>' 
          offset=offset+size_int+nc*size_int
          write(qunit,'(a)')'</Cells>'
          write(qunit,'(a)')'</Piece>'
        end if
       end do
     end if
    end do
 
    write(qunit,'(a)')'</UnstructuredGrid>'
    write(qunit,'(a)')'<AppendedData encoding="raw">'
    close(qunit)
    ! next to make gfortran compiler happy, as it does not know
    ! form='binary' and produces error on compilation
    !bufform='binary'
    !open(qunit,file=pfilename,form=bufform,position='append')
    !This should in principle do also for gfortran (tested with gfortran 4.6.0 and Intel 11.1):
    open(qunit,file=pfilename,access='stream',form='unformatted',position='append')
    buf='_'
    write(qunit) trim(buf)
 
    do level=levmin,levmax
     if (writelevel(level)) then
       do morton_no=morton_start(mype),morton_stop(mype)
        igrid=sfc_to_igrid(morton_no)
        if (node(plevel_,igrid)/=level) cycle
          ! only output a grid when fully within clipped region selected
          ! by writespshift array
          if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
              *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
             <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
            call calc_x(igrid,xc,xcc)
            call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
                           ixc^l,ixcc^l,.true.)
            do iw=1,nw
              if(.not.w_write(iw))cycle
              select case(convert_type)
                case('pvtuBmpi')
                  write(qunit) length
                  write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
                case('pvtuBCCmpi')
                  write(qunit) lengthcc
                  write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
              end select 
            enddo
            do iw=nw+1,nw+nwauxio
              select case(convert_type)
                case('pvtuBmpi')
                  write(qunit) length
                  write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
                case('pvtuBCCmpi')
                  write(qunit) lengthcc
                  write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
              end select 
            enddo
            write(qunit) length_coords
            {do ix^db=ixcmin^db,ixcmax^db \}
               x_vtk(1:3)=zero;
               x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
               do k=1,3
                 write(qunit) real(x_vtk(k))
               end do
            {end do \}
            write(qunit) length_conn
            call write_connvtk_binary(qunit)
            write(qunit) length_offsets
            do icel=1,nc
              write(qunit) icel*(2**^nd)
            end do
            vtk_type=vtk_cell_type()
            write(qunit) size_int*nc
            do icel=1,nc
              write(qunit) vtk_type
            end do
        end if
      end do
     end if
    end do
 
    close(qunit)
    open(qunit,file=pfilename,status='unknown',form='formatted',position='append')
    write(qunit,'(a)')'</AppendedData>'
    write(qunit,'(a)')'</VTKFile>'
    close(qunit)
 
  end subroutine punstructuredvtkb_mpi
  {^iftwod
  ! subroutines to convert 2.5D data to 3D data
  subroutine unstructuredvtkb23(qunit)
    ! output for vtu format to paraview, binary version output
    ! not parallel, uses calc_grid to compute nwauxio variables
    use mod_global_parameters
    use mod_physics
    use mod_calculate_xw
 
    integer, intent(in) ::    qunit
 
    double precision ::  x_VTK(1:3)
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1&
       -1:ixMhi1,3) :: xC_TMP
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,&
       3)   :: xCC_TMP
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1&
       -1:ixMhi1,nw+nwauxio)   :: wC_TMP
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw&
       +nwauxio)     :: wCC_TMP
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:nw&
       +nwauxio)   :: w
    double precision :: normconv(0:nw+nwauxio)
    double precision :: zlength
    double precision ::d3grid,zlengsc,zgridsc
    integer*8 :: offset
    integer::               igrid,iigrid,level,igonlevel,icel,ixCmin1,ixCmin2,&
       ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,&
       ixCCmax2,ixCCmax3
    integer::               NumGridsOnLevel(1:nlevelshi)
    integer :: nx1,nx2,nx3,nxC1,nxC2,nxC3,nodesonlevel,elemsonlevel,nc,np,&
       VTK_type,ix1,ix2,ix3
    integer :: size_length,recsep,k,iw
    integer :: length,lengthcc,offset_points,offset_cells, length_coords,&
       length_conn,length_offsets
    integer :: i3grid,n3grid
    logical ::   fileopen
    character::  buffer
    character(len=6)::  bufform
    character(len=80)::  filename
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:3+nw+nwauxio)
    character(len=1024) :: outfilehead
 
    if(npe>1)then
      if(mype==0) print *,'unstructuredvtkB23 not parallel, use vtumpi'
      call mpistop('npe>1, unstructuredvtkB23')
    end if
 
    offset=0
    recsep=4
    size_length=4
    inquire(qunit,opened=fileopen)
    if(.not.fileopen)then
      ! generate filename 
      write(filename,'(a,a,i4.4,a)') trim(base_filename),"3D",snapshotini,".vtu"
      ! Open the file for the header part
      open(qunit,file=filename,status='replace')
    endif
    call getheadernames(wnamei,xandwnamei,outfilehead)
    ! generate xml header
    write(qunit,'(a)')'<?xml version="1.0"?>'
    write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
    write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
    write(qunit,'(a)')'<UnstructuredGrid>'
    write(qunit,'(a)')'<FieldData>'
    write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                       'NumberOfTuples="1" format="ascii">'
    write(qunit,'(f10.2)') real(global_time*time_convert_factor)
    write(qunit,'(a)')'</DataArray>'
    write(qunit,'(a)')'</FieldData>'
 
    ! number of cells, number of corner points, per grid.
    nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
    nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
    nc=nx1*nx2*nx3
    np=nxc1*nxc2*nxc3
 
    length=np*size_real
    lengthcc=nc*size_real
 
    length_coords=3*length
    length_conn=2**3*size_int*nc
    length_offsets=nc*size_int
 
    ! Note: using the w_write, writelevel, writespshift
    ! we can clip parts of the grid away, select variables, levels etc.
    zgridsc=2.d0
    zlengsc=2.d0*zgridsc
    zlength=zlengsc*(xprobmax1-xprobmin1)
    do level=levmin,levmax
     if (writelevel(level)) then
       do iigrid=1,igridstail; igrid=igrids(iigrid);
        if (node(plevel_,igrid)/=level) cycle
        block=>ps(igrid)
        ! only output a grid when fully within clipped region selected
        ! by writespshift array
        if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
           *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
           +(xprobmax2-xprobmin2)*writespshift(2,1))&
           .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
           *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
           -xprobmin2)*writespshift(2,2))) then
            d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
          n3grid=nint(zlength/d3grid)
          do i3grid=1,n3grid !subcycles
          select case(convert_type)
           case('vtuB23')
             ! we write out every grid as one VTK PIECE
             write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
                '" NumberOfCells="',nc,'">'
             write(qunit,'(a)')'<PointData>'
             do iw=1,nw
                if(.not.w_write(iw))cycle
                write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                   trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+length+size_int
             enddo
             if(nwauxio>0)then
              do iw=nw+1,nw+nwauxio
                 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                    trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                 write(qunit,'(a)')'</DataArray>'
                 offset=offset+length+size_int
              enddo
             endif
             write(qunit,'(a)')'</PointData>'
    
             write(qunit,'(a)')'<Points>'
             write(qunit,'(a,i16,a)') &
    '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
                offset,'"/>'
             ! write cell corner coordinates in a backward dimensional loop, always 3D output
             offset=offset+length_coords+size_int
             write(qunit,'(a)')'</Points>'
           case('vtuBCC23')
             ! we write out every grid as one VTK PIECE
             write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
                '" NumberOfCells="',nc,'">'
             write(qunit,'(a)')'<CellData>'
             do iw=1,nw
                if(.not.w_write(iw))cycle
                write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                   trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+lengthcc+size_int
             enddo
             if(nwauxio>0)then
              do iw=nw+1,nw+nwauxio
                 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                    trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                 write(qunit,'(a)')'</DataArray>'
                 offset=offset+lengthcc+size_int
              enddo
             endif
             write(qunit,'(a)')'</CellData>'
             write(qunit,'(a)')'<Points>'
             write(qunit,'(a,i16,a)') &
    '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
                offset,'"/>'
             ! write cell corner coordinates in a backward dimensional loop, always 3D output
             offset=offset+length_coords+size_int
             write(qunit,'(a)')'</Points>'
          end select
          write(qunit,'(a)')'<Cells>'
          ! connectivity part
          write(qunit,'(a,i16,a)')&
    '<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
             offset,'"/>'
          offset=offset+length_conn+size_int
          ! offsets data array
          write(qunit,'(a,i16,a)') &
             '<DataArray type="Int32" Name="offsets" format="appended" offset="',&
             offset,'"/>'
          offset=offset+length_offsets+size_int
          ! VTK cell type data array
          write(qunit,'(a,i16,a)') &
             '<DataArray type="Int32" Name="types" format="appended" offset="',&
             offset,'"/>'
          offset=offset+size_length+nc*size_int
          write(qunit,'(a)')'</Cells>'
          write(qunit,'(a)')'</Piece>'
          end do !subcycles
        end if
       end do
     end if
    end do
 
    write(qunit,'(a)')'</UnstructuredGrid>'
    write(qunit,'(a)')'<AppendedData encoding="raw">'
    close(qunit)
    open(qunit,file=filename,form='unformatted',access='stream',status='old',position='append')
    buffer='_'
    write(qunit) trim(buffer)
 
    do level=levmin,levmax
     if (writelevel(level)) then
       do iigrid=1,igridstail; igrid=igrids(iigrid);
          if (node(plevel_,igrid)/=level) cycle
          block=>ps(igrid)
          ! only output a grid when fully within clipped region selected
          ! by writespshift array
          if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
             *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
             +(xprobmax2-xprobmin2)*writespshift(2,1))&
             .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
             *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
             -xprobmin2)*writespshift(2,2))) then
              d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
           n3grid=nint(zlength/d3grid)
           ! In case primitives to be saved: use primitive subroutine
           !  extra layer around mesh only needed when storing corner values and averaging
           if(saveprim) then
             call phys_to_primitive(ixglo1,ixglo2,ixghi1,ixghi2,&
             ixglo1,ixglo2,ixghi1,ixghi2,ps(igrid)%w,ps(igrid)%x)
           endif
           ! using array w so that new output auxiliaries can be calculated by the user
           ! extend 2D data to 3D insuring variables are independent on the third coordinate
           do ix3=ixglo1,ixghi1
            w(ixglo1:ixghi1,ixglo2:ixghi2,ix3,1:nw)=ps(igrid)%w(ixglo1:ixghi1,&
              ixglo2:ixghi2,1:nw)
           end do
           do i3grid=1,n3grid !subcycles
            call calc_grid23(qunit,igrid,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
               ixcmin1,ixcmin2,ixcmin3,ixcmax1,ixcmax2,ixcmax3,ixccmin1,ixccmin2,&
               ixccmin3,ixccmax1,ixccmax2,ixccmax3,.true.,i3grid,d3grid,w,zlength,zgridsc)
            do iw=1,nw
              if(.not.w_write(iw))cycle
              select case(convert_type)
                case('vtuB23')
                  write(qunit) length
                  write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                     =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
                case('vtuBCC23')
                  write(qunit) lengthcc
                  write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                     =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
                     =ixccmin3,ixccmax3)
              end select
            enddo
            if(nwauxio>0)then
             do iw=nw+1,nw+nwauxio
               select case(convert_type)
                 case('vtuB23')
                   write(qunit) length
                   write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                      =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
                 case('vtuBCC23')
                   write(qunit) lengthcc
                   write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                      =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
                      =ixccmin3,ixccmax3)
               end select
             end do
            end if
            write(qunit) length_coords
            do ix3=ixcmin3,ixcmax3
            do ix2=ixcmin2,ixcmax2
            do ix1=ixcmin1,ixcmax1
              x_vtk(1:3)=zero;
              x_vtk(1:3)=xc_tmp(ix1,ix2,ix3,1:3)*normconv(0);
              do k=1,3
               write(qunit) real(x_vtk(k))
              end do
            end do
            end do
            end do
            write(qunit) length_conn
            do ix3=1,nx3
            do ix2=1,nx2
            do ix1=1,nx1
              write(qunit)&
              (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
              (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
              (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1-1,&
              (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1,&
               ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
               ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
               ix3*nxc2*nxc1+    ix2*nxc1+ix1-1,&
               ix3*nxc2*nxc1+    ix2*nxc1+ix1
            end do
            end do
            end do
            write(qunit) length_offsets
            do icel=1,nc
              write(qunit) icel*(2**3)
            end do
            vtk_type=11
            write(qunit) size_int*nc
            do icel=1,nc
              write(qunit) vtk_type
            end do
           end do !subcycles
          end if
       end do
     end if
    end do
 
    close(qunit)
    open(qunit,file=filename,status='unknown',form='formatted',position='append')
 
    write(qunit,'(a)')'</AppendedData>'
    write(qunit,'(a)')'</VTKFile>'
    close(qunit)
 
  end subroutine unstructuredvtkb23
 
  subroutine unstructuredvtkbsym23(qunit)
    ! output for vtu format to paraview, binary version output
    ! not parallel, uses calc_grid to compute nwauxio variables
    ! use this subroutine  when the physical domain is symmetric/asymmetric about (0,y,z) 
    ! plane, xprobmin1=0 and the computational domain is a half of the physical domain
    use mod_global_parameters
    use mod_physics
    use mod_calculate_xw
 
    integer, intent(in) ::    qunit
 
    double precision ::  x_VTK(1:3)
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1&
       -1:ixMhi1,3) :: xC_TMP
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,&
       3)   :: xCC_TMP
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1&
       -1:ixMhi1,nw+nwauxio)   :: wC_TMP
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw&
       +nwauxio)     :: wCC_TMP
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:nw&
       +nwauxio)   :: w
    double precision :: normconv(0:nw+nwauxio)
    double precision ::d3grid,zlengsc,zgridsc
    double precision :: zlength
    integer*8 :: offset
    integer::               igrid,iigrid,level,igonlevel,icel,ixCmin1,ixCmin2,&
       ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,&
       ixCCmax2,ixCCmax3
    integer::               NumGridsOnLevel(1:nlevelshi)
    integer :: nx1,nx2,nx3,nxC1,nxC2,nxC3,nodesonlevel,elemsonlevel,nc,np,&
       VTK_type,ix1,ix2,ix3
    integer :: size_length,recsep,k,iw
    integer :: length,lengthcc,offset_points,offset_cells, length_coords,&
       length_conn,length_offsets
    integer :: i3grid,n3grid
    logical ::   fileopen
    character(len=80)::  filename
    character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:3+nw+nwauxio)
    character(len=1024) :: outfilehead
    character::  buffer
    character(len=6)::  bufform
 
    if(npe>1)then
      if(mype==0) print *,'unstructuredvtkBsym23 not parallel, use vtumpi'
      call mpistop('npe>1, unstructuredvtkBsym23')
    end if
 
    offset=0
    recsep=4
    size_length=4
 
    inquire(qunit,opened=fileopen)
    if(.not.fileopen)then
      ! generate filename 
      write(filename,'(a,a,i4.4,a)') trim(base_filename),"3D",snapshotini,".vtu"
      ! Open the file for the header part
      open(qunit,file=filename,status='unknown')
    end if
 
    call getheadernames(wnamei,xandwnamei,outfilehead)
    ! generate xml header
    write(qunit,'(a)')'<?xml version="1.0"?>'
    write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
    write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
    write(qunit,'(a)')'<UnstructuredGrid>'
    write(qunit,'(a)')'<FieldData>'
    write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                       'NumberOfTuples="1" format="ascii">'
    write(qunit,'(f10.2)') real(global_time*time_convert_factor)
    write(qunit,'(a)')'</DataArray>'
    write(qunit,'(a)')'</FieldData>'
 
    ! number of cells, number of corner points, per grid.
    nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
    nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
    nc=nx1*nx2*nx3
    np=nxc1*nxc2*nxc3
 
    length=np*size_real
    lengthcc=nc*size_real
 
    length_coords=3*length
    length_conn=2**3*size_int*nc
    length_offsets=nc*size_int
 
    ! Note: using the w_write, writelevel, writespshift
    ! we can clip parts of the grid away, select variables, levels etc.
    zlengsc=4.d0
    zgridsc=2.d0
    zlength=zlengsc*(xprobmax1-xprobmin1)
    do level=levmin,levmax
     if (writelevel(level)) then
       do iigrid=1,igridstail; igrid=igrids(iigrid);
        if (node(plevel_,igrid)/=level) cycle
        block=>ps(igrid)
        ! only output a grid when fully within clipped region selected
        ! by writespshift array
        if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
           *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
           +(xprobmax2-xprobmin2)*writespshift(2,1))&
           .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
           *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
           -xprobmin2)*writespshift(2,2))) then
         d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
         n3grid=nint(zlength/d3grid)
         do i3grid=1,n3grid !subcycles
          !! original domain ----------------------------------start 
          select case(convert_type)
           case('vtuBsym23')
             ! we write out every grid as one VTK PIECE
             write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
                '" NumberOfCells="',nc,'">'
             write(qunit,'(a)')'<PointData>'
             do iw=1,nw
                if(.not.w_write(iw))cycle
                write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                   trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+length+size_length
             enddo
             if(nwauxio>0)then
              do iw=nw+1,nw+nwauxio
                 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                    trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                 write(qunit,'(a)')'</DataArray>'
                 offset=offset+length+size_length
              enddo
             endif
             write(qunit,'(a)')'</PointData>'
             write(qunit,'(a)')'<Points>'
             write(qunit,'(a,i16,a)') &
    '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
                offset,'"/>'
             ! write cell corner coordinates in a backward dimensional loop, always 3D output
             offset=offset+length_coords+size_length
             write(qunit,'(a)')'</Points>'
           case('vtuBCCsym23')
             ! we write out every grid as one VTK PIECE
             write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
                '" NumberOfCells="',nc,'">'
             write(qunit,'(a)')'<CellData>'
             do iw=1,nw
                if(.not.w_write(iw))cycle
                write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                   trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+lengthcc+size_length
             enddo
             if(nwauxio>0)then
              do iw=nw+1,nw+nwauxio
                 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                    trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                 write(qunit,'(a)')'</DataArray>'
                 offset=offset+lengthcc+size_length
              enddo
             endif
             write(qunit,'(a)')'</CellData>'
    
             write(qunit,'(a)')'<Points>'
             write(qunit,'(a,i16,a)') &
    '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
                offset,'"/>'
             ! write cell corner coordinates in a backward dimensional loop, always 3D output
             offset=offset+length_coords+size_length
             write(qunit,'(a)')'</Points>'
          end select
          write(qunit,'(a)')'<Cells>'
          ! connectivity part
          write(qunit,'(a,i16,a)')&
    '<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
             offset,'"/>'
          offset=offset+length_conn+size_length
          ! offsets data array
          write(qunit,'(a,i16,a)') &
             '<DataArray type="Int32" Name="offsets" format="appended" offset="',&
             offset,'"/>'
          offset=offset+length_offsets+size_length
          ! VTK cell type data array
          write(qunit,'(a,i16,a)') &
             '<DataArray type="Int32" Name="types" format="appended" offset="',&
             offset,'"/>'
          offset=offset+size_length+nc*size_int
          write(qunit,'(a)')'</Cells>'
          write(qunit,'(a)')'</Piece>'
          !! original domain ----------------------------------end 
          !! symetric/asymetric mirror domain -----------------start
          select case(convert_type)
           case('vtuBsym23')
             ! we write out every grid as one VTK PIECE
             write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
                '" NumberOfCells="',nc,'">'
             write(qunit,'(a)')'<PointData>'
             do iw=1,nw
                if(.not.w_write(iw))cycle
                write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                   trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+length+size_length
             enddo
             if(nwauxio>0)then
              do iw=nw+1,nw+nwauxio
                 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                    trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                 write(qunit,'(a)')'</DataArray>'
                 offset=offset+length+size_length
              enddo
             endif
             write(qunit,'(a)')'</PointData>'
             write(qunit,'(a)')'<Points>'
             write(qunit,'(a,i16,a)') &
    '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
                offset,'"/>'
             ! write cell corner coordinates in a backward dimensional loop, always 3D output
             offset=offset+length_coords+size_length
             write(qunit,'(a)')'</Points>'
           case('vtuBCCsym23')
             ! we write out every grid as one VTK PIECE
             write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
                '" NumberOfCells="',nc,'">'
             write(qunit,'(a)')'<CellData>'
             do iw=1,nw
                if(.not.w_write(iw))cycle
                write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                   trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                write(qunit,'(a)')'</DataArray>'
                offset=offset+lengthcc+size_length
             enddo
             if(nwauxio>0)then
              do iw=nw+1,nw+nwauxio
                 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
                    trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
                 write(qunit,'(a)')'</DataArray>'
                 offset=offset+lengthcc+size_length
              enddo
             endif
             write(qunit,'(a)')'</CellData>'
             write(qunit,'(a)')'<Points>'
             write(qunit,'(a,i16,a)') &
    '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
                offset,'"/>'
             ! write cell corner coordinates in a backward dimensional loop, always 3D output
             offset=offset+length_coords+size_length
             write(qunit,'(a)')'</Points>'
          end select
          write(qunit,'(a)')'<Cells>'
          ! connectivity part
          write(qunit,'(a,i16,a)')&
    '<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
             offset,'"/>'
          offset=offset+length_conn+size_length
          ! offsets data array
          write(qunit,'(a,i16,a)') &
             '<DataArray type="Int32" Name="offsets" format="appended" offset="',&
             offset,'"/>'
          offset=offset+length_offsets+size_length
          ! VTK cell type data array
          write(qunit,'(a,i16,a)') &
             '<DataArray type="Int32" Name="types" format="appended" offset="',&
             offset,'"/>'
          offset=offset+size_length+nc*size_int
          write(qunit,'(a)')'</Cells>'
          write(qunit,'(a)')'</Piece>'
          !! symetric/asymetric mirror domain -----------------end
         end do !subcycles
        end if
       end do
     end if
    end do
 
    write(qunit,'(a)')'</UnstructuredGrid>'
    write(qunit,'(a)')'<AppendedData encoding="raw">'
    close(qunit)
    open(qunit,file=filename,form='unformatted',access='stream',status='old',position='append')
    buffer='_'
    write(qunit) trim(buffer)
    do level=levmin,levmax
     if (writelevel(level)) then
       do iigrid=1,igridstail; igrid=igrids(iigrid);
          if (node(plevel_,igrid)/=level) cycle
          block=>ps(igrid)
          ! only output a grid when fully within clipped region selected
          ! by writespshift array
          if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
             *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
             +(xprobmax2-xprobmin2)*writespshift(2,1))&
             .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
             *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
             -xprobmin2)*writespshift(2,2))) then
           d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
           n3grid=nint(zlength/d3grid)
           ! In case primitives to be saved: use primitive subroutine
           !  extra layer around mesh only needed when storing corner values and averaging
           if(saveprim) then
             call phys_to_primitive(ixglo1,ixglo2,ixghi1,ixghi2,&
             ixglo1,ixglo2,ixghi1,ixghi2,ps(igrid)%w,ps(igrid)%x)
           endif
           ! using array w so that new output auxiliaries can be calculated by the user
           ! extend 2D data to 3D insuring variables are independent on the third coordinate
           do ix3=ixglo1,ixghi1
            w(ixglo1:ixghi1,ixglo2:ixghi2,ix3,1:nw)=ps(igrid)%w(ixglo1:ixghi1,&
              ixglo2:ixghi2,1:nw)
           end do
           do i3grid=1,n3grid !subcycles
            call calc_grid23(qunit,igrid,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
               ixcmin1,ixcmin2,ixcmin3,ixcmax1,ixcmax2,ixcmax3,ixccmin1,ixccmin2,&
               ixccmin3,ixccmax1,ixccmax2,ixccmax3,.true.,i3grid,d3grid,w,zlength,zgridsc)
            !! original domain ----------------------------------start 
            do iw=1,nw
              if(.not.w_write(iw))cycle
              select case(convert_type)
                case('vtuBsym23')
                  write(qunit) length
                  write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                     =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
                case('vtuBCCsym23')
                  write(qunit) lengthcc
                  write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                     =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
                     =ixccmin3,ixccmax3)
              end select
            enddo
            if(nwauxio>0)then
             do iw=nw+1,nw+nwauxio
               select case(convert_type)
                 case('vtuBsym23')
                   write(qunit) length
                   write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                      =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
                 case('vtuBCCsym23')
                   write(qunit) lengthcc
                   write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                      =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
                      =ixccmin3,ixccmax3)
               end select
             enddo
            endif
            write(qunit) length_coords
            do ix3=ixcmin3,ixcmax3
            do ix2=ixcmin2,ixcmax2
            do ix1=ixcmin1,ixcmax1
              x_vtk(1:3)=zero;
              x_vtk(1:3)=xc_tmp(ix1,ix2,ix3,1:3)*normconv(0);
              do k=1,3
               write(qunit) real(x_vtk(k))
              end do
            end do
            end do
            end do
            write(qunit) length_conn
            do ix3=1,nx3
            do ix2=1,nx2
            do ix1=1,nx1
              write(qunit)&
              (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
              (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
              (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1-1,&
              (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1,&
               ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
               ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
               ix3*nxc2*nxc1+    ix2*nxc1+ix1-1,&
               ix3*nxc2*nxc1+    ix2*nxc1+ix1
            end do
            end do
            end do
            write(qunit) length_offsets
            do icel=1,nc
               write(qunit) icel*(2**3)
            end do
            vtk_type=11
            write(qunit) size_int*nc
            do icel=1,nc
             write(qunit) vtk_type
            end do
            !! original domain ----------------------------------end 
            !! symetric/asymetric mirror domain -----------------start
            do iw=1,nw
              if(.not.w_write(iw))cycle
              if(iw==2 .or. iw==4 .or. iw==7) then
               wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,iw)=&
               -wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,iw)
               wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,iw)=&
               -wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,iw)
              end if
              select case(convert_type)
                case('vtuBsym23')
                  write(qunit) length
                  write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                     =ixcmax1,ixcmin1,-1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
                case('vtuBCCsym23')
                  write(qunit) lengthcc
                  write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                     =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
                     =ixccmin3,ixccmax3)
              end select
            enddo
            if(nwauxio>0)then
             do iw=nw+1,nw+nwauxio
               select case(convert_type)
                 case('vtuBsym23')
                   write(qunit) length
                   write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                      =ixcmax1,ixcmin1,-1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
                 case('vtuBCCsym23')
                   write(qunit) lengthcc
                   write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
                      =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
                      =ixccmin3,ixccmax3)
               end select
             end do
            end if
            write(qunit) length_coords
            do ix3=ixcmin3,ixcmax3
            do ix2=ixcmin2,ixcmax2
            do ix1=ixcmax1,ixcmin1,-1
              x_vtk(1:3)=zero;
              x_vtk(1:3)=xc_tmp(ix1,ix2,ix3,1:3)*normconv(0);
              x_vtk(1)=-x_vtk(1)
              do k=1,3
               write(qunit) real(x_vtk(k))
              end do
            end do
            end do
            end do
            write(qunit) length_conn
            do ix3=1,nx3
            do ix2=1,nx2
            do ix1=nx1,1,-1
              write(qunit)&
              (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
              (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
              (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1-1,&
              (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1,&
               ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
               ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
               ix3*nxc2*nxc1+    ix2*nxc1+ix1-1,&
               ix3*nxc2*nxc1+    ix2*nxc1+ix1
            end do
            end do
            end do
            write(qunit) length_offsets
            do icel=1,nc
               write(qunit) icel*(2**3)
            end do
            vtk_type=11
            write(qunit) size_int*nc
            do icel=1,nc
             write(qunit) vtk_type
            end do
            !! symetric/asymetric mirror domain -----------------end
           end do !subcycles
          end if
       end do
     end if
    end do
    close(qunit)
    open(qunit,file=filename,status='unknown',form='formatted',position='append')
    write(qunit,'(a)')'</AppendedData>'
    write(qunit,'(a)')'</VTKFile>'
    close(qunit)
 
  end subroutine unstructuredvtkbsym23
 
  subroutine calc_grid23(qunit,igrid,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,normconv,&
  ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,&
   ixCCmax1,ixCCmax2,ixCCmax3,first,i3grid,d3grid,w,zlength,zgridsc)
    ! this subroutine computes both corner as well as cell-centered values
    ! it handles how we do the center to corner averaging, as well as 
    ! whether we switch to cartesian or want primitive or conservative output,
    ! handling the addition of B0 in B0+B1 cases, ...
    ! the normconv is passed on to specialvar_output for extending with
    ! possible normalization values for the nw+1:nw+nwauxio entries
    use mod_global_parameters
    integer, intent(in) :: qunit, igrid,i3grid
    logical, intent(in) :: first
 
    double precision :: dx1,dx2,dx3,d3grid,zlength,zgridsc
    double precision :: ldw(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1),&
        dwC(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1)
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1&
       -1:ixMhi1,3) :: xC
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,&
       3)   :: xCC
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1&
       -1:ixMhi1,nw+nwauxio)   :: wC
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw&
       +nwauxio)     :: wCC
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1&
       -1:ixMhi1,3) :: xC_TMP
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,&
       3)   :: xCC_TMP
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1&
       -1:ixMhi1,nw+nwauxio)   :: wC_TMP
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw&
       +nwauxio)     :: wCC_TMP
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:nw&
       +nwauxio)   :: w
    double precision,dimension(0:nw+nwauxio)       :: normconv
    integer :: nx1,nx2,nx3, nxC1,nxC2,nxC3, ix1,ix2,ix3, ix, iw, level, idir
    integer :: ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,&
       ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3,nxCC1,nxCC2,nxCC3
    integer :: idims,jxCmin1,jxCmin2,jxCmin3,jxCmax1,jxCmax2,jxCmax3
    logical, save :: subfirst=.true.
 
    ! following only for allowing compiler to go through with debug on
    nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
    level=node(plevel_,igrid)
    dx1=dx(1,level);dx2=dx(2,level);dx3=zgridsc*dx(1,level);
    ! for normalization within the code
    if(saveprim) then
      normconv(0) = length_convert_factor
      normconv(1:nw) = w_convert_factor
    else
      normconv(0)=length_convert_factor
      ! assuming density
      normconv(1)=w_convert_factor(1)
      ! assuming momentum=density*velocity
      if (nw>=2) normconv(2:2+3)=w_convert_factor(1)*w_convert_factor(2:2+3)
      ! assuming energy/pressure and magnetic field
      if (nw>=2+3) normconv(2+3:nw)=w_convert_factor(2+3:nw)
    end if
    ! coordinates of cell centers
    nxcc1=nx1;nxcc2=nx2;nxcc3=nx3;
    ixccmin1=ixmlo1;ixccmin2=ixmlo2;ixccmin3=ixmlo1; ixccmax1=ixmhi1
    ixccmax2=ixmhi2;ixccmax3=ixmhi1;
    do ix=ixccmin1,ixccmax1
      xcc(ix,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1)=rnode(rpxmin1_,igrid)&
          +(dble(ix-ixccmin1)+half)*dx1
    end do
    do ix=ixccmin2,ixccmax2
      xcc(ixccmin1:ixccmax1,ix,ixccmin3:ixccmax3,2)=rnode(rpxmin2_,igrid)&
          +(dble(ix-ixccmin2)+half)*dx2
    end do
    do ix=ixccmin3,ixccmax3
      xcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ix,3)=-zlength/two+&
         dble(i3grid-1)*d3grid+(dble(ix-ixccmin3)+half)*dx3
    end do
 
    ! coordinates of cell corners
    nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
    ixcmin1=ixmlo1-1;ixcmin2=ixmlo2-1;ixcmin3=ixmlo1-1; ixcmax1=ixmhi1
    ixcmax2=ixmhi2;ixcmax3=ixmhi1;
    do ix=ixcmin1,ixcmax1
      xc(ix,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1)=rnode(rpxmin1_,igrid)&
          +dble(ix-ixcmin1)*dx1
    end do
    do ix=ixcmin2,ixcmax2
      xc(ixcmin1:ixcmax1,ix,ixcmin3:ixcmax3,2)=rnode(rpxmin2_,igrid)&
          +dble(ix-ixcmin2)*dx2
    end do
    do ix=ixcmin3,ixcmax3
      xc(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ix,3)=-zlength/two+&
          dble(i3grid-1)*d3grid+dble(ix-ixcmin3)*dx3
    end do
 
    if (nwextra>0) then
     ! here we actually fill the ghost layers for the nwextra variables using 
     ! continuous extrapolation (as these values do not exist normally in ghost
     ! cells)
     do idims=1,3
      select case(idims)
       case(1)
        jxcmin1=ixghi1+1-nghostcells;jxcmin2=ixglo2;jxcmin3=ixglo1;
        jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixghi1;
         do ix1=jxcmin1,jxcmax1
             w(ix1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw) = w(jxcmin1&
                -1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw)
         end do
         jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixglo1;
         jxcmax1=ixglo1-1+nghostcells;jxcmax2=ixghi2;jxcmax3=ixghi1;
         do ix1=jxcmin1,jxcmax1
             w(ix1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw) = w(jxcmax1&
                +1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw)
         end do
       case(2)
         jxcmin1=ixglo1;jxcmin2=ixghi2+1-nghostcells;jxcmin3=ixglo1;
         jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixghi1;
         do ix2=jxcmin2,jxcmax2
             w(jxcmin1:jxcmax1,ix2,jxcmin3:jxcmax3,nw-nwextra+1:nw) &
                = w(jxcmin1:jxcmax1,jxcmin2-1,jxcmin3:jxcmax3,nw-nwextra+1:nw)
         end do
         jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixglo1;
         jxcmax1=ixghi1;jxcmax2=ixglo2-1+nghostcells;jxcmax3=ixghi1;
         do ix2=jxcmin2,jxcmax2
             w(jxcmin1:jxcmax1,ix2,jxcmin3:jxcmax3,nw-nwextra+1:nw) &
                = w(jxcmin1:jxcmax1,jxcmax2+1,jxcmin3:jxcmax3,nw-nwextra+1:nw)
         end do
       case(3)
         jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixghi1+1-nghostcells;
         jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixghi1;
         do ix3=jxcmin3,jxcmax3
             w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,ix3,nw-nwextra+1:nw) &
                = w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,jxcmin3-1,nw-nwextra+1:nw)
         end do
         jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixglo1;
         jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixglo1-1+nghostcells;
         do ix3=jxcmin3,jxcmax3
             w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,ix3,nw-nwextra+1:nw) &
                = w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,jxcmax3+1,nw-nwextra+1:nw)
         end do
      end select
     end do
    end if
    ! next lines needed when specialvar_output uses gradients
    ! and later on when dwlimiter2 is used 
    if(nwauxio>0)then
      ! auxiliary io variables can be computed and added by user
      ! next few lines ensure correct usage of routines like divvector etc
      dxlevel(1)=rnode(rpdx1_,igrid);dxlevel(2)=rnode(rpdx2_,igrid)
      ! default (no) normalization for auxiliary variables
      normconv(nw+1:nw+nwauxio)=one
      ! maybe need for restriction to ixG^LL^LSUB1 
      call specialvar_output23(ixglo1,ixglo2,ixglo1,ixghi1,ixghi2,ixghi1,ixglo1&
         +1,ixglo2+1,ixglo1+1,ixghi1-1,ixghi2-1,ixghi1-1,w,xcc,normconv)
    endif
    ! compute the cell-center values for w first
    !===========================================
    ! cell center values obtained from mere copy, while B0+B1 split handled here
    wcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,:)=w(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,:)
    if(b0field) then
      do ix3=ixccmin3,ixccmax3
        do ix2=ixccmin2,ixccmax2
          do ix1=ixccmin1,ixccmax1
            wcc(ix1,ix2,ix3,iw_mag(:))=wcc(ix1,ix2,ix3,iw_mag(:))+ps(igrid)%B0(ix1,ix2,&
                :,0)
          end do
        end do
      end do
    end if
    if(.not.saveprim .and. b0field .and. iw_e>0) then
       do ix3=ixccmin3,ixccmax3
       do ix2=ixccmin2,ixccmax2
       do ix1=ixccmin1,ixccmax1
           wcc(ix1,ix2,ix3,iw_e)=w(ix1,ix2,ix3,iw_e) +half*sum(ps(igrid)%B0(ix1,&
              ix2,:,0)**2 ) + sum(w(ix1,ix2,ix3,&
              iw_mag(:))*ps(igrid)%B0(ix1,ix2,:,0))
       end do
       end do
       end do
    end if
    ! compute the corner values for w now by averaging
    !=================================================
    if(slab_uniform)then
       ! for slab symmetry: no geometrical info required
       do iw=1,nw+nwauxio
          if (b0field.and.iw>iw_mag(1)-1.and.iw<=iw_mag(ndir)) then
             idir=iw-iw_mag(1)+1
             do ix3=ixcmin3,ixcmax3
             do ix2=ixcmin2,ixcmax2
             do ix1=ixcmin1,ixcmax1
               wc(ix1,ix2,ix3,iw)=sum(w(ix1:ix1+1,ix2:ix2+1,ix3,iw) &
                  +ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
                  ,idir,0))/dble(2**3)+&
                  sum(w(ix1:ix1+1,ix2:ix2+1,ix3+1,iw) &
                  +ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
                  ,idir,0))/dble(2**3)
             end do
             end do
             end do
          else
             do ix3=ixcmin3,ixcmax3
             do ix2=ixcmin2,ixcmax2
             do ix1=ixcmin1,ixcmax1
                wc(ix1,ix2,ix3,iw)=sum(w(ix1:ix1+1,ix2:ix2+1,ix3:ix3&
                   +1,iw))/dble(2**3)
             end do
             end do
             end do
          end if
       end do
       if(.not.saveprim .and. b0field .and. iw_e>0) then
          do ix3=ixcmin3,ixcmax3
          do ix2=ixcmin2,ixcmax2
          do ix1=ixcmin1,ixcmax1
             wc(ix1,ix2,ix3,iw_e)=sum( w(ix1:ix1+1,ix2:ix2+1,ix3,iw_e) &
                +half*sum(ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
                ,:,0)**2,dim=ndim+1) + sum( w(ix1:ix1+1,ix2:ix2+1,ix3&
                ,iw_mag(:))*ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
                ,:,0),dim=ndim+1) ) /dble(2**3)+&
                sum( w(ix1:ix1+1,ix2:ix2+1,ix3+1,iw_e) &
                +half*sum(ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
                ,:,0)**2,dim=ndim+1) + sum( w(ix1:ix1+1,ix2:ix2+1,ix3&
                +1,iw_mag(:))*ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
                ,:,0),dim=ndim+1) ) /dble(2**3)
          end do
          end do
          end do
       end if
    end if
    ! keep the coordinate and vector components
    xc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:3)          &
       = xc(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:3)
    wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw&
       +nwauxio)    = wc(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw&
       +nwauxio)
    xcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,&
       1:3)        = xcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,&
       ixccmin3:ixccmax3,1:3)
    wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:nw&
       +nwauxio)  = wcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,&
       1:nw+nwauxio)
  end subroutine calc_grid23
 
  subroutine save_connvtk23(qunit,igrid)
    ! this saves the basic line, pixel and voxel connectivity,
    ! as used by VTK file outputs for unstructured grid
    use mod_global_parameters
 
    integer, intent(in) :: qunit, igrid
 
    integer :: nx1,nx2,nx3, nxC1,nxC2,nxC3, ix1,ix2,ix3
 
    nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
    nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
    do ix3=1,nx3
    do ix2=1,nx2
    do ix1=1,nx1
          write(qunit,'(8(i7,1x))')&
                     (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
                     (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
                     (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1-1,&
                     (ix3-1)*nxc2*nxc1+    ix2*nxc1+ix1,&
                         ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
                         ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
                         ix3*nxc2*nxc1+    ix2*nxc1+ix1-1,&
                         ix3*nxc2*nxc1+    ix2*nxc1+ix1
    
    end do
    end do
    end do
  end subroutine save_connvtk23
 
  subroutine specialvar_output23(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
   ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,w,x,normconv)
    ! this subroutine can be used in convert, to add auxiliary variables to the
    ! converted output file, for further analysis using tecplot, paraview, ....
    ! these auxiliary values need to be stored in the nw+1:nw+nwauxio slots
    ! the array normconv can be filled in the (nw+1:nw+nwauxio) range with
    ! corresponding normalization values (default value 1)
    use mod_global_parameters
 
    integer, intent(in)           :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
       ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
    double precision, intent(in)  :: x(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3,1:3)
    double precision              :: w(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3,nw+nwauxio)
    double precision              :: normconv(0:nw+nwauxio)
    
    double precision :: qvec(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:ndir),&
       curlvec(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:ndir)
    integer                       :: idirmin
 
    ! output Te
    !if(saveprim)then
    !  w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+1)=w(ixOmin1:ixOmax1,&
    !   ixOmin2:ixOmax2,ixOmin3:ixOmax3,iw_e)/w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
    !   ixOmin3:ixOmax3,iw_rho)
    !endif
    !!! store current
    ! qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,1)=w(ixImin1:ixImax1,&
    !    ixImin2:ixImax2,ixImin3:ixImax3,mag(1))
    ! qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,2)=w(ixImin1:ixImax1,&
    !    ixImin2:ixImax2,ixImin3:ixImax3,mag(2))
    ! qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,3)=w(ixImin1:ixImax1,&
    !    ixImin2:ixImax2,ixImin3:ixImax3,mag(3));
    !call curlvector3D(qvec,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,ixOmin1,&
    !   ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,curlvec,idirmin,1,ndir)
    !w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+2)=curlvec&
    !   (ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,1)
    !w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+3)=curlvec&
    !   (ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,2)
    !w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+4)=curlvec&
    !   (ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,3);
  end subroutine specialvar_output23 
  \}
 
end module mod_convert_files