module initial
use precision 
implicit none

public gtc_initial, timer, gtc_final
save
contains
!=============================================================================
subroutine init_MPI
    use global_parameters

    integer ierror
    character(len=16) cdate(4)

! MPI initialize, total # of PE, and rank of PE
    call mpi_init(ierror)
    call mpi_comm_size(mpi_comm_world,numberpe,ierror)
    call mpi_comm_rank(mpi_comm_world,mype,ierror)

    if(mype==0)then
        open(gtcout,file='gtc.out',form='formatted',status='replace')
        call date_and_time(cdate(1),cdate(2))
        write(gtcout,*)"GTC Version xxx"
        write(gtcout,*)"Compiled on ",cdate(1)," at ",cdate(2)
#ifdef _FRC
        write(gtcout,'(/,"===================================================")')
        write(gtcout,*)' compiled for field-reversed configuration equilibrium.'
        write(gtcout,'("===================================================",/)')
#elif _TOKAMAK
        write(gtcout,'(/,"===================================================")')
        write(gtcout,*)' compiled for tokamak equilibrium.'
        write(gtcout,'("===================================================",/)')
#endif
 
        write(gtcout,'(/,"===================================================")')
        write(gtcout,*)' Compiled with PETSc.'
        write(gtcout,'("===================================================",/)')

#ifdef _OPENMP
!$omp parallel
!$omp master
        nthreads=omp_get_num_threads()
        write(gtcout,'(/,"===================================================")')
        write(gtcout,*)' Number of OpenMP threads = ',nthreads
        write(gtcout,'("===================================================",/)')
!$omp end master
!$omp end parallel
#else
        write(gtcout,'(/,"===================================================")')
        write(gtcout,*)' Run without OpenMP threads'
        write(gtcout,'("===================================================",/)')
#endif
    endif

end subroutine init_MPI
!=============================================================================

subroutine gtc_initial(cdate,timecpu)
    use global_parameters
    use particle_array
#ifdef _PETSc
    use petsc_gtc, only: petsc_init
#endif
    use read_in_parameters, only: irun
    implicit none

    integer ierror,i,m,nsize
    real(doubleprec) timecpu(10)
    integer,dimension(:),allocatable :: mget,mput
    character(len=16) cdate(4)
#ifdef _OPENACC
    integer :: numDevice,myDevice,ierr
#endif

    call date_and_time(cdate(1),cdate(2))
    call init_MPI

! initialize GPU
#ifdef _OPENACC
    call acc_init(acc_device_nvidia)
    numDevice=acc_get_num_devices(acc_device_nvidia)
    if(numDevice/=0)then
        myDevice=mod(mype,numDevice)
        call acc_set_device_num(myDevice,acc_device_nvidia)
        myDevice=acc_get_device_num(acc_device_nvidia)
    ! print *, 'pe=',mype,'using ',idevice
    !call setupgpu(1)
        if(mype==0)then
            write(gtcout,'(/,"===================================")')
            write(gtcout,*)' Number of OpenACC Devices = ',numDevice 
            write(gtcout,'("===================================",/)')
        endif
    else
        call acc_set_device_type(acc_device_host)
        if(mype==0)then
            write(gtcout,'(/,"===================================")')
            write(gtcout,*)' Run without OpenACC Devices'
            write(gtcout,'("===================================",/)')
        endif
    endif
#endif

