program slices
      implicit none
! Converts 2D slices from var.dat snapshot to the Dx readable 
! binary data files.
! The data are read in again by .dx routines on regular rectangular
! mesh. The dx_slices.f90 should be copied to the 
! directory of the run. Compile dx_slices.f90 and run it. 
! .net files should be linked to run directory (from pencil-code/dx).
! Then, Data Explorer is ready to be started from run directory.
! This slices routine is better to be integrated into the main code - to
! do in the future.

! Include grid sizes
      include 'src/cparam.local'
      include 'src/cparam.inc'

! This was taken from cparam.f90 file
  integer, parameter :: nghost=3
  integer, parameter :: nx=nxgrid,ny=nygrid/nprocy,nz=nzgrid/nprocz,nw=nx*ny*nz
  integer, parameter :: mx=nx+2*nghost,l1=1+nghost,l2=mx-nghost
  integer, parameter :: my=ny+2*nghost,m1=1+nghost,m2=my-nghost
  integer, parameter :: mz=nz+2*nghost,n1=1+nghost,n2=mz-nghost
  integer, parameter :: mw=mx*my*mz
  integer, parameter :: l1i=l1+nghost-1,l2i=l2-nghost+1
  integer, parameter :: m1i=m1+nghost-1,m2i=m2-nghost+1
  integer, parameter :: n1i=n1+nghost-1,n2i=n2-nghost+1

! These are definitions intrincis for todx_magn.f90
      real, dimension(mx,my,mz,mvar) :: f
      real, dimension(mx) :: x
      real, dimension(my) :: y
      real, dimension(mz) :: z    
      real, dimension(mx,my,mz,3) :: bbb
      real :: x00,y00,z00,t,dx,dy,dz
      integer :: i,j,k,n
      integer :: iux,iuy,iuz,ilnrho,iss,iax,iay,iaz
      character (len=2) :: sltype
      integer :: isltype,isl,nfile

! Read snapshot file
      open(unit=17, file='./tmp/proc0/var.dat', form='unformatted')
      read(17) f
      read(17) t,x,y,z,dx,dy,dz
      close(17)
      print*,'t=',t,' dx=',dx,' dy=',dy,' dz=',dz
      print*,'mx=',mx,' my=',my,' mz=',mz

! Define here the position of the corner of the box formed by only 
! true grid points
      x00=x(nghost+1)
      y00=y(nghost+1)
      z00=z(nghost+1)
      print*, 'x0=',x00,' y0=',y00,' z0=',z00
      print*,'mvar=',mvar

! Read in the type and position of the slice
      print*,'type of slice, xy, yz, or xz ?'
      read*, sltype
      print*,'number of the slice files (must be from 1 to 9) ?'
      read*, nfile
      select case(sltype)
         case('xy') 
            isltype=1
            print*,'position of the slice on z-axis of true grid, isl ?'
            print*,'must be 1 .le. isl .le.',nzgrid
            read*,isl
            if(isl.lt.1.or.isl.gt.nzgrid) then
               print*,'isl is out of grid limits'
               stop
            endif   

         case('yz')
            isltype=2
            print*,'position of the slice on x-axis of true grid, isl ?'
            print*,'must be 1 .le. isl .le.',nxgrid
            read*,isl
            if(isl.lt.1.or.isl.gt.nxgrid) then
               print*,'isl is out of grid limits'
               stop
            endif           

         case('xz')
            isltype=3
            print*,'position of the slice on y-axis of true grid, isl ?'
            print*,'must be 1 .le. isl .le.',nygrid
            read*,isl
            if(isl.lt.1.or.isl.gt.nygrid) then
               print*,'isl is out of grid limits'
               stop
            endif

         case default
        !
        !  Catch unknown values
        !
        print*,'No such value for sltype: ', trim(sltype)
        stop
      endselect     

!!
!!! Write out slices -- global select 
!!
      select case(sltype)

!!!!!!!
!!! Case of xy-slice
!!!!!!!
      case('xy')

! Write out parallel velocity field in the slice
      iux=1
      iuy=2
      iuz=3
      open(unit=18, file='./tmp/uuu_dx_sl'//char(nfile+48)//'.bin', &
        form='unformatted')
      write(18) (((real(f(i,j,isl,n)),n=iux,iuy),j=nghost+1,my-nghost), &
        i=nghost+1,mx-nghost)
      close(18)

