! $Id$
!
!  This module produces slices for animation purposes.
!
module Slices
!
! 16-nov-11/MR: I/O error handling generally introduced
!
  use Cdata
  use Messages
  use Mpicomm, only: mpiallreduce_or
  use Sub, only: xlocation, zlocation, update_snaptime, read_snaptime, position
!
  implicit none
!
  public :: wvid, wvid_prepare, setup_slices
!
  real, public :: tvid
  integer, public :: nvid
!
  private
!
  real :: tslice=0.
  real, target, dimension(:,:), allocatable :: slice_xy,slice_xy2,slice_xy3,slice_xy4
  real, target, dimension(:,:), allocatable :: slice_xz,slice_xz2,slice_yz,slice_r
!
contains
!***********************************************************************
    subroutine wvid_prepare
!
!  Prepare lvideo for writing slices into video file
!  This is useful for visualization of scalar field (or one component
!  of a vector field) on the periphery of a box.
!  Can be visualized in idl using rvid_box.pro
!
!  20-oct-97/axel: coded
!  08-oct-02/tony: increased size of file to handle datadir//'/tvid.dat'
!  13-nov-02/axel: added more fields
!  18-mar-03/axel: added dust velocity
!
      logical, save :: lfirst_call=.true.
      character (len=fnlen) :: file
!
!  Output vid-data in 'dvid' time intervals
!
      file = trim(datadir)//'/tvid.dat'
      if (lfirst_call) then
        call read_snaptime(file,tvid,nvid,dvid,t)
        lfirst_call=.false.
      endif
!
!  This routine sets lvideo=T whenever its time to write a slice
!
      call update_snaptime(file,tvid,nvid,dvid,t,lvideo)
!
!  Save current time so that the time that is written out in
!  output_slice() is not from the next time step
!
      if (lvideo) tslice = t
!
    endsubroutine wvid_prepare
!***********************************************************************
    subroutine wvid(f)
!
!  Write slices for animation of scalar field
!  (or one component of a vector field) on the perifery of a box.
!  Can be visualized in idl using rvid_box.pro.
!
!  13-nov-02/axel: added more fields
!  22-sep-07/axel: changed Xy to xy2, to be compatible with Mac
!  28-oct-18/PABourdin: moved output to IO modules 'output_slice'
!
      use General,         only: itoa
      use IO,              only: output_slice
      use Slices_methods,  only: assign_slices_scal
      use Chemistry,       only: get_slices_chemistry
      use Chiral,          only: get_slices_chiral
      use Cosmicray,       only: get_slices_cosmicray
      use Density,         only: get_slices_density,get_slices_pressure
      use Dustdensity,     only: get_slices_dustdensity
      use Dustvelocity,    only: get_slices_dustvelocity
      use EquationOfState, only: get_slices_eos
      use Energy,          only: get_slices_energy
      use Heatflux,        only: get_slices_heatflux
      use Hydro,           only: get_slices_hydro
      use Interstellar,    only: get_slices_interstellar
      use Magnetic,        only: get_slices_magnetic
      use Particles_main,  only: get_slices_particles
      use Pscalar,         only: get_slices_pscalar
      use Radiation,       only: get_slices_radiation
      use Shock,           only: get_slices_shock
      use Special,         only: get_slices_special
      use Testfield,       only: get_slices_testfield
      use Testflow,        only: get_slices_testflow
      use Testscalar,      only: get_slices_testscalar
      use Mpicomm,         only: mpiwtime
!
      real, dimension (mx,my,mz,mfarray), intent(IN) :: f
!
      logical :: lslices_legacy=.true.
      integer :: inamev
!
      type (slice_data) :: slices
      character (LEN=labellen) :: sname
      real :: time1
!
      slices%index=0
!
!  Loop over slices.
!
      inamev=1
      if (ip<=12.and.lroot) time1=mpiwtime()
      do while (inamev <= nnamev)
!
        if (trim(cformv(inamev))=='') then
          inamev=inamev+1
          cycle           ! skip undefined slices
        endif

        sname=trim(cnamev(inamev))

        slices%name=sname
        call assign_slices_scal(slices,slice_xy,slice_xz,slice_yz,slice_xy2,slice_xy3,slice_xy4,slice_xz2)
        if (lwrite_slice_r) slices%r => slice_r   ! no interpolation
        slices%ready=.false.
