! $Id: grid2grid.f90,v 1.11 2007/12/03 06:05:45 brandenb Exp $

!-------------------------------------------------------------------------------
PROGRAM grid2grid
!
! Regrid 3D-data, to change resolution and possibly staggering.
! This version has been hardwired for pencil code usage only.
!
  USE params
  USE io
  implicit none
  integer m,mxin,myin,mzin,mx,my,mz,mv,var,iz,vx(2),i
  integer gx,gy,gz,nx,ny,nz
  real, allocatable, dimension(:,:,:):: vin,vout,bx,by,bz
  real, allocatable, dimension(:,:,:):: lr,ux,uy,uz
  real, allocatable, dimension(:,:,:):: lr0,ux0,uy0,uz0,bx0,by0,bz0
  real, allocatable, dimension(:):: x,y,z
  real, allocatable, dimension(:):: xx,yy,zz
  character(len=mfile) file,infile,outfile,arg
  real offset_in,sx,sy,sz,average
  real, dimension(3):: off_in,off_out
  real :: t=0.,dx,dy,dz
  logical do_smooth,do_sink,do_seq
  namelist /in/file,vx,offset,lio,do_exp,do_seq,do_sink
  namelist /out/file,m,mv,offset,do_smooth,do_exp,do_seq
  character(len=80),save:: id="$Id: grid2grid.f90,v 1.11 2007/12/03 06:05:45 brandenb Exp $"
  character(len=5) :: chproc
! integer :: iproc, ipy, ipz, nprocy=8, nprocz=64
  integer :: iproc, ipy, ipz, nprocy=8, nprocz=4
! integer :: iproc, ipy, ipz, nprocy=4, nprocz=2
  integer :: m1, m2, n1, n2
  logical :: lread=.true.
!-------------------------------------------------------------------------------
  print'(1x,a)',hl,id,hl

  call get_arg(arg)                                                             ! optional param file name
  open(txt_unit,file=trim(arg),form='formatted',status='old')                   ! open param file

  mv=4                                                                          ! number of variables (hd=4, mhd=7)
  vx=(/2,5/)                                                                    ! indices to vector x-components
  offset=-0.5                                                                   ! staggering (-0.5 for left face)
  do_smooth=.true.                                                              ! mimic actual resolution
  do_exp=.false.                                                                ! input file has ln(rho)
  do_sink=.false.                                                               ! add sink particles
  infile='test.dat'

  do_seq=.false.
  write(*,*) ' '; read(txt_unit,in)                                             ! infile parameters
  infile=file                                                                   ! remember
  do_seq_read=do_seq
  if (index(infile,'.seq') > 0) do_seq_read=.true.                              ! *.seq files are sequential
  print *,'index:',index(infile,'.seq')
  write(*,in); write(*,*) ' '

  call openr_grid(in_unit,infile,mx,my,mz,mv,sx,sy,sz)                          ! open input file
  mxin=mx; myin=my; mzin=mz; offset_in=offset

  write(*,*) ' '
  outfile='grid.dat'; m=50                                                      ! outfile parameters
  do_seq=.false.
  offset=0.0                                                                    ! staggering (0.0 for cell center)
  read(txt_unit,out); write(*,out)
  outfile=file; mx=m; my=m; mz=m
  do_seq_write=do_seq
  if (index(outfile,'.seq') > 0) do_seq_write=.true.                            ! *.seq files are sequential
  write(*,*) ' '
!
!  define size of target grid
!  gx,gy,gz is the mesh size with ghost zones
!  nx,ny,nz is the mesh size without ghost zones
!
  nx=mx; ny=my/nprocy; nz=mz/nprocz
  gx=6+nx; gy=6+ny; gz=6+nz
  print*,'nx,ny,nz=',nx,ny,nz
  print*,'gx,gy,gz=',gx,gy,gz
!
!  x,y,z is the full mesh, xx,yy,zz is the submesh per processor
!
  dx=1./mx; dy=1./my; dz=1./mz
  allocate (x(mx),y(my),z(mz))
  allocate (xx(gx),yy(gy),zz(gz))
!
!  mesh values of full mesh
!
  do i=1,mx; x(i)=dx*(i-1); enddo
  do i=1,my; y(i)=dy*(i-1); enddo
  do i=1,mz; z(i)=dz*(i-1); enddo

! call openw_grid(out_unit,outfile,mx,my,mz,mv,sx,sy,sz)                        ! open output file

  allocate (vin(mxin,myin,mzin), vout(mx,my,mz))
  allocate (lr0(mx,my,mz),ux0(mx,my,mz),uy0(mx,my,mz),uz0(mx,my,mz))
  allocate (lr(gx,gy,gz),ux(gx,gy,gz),uy(gx,gy,gz),uz(gx,gy,gz))
  allocate (bx0(mx,my,mz),by0(mx,my,mz),bz0(mx,my,mz))
  allocate (bx(gx,gy,gz),by(gx,gy,gz),bz(gx,gy,gz))
