! $Id$
!
!  This module tracks the evolution of particles in terms of
!  user defined "states".  Currently the only implemented state is in terms
!  of the local gas velocity.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
!
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! MPVAR CONTRIBUTION 5
! MAUX CONTRIBUTION 0
!
! CPARAM logical, parameter :: lparticles_diagnos_state=.true.
!
!
!***************************************************************
module Particles_diagnos_state
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
  use Particles_cdata
  use Particles_mpicomm
  use Particles_sub
!
  implicit none
!
  include 'particles_diagnos_state.h'
!
  character (len=labellen) :: state_def ='velz'
  real :: slow_vel_trigger=0.,slow_vel_trigger2=0.
!
  namelist /particles_diagnos_state_run_pars/ &
      slow_vel_trigger, state_def
!
  integer :: idiag_upm=0, idiag_vprms=0, idiag_uzplus=0, idiag_uzminus=0
  integer :: idiag_uzminuscount=0, idiag_uzpluscount=0
!
  contains
!***********************************************************************
    subroutine register_pars_diagnos_state()
!
!  Set up indices for access to the fp and dfp arrays
!
      if (lroot) call svn_id( &
           "$Id$")
!
      call append_npvar('ipss',ipss)  ! particle state
      call append_npvar('ipst',ipst)  ! time particle entered current state
      call append_npvar('ipxx',ipxx)  ! x position particle entered current state
      call append_npvar('ipyy',ipyy)  ! y position particle entered current state
      call append_npvar('ipzz',ipzz)  ! z position particle entered current state
!
    endsubroutine register_pars_diagnos_state
!***********************************************************************
    subroutine initialize_pars_diagnos_state(f)
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      slow_vel_trigger2=slow_vel_trigger**2
!
      call keep_compiler_quiet(f)
!
    endsubroutine initialize_pars_diagnos_state
!***********************************************************************
    subroutine init_particles_diagnos_state(fp)
!
      real, dimension (mpar_loc,mparray) :: fp
!
      integer :: k
!
      intent (out) :: fp
!
      do k=1,npar_loc
        fp(k,ipst)=t
        fp(k,ipss)=-1.
        fp(k,ipxx:ipzz)=fp(k,ixp:izp)
      enddo
!
    endsubroutine init_particles_diagnos_state
!***********************************************************************
    subroutine insert_particles_diagnos_state(fp, npar_loc_old)
!
! Allow for the insertion of particles
!
      real, dimension (mpar_loc,mparray) :: fp
      integer :: npar_loc_old
!
      integer :: k
!
      intent (inout) :: fp
!
      do k=npar_loc_old+1,npar_loc
        fp(k,ipst)=t
        fp(k,ipss)=-1.
        fp(k,ipxx:ipzz)=fp(k,ixp:izp)
      enddo
!
    endsubroutine insert_particles_diagnos_state
!***********************************************************************
    subroutine read_pars_diag_state_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=particles_diagnos_state_run_pars, IOSTAT=iostat)
!
    endsubroutine read_pars_diag_state_run_pars
!***********************************************************************
    subroutine write_pars_diag_state_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=particles_diagnos_state_run_pars)
!
    endsubroutine write_pars_diag_state_run_pars
!***********************************************************************
    subroutine rprint_particles_diagnos_state(lreset,lwrite)
!
      use Diagnostics
!
      logical :: lreset
      logical, optional :: lwrite
!
      integer :: iname,inamez,inamey,inamex,inamexy,inamexz,inamer,inamerz
!
!  Reset everything in case of reset.
!
      if (lreset) then
        idiag_upm=0; idiag_uzminus=0; idiag_uzplus=0
        idiag_uzminuscount=0; idiag_uzpluscount=0
      endif
!
!  Run through all possible names that may be listed in print.in.
!
      if (lroot .and. ip<14) print*,'rprint_particles: run through parse list'
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'upm',idiag_upm)
        call parse_name(iname,cname(iname),cform(iname),'uzminus',idiag_uzminus)
        call parse_name(iname,cname(iname),cform(iname),'uzplus',idiag_uzplus)
        call parse_name(iname,cname(iname),cform(iname),'uzminuscount',idiag_uzminuscount)
        call parse_name(iname,cname(iname),cform(iname),'uzpluscount',idiag_uzpluscount)
      enddo
!
    endsubroutine rprint_particles_diagnos_state
!***********************************************************************
    subroutine persistence_check(fp, k, uup)
