! This file is part of GTC version 4.4 
! GTC version 4.4 is released under the 3-Clause BSD license:

! Copyright (c) 2002,2010,2016, GTC Team (team leader: Zhihong Lin, zhihongl@uci.edu)
! All rights reserved.

! Redistribution and use in source and binary forms, with or without 
! modification, are permitted provided that the following conditions are met:

! 1. Redistributions of source code must retain the above copyright notice, 
!    this list of conditions and the following disclaimer.

! 2. Redistributions in binary form must reproduce the above copyright notice, 
!    this list of conditions and the following disclaimer in the documentation 
!    and/or other materials provided with the distribution.

! 3. Neither the name of the GTC Team nor the names of its contributors may be 
!    used to endorse or promote products derived from this software without 
!    specific prior written permission.
! ==============================================================================


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                            !
!                 Gyrokinetic Toroidal Code (GTC)                            !
!                          Version 4.4                                       !
!                               2020                                         !
!                             GTC Team                                       !
!                                                                            !
! Developers:                                                                !
! Z. Lin (zhihongl@uci.edu), I. Holod, Y. Xiao, W.L. Zhang, S. Klasky,       !
! S. Ethier, W.J. Deng, H.S. Zhang, Z.X. Wang, D. Fulton, J. McClenaghan,    !
! A. Kuley, P. Jiang, J. Bao, G. Dong, C.K. Lau, D.J. Liu, D. Spong,         !
! S. Taimourzadeh, X. Cheng, X.S. Wei, Y.Q. Liu, P. Wang, W. Joubert,        ! 
! J.Y. Cheng, B. Wang, Y. Zhao, X.F. Meng, S.H. Wu, Q. Zhang,                !
! W.M. Tang, W.W. Lee, R.B. White, L. Shi, K. S. Fang, H. Y. Wang            !
! P. Liu, G.J. Choi, J.H. Nicolau, W. Wang, J.Y. Fu
!                                                                            !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program GTC
  use system_env
  use precision, only: doubleprec
  use global_parameters, only: nhybrid, istep, mstep, irk, idiag,&
    track_particles, ndata3d, magnetic, ndiag, irun, msnap, ncyclefe, nhybrid,&
    ncyclee, ihybrid, deltabpara, mype, gtcout, ipoincare, antenna_on, numereq
  use particle_array, only: fload, feload, iload, icoll, ecoll
  use interfaces,only:push
#ifdef ADIOS2
  use adios2_common
#endif
  !use equilibrium
  implicit none

  integer i, radial_bin_freq
  real(doubleprec) timecpu(13)
  character(len=16) cdate(4)
  radial_bin_freq = 5

! Initialize MPI, OpenMP, Adios, PETSc, Random number generator etc
  CALL INITIAL(cdate,timecpu)


  if(mype==0) then
    write(gtcout, *) "============================================="
    write(gtcout, *) "Setup Starts:"
    write(gtcout, *) "" !this is a new line
  endif

! input parameters, equilibrium & profile data, coefficients for gyroaveraging & interpolation
  CALL SETUP

  if(mype==0) then
    write(gtcout, *) "============================================="
    write(gtcout, *) "Load Starts:"
    write(gtcout, *) "" !this is a new line
  endif
! initialize particle position and velocity for all species
  call load


  if(mype==0) then
    write(gtcout, *) "============================================="
    write(gtcout, *) "Initial Particle Locating Starts:"
    write(gtcout, *) "" !this is a new line
  endif
! calculate ion gather-scatter coefficients for chargei and pushi
  call locate("thermal-ion")
  if(nhybrid>0)call locate("thermal-electron")
  if(fload>0)call locate("fast-ion")
  if(feload>0)call locate("fast-electron")


  if(mype==0) then
    write(gtcout, *) "============================================="
    write(gtcout, *) "Initial snapshot Starts:"
    write(gtcout, *) "" !this is a new line
  endif
! write initialization snapshot at irun==0
  if(irun==0) call snapshot
  call timer(timecpu(13),timecpu(10))