!
!  By default assume we're not using module hooks to get the slice contents
!
        lslices_legacy=.true.
!
!  Get slice information from the modules.
!
        lslices_legacy=.false.
        if (lchemistry)    call get_slices_chemistry   (f,slices)
        if (lchiral)       call get_slices_chiral      (f,slices)
        if (lcosmicray)    call get_slices_cosmicray   (f,slices)
        if (ldensity .or. lanelastic) &
                           call get_slices_density     (f,slices)
        if (lanelastic)    call get_slices_pressure    (f,slices)
        if (ldustdensity)  call get_slices_dustdensity (f,slices)
        if (ldustvelocity) call get_slices_dustvelocity(f,slices)
        if (lenergy)       call get_slices_energy      (f,slices)
        if (leos)          call get_slices_eos         (f,slices)
        if (lheatflux)     call get_slices_heatflux    (f,slices)
        if (lhydro.or.lhydro_kinematic) &
                           call get_slices_hydro       (f,slices)
        if (linterstellar) call get_slices_interstellar(f,slices)
        if (lmagnetic)     call get_slices_magnetic    (f,slices)
        if (lparticles)    call get_slices_particles   (f,slices)
        if (lpscalar)      call get_slices_pscalar     (f,slices)
        if (lradiation)    call get_slices_radiation   (f,slices)
        if (lshock)        call get_slices_shock       (f,slices)
        if (lspecial)      call get_slices_special     (f,slices)
        if (ltestfield)    call get_slices_testfield   (f,slices)
        if (ltestflow)     call get_slices_testflow    (f,slices)
        if (ltestscalar)   call get_slices_testscalar  (f,slices)
!
        if (lslices_legacy) then
          inamev=inamev+1
          cycle
        endif
!
!  Slice, or component of slice, ready for saving.
!
        if (slices%ready) then
          if (slices%index==0) then    ! If this wasn't a multi slice...
            inamev=inamev+1
          else
            sname=trim(sname)//trim(itoa(slices%index))
          endif
          call output_slice(lwrite_slice_yz , tslice, sname, 'yz',  x(ix_loc), ix, slices%yz)
          call output_slice(lwrite_slice_xz , tslice, sname, 'xz',  y(iy_loc), iy, slices%xz)
          call output_slice(lwrite_slice_xz2, tslice, sname, 'xz2', y(iy2_loc), iy2, slices%xz2)
          call output_slice(lwrite_slice_xy , tslice, sname, 'xy',  z(iz_loc), iz, slices%xy)
          call output_slice(lwrite_slice_xy2, tslice, sname, 'xy2', z(iz2_loc), iz2, slices%xy2)
          call output_slice(lwrite_slice_xy3, tslice, sname, 'xy3', z(iz3_loc), iz3, slices%xy3)
          call output_slice(lwrite_slice_xy4, tslice, sname, 'xy4', z(iz4_loc), iz4, slices%xy4)
          call output_slice(lwrite_slice_r  , tslice, sname, 'r',   r_rslice, 0, slices%r)
        else
          if (slices%index/=0) slices%index=0
          inamev=inamev+1
        endif
      enddo
      if (ip<=12.and.lroot) print*,'wvid: written slices in ',&
                                   mpiwtime()-time1,' seconds'
!
    endsubroutine wvid
!***********************************************************************
    subroutine setup_slices