! Write out ln density field for Dx
      if(mvar.ge.4) then
      ilnrho=4
      open(unit=19, file='./tmp/lnrho_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(19) ((real(f(i,j,isl,ilnrho)),j=nghost+1,my-nghost), & 
        i=nghost+1,mx-nghost)
      close(19)
      endif

! Write out entropy field for Dx
      if(mvar.ge.5) then
      iss=5
      open(unit=20, file='./tmp/ss_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(20) ((real(f(i,j,isl,iss)),j=nghost+1,my-nghost), & 
        i=nghost+1,mx-nghost)
      close(20)
      endif

! Create .dx scripts describing the data in Data Explorer format
      open(unit=21, file='./tmp/uuu_sl'//char(nfile+48)//'.dx')
      write(21,*) 'object 1 class gridpositions counts ', mx-2*nghost, &
        my-2*nghost
      write(21,*) 'origin  ', x00, y00 
      write(21,*) 'delta  ', dx,'  0   '
      write(21,*) 'delta  ', '  0  ',dy
      write(21,*) '  '
      write(21,*) 'object 2 class gridconnections counts ', mx-2*nghost, &
        my-2*nghost
      write(21,*) '  '
      write(21,*) 'object 3 class array type float rank 1 shape 2'
      write(21,*) 'items ', (mx-2*nghost)*(my-2*nghost), ' lsb binary'
      write(21,*) 'data file "uuu_dx_sl'//char(nfile+48)//'.bin", 4'
      write(21,*) '  '
      write(21,*) 'attribute "dep" string "positions" '
      write(21,*) '  '
      write(21,*) 'object "regular positions regular connections" class field'
      write(21,*) 'component "positions" value 1'
      write(21,*) 'component "connections" value 2'
      write(21,*) 'component "data" value 3'
      write(21,*) 'end'
      close(21)

      if(mvar.ge.4) then         
      open(unit=22, file='./tmp/lnrho_sl'//char(nfile+48)//'.dx')
      write(22,*) 'object 1 class gridpositions counts ', mx-2*nghost,  &
        my-2*nghost
      write(22,*) 'origin  ', x00, y00 
      write(22,*) 'delta  ', dx,'  0   '
      write(22,*) 'delta  ', '  0  ',dy
      write(22,*) '  '
      write(22,*) 'object 2 class gridconnections counts ', mx-2*nghost,  &
        my-2*nghost
      write(22,*) '  '
      write(22,*) 'object 3 class array type float rank 0'
      write(22,*) 'items ', (mx-2*nghost)*(my-2*nghost), ' lsb binary'
      write(22,*) 'data file "lnrho_dx_sl'//char(nfile+48)//'.bin", 4'
      write(22,*) '  '
      write(22,*) 'attribute "dep" string "positions" '
      write(22,*) '  '
      write(22,*) 'object "regular positions regular connections" class field'
      write(22,*) 'component "positions" value 1'
      write(22,*) 'component "connections" value 2'
      write(22,*) 'component "data" value 3'
      write(22,*) 'end'
      close(22)
      endif

      if(mvar.ge.5) then
      open(unit=23, file='./tmp/ss_sl'//char(nfile+48)//'.dx')
      write(23,*) 'object 1 class gridpositions counts ', mx-2*nghost, &
        my-2*nghost
      write(23,*) 'origin  ', x00, y00 
      write(23,*) 'delta  ', dx,'  0   '
      write(23,*) 'delta  ', '  0  ',dy
      write(23,*) '  '
      write(23,*) 'object 2 class gridconnections counts ', mx-2*nghost, &
        my-2*nghost
      write(23,*) '  '
      write(23,*) 'object 3 class array type float rank 0'
      write(23,*) 'items ', (mx-2*nghost)*(my-2*nghost), ' lsb binary'
      write(23,*) 'data file "ss_dx_sl'//char(nfile+48)//'.bin", 4'
      write(23,*) '  '
      write(23,*) 'attribute "dep" string "positions" '
      write(23,*) '  '
      write(23,*) 'object "regular positions regular connections" class field'
      write(23,*) 'component "positions" value 1'
      write(23,*) 'component "connections" value 2'
      write(23,*) 'component "data" value 3'
      write(23,*) 'end'
      close(23)
      endif



      if(mvar.ge.6) then
! Process magnetic part if mvar>5: no kinematic dynamo yet, needs to be added
! Write out 2D slice of parallel vector potential for Dx
      iax=6
      iay=7
      iaz=8
      open(unit=18, file='./tmp/aaa_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(18) (((real(f(i,j,isl,n)),n=iax,iay),j=nghost+1,my-nghost), & 
        i=nghost+1,mx-nghost)
      close(18)
      
! Calculate magnetic field, using central differences.
      do i=2,mx-1
         do j=2,my-1
            do k=2,mz-1
               bbb(i,j,k,1)=(f(i,j+1,k,iaz)-f(i,j-1,k,iaz))/2./dy-  &
       (f(i,j,k+1,iay)-f(i,j,k-1,iay))/2./dz
               bbb(i,j,k,2)=(f(i,j,k+1,iax)-f(i,j,k-1,iax))/2./dz-  &
       (f(i+1,j,k,iaz)-f(i-1,j,k,iaz))/2./dx
               bbb(i,j,k,3)=(f(i+1,j,k,iay)-f(i-1,j,k,iay))/2./dx-  &
       (f(i,j+1,k,iax)-f(i,j-1,k,iax))/2./dy               
            enddo
         enddo
      enddo   

! Write out parallel magnetic field for Dx
      open(unit=18, file='./tmp/bbb_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(18) (((real(bbb(i,j,isl,n)),n=1,2),j=nghost+1,my-nghost), & 
        i=nghost+1,mx-nghost)
      close(18)

! Create  .dx scripts describing magnetic data 
      open(unit=24, file='./tmp/aaa_sl'//char(nfile+48)//'.dx')
      write(24,*) 'object 1 class gridpositions counts ', mx-2*nghost,  &
        my-2*nghost
      write(24,*) 'origin  ', x00, y00 
      write(24,*) 'delta  ', dx,'  0   '
      write(24,*) 'delta  ', '  0  ',dy
      write(24,*) '  '
      write(24,*) 'object 2 class gridconnections counts ', mx-2*nghost,  &
        my-2*nghost
      write(24,*) '  '
      write(24,*) 'object 3 class array type float rank 1 shape 2'
      write(24,*) 'items ', (mx-2*nghost)*(my-2*nghost), ' lsb binary'
      write(24,*) 'data file "aaa_dx_sl'//char(nfile+48)//'.bin", 4'
      write(24,*) '  '
      write(24,*) 'attribute "dep" string "positions" '
      write(24,*) '  '
      write(24,*) 'object "regular positions regular connections" class field'
      write(24,*) 'component "positions" value 1'
      write(24,*) 'component "connections" value 2'
      write(24,*) 'component "data" value 3'
      write(24,*) 'end'
      close(24)

      open(unit=25, file='./tmp/bbb_sl'//char(nfile+48)//'.dx')
      write(25,*) 'object 1 class gridpositions counts ', mx-2*nghost, &
        my-2*nghost
      write(25,*) 'origin  ', x00, y00 
      write(25,*) 'delta  ', dx,'  0   '
      write(25,*) 'delta  ', '  0  ',dy
      write(25,*) '  '
      write(25,*) 'object 2 class gridconnections counts ', mx-2*nghost, &
        my-2*nghost
      write(25,*) '  '
      write(25,*) 'object 3 class array type float rank 1 shape 2'
      write(25,*) 'items ', (mx-2*nghost)*(my-2*nghost), ' lsb binary'
      write(25,*) 'data file "bbb_dx_sl'//char(nfile+48)//'.bin", 4'
      write(25,*) '  '
      write(25,*) 'attribute "dep" string "positions" '
      write(25,*) '  '
      write(25,*) 'object "regular positions regular connections" class field'
      write(25,*) 'component "positions" value 1'
      write(25,*) 'component "connections" value 2'
      write(25,*) 'component "data" value 3'
      write(25,*) 'end'
      close(25)


      endif

!!!!!!!!
!!!  Case of yz-slice
!!!!!!!!
      case('yz')

! Write out parallel velocity field in the slice
      iux=1
      iuy=2
      iuz=3
      open(unit=18, file='./tmp/uuu_dx_sl'//char(nfile+48)//'.bin', &
        form='unformatted')
      write(18) (((real(f(isl,j,k,n)),n=iuy,iuz),k=nghost+1,mz-nghost), &
        j=nghost+1,my-nghost)
      close(18)

! Write out ln density field for Dx
      if(mvar.ge.4) then
      ilnrho=4
      open(unit=19, file='./tmp/lnrho_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(19) ((real(f(isl,j,k,ilnrho)),k=nghost+1,mz-nghost), & 
        j=nghost+1,my-nghost)
      close(19)
      endif

! Write out entropy field for Dx
      if(mvar.ge.5) then
      iss=5
      open(unit=20, file='./tmp/ss_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(20) ((real(f(isl,j,k,iss)),k=nghost+1,mz-nghost), & 
        j=nghost+1,my-nghost)
      close(20)
      endif

! Create .dx scripts describing the data in Data Explorer format
      open(unit=21, file='./tmp/uuu_sl'//char(nfile+48)//'.dx')
      write(21,*) 'object 1 class gridpositions counts ', my-2*nghost, &
        mz-2*nghost
      write(21,*) 'origin  ', y00, z00 
      write(21,*) 'delta  ', dy,'  0   '
      write(21,*) 'delta  ', '  0  ',dz
      write(21,*) '  '
      write(21,*) 'object 2 class gridconnections counts ', my-2*nghost, &
        mz-2*nghost
      write(21,*) '  '
      write(21,*) 'object 3 class array type float rank 1 shape 2'
      write(21,*) 'items ', (my-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(21,*) 'data file "uuu_dx_sl'//char(nfile+48)//'.bin", 4'
      write(21,*) '  '
      write(21,*) 'attribute "dep" string "positions" '
      write(21,*) '  '
      write(21,*) 'object "regular positions regular connections" class field'
      write(21,*) 'component "positions" value 1'
      write(21,*) 'component "connections" value 2'
      write(21,*) 'component "data" value 3'
      write(21,*) 'end'
      close(21)

      if(mvar.ge.4) then         
      open(unit=22, file='./tmp/lnrho_sl'//char(nfile+48)//'.dx')
      write(22,*) 'object 1 class gridpositions counts ', my-2*nghost,  &
        mz-2*nghost
      write(22,*) 'origin  ', y00, z00 
      write(22,*) 'delta  ', dy,'  0   '
      write(22,*) 'delta  ', '  0  ',dz
      write(22,*) '  '
      write(22,*) 'object 2 class gridconnections counts ', my-2*nghost,  &
        mz-2*nghost
      write(22,*) '  '
      write(22,*) 'object 3 class array type float rank 0'
      write(22,*) 'items ', (my-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(22,*) 'data file "lnrho_dx_sl'//char(nfile+48)//'.bin", 4'
      write(22,*) '  '
      write(22,*) 'attribute "dep" string "positions" '
      write(22,*) '  '
      write(22,*) 'object "regular positions regular connections" class field'
      write(22,*) 'component "positions" value 1'
      write(22,*) 'component "connections" value 2'
      write(22,*) 'component "data" value 3'
      write(22,*) 'end'
      close(22)
      endif

      if(mvar.ge.5) then
      open(unit=23, file='./tmp/ss_sl'//char(nfile+48)//'.dx')
      write(23,*) 'object 1 class gridpositions counts ', my-2*nghost, &
        mz-2*nghost
      write(23,*) 'origin  ', y00, z00 
      write(23,*) 'delta  ', dy,'  0   '
      write(23,*) 'delta  ', '  0  ',dz
      write(23,*) '  '
      write(23,*) 'object 2 class gridconnections counts ', my-2*nghost, &
        mz-2*nghost
      write(23,*) '  '
      write(23,*) 'object 3 class array type float rank 0'
      write(23,*) 'items ', (my-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(23,*) 'data file "ss_dx_sl'//char(nfile+48)//'.bin", 4'
      write(23,*) '  '
      write(23,*) 'attribute "dep" string "positions" '
      write(23,*) '  '
      write(23,*) 'object "regular positions regular connections" class field'
      write(23,*) 'component "positions" value 1'
      write(23,*) 'component "connections" value 2'
      write(23,*) 'component "data" value 3'
      write(23,*) 'end'
      close(23)
      endif



      if(mvar.ge.6) then
! Process magnetic part if mvar>5: no kinematic dynamo yet, needs to be added
! Write out 2D slice of parallel vector potential for Dx
      iax=6
      iay=7
      iaz=8
      open(unit=18, file='./tmp/aaa_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(18) (((real(f(isl,j,k,n)),n=iay,iaz),k=nghost+1,mz-nghost), & 
        j=nghost+1,my-nghost)
      close(18)

! Calculate magnetic field, using central differences.
      do i=2,mx-1
         do j=2,my-1
            do k=2,mz-1
               bbb(i,j,k,1)=(f(i,j+1,k,iaz)-f(i,j-1,k,iaz))/2./dy-  &
       (f(i,j,k+1,iay)-f(i,j,k-1,iay))/2./dz
               bbb(i,j,k,2)=(f(i,j,k+1,iax)-f(i,j,k-1,iax))/2./dz-  &
       (f(i+1,j,k,iaz)-f(i-1,j,k,iaz))/2./dx
               bbb(i,j,k,3)=(f(i+1,j,k,iay)-f(i-1,j,k,iay))/2./dx-  &
       (f(i,j+1,k,iax)-f(i,j-1,k,iax))/2./dy               
            enddo
         enddo
      enddo   

! Write out parallel magnetic field for Dx
      open(unit=18, file='./tmp/bbb_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(18) (((real(bbb(isl,j,k,n)),n=2,3),k=nghost+1,mz-nghost), & 
        j=nghost+1,my-nghost)
      close(18)

! Create  .dx scripts describing magnetic data 
      open(unit=24, file='./tmp/aaa_sl'//char(nfile+48)//'.dx')
      write(24,*) 'object 1 class gridpositions counts ', my-2*nghost,  &
        mz-2*nghost
      write(24,*) 'origin  ', y00, z00 
      write(24,*) 'delta  ', dy,'  0   '
      write(24,*) 'delta  ', '  0  ',dz
      write(24,*) '  '
      write(24,*) 'object 2 class gridconnections counts ', my-2*nghost,  &
        mz-2*nghost
      write(24,*) '  '
      write(24,*) 'object 3 class array type float rank 1 shape 2'
      write(24,*) 'items ', (my-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(24,*) 'data file "aaa_dx_sl'//char(nfile+48)//'.bin", 4'
      write(24,*) '  '
      write(24,*) 'attribute "dep" string "positions" '
      write(24,*) '  '
      write(24,*) 'object "regular positions regular connections" class field'
      write(24,*) 'component "positions" value 1'
      write(24,*) 'component "connections" value 2'
      write(24,*) 'component "data" value 3'
      write(24,*) 'end'
      close(24)

      open(unit=25, file='./tmp/bbb_sl'//char(nfile+48)//'.dx')
      write(25,*) 'object 1 class gridpositions counts ', my-2*nghost, &
        mz-2*nghost
      write(25,*) 'origin  ', y00, z00 
      write(25,*) 'delta  ', dy,'  0   '
      write(25,*) 'delta  ', '  0  ',dz
      write(25,*) '  '
      write(25,*) 'object 2 class gridconnections counts ', my-2*nghost, &
        mz-2*nghost
      write(25,*) '  '
      write(25,*) 'object 3 class array type float rank 1 shape 2'
      write(25,*) 'items ', (my-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(25,*) 'data file "bbb_dx_sl'//char(nfile+48)//'.bin", 4'
      write(25,*) '  '
      write(25,*) 'attribute "dep" string "positions" '
      write(25,*) '  '
      write(25,*) 'object "regular positions regular connections" class field'
      write(25,*) 'component "positions" value 1'
      write(25,*) 'component "connections" value 2'
      write(25,*) 'component "data" value 3'
      write(25,*) 'end'
      close(25)


      endif


!!!!!!!!
!!!  Case of xz-slice
!!!!!!!!
      case('xz')

! Write out parallel velocity field in the slice
      iux=1
      iuy=2
      iuz=3
      open(unit=18, file='./tmp/uuu_dx_sl'//char(nfile+48)//'.bin', &
        form='unformatted')
      write(18) (((real(f(i,isl,k,n)),n=iux,iuz,2),k=nghost+1,mz-nghost), &
        i=nghost+1,mx-nghost)
      close(18)

! Write out ln density field for Dx
      if(mvar.ge.4) then
      ilnrho=4
      open(unit=19, file='./tmp/lnrho_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(19) ((real(f(i,isl,k,ilnrho)),k=nghost+1,mz-nghost), & 
        i=nghost+1,mx-nghost)
      close(19)
      endif

! Write out entropy field for Dx
      if(mvar.ge.5) then
      iss=5
      open(unit=20, file='./tmp/ss_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(20) ((real(f(i,isl,k,iss)),k=nghost+1,mz-nghost), & 
        i=nghost+1,mx-nghost)
      close(20)
      endif

! Create .dx scripts describing the data in Data Explorer format
      open(unit=21, file='./tmp/uuu_sl'//char(nfile+48)//'.dx')
      write(21,*) 'object 1 class gridpositions counts ', mx-2*nghost, &
        mz-2*nghost
      write(21,*) 'origin  ', x00, z00 
      write(21,*) 'delta  ', dx,'  0   '
      write(21,*) 'delta  ', '  0  ',dz
      write(21,*) '  '
      write(21,*) 'object 2 class gridconnections counts ', mx-2*nghost, &
        mz-2*nghost
      write(21,*) '  '
      write(21,*) 'object 3 class array type float rank 1 shape 2'
      write(21,*) 'items ', (mx-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(21,*) 'data file "uuu_dx_sl'//char(nfile+48)//'.bin", 4'
      write(21,*) '  '
      write(21,*) 'attribute "dep" string "positions" '
      write(21,*) '  '
      write(21,*) 'object "regular positions regular connections" class field'
      write(21,*) 'component "positions" value 1'
      write(21,*) 'component "connections" value 2'
      write(21,*) 'component "data" value 3'
      write(21,*) 'end'
      close(21)

      if(mvar.ge.4) then         
      open(unit=22, file='./tmp/lnrho_sl'//char(nfile+48)//'.dx')
      write(22,*) 'object 1 class gridpositions counts ', mx-2*nghost,  &
        mz-2*nghost
      write(22,*) 'origin  ', x00, z00 
      write(22,*) 'delta  ', dx,'  0   '
      write(22,*) 'delta  ', '  0  ',dz
      write(22,*) '  '
      write(22,*) 'object 2 class gridconnections counts ', mx-2*nghost,  &
        mz-2*nghost
      write(22,*) '  '
      write(22,*) 'object 3 class array type float rank 0'
      write(22,*) 'items ', (mx-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(22,*) 'data file "lnrho_dx_sl'//char(nfile+48)//'.bin", 4'
      write(22,*) '  '
      write(22,*) 'attribute "dep" string "positions" '
      write(22,*) '  '
      write(22,*) 'object "regular positions regular connections" class field'
      write(22,*) 'component "positions" value 1'
      write(22,*) 'component "connections" value 2'
      write(22,*) 'component "data" value 3'
      write(22,*) 'end'
      close(22)
      endif

      if(mvar.ge.5) then
      open(unit=23, file='./tmp/ss_sl'//char(nfile+48)//'.dx')
      write(23,*) 'object 1 class gridpositions counts ', mx-2*nghost, &
        mz-2*nghost
      write(23,*) 'origin  ', x00, z00 
      write(23,*) 'delta  ', dx,'  0   '
      write(23,*) 'delta  ', '  0  ',dz
      write(23,*) '  '
      write(23,*) 'object 2 class gridconnections counts ', mx-2*nghost, &
        mz-2*nghost
      write(23,*) '  '
      write(23,*) 'object 3 class array type float rank 0'
      write(23,*) 'items ', (mx-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(23,*) 'data file "ss_dx_sl'//char(nfile+48)//'.bin", 4'
      write(23,*) '  '
      write(23,*) 'attribute "dep" string "positions" '
      write(23,*) '  '
      write(23,*) 'object "regular positions regular connections" class field'
      write(23,*) 'component "positions" value 1'
      write(23,*) 'component "connections" value 2'
      write(23,*) 'component "data" value 3'
      write(23,*) 'end'
      close(23)
      endif


      if(mvar.ge.6) then
! Process magnetic part if mvar>5: no kinematic dynamo yet, needs to be added
! Write out 2D slice of parallel vector potential for Dx
      iax=6
      iay=7
      iaz=8
      open(unit=18, file='./tmp/aaa_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(18) (((real(f(i,isl,k,n)),n=iax,iaz,2),k=nghost+1,mz-nghost), & 
        i=nghost+1,mx-nghost)
      close(18)
      
! Calculate magnetic field, using central differences.
      do i=2,mx-1
         do j=2,my-1
            do k=2,mz-1
               bbb(i,j,k,1)=(f(i,j+1,k,iaz)-f(i,j-1,k,iaz))/2./dy-  &
       (f(i,j,k+1,iay)-f(i,j,k-1,iay))/2./dz
               bbb(i,j,k,2)=(f(i,j,k+1,iax)-f(i,j,k-1,iax))/2./dz-  &
       (f(i+1,j,k,iaz)-f(i-1,j,k,iaz))/2./dx
               bbb(i,j,k,3)=(f(i+1,j,k,iay)-f(i-1,j,k,iay))/2./dx-  &
       (f(i,j+1,k,iax)-f(i,j-1,k,iax))/2./dy               
            enddo
         enddo
      enddo   

! Write out parallel magnetic field for Dx
      open(unit=18, file='./tmp/bbb_dx_sl'//char(nfile+48)//'.bin', & 
        form='unformatted')
      write(18) (((real(bbb(i,isl,k,n)),n=1,3,2),k=nghost+1,mz-nghost), & 
        i=nghost+1,mx-nghost)
      close(18)

! Create  .dx scripts describing magnetic data 
      open(unit=24, file='./tmp/aaa_sl'//char(nfile+48)//'.dx')
      write(24,*) 'object 1 class gridpositions counts ', mx-2*nghost,  &
        mz-2*nghost
      write(24,*) 'origin  ', x00, z00 
      write(24,*) 'delta  ', dx,'  0   '
      write(24,*) 'delta  ', '  0  ',dz
      write(24,*) '  '
      write(24,*) 'object 2 class gridconnections counts ', mx-2*nghost,  &
        mz-2*nghost
      write(24,*) '  '
      write(24,*) 'object 3 class array type float rank 1 shape 2'
      write(24,*) 'items ', (mx-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(24,*) 'data file "aaa_dx_sl'//char(nfile+48)//'.bin", 4'
      write(24,*) '  '
      write(24,*) 'attribute "dep" string "positions" '
      write(24,*) '  '
      write(24,*) 'object "regular positions regular connections" class field'
      write(24,*) 'component "positions" value 1'
      write(24,*) 'component "connections" value 2'
      write(24,*) 'component "data" value 3'
      write(24,*) 'end'
      close(24)

      open(unit=25, file='./tmp/bbb_sl'//char(nfile+48)//'.dx')
      write(25,*) 'object 1 class gridpositions counts ', mx-2*nghost, &
        mz-2*nghost
      write(25,*) 'origin  ', x00, z00 
      write(25,*) 'delta  ', dx,'  0   '
      write(25,*) 'delta  ', '  0  ',dz
      write(25,*) '  '
      write(25,*) 'object 2 class gridconnections counts ', mx-2*nghost, &
        mz-2*nghost
      write(25,*) '  '
      write(25,*) 'object 3 class array type float rank 1 shape 2'
      write(25,*) 'items ', (mx-2*nghost)*(mz-2*nghost), ' lsb binary'
      write(25,*) 'data file "bbb_dx_sl'//char(nfile+48)//'.bin", 4'
      write(25,*) '  '
      write(25,*) 'attribute "dep" string "positions" '
      write(25,*) '  '
      write(25,*) 'object "regular positions regular connections" class field'
      write(25,*) 'component "positions" value 1'
      write(25,*) 'component "connections" value 2'
      write(25,*) 'component "data" value 3'
      write(25,*) 'end'
      close(25)


      endif


!!!!! Finishing choosing the slices
      case default
        !
        !  Catch unknown values
        !
        print*,'No such value for sltype: ', trim(sltype)
        stop
      endselect     


      stop
      end