program GTC
    use global_parameters
    use read_in_parameters, only: mstep, ndiag, msnap, irun, track_particles, &
                                  fluid_on, antenna_on, global_solver_on, ifullk, kinetice
    use initial, only: gtc_initial, timer, gtc_final
    use tracking, only: tag, partout
    use data_output, only: diagnosis, snapshot
    use restart, only: restart_io
    use fluid, only: fluid_initial, pushfluid
    use antenna, only: antenna_initial, antenna_profile
    use particle_array, only: zfki
    use equilibrium, only: psirz_Eq
    use spline_function, only: spline2d
    use equilibrium, only: lsr, lsz, spdR, spdZ, zmax
    use pushe_interface, only: pushe
    implicit none
 
    integer i,ierror,isub
    real(doubleprec) timecpu(10)
    character(len=16) cdate(4)
 
!===================================================================
!   Initialize MPI, OpenMP, Adios, PETSc, Random number generator etc
!===================================================================
    call gtc_initial(cdate,timecpu)

!===================================================================
!   Input parameters, equilibrium & profile data, 
!   coefficients for gyroaveraging & interpolation
!===================================================================
    call setup

!===================================================================
!   Initialize particle quantities on simulation grids
!   Load kinetic markers
!===================================================================
    call load_particle
    if(ifullk==0)then
        call locate('gki-thermal')
    else
        call locate('fki-thermal')
    endif
    call locate('electron')
    if(track_particles==1)call tag

!=================================================================== 
!   Initialize fluid quantities on simulation grids 
!    (not be controlled by fluid_on)
!=================================================================== 
    call fluid_initial

!=================================================================== 
!   Initialize external antenna quantities on simulation grids
!    (not be controlled by antenna_on)
!=================================================================== 
    call antenna_initial

!=================================================================== 
!   Start time integrator
!=================================================================== 
    if(mype==0) then
        write(gtcout, *) "============================================="
        write(gtcout, *) "Main Loop Starts:"
        write(gtcout, *) "============================================="
        write(gtcout, *) "" !this is a new line
        flush(gtcout)
    endif

    call timer(timecpu(10),timecpu(7))
    do istep=1,mstep
        do irk=1,2
! idiag=0: do time history diagnosis at irk=1 & istep=ndiag*N
            idiag=mod(irk+1,2)+mod(istep,ndiag)
            if(idiag==0 .and. track_particles==1)call partout
            call timer(timecpu(10),timecpu(8))

            call field_gradient
            if(fluid_on==1)call pushfluid
            if(ifullk==0)then
                call pushi
                call timer(timecpu(10),timecpu(1))
                call shift("gki-thermal")
                call timer(timecpu(10),timecpu(2))
                call locate('gki-thermal')
                call charge('gki-thermal')
                call timer(timecpu(10),timecpu(3))
            else
                call pushfki
                call timer(timecpu(10),timecpu(1))
                call shift("fki-thermal")
                call timer(timecpu(10),timecpu(2))
                call locate('fki-thermal')
                call charge('fki-thermal')
                call timer(timecpu(10),timecpu(3))
            endif
            if(kinetice==1)then
                call pushe
                call shift("electron")
                call locate("electron")
                call charge("electron")
                call timer(timecpu(10),timecpu(4))
            endif
                       
#ifdef _PETSc  
            if(global_solver_on==0)then
#endif
                call local_pssolver
                call timer(timecpu(10),timecpu(6))
#ifdef _PETSc
            elseif(global_solver_on==1)then
                call global_pssolver
                call timer(timecpu(10),timecpu(6))
            else
                if(mype==0)write(gtcout,*)"wrong global_solver_on value"
                stop
            endif
#endif
            if(antenna_on==1)call antenna_profile

            if(idiag==0)call diagnosis
            call timer(timecpu(10),timecpu(8))
        enddo
      !! profile snapshots
        if(msnap/=0)then
            if(mod(istep,mstep/msnap)==0 .or. istep+irun==1)call snapshot
            if(mod(istep,mstep/msnap)==0)call restart_io("write")
            call timer(timecpu(10),timecpu(8))
        endif
    enddo
! Program finalize
    call gtc_final(cdate,timecpu)

end program GTC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