!
!  Determine slice positions and whether slices are to be written on this
!  processor.
!
!  29-may-06/tobi: wrapped code from param_io.f90 into this subroutine
!  21-apr-15/MR: corrected i[xyz]_loc determination, see subroutine position
!
!  set slice position. The default for slice_position is 'p' for periphery,
!  although setting ix, iy, iz, iz2 by hand will overwrite this.
!  If slice_position is not 'p', then ix, iy, iz, iz2 are overwritten.
!
      use Slices_methods, only: alloc_slice_buffers, alloc_rslice, prep_rslice
      use General, only: itoa
      use IO, only: output_slice_position
      use Mpicomm, only: set_rslice_communicator
    
      character(LEN=80) :: text, data

      lwrite_slice_xy=.false. 
      lwrite_slice_xz=.false. 
      lwrite_slice_yz=.false. 
      lwrite_slice_xy2=.false.
      lwrite_slice_xy3=.false.
      lwrite_slice_xy4=.false.
      lwrite_slice_xz2=.false.
      lwrite_slice_r=.false.

      if (slice_position=='p' .or. slice_position=='S') then

        lwrite_slice_xy2=llast_proc_z; if (lwrite_slice_xy2)iz2_loc=n2
        lwrite_slice_xy=lfirst_proc_z; if (lwrite_slice_xy) iz_loc=n1
        lwrite_slice_xz=lfirst_proc_y; if (lwrite_slice_xz) iy_loc=m1
        lwrite_slice_yz=lfirst_proc_x; if (lwrite_slice_yz) ix_loc=l1
!
!  slice position in middle of the box in dependence of nprocy,nprocz
!  second horizontal slice is the uppermost layer
!
      elseif (slice_position=='m') then
!
!  xy2 is top layer as default.
!  Please set iz2 in run.in to select a different layer
!  where nghost+1 <= iz2 <= mzgrid-nghost
!
        lwrite_slice_yz=(ipx==nprocx/2)
        if (lwrite_slice_yz) then
          if (mod(nprocx,2)==0) then; ix_loc=l1; else; ix_loc=(l1+l2)/2; endif
        endif
        lwrite_slice_xz=(ipy==nprocy/2)
        if (lwrite_slice_xz) then
          if (mod(nprocy,2)==0) then; iy_loc=m1; else; iy_loc=(m1+m2)/2; endif
        endif
        lwrite_slice_xy=(ipz==nprocz/2)
        if (lwrite_slice_xy) then
          if (mod(nprocz,2)==0) then; iz_loc=n1; else; iz_loc=(n1+n2)/2; endif
        endif
        lwrite_slice_xy2=llast_proc_z; if (lwrite_slice_xy2) iz2_loc=n2
!
!  slice positions for spherical coordinates
!  w is for "wedges" since the outputs are
!  the midplane (rphi,theta=y(mpoint)) and four
!  wedges in rtheta (xy)
!
      elseif (slice_position=='w') then

        if (nprocx>1) call warning('setup_slice', &
            'slice_position=w may be wrong for nprocx>1')
        !midplane slices
        !ix_loc=nxgrid/2+nghost
        iy = nygrid/2+nghost           !MR: nghost not tb added!
        !meridional wedges, at 4 different
        !equally spaced azimuthal locations
        iz =  0*nzgrid/4+1+nghost
        iz2=  1*nzgrid/4+1+nghost
        iz3=  2*nzgrid/4+1+nghost
        iz4=  3*nzgrid/4+1+nghost
        ix_loc=0; iy_loc=0
!
!  Another slice position for spherical coordinates
!  s is for "surface" meaning theta-phi sections
!  keep iz_loc=n1, corresponding to a meridional slice on n=n1.
!
      elseif (slice_position=='s') then

        if (nprocx>1) call warning('setup_slice', &
            'slice_position=s may be wrong for nprocx>1')

        call xlocation(xtop_slice,ix_loc,lwrite_slice_yz)
        lwrite_slice_xy2=(ipz==nprocz/4); if (lwrite_slice_xy2) iz2_loc=n2
        lwrite_slice_xy=lfirst_proc_z; if (lwrite_slice_xy) iz_loc=n1
        lwrite_slice_xz=.false.; iy_loc=0
!
! Another slice position for spherical coordinates, for global disks with
! buffer zones. It will read the midplane (xz2), and three other surfaces:
! the theta-phi wall at constant radius (yz); the meridional plane at constant
! azimuth (xy); and the upper "lid", a radius-azimuth surface at constant
! theta, in the upper disk (xz). Both xy and yz are taken 10 grid cells away from
! the beginning of the grid. This is to avoid the boundary.
!
      elseif (slice_position=='d') then

        lwrite_slice_xy=lfirst_proc_z; if (lwrite_slice_xy) iz_loc=n1
        lwrite_slice_xz=lfirst_proc_y; if (lwrite_slice_xz) iy_loc=min(m1+10,m2)
        lwrite_slice_yz=lfirst_proc_x; if (lwrite_slice_yz) ix_loc=min(l1+10,l2)
        lwrite_slice_xz2=(ipy==nprocy/2)
        if (lwrite_slice_xz2) then
          if (mod(nprocy,2)==0) then; iy2_loc=m1; else; iy2_loc=(m1+m2)/2; endif
        endif