! main time loop
  if(mype==0) then
    write(gtcout, *) "============================================="
    write(gtcout, *) "Main Loop Starts:"
    write(gtcout, *) "" !this is a new line
  endif
  do istep=1,mstep
     call nvtxStart("step")
     if (mype.eq.0) print *, 'istep=',istep
     flush(output_unit)
     do irk=1,2
! idiag=0: do time history diagnosis at irk=1 & istep=ndiag*N
        idiag=mod(irk+1,2)+mod(istep,ndiag)

        ! gradients of phi, phiind, fluidue, apara; and zonal flow
        call nvtxStart("field_gradient")
        CALL FIELD_GRADIENT
        call nvtxEnd()
        call timer(timecpu(13),timecpu(8))
! push ion
        if(iload/=0)then
           if(iload==9)then
             call pushifk_boris
           else
             call nvtxStart("pushi")
             call push("thermal-ion")
             call nvtxEnd()
           endif
           call timer(timecpu(13),timecpu(1))
           ! redistribute ion across PEs
           call nvtxStart("shifti")
           call shift("thermal-ion")
           call nvtxEnd()
           call timer(timecpu(13),timecpu(2))
           if (mod(istep,radial_bin_freq).eq.0.and.irk.eq.2) then
              call radial_bin("thermal-ion")
           endif
           ! ion perturbed density
           call nvtxStart("locatei")
           call locate("thermal-ion")
           call nvtxEnd()
           call timer(timecpu(13),timecpu(1))
           call nvtxStart("chargei")
           call charge("thermal-ion")
           call nvtxEnd()
           call timer(timecpu(13),timecpu(3))
        endif
! push fast ion
        if(fload>0)then
           call nvtxStart("pushf")
           call push("fast-ion")
           call nvtxEnd()
           call nvtxStart("shiftf")
           call shift("fast-ion")
           call nvtxEnd()
           if (mod(istep,radial_bin_freq).eq.0.and.irk.eq.2) then
              call radial_bin("fast-ion")
           endif
           call nvtxStart("locatef")
           call locate("fast-ion")
           call nvtxEnd()
           call nvtxStart("chargef")
           call charge("fast-ion")
           call nvtxEnd()
           call timer(timecpu(13),timecpu(7))
        endif
! push fast electron
        if(feload>0)then
#ifdef _FRC
          call push("fast-electron")
          call timer(timecpu(13),timecpu(4))
          call shift("fast-electron")
          call timer(timecpu(13),timecpu(5))
          if (mod(istep,radial_bin_freq).eq.0.and.irk.eq.2) then
             call radial_bin("fast-electron")
          endif
          call locate("fast-electron")
          call timer(timecpu(13),timecpu(4))
          call charge("fast-electron")
          call timer(timecpu(13),timecpu(6))
#else
          do i=1,ncyclefe*irk
             call push("fast-electron",i,1)
             call timer(timecpu(13),timecpu(4))
             call shift("fast-electron")
             call timer(timecpu(13),timecpu(5))
             call locate("fast-electron")
             call push("fast-electron",i,2)
             call timer(timecpu(13),timecpu(4))
             call shift("fast-electron")
             call timer(timecpu(13),timecpu(5))
             if (mod(istep,radial_bin_freq).eq.0.and.irk.eq.2) then
                call radial_bin("fast-electron")
             endif
             call locate("fast-electron")
             call timer(timecpu(13),timecpu(4))
          enddo
          call charge("fast-electron")
          call timer(timecpu(13),timecpu(6))
#endif
        endif

        if(magnetic==0)then
! solve GK Poisson equation for phi using adiabatic electron
#ifdef _FRC
           CALL POISSON_FRC("ES-DKE")
#else
! don't solve poisson eq when doing poincare     
        if(ipoincare==0)CALL POISSON_SOLVER("ES-adiabatic")
#endif
        else
           ! push fields of apara, fluidne
           call nvtxStart("pushfield")
           CALL PUSHFIELD
           call nvtxEnd()
           call timer(timecpu(13),timecpu(9))
           ! solver fields of phi, phiind, fluidue
           call nvtxStart("poisson_solver")
           CALL POISSON_SOLVER("EM-hybrid")
           call nvtxEnd()
           call nvtxStart("field_solver")
           CALL FIELD_SOLVER
           call nvtxEnd()
