program todx_magn
      implicit none
! Converts var.dat snapshot to the Dx readable binary data files.
! The data are read in again by .dx routines on regular rectangular
! mesh. The todx_magn.f90 should be copied to the 
! directory of the run. Compile todx_magn.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.

! 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,ilncc
      integer :: i_t,i_it,i_dt,i_dtc,nname,i_u2m,i_um2,i_o2m
      integer :: i_oum,i_urms,i_umax,i_orms,i_omax,i_uxm
      integer :: i_uym,i_uzm,iuu,i_eth,i_ekin,i_rhom
      integer :: i_ssm,i_abm,i_jbm,i_b2m,i_bm2,i_j2m,i_jm2
      integer :: iaa,nnamez,i_bxmz,i_bymz,i_bzmz,i_bmx,i_bmy
      integer :: i_bmz,nnamexy,i_bxmxy,i_bymxy,i_bzmxy,i_frms
      integer :: i_fmax,i_Erad_rms,i_Erad_max,i_Egas_rms
      integer :: i_Egas_max,ie,ifx,ify,ifz,i_ccm,i_ccmax,i_rhoccm

! Read file of integer parameters, index.pro
      include 'tmp/index.pro'
      
      print*,'iux=',iux,' iuy=',iuy,' iuz=',iuz,' ilnrho=',  &
       ilnrho,' iss=',iss,' iax=',iax,' iay=', iay,' iaz=', &
       iaz,' ilncc=',ilncc

! 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
      
! Write out velocity field for Dx
      if(iux.ne.0) then
         if(iuy-iux.ne.1.or.iuz-iuy.ne.1) then
            print*,'iux,iuy,iuz are not in order, exiting'
            stop
         else   
      open(unit=18, file='./tmp/uuu_dx.bin', form='unformatted')
      write(18) ((((real(f(i,j,k,n)),n=iux,iuz),k=nghost+1,mz-nghost), &
       j=nghost+1,my-nghost),i=nghost+1,mx-nghost)
      close(18)

      open(unit=21, file='./tmp/uuu.dx')
      write(21,*) 'object 1 class gridpositions counts ', mx-2*nghost, &
        my-2*nghost, mz-2*nghost
      write(21,*) 'origin  ', x00, y00, z00 
      write(21,*) 'delta  ', dx,'  0   ', '   0    '
      write(21,*) 'delta  ', '  0  ',dy,'  0  '
      write(21,*) 'delta','  0  ','  0  ',dz
      write(21,*) '  '
      write(21,*) 'object 2 class gridconnections counts ', mx-2*nghost, &
        my-2*nghost, mz-2*nghost
      write(21,*) '  '
      write(21,*) 'object 3 class array type float rank 1 shape 3'
      write(21,*) 'items ', (mx-2*nghost)*(my-2*nghost)*(mz-2*nghost),  &
        ' lsb binary'
      write(21,*) 'data file "uuu_dx.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)

      print*,'velocity field has been written out'

      endif
      endif


! Write out ln density field for Dx
      if(ilnrho.ne.0) then

      open(unit=19, file='./tmp/lnrho_dx.bin', form='unformatted')
      write(19) (((real(f(i,j,k,ilnrho)),k=nghost+1,mz-nghost), &
       j=nghost+1,my-nghost),i=nghost+1,mx-nghost)
      close(19)

      open(unit=22, file='./tmp/lnrho.dx')
      write(22,*) 'object 1 class gridpositions counts ', mx-2*nghost,  &
        my-2*nghost, mz-2*nghost
      write(22,*) 'origin  ', x00, y00, z00 
      write(22,*) 'delta  ', dx,'  0   ', '   0    '
      write(22,*) 'delta  ', '  0  ',dy,'  0  '
      write(22,*) 'delta','  0  ','  0  ',dz
      write(22,*) '  '
      write(22,*) 'object 2 class gridconnections counts ', mx-2*nghost,  &
        my-2*nghost, mz-2*nghost
      write(22,*) '  '
      write(22,*) 'object 3 class array type float rank 0'
      write(22,*) 'items ', (mx-2*nghost)*(my-2*nghost)*(mz-2*nghost),  &
        ' lsb binary'
      write(22,*) 'data file "lnrho_dx.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)

      print*,'density field has been written out'

      endif