!
!  later we may also want to write other slices
!
        !call xlocation(xtop_slice,ix2_loc,lwrite_slice_yz2)
        !call ylocation(ytop_slice,iy2_loc,lwrite_slice_xz2)
!
!  slice position when the first meshpoint in z is the equator (sphere)
!  For one z-processor, iz remains n1, but iz2 is set to the middle.
!
!  TH: A certain processor layout is implied here
!
      elseif (slice_position=='e') then

        if (nprocx>1) call warning('setup_slice', &
            'slice_position=e may be wrong for nprocx>1')

        lwrite_slice_xy=lfirst_proc_z; if (lwrite_slice_xy) iz_loc=n1
        lwrite_slice_yz=(ipx==nprocx/2); if (lwrite_slice_yz) ix_loc=(l1+l2)/2

        lwrite_slice_xy2=(ipz==nprocz/4)
        if (lwrite_slice_xy2) then
          if (nprocz==1) then
            iz2_loc=(n1+n2)/2   !MR: not iz2_loc=(iz+n2)/2!  
          else
            iz2_loc=n2
          endif
        endif

        lwrite_slice_xz=(ipy==nprocy/2)
        if (lwrite_slice_xz) then
          if (nprocy==1) then
            iy_loc=(m1+m2)/2
          else
            iy_loc=m1
          endif
        endif
!
!  slice position similar to periphery (p), but the two z positions
!  can now be given in terms of z (zbot_slice, ztop_slice).
!
      elseif (slice_position=='c') then

        call zlocation(zbot_slice,iz_loc,lwrite_slice_xy)
        call zlocation(ztop_slice,iz2_loc,lwrite_slice_xy2)
        lwrite_slice_xz=lfirst_proc_y; if (lwrite_slice_xz) iy_loc=m1
        lwrite_slice_yz=lfirst_proc_x; if (lwrite_slice_yz) ix_loc=l1