!
      real, dimension (mpar_loc,mparray) :: fp
      real :: r_new
      integer :: k
      real, dimension(3) :: uup
      logical :: l_swap
!
!  check for state change
!
      l_swap=.false.
      call evaluate_state(fp, k, uup, l_swap, r_new)
!
!  execute state change, write state change to file
!
      if (l_swap) then
        if (.not. (fp(k,ipss) < 0)) call &
            data_store_persistence(fp, k, r_new, .false.)
        fp(k,ipss)=r_new
        fp(k,ipst)=t
        fp(k,ipzz)=fp(k,izp)
      endif
!
!  if final particle in pencil, close file
!
      if (k==k2_imn(imn)) &
          call data_store_persistence(fp, -1, -1., .true.)
!
!  some state diagnostics
!
      if (ldiagnos) then
        if (idiag_upm/=0) &
            call sum_par_name((/sum((fp(k,ivpx:ivpz)-uup)**2)/),idiag_upm,lsqrt=.true.)
        if (idiag_uzminus/=0) then
          if (uup(3) < 0) then
            call sum_par_name((/ uup(3) /),idiag_uzminus)
            call sum_par_name((/ 1. /),idiag_uzminuscount)
          else
            call sum_par_name((/ 0. /),idiag_uzminus)
            call sum_par_name((/ 0. /),idiag_uzminuscount)
          endif
        endif
        if (idiag_uzplus/=0) then
          if (uup(3) >= 0) then
            call sum_par_name((/ uup(3) /),idiag_uzplus)
            call sum_par_name((/ 1. /),idiag_uzpluscount)
          else
            call sum_par_name((/ 0. /),idiag_uzplus)
            call sum_par_name((/ 0. /),idiag_uzpluscount)
          endif
        endif
      endif
!
     endsubroutine persistence_check
!***********************************************************************
    subroutine evaluate_state(fp, k, uup, l_swap, r_new)
!
      real, dimension (mpar_loc,mparray) :: fp
      real :: r_new
      integer :: k
      real, dimension(3) :: uup
      logical :: l_swap
!
! as more cases are added, the call to
! evaluate_state may need to be moved/changed/elaborated on
!
      select case(state_def)
      case ('velz') !gas z velocity up or down in lab frame (particle settling)
        if (uup(3) > 0) then
          r_new=2
        else
          r_new=1
        endif
      case default
        call fatal_error('evaluate_state','No such such value for '// &
          'state_def: '//trim(state_def))
      endselect
!
      if (abs(fp(k, ipss)-r_new) > 0.1) then
        l_swap=.true.
      elseif (abs(fp(k, ipss)-r_new) <= 0.1) then
        l_swap=.false.
      endif
!
    endsubroutine evaluate_state
!***********************************************************************
    subroutine data_store_persistence(fp, k, r_new, l_close)
!
      use General, only: safe_character_assign, itoa
!
      real, dimension (mpar_loc,mparray) :: fp
      real :: r_new
      integer :: k
      logical :: l_close                   !do we close the output file?
!
      logical, save       :: l_opened=.false. !is state output file open
!
      integer             :: lun_input=91.
      real                :: time_span
      real, dimension(3)  :: travel
      character (len=128) :: filename
!
!  Close file if told to and it is opened
!  Open file if it is closed, and not told to close it
!
      if (l_close .and. l_opened) then
        close(UNIT=lun_input)
        l_opened=.false.
      elseif ((.not. l_close) .and. (.not. l_opened)) then
        call safe_character_assign(filename,trim(datadir)//&
            '/persistence'//trim(itoa(iproc))//'.dat')
        open(UNIT=lun_input,FILE=trim(filename),&
            POSITION='append',FORM='formatted')
        l_opened=.true.
      endif
!
      if (l_opened) then
        time_span=t-fp(k,ipst)
        travel=fp(k,ixp:izp)-fp(k,ipxx:ipzz)
        write(UNIT=lun_input,FMT= &
            '(f12.5, i9, f4.1, f12.5, f4.1, f12.5, f12.5, f12.5)') &
            t, ipar(k), fp(k,ipss),&
            time_span, r_new, travel(1), travel(2), travel(3)
!
!  ie: current time, particle label, old state, time spend in old state,
!      new state, distance traveled in old state (x,y,z)
!
      endif
!
    endsubroutine data_store_persistence
!***********************************************************************
endmodule Particles_diagnos_state