! Write out entropy field for Dx
      if(iss.ne.0) then

      open(unit=20, file='./tmp/ss_dx.bin', form='unformatted')
      write(20) (((real(f(i,j,k,iss)),k=nghost+1,mz-nghost),  &
       j=nghost+1,my-nghost),i=nghost+1,mx-nghost)
      close(20)

      open(unit=23, file='./tmp/ss.dx')
      write(23,*) 'object 1 class gridpositions counts ', mx-2*nghost, &
        my-2*nghost, mz-2*nghost
      write(23,*) 'origin  ', x00, y00, z00 
      write(23,*) 'delta  ', dx,'  0   ', '   0    '
      write(23,*) 'delta  ', '  0  ',dy,'  0  '
      write(23,*) 'delta','  0  ','  0  ',dz
      write(23,*) '  '
      write(23,*) 'object 2 class gridconnections counts ', mx-2*nghost, &
        my-2*nghost, mz-2*nghost
      write(23,*) '  '
      write(23,*) 'object 3 class array type float rank 0'
      write(23,*) 'items ', (mx-2*nghost)*(my-2*nghost)*(mz-2*nghost),  &
        ' lsb binary'
      write(23,*) 'data file "ss_dx.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)

      print*,'entropy field has been written out'

      endif


! Write out vector potential and magnetic field for Dx
      if(iax.ne.0) then
         if(iay-iax.ne.1.or.iaz-iay.ne.1) then
            print*,'iax,iay,iaz are not in order, exiting'
            stop
         else   

      open(unit=18, file='./tmp/aaa_dx.bin', form='unformatted')
      write(18) ((((real(f(i,j,k,n)),n=iax,iaz),k=nghost+1,mz-nghost), &
        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 magnetic field for Dx
      open(unit=18, file='./tmp/bbb_dx.bin', form='unformatted')
      write(18) ((((real(bbb(i,j,k,n)),n=1,3),k=nghost+1,mz-nghost), &
       j=nghost+1,my-nghost),i=nghost+1,mx-nghost)
      close(18)

! Create  .dx scripts describing magnetic data 
      open(unit=24, file='./tmp/aaa.dx')
      write(24,*) 'object 1 class gridpositions counts ', mx-2*nghost,  &
        my-2*nghost, mz-2*nghost
      write(24,*) 'origin  ', x00, y00, z00 
      write(24,*) 'delta  ', dx,'  0   ', '   0    '
      write(24,*) 'delta  ', '  0  ',dy,'  0  '
      write(24,*) 'delta','  0  ','  0  ',dz
      write(24,*) '  '
      write(24,*) 'object 2 class gridconnections counts ', mx-2*nghost,  &
        my-2*nghost, mz-2*nghost
      write(24,*) '  '
      write(24,*) 'object 3 class array type float rank 1 shape 3'
      write(24,*) 'items ', (mx-2*nghost)*(my-2*nghost)*(mz-2*nghost),  &
        ' lsb binary'
      write(24,*) 'data file "aaa_dx.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.dx')
      write(25,*) 'object 1 class gridpositions counts ', mx-2*nghost, &
        my-2*nghost, mz-2*nghost
      write(25,*) 'origin  ', x00, y00, z00 
      write(25,*) 'delta  ', dx,'  0   ', '   0    '
      write(25,*) 'delta  ', '  0  ',dy,'  0  '
      write(25,*) 'delta','  0  ','  0  ',dz
      write(25,*) '  '
      write(25,*) 'object 2 class gridconnections counts ', mx-2*nghost, &
        my-2*nghost, mz-2*nghost
      write(25,*) '  '
      write(25,*) 'object 3 class array type float rank 1 shape 3'
      write(25,*) 'items ', (mx-2*nghost)*(my-2*nghost)*(mz-2*nghost),  &
        ' lsb binary'
      write(25,*) 'data file "bbb_dx.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)

      print*,'vector potential and magnetic field have been written out'

      endif
      endif

! Write out passive scalar field for Dx
      if(ilncc.ne.0) then

      open(unit=19, file='./tmp/lncc_dx.bin', form='unformatted')
      write(19) (((real(f(i,j,k,ilncc)),k=nghost+1,mz-nghost), &
       j=nghost+1,my-nghost),i=nghost+1,mx-nghost)
      close(19)

      open(unit=22, file='./tmp/lncc.dx')
      write(22,*) 'object 1 class gridpositions counts ', mx-2*nghost,  &
        my-2*nghost, mz-2*nghost
      write(22,*) 'origin  ', x00, y00, z00 
      write(22,*) 'delta  ', dx,'  0   ', '   0    '
      write(22,*) 'delta  ', '  0  ',dy,'  0  '
      write(22,*) 'delta','  0  ','  0  ',dz
      write(22,*) '  '
      write(22,*) 'object 2 class gridconnections counts ', mx-2*nghost,  &
        my-2*nghost, mz-2*nghost
      write(22,*) '  '
      write(22,*) 'object 3 class array type float rank 0'
      write(22,*) 'items ', (mx-2*nghost)*(my-2*nghost)*(mz-2*nghost),  &
        ' lsb binary'
      write(22,*) 'data file "lncc_dx.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)

      print*,'passive scalar field has been written out'

      endif

      stop
      end