! initiate adios
#if ADIOS
    CALL adios_init ("config.xml"//char(0), MPI_COMM_WORLD, MPI_COMM_SELF, MPI_INFO_NULL)
#endif

!!XY rolling restart, initialize restart directory and control var
    irest=0
    FileExit=0

! initialize timer
    istep=0
    timecpu=0.0
    timecpu(9)=mpi_wtime()
    timecpu(10)=timecpu(9)

! numerical constant
    pi=4.0_lk*atan(1.0_lk)
    pi2=2.0_lk*pi
    pi2_inv=1.0_lk/pi2

! **** Use the intrinsic F90 random number generator *****
! initialize f90 random number generator
    ! call random_seed !comemt out this line for identitical random number
    call random_seed(size=nsize)
    allocate(mget(nsize),mput(nsize))
    call random_seed(get=mget)
    do i=1,nsize
     !TODO: system_clock causes problem on Titan with pgi compiler

     !call system_clock(m)   !random initialization for collision

        if(irun==0)m=1         !same initialization
         mput(i)=111111*(mype+1)+m+mget(i)
    enddo
    call random_seed(put=mput)
    deallocate(mget,mput)

! Read the input file that contains the run parameters
    call read_input_params

! Set up the particle decomposition within each toroidal domain
    call set_particle_decomp

! PETSc initialize
#ifdef _PETSc
    call petsc_init
#endif

end subroutine gtc_initial
!=====================================================================
subroutine read_input_params
!=====================================================================

    use global_parameters
    use read_in_parameters
    use precision
    implicit none

    logical file_exist
    integer ierror
  
    namelist /control_input/ irun, mstep, msnap, ndiag, nonlinear, irestore, toroidaln, tstep, &
                             track_particles, numerfrc, gamma0, frcfullf, global_solver_on, &
                             local_solver_model, nfilter, sol_pfilter, uni_profile, diag_region, &
                             sol_para_periodicity, antenna_on, fluid_on, min_r, ifullk, ifluxtube, izonal
    namelist /boundary_decay/ BC_type,BC_nbound0,BC_nbound1
    namelist /modes_input/ n_modes, m_modes
    namelist /simulation_domain_input/ Zsol_domain0, Zsol_domain1, Rdomain0, Rdomain1
    namelist /simulation_grid_input/ mpsisol_input, mpsicore_input, mparasol_input, mparacore_input, mtoroidal_input, iflux_core, iflux_sol, para_diag_location_core, para_diag_location_sol
    namelist /particle_input/ micell, aion, qion, ngyroi, &
                              mecell, aelectron, qelectron, ngyroe, kinetice
    namelist /equilibrium_input/ etemp0_input, eden0_input, r0_input, b0_input
    namelist /constructed_profile/ construct_nT,n_coef,t_coef,n0_1d,t0_1d
 
  
  ! Since it is preferable to have only one MPI process reading the input file,
  ! we choose the master process to set the default run parameters and to read
  ! the input file. The parameters will then be broadcast to the other processes.
  
    
    if(mype==0) then
        nonlinear=0
        irestore=0
        gamma0=0.0_lk
    !default values for recently added parameters. 
        n_modes=(/1,1,1,1,1,1,1,1/)
        m_modes=(/1,2,3,4,5,6,7,8/)
        kinetice=0
        ifluxtube=0
        izonal=1
!!! default parameters for boundary decay
        BC_type=0
        BC_nbound0=0
        BC_nbound1=0
!!! default coefficients for temperature/density 1d profile
        construct_nT = 0
        n_coef = (/ 0.25_lk, 1.5_lk, 0.75_lk /)
        t_coef = (/ 0.42_lk, 1.5_lk, 0.75_lk /)
        t0_1d = 0.00412195_lk*1.0_lk
        n0_1d = 0.75_lk
  ! Test if the input file gtc.in exists
        inquire(file='gtc.in',exist=file_exist)
        if (file_exist) then
            open(55,file='gtc.in',status='old')
            read(55,nml=control_input,iostat=ierror)
            read(55,nml=boundary_decay,iostat=ierror)
            read(55,nml=modes_input,iostat=ierror)
            read(55,nml=simulation_domain_input,iostat=ierror)
            read(55,nml=simulation_grid_input,iostat=ierror)
            read(55,nml=particle_input,iostat=ierror)
            read(55,nml=equilibrium_input,iostat=ierror)
            read(55,nml=constructed_profile,iostat=ierror)
            if(ierror > 0)write(*,*)'WARNING: check input_parameters'
            close(55)
        else
            write(*,*)'Cannot find file gtc.in !!!'
            stop
        endif
        if(nonlinear==0) izonal=0
        if(nonlinear==0) irestore=0
  ! treat tearing mode
  
        mstep=max(ndiag,mstep)
        msnap=min(msnap,mstep/ndiag)
  
      
  ! write parameters to standard output
        write(gtcout,control_input)
        write(gtcout,boundary_decay)
        write(gtcout,modes_input)
        write(gtcout,simulation_domain_input)
        write(gtcout,simulation_grid_input)
        write(gtcout,particle_input)
        write(gtcout,equilibrium_input)
        write(gtcout,constructed_profile)
  
      ! approximate # of particle per cell
        mi=micell
        me=mecell
        mtoroidal=mtoroidal_input

        if((iflux_core > mpsicore_input) .or. (iflux_sol > mpsisol_input))then
            write(gtcout,*)"wrong iflux_core or iflux_sol value, which is out of 1d array bound"
            stop
        endif


    endif
  
  ! Now send the parameter values to all the other MPI processes
    call broadcast_input_params

end subroutine read_input_params

!=============================================================================
subroutine broadcast_input_params
!=============================================================================
    use global_parameters
    use read_in_parameters
    implicit none
  
    integer ierror, icount
  
! The master process, mype=0, holds all the input parameters. We need
! to broadcast their values to the other processes. Instead of issuing
! an expensive MPI_BCAST() for each parameter, it is better to pack
! everything in a single vector, broadcast it, and unpack it.
! mpi broadcast example
!      call MPI_BCAST(buffer,count,datatype,root,comm,ierror)
!      [ INOUT buffer] starting address of buffer (choice)
!      [ IN count] number of entries in buffer (integer)
!      [ IN datatype] data type of buffer (handle)
!      [ IN root] rank of broadcast root (integer)
!      [ IN comm] communicator (handle)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    icount=1
! control_input
! integer
    call MPI_BCAST(irun,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(mstep,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(msnap,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)

    call MPI_BCAST(ndiag,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(nonlinear,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(irestore,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(toroidaln,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)

    call MPI_BCAST(track_particles,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(numerfrc,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(gamma0,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)

    call MPI_BCAST(frcfullf,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(global_solver_on,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(local_solver_model,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)

    call MPI_BCAST(nfilter,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(uni_profile,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(diag_region,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)

    call MPI_BCAST(sol_pfilter,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(sol_para_periodicity,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(antenna_on,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)

    call MPI_BCAST(fluid_on,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(ifullk,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(ifluxtube,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(izonal,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)

! real
    call MPI_BCAST(tstep,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(min_r,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)

! boundary decay control parameters
    icount=1
    call MPI_BCAST(BC_type,icount,MPI_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(BC_nbound0,icount,MPI_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(BC_nbound1,icount,MPI_Rsize,0,MPI_COMM_WORLD,ierror)
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    icount=num_modes   ! 8 as default
! modes_input
    call MPI_BCAST(n_modes,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(m_modes,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    icount=1
! simulation_domain_input
    call MPI_BCAST(Zsol_domain0,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(Zsol_domain1,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(Rdomain0,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(Rdomain1,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    icount=1
! simulation_grid_input
    call MPI_BCAST(mpsisol_input,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(mpsicore_input,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(mparasol_input,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(mparacore_input,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(mtoroidal_input,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(iflux_core,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(iflux_sol,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(para_diag_location_core,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(para_diag_location_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
! notice that we need to broadcast mtoroidal (IMPORTANT)
    call MPI_BCAST(mtoroidal,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)

! particle_input
    call MPI_BCAST(mi,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(aion,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(qion,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(ngyroi,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(kinetice, icount, MPI_INTEGER, 0, MPI_COMM_WORLD, ierror)

    call MPI_BCAST(me,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(aelectron,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(qelectron,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(ngyroe,icount,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)

! equilibrium_input
    call MPI_BCAST(etemp0_input,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(eden0_input,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(r0_input,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(b0_input,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)

! constructed temperature/density profile input
!    call MPI_BCAST(b0_input,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    icount=3
    call MPI_BCAST(n_coef,icount,MPI_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(t_coef,icount,MPI_Rsize,0,MPI_COMM_WORLD,ierror)
    icount=1
    call MPI_BCAST(construct_nT,icount,MPI_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(n0_1d,icount,MPI_Rsize,0,MPI_COMM_WORLD,ierror)
    call MPI_BCAST(t0_1d,icount,MPI_Rsize,0,MPI_COMM_WORLD,ierror)
!$acc update device(n_modes)
   
end subroutine broadcast_input_params

!=============================================================================
subroutine set_particle_decomp
!=============================================================================
    use global_parameters
    use read_in_parameters, only: toroidaln
    implicit none
  
    integer  :: i,j,k,pe_number,mtest,ierror
    logical file_exist 
  
  
  ! ----- First we verify the consistency of mtoroidal and npartdom -------
  ! The number of toroidal domains (mtoroidal) times the number of particle
  ! "domains" (npartdom) needs to be equal to the number of processor "numberpe".
  ! numberpe cannot be changed since it is given on the command line.
  
  ! numberpe must be a multiple of mtoroidal
    if(mod(numberpe,mtoroidal) /= 0)then
        write(gtcout,*)'Wrong PE #; PE=',numberpe,'Mtoroidal=',mtoroidal
        stop
    endif
  
  !number of particle decomposition.
    npartdom=numberpe/mtoroidal
    if(mype==0)then
        write(gtcout,*)'*******************************************************'
        write(gtcout,*)'  Using npartdom=',npartdom,' and mtoroidal=',mtoroidal
        write(gtcout,*)'*******************************************************'
        write(gtcout,*)
    endif
  
  ! We now give each PE (task) a unique domain identified by 2 numbers: the
  ! particle and toroidal domain numbers.
  !    particle_domain_location = rank of the particle domain holding mype
  !    toroidal_domain_location = rank of the toroidal domain holding mype
  !
  ! On the IBM SP, the MPI tasks are distributed in an orderly fashion to each
  ! node unless the LoadLeveler instruction "#@ blocking = unlimited" is used.
  ! On Seaborg for example, the first 16 tasks (mype=0-15) will be assigned to
  ! the first node that has been allocated to the job, then the next 16
  ! (mype=16-31) will be assigned to the second node, etc. When not using the
  ! OpenMP, we want the particle domains to sit on the same node because
  ! communication is more intensive. To achieve this, successive PE numbers are
  ! assigned to the particle domains first.
  ! It is easy to achieve this ordering by simply using mype/npartdom for
  ! the toroidal domain and mod(mype,npartdom) for the particle domain.
  !
  !  pe_number=0
  !  do j=0,mtoroidal-1
  !     do i=0,npartdom-1
  !        pe_grid(i,j)=pe_number
  !        particle_domain_location(pe_number)=i
  !        toroidal_domain_location(pe_number)=j
  !        pe_number=pe_number+1
  !     enddo
  !  enddo
  
    particle_domain_location=mod(mype,npartdom)
    toroidal_domain_location=mype/npartdom
  
  !  write(0,*)'mype=',mype,"  particle_domain_location =",&
  !            particle_domain_location,' toroidal_domain_location =',&
  !            toroidal_domain_location,' pi=',pi
  
  ! Domain decomposition in toroidal direction.
    if(toroidaln<1)then
        write(gtcout,*)'Please check gtc.in for toroidaln (must be an integer equal to or larger than 1)!!!',&
                      'You put in toroidaln=',toroidaln,' but now it is toroidaln=1.'
        toroidaln=1
    endif     
    torbound=real(pi2)/real(toroidaln)
    nmode=toroidaln
    zeta0=torbound*real(toroidal_domain_location)/real(mtoroidal)
    zeta1=torbound*real(toroidal_domain_location+1)/real(mtoroidal)
  
  ! grid spacing in the toroidal direction, one toroidal grip per toroidal domain
    deltaz=zeta1-zeta0
  
    if(mype==0)then
        write(gtcout,*)'Toroidal domain : [0,2pi/',toroidaln,'], deltaz=',deltaz
    endif
  
  ! ---- Create particle domain communicator and toroidal communicator -----
  ! We now need to create a new communicator which will include only the
  ! processes located in the same toroidal domain. The particles inside
  ! each toroidal domain are divided equally between "npartdom" processes.
  ! Each one of these processes will do a charge deposition on a copy of
  ! the same grid, requiring a toroidal-domain-wide reduction after the
  ! deposition. The new communicator will allow the reduction to be done
  ! only between those processes located in the same toroidal domain.
  !
  ! We also need to create a purely toroidal communicator so that the
  ! particles with the same particle domain id can exchange with their
  ! toroidal neighbors.
  !
  ! Form 2 subcommunicators: one that includes all the processes located in
  ! the same toroidal domain (partd_comm), and one that includes all the
  ! processes part of the same particle domain (toroidal_comm).
  ! Here is how to create a new communicator from an old one by using
  ! the MPI call "MPI_COMM_SPLIT()".
  ! All the processes passing the same value of "color" will be placed in
  ! the same communicator. The "rank_in_new_comm" value will be used to
  ! set the rank of that process on the communicator.
  !  call MPI_COMM_SPLIT(old_comm,color,rank_in_new_comm,new_comm,ierror)
  
  ! particle domain communicator (for communications between the particle
  ! domains WITHIN the same toroidal domain)
    call MPI_COMM_SPLIT(MPI_COMM_WORLD,toroidal_domain_location,&
                        particle_domain_location,partd_comm,ierror)
  
  ! toroidal communicator (for communications BETWEEN toroidal domains of same
  ! particle domain number)
    call MPI_COMM_SPLIT(MPI_COMM_WORLD,particle_domain_location,&
                        toroidal_domain_location,toroidal_comm,ierror)
  
    call mpi_comm_size(partd_comm,nproc_partd,ierror)
    call mpi_comm_rank(partd_comm,myrank_partd,ierror)
  
    call mpi_comm_size(toroidal_comm,nproc_toroidal,ierror)
    call mpi_comm_rank(toroidal_comm,myrank_toroidal,ierror)
  
  !  write(0,*)'mype=',mype,'  nproc_toroidal=',nproc_toroidal,&
  !       ' myrank_toroidal=',myrank_toroidal,'  nproc_partd=',nproc_partd,&
  !       ' myrank_partd=',myrank_partd
  
    if(nproc_partd/=npartdom)then
        write(0,*)'*** nproc_partd=',nproc_partd,' NOT EQUAL to npartdom=',npartdom
        call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    endif
  
    if(nproc_toroidal/=mtoroidal)then
        write(*,*)'*** nproc_toroidal=',nproc_toroidal,' NOT EQUAL to mtoroidal=',&
                mtoroidal
        call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    endif
  
  ! We now find the toroidal neighbors of the current toroidal domain and
  ! store that information in 2 easily accessible variables. This information
  ! is needed several times inside the code, such as when particles cross
  ! the domain boundaries. We will use the toroidal communicator to do these
  ! transfers so we don't need to worry about the value of myrank_partd.
  ! We have periodic boundary conditions in the toroidal direction so the
  ! neighbor to the left of myrank_toroidal=0 is (mtoroidal-1).
  
    left_pe=mod(myrank_toroidal-1+mtoroidal,mtoroidal)
    right_pe=mod(myrank_toroidal+1,mtoroidal)
  
end subroutine set_particle_decomp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine gtc_final(cdate,timecpu)
    use global_parameters
#ifdef _PETSc
#include "petscversion.h"
#if PETSC_VERSION_GT(3,8,0)
    use petsc_gtc, only: petsc_final, petsc_ps
#else
    use petsc_gtc, only: petsc_final
#endif
#endif
    implicit none
 
    integer ierror,iorestart,i
    real(doubleprec) timecpu(10)
    character(len=16) ic(9)
    character(len=16) cdate(4)
 
!total cpu and wall clock time
    timecpu(10)=mpi_wtime()
    timecpu(9)=timecpu(10)-timecpu(9)
    ic(1)='pusher:'
    ic(2)='shift:'
    ic(3)='charge:'
    ic(4)='electron:'
    ic(5)='fast:'
    ic(6)='Poisson:'
    ic(7)='initial:'
    ic(8)='IO:'
    ic(9)='total:'
    if(mype==0)then
        write(gtcout,*)'CPU TIME USAGE (in SEC):'
        !write(gtcout,*)ic
        !write(gtcout,'(8(1pe8.1),/)')timecpu(1:8)
        write(gtcout,'(A16,1pe12.5,1pe16.5)')&
              (ic(i),timecpu(i),timecpu(i)/timecpu(9), i=1,9)
 
! Restart file info
        FileExit=1
        iorestart=345
        open(iorestart,file="FileExit.out",status="replace")
        write(iorestart,"(A9,i1)")"FileExit=",FileExit
        write(iorestart,"(A9,i5)")"irest   =",irest
        if(mod(irest+1,2)==0)then
            write(iorestart,"(A12)")"restart_dir1"
        else
            write(iorestart,"(A12)")"restart_dir2"
        endif
        close(iorestart)
 
! record program end time
        call date_and_time(cdate(3),cdate(4))
        write(gtcout,*)'Program starts at DATE=', cdate(1), 'TIME=', cdate(2)
        write(gtcout,*)'Program ends at   DATE=', cdate(3), 'TIME=', cdate(4)
 
! close standard output file
        if(gtcout /= 6 .and. gtcout /= 0)close(gtcout)
   endif
 
! PETSc finalize
#ifdef _PETSc
#if PETSC_VERSION_GT(3,8,0)
    call petsc_final(petsc_ps)
#else
    call petsc_final
#endif
#endif
 
#if ADIOS
    call adios_finalize (mype)
#endif
 
! MPI finalize
    call mpi_finalize(ierror)
 
end subroutine gtc_final
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!=========================================
subroutine timer(t0,timecpu)
!=========================================
    use precision
    implicit none
    real(doubleprec) t0,t1,dt,timecpu

! Get cpu usage time since the begine of run and subtract value from the previous call
    t1=mpi_wtime()
    dt=t1-t0
    t0=t1
    timecpu=timecpu+dt

! Get wall clock time and subtract value from the previous call
!  t1wc=mpi_wtime()
!  dtwc=t1wc-t0wc
!  t0wc=t1wc
!  timewc=timewc+dtwc

end subroutine timer

end module initial