!
!  read in everything
!  loop over all variables
!
print*,'mv=',mv
if (lread) then
  do var=1,mv
    call read_grid(in_unit,mxin,myin,mzin,var,vin)
    print '(i3,1x,a,g12.4)',var,'interpolating, average =',average(mxin,myin,mzin,vin)
    off_in = 0.
    off_out = 0.
    do i=1,2
      if (var == vx(i)  ) off_in=(/offset_in,0.,0./)
      if (var == vx(i)+1) off_in=(/.0,offset_in,.0/)
      if (var == vx(i)+2) off_in=(/.0,.0,offset_in/)
      if (var == vx(i)  ) off_out=(/offset,.0,.0/)
      if (var == vx(i)+1) off_out=(/.0,offset,.0/)
      if (var == vx(i)+2) off_out=(/.0,.0,offset/)
    end do
    call interpolate (vin,mxin,myin,mzin,vout,mx,my,mz,off_in,off_out)
    if (do_smooth) call smooth3t(mx,my,mz,vout)
    print '(i3,1x,a,g12.4)',var,'average =',average(mx,my,mz,vout)

    if (var < vx(2)) then
      if (var == 1) then
        lr0=vout
      elseif (var == 2) then
        ux0=vout
      elseif (var == 3) then
        uy0=vout
      elseif (var == 4) then
        uz0=vout
      endif
    else if (var == vx(2)) then
      bx0=vout
    else if (var == vx(2)+1) then
      by0=vout
    else if (var == vx(2)+2) then
      bz0=vout
    end if
  end do
  deallocate (vin,vout)
  close(txt_unit)
  close(in_unit)
  close(out_unit)
endif
!
!  loop over all processors to write out the data
!
  print*,'loop over all processors to write out the data'
  iproc=0
  do ipz=0,nprocz-1
  do ipy=0,nprocy-1
!
!  processor names
!
    call chn(iproc,chproc,'directory_names')
    print*,'data/proc'//trim(chproc)//'/var.dat'
    open(8,file='data/proc'//trim(chproc)//'/var.dat',form='unformatted')
!
!  first and last mesh point of full mesh
!
!-- l1=ipx*nx+1; l2=l1+nx-1
    m1=ipy*ny+1; m2=m1+ny-1
    n1=ipz*nz+1; n2=n1+nz-1
!
!  print range
!
    print*,'m1,m2=',m1,m2
    print*,'n1,n2=',n1,n2
!
!  loop over all variables
!
    xx(4:nx+3)=x
    yy(4:ny+3)=y(m1:m2)
    print*,'yy=',yy(4),yy(ny+3)
    zz(4:nz+3)=y(n1:n2)
    print*,'zz=',zz(4),zz(nz+3)
    ux(4:nx+3,4:ny+3,4:nz+3)=ux0(:,m1:m2,n1:n2)
    uy(4:nx+3,4:ny+3,4:nz+3)=uy0(:,m1:m2,n1:n2)
    uz(4:nx+3,4:ny+3,4:nz+3)=uz0(:,m1:m2,n1:n2)
    lr(4:nx+3,4:ny+3,4:nz+3)=lr0(:,m1:m2,n1:n2)
    bx(4:nx+3,4:ny+3,4:nz+3)=bx0(:,m1:m2,n1:n2)
    by(4:nx+3,4:ny+3,4:nz+3)=by0(:,m1:m2,n1:n2)
    bz(4:nx+3,4:ny+3,4:nz+3)=bz0(:,m1:m2,n1:n2)
    write(8) ux,uy,uz,lr,bx,by,bz
    write(8) t,xx,yy,zz,dx,dy,dz
    close(8)
    iproc=iproc+1
  enddo
  enddo
  deallocate (lr,ux,uy,uz)
  deallocate (lr0,ux0,uy0,uz0)
!
END
!***********************************************************************
    subroutine chn(n,ch,label)
!
!  make a character out of a number
!  take care of numbers that have less than 4 digits
!  30-sep-97/axel: coded
!
      character (len=5) :: ch
      character (len=*), optional :: label
      integer :: n
!
      intent(in) :: n
!
      ch='     '
      if (n<0) stop 'chn: lt1'
      if (n<10) then
        write(ch(1:1),'(i1)') n
      elseif (n<100) then
        write(ch(1:2),'(i2)') n
      elseif (n<1000) then
        write(ch(1:3),'(i3)') n
      elseif (n<10000) then
        write(ch(1:4),'(i4)') n
      elseif (n<100000) then
        write(ch(1:5),'(i5)') n
      else
        if (present(label)) print*, 'CHN: <', label, '>'
        print*,'CHN: n=',n
        stop "CHN: n too large"
      endif
!
    endsubroutine chn