! push field of bpara
           if(deltabpara==1 .and. nhybrid==0)then  
              call bparasolver("adiabatic")
           endif
        endif
        call timer(timecpu(13),timecpu(8))

        do ihybrid=1,nhybrid
! push field of bpara
          if(deltabpara==1)then
            call bparasolver("kinetic")
          endif
          ! time derivative of effective potential, dphieff needs to be re-written for EM with KE
          call nvtxStart("dphieff")
          CALL DPHIEFF
          call nvtxEnd()
          call timer(timecpu(13),timecpu(8))           

           !push electron, sub-cycling
           do i=1,ncyclee*irk
             ! 1st RK step
             call nvtxStart("pushe1")
             call push("thermal-electron",icycle=i,irksub=1,ihybrid=ihybrid)
             call nvtxEnd()
             call timer(timecpu(13),timecpu(4))
             call nvtxStart("shifte1")
             call shift("thermal-electron")
             call nvtxEnd()
             call timer(timecpu(13),timecpu(5))
             call nvtxStart("locatee1")
             call locate("thermal-electron")
             call nvtxEnd()
             ! 2nd RK step
             call nvtxStart("pushe2")
             call push("thermal-electron",icycle=i,irksub=2,ihybrid=ihybrid)
             call nvtxEnd()
             call timer(timecpu(13),timecpu(4))
             call nvtxStart("shifte2")
             call shift("thermal-electron")
             call nvtxEnd()
              call timer(timecpu(13),timecpu(5))

              if (mod(i,radial_bin_freq).eq.0) then
                 call nvtxStart("radial_bin")
                 call radial_bin("thermal-electron")
                 call nvtxEnd()
              endif
              call nvtxStart("locatee2")
              call locate("thermal-electron")
              call nvtxEnd()
              call timer(timecpu(13),timecpu(4))
           enddo

           ! nonadiabatic electron charge density
           call nvtxStart("chargee")
           call charge("thermal-electron")
           call nvtxEnd()
           call timer(timecpu(13),timecpu(6))
! solve GK Poisson equation using non-adiabatic electron
           if (magnetic==0) then
               ! CALL POISSON_SOLVER("ES-hybrid-subtractKineticPart")
                CALL POISSON_SOLVER("ES-hybrid-resolvePoisson")
           endif
           call timer(timecpu(13),timecpu(8))
        enddo
! write diagnostics and 1D data
! when doing poincare,diagnosis is of no use
! and represents a large portion of total computational time 
        if(idiag==0.and.ipoincare==0) then
           call nvtxStart("diagnosis")
           CALL DIAGNOSIS
           call nvtxEnd()
          if(track_particles==1)then
            CALL PARTOUT !write particle data
            call timer(timecpu(13),timecpu(11))
          endif
        endif
        call timer(timecpu(13),timecpu(11))
     enddo

! lorentz and like species collisions
! ions
      if(icoll>0 .and. mod(istep,icoll)==0)then
        call fokkerplanck('thermal-ion')
        call timer(timecpu(13),timecpu(1))
      endif
#ifdef _FRC
! DKE
      if(ecoll>0 .and. mod(istep,ecoll))then
        call lorentz('fast-electron')
        call timer(timecpu(13),timecpu(4))
       !call fokkerplanck('fast-electron')
      endif
#else
! hybrid electrons
     if(ecoll>0 .and. mod(istep,ecoll)==0)then
        call lorentz('thermal-electron')
        call timer(timecpu(13),timecpu(4))
     !  call fokkerplanck('thermal-electron')
     endif
#endif     
! profile snapshots, write particle information to restart file
     if(msnap/=0)then
       if(mod(istep,mstep/msnap)==0 .or. istep+irun==1) then
           CALL SNAPSHOT
          !CALL PHISHOT
           if(ndata3d==1) call DATAOUT3D !write 3D data
       endif
       
       if(mod(istep,mstep/msnap)==0)CALL RESTART_IO("write")
       call timer(timecpu(13),timecpu(11))
    endif
    call nvtxEnd()
  enddo

! Program finalize
#ifdef ADIOS2
  call adios2_finalize_common()
#endif
 CALL FINAL(cdate,timecpu)

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