!
!  periphery of the box, but the other way around
!
      elseif (slice_position=='q') then

        lwrite_slice_xy2=lfirst_proc_z; if (lwrite_slice_xy2) iz2_loc=n1
        lwrite_slice_xy=llast_proc_z; if (lwrite_slice_xy) iz_loc=n2
        lwrite_slice_xz=llast_proc_y; if (lwrite_slice_xz) iy_loc=m2
        lwrite_slice_yz=llast_proc_x; if (lwrite_slice_yz) ix_loc=l2

      else
        if (lroot) &
          call fatal_error('setup_slices', &
                           'No such value for slice_position: '//slice_position)
      endif
!
!  Spherical admits only position 'w', 's', or 'd'. Break if this is not met.
!  Also, turn extra r-theta slices to false in case of
!  non-spherical coordinates
!
      if (slice_position/='w'.and.slice_position/='s') then
        lwrite_slice_xy3=.false.
        lwrite_slice_xy4=.false.
      endif
      if (slice_position/='d') then
        lwrite_slice_xz2=.false.
      endif
!
!  Overwrite slice positions if any ix,iy,iz,iz2,iz3,iz4 is greater then zero
!  position sets lwrite_slice_* according to whether or not the executing proc 
!  contributes to the slice data
!
      call position(ix,ipx,nx,ix_loc,lwrite_slice_yz)
      call position(iy,ipy,ny,iy_loc,lwrite_slice_xz)
      call position(iy2,ipy,ny,iy2_loc,lwrite_slice_xz2)
      call position(iz,ipz,nz,iz_loc,lwrite_slice_xy)
      call position(iz2,ipz,nz,iz2_loc,lwrite_slice_xy2)
      call position(iz3,ipz,nz,iz3_loc,lwrite_slice_xy3)
      call position(iz4,ipz,nz,iz4_loc,lwrite_slice_xy4)
!
!  Slices for star-in-box runs.
!
      lwrite_slice_r = (r_rslice>0.) .and. (nth_rslice>0) .and. (nph_rslice>0)
      if (lwrite_slice_r) then
        if (coord_system/='cartesian') then
          call warning('setup_slices','r-slices only meaningful in Cartesian runs')
          lwrite_slice_r=.false.
        elseif (r_rslice>min(xyz1(1),-xyz0(1)) .or. r_rslice>min(xyz1(2),-xyz0(2)) &
                                               .or. r_rslice>min(xyz1(3),-xyz0(3)) ) then
          call warning('setup_slices','rslice partially outside box')
          lwrite_slice_r=.false.
        elseif (dimensionality<3) then
          call warning('setup_slices','r-slices only meaningful in 3D- runs')
          lwrite_slice_r=.false.
        elseif (any(.not.lequidist)) then
          call not_implemented('setup_slices','r-slices in non-equidistannt grids')
        elseif (nprocx>1.and.mod(nprocx,2)/=0 .or. nprocy>1.and.mod(nprocy,2)/=0) then
          call not_implemented('setup_slices','r-slices for odd nprocx or nprocy')
        else
          call prep_rslice
          call set_rslice_communicator
        endif
      endif

      if (.not.lactive_dimension(3)) then
        lwrite_slice_xy2=.false.; lwrite_slice_xy3=.false.; lwrite_slice_xy4=.false.
        iz2_loc=1; iz3_loc=1; iz4_loc=1
      endif 

      if (.not.lactive_dimension(2)) then
        lwrite_slice_xz2=.false.; iy2_loc=1
      endif

      call output_slice_position
!
      if (lroot.and.dvid>0.) then  !MR: tbdone: write global position from all procs

        write (*,*)'setup_slices: slice_position = '//slice_position
        text=''; data=''
        if (lwrite_slice_yz) then
          text='ix_loc,'; data=itoa(ix_loc)
        endif
        if (lwrite_slice_xz) then
          text=trim(text)//'iy_loc,'; data=trim(data)//' '//itoa(iy_loc)
        endif
        if (lwrite_slice_xy) then
          text=trim(text)//'iz_loc,'; data=trim(data)//' '//itoa(iz_loc)
        endif
        if (lwrite_slice_xy2) then
          text=trim(text)//'iz2_loc,'; data=trim(data)//' '//itoa(iz2_loc)
        endif
        if (lwrite_slice_xy3) then
          text=trim(text)//'iz3_loc,'; data=trim(data)//' '//itoa(iz3_loc)
        endif
        if (lwrite_slice_xy4) then
          text=trim(text)//'iz4_loc,'; data=trim(data)//' '//itoa(iz4_loc)
        endif
        if (lwrite_slice_xz2) then
          text=trim(text)//'iy2_loc,'; data=trim(data)//' '//itoa(iy2_loc)
        endif
        write (*,*) 'setup_slices: '//trim(text)//' (video files) = '//data
      endif
!
      call alloc_slice_buffers(slice_xy,slice_xz,slice_yz,slice_xy2,slice_xy3,slice_xy4,slice_xz2)
      if (lwrite_slice_r .and. .not.allocated(slice_r)) call alloc_rslice(slice_r)

    endsubroutine setup_slices
!***********************************************************************
    subroutine prep_xy_slice(izloc)

      use General, only: indgen

      integer, intent(IN) :: izloc

      real, dimension(nygrid/2) :: yloc
      !real, dimension(nygrid/2,1) :: thphprime

      if (ipz<=nprocz/3) then
!
! Line trough Southern cap
!
        yloc=xyz1(2)+indgen(nygrid/2)*dy

        !call yy_transform_strip_other(yloc,(/z(izloc)/),thphprime)
        !nok=prep_interp(thphprime,intcoeffs)

      elseif (ipz>=2*nprocz/3) then
!
! Line trough Nouthern cap
!
        yloc=xyz0(2)-(nygrid/2+1-indgen(nygrid/2))*dy 

        !call yy_transform_strip_other(yloc,(/z(izloc)/),thphprime)
        !nok=prep_interp(thphprime,intcoeffs,iyinyang_intpol_type,thrange_cap)

      endif 

    endsubroutine prep_xy_slice
!***********************************************************************
endmodule Slices