!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine setup
    use precision
    use global_parameters
    use particle_array
    use read_in_parameters, only: track_particles,ngyroi
    use equilibrium, only: eqdata
    use grids_setup, only: mgrid_sol, mgrid_core
    implicit none

    integer ierror
    real(lk) randNumCheck
    namelist /key_parameters/ numberpe,nspecies,nparami,nparame
    if(mype==0) then
        write(gtcout, *) "============================================="
        write(gtcout, *) "Setup Starts:"
        write(gtcout, *) "============================================="
        write(gtcout, *) "" !this is a new line
        flush(gtcout)
    endif
    
    call eqdata
    call grids_initial
    call grids_mpi_initial

! number of particle  
    misol=int(real(mi)*real(mgrid_sol)/real(npartdom))
    mesol=int(real(me)*real(mgrid_sol)/real(npartdom))
    micore=int(real(mi)*real(mgrid_core)/real(npartdom))
    mecore=int(real(me)*real(mgrid_core)/real(npartdom))
    mi=micore+misol
    mfki=mi
    me=mecore+mesol
    nspecies=2
! particle attributes: 1: R; 2: \zeta; 3: Z; 4: v_{||}; 5: weight=\delta f/g; 6: \sqrt\mu;
!                      7: \psi; 8: energy; 9: actual zeta in partial torus simulation
    nparami=7
    nparame=7

! Particles is tracked by tagging them with a number with an extra element to the particle array
    if(track_particles == 1)then
        nparami=nparami+2
        nparame=nparame+2
    endif

!write test random number to check RNG initialization
    call random_number(randNumCheck)
    if(mype==0)write(gtcout,*)'Rand Num Check = ',randNumCheck

! write out key parameter
    if(mype == 0) then
        write(gtcout,key_parameters)
        call FLUSH(gtcout)
    endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! test gyro-averaging initialization
! gyro-averaging for sqrt(mu)=rho0 on grid of magnetic coordinates
    call gyroinitial(ngyroi)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    call MPI_BARRIER(MPI_COMM_WORLD,ierror)

end subroutine setup
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine grids_initial
    use precision
    use global_parameters, only: shear_rate,mype,gtcout, min_b
    use grids_setup, only: frc_mesh_create
    use grids_eqdata, only: frc_mesh_equilibrium_field_initial
    use grids_fields, only: frc_mesh_perturb_field_initial
    use mapping_splines, only: frc_coordinate_spline, er_init
    use simulation_domain, only: frc_domain_setup
    use laplacian, only: lapperp_initial
    implicit none

    min_b=0.01_lk
    shear_rate=0.0_lk!-7.0e-4_lk
    if(mype==0)then
        write(gtcout,*)"********************************************"
        write(gtcout,*)"domain_setup starts"
        write(gtcout,*)""
    endif
!Determine the simulation domain (SOL, Core and SOLCore coupled)
    call frc_domain_setup

    if(mype==0)then
        write(gtcout,*)"********************************************"
        write(gtcout,*)"coordinate_spline starts"
        write(gtcout,*)""
    endif
! S(\psi,Z) and S(\psi,\theta) are created based on simulation domain
    call frc_coordinate_spline   

    if(mype==0)then
        write(gtcout,*)"********************************************"
        write(gtcout,*)"er_init starts"
        write(gtcout,*)""
    endif
    call er_init

    if(mype==0)then
        write(gtcout,*)"********************************************"
        write(gtcout,*)"mesh_create starts"
        write(gtcout,*)""
    endif
! Field aligned mesh_create is created
    call frc_mesh_create

    if(mype==0)then
        write(gtcout,*)"********************************************"
        write(gtcout,*)"mesh_equilibrium_field_initial starts"
        write(gtcout,*)""
    endif
    call frc_mesh_equilibrium_field_initial

    if(mype==0)then
        write(gtcout,*)"********************************************"
        write(gtcout,*)"mesh_perturb_field_initial starts"
        write(gtcout,*)""
    endif
    call frc_mesh_perturb_field_initial

    if(mype==0)then
        write(gtcout,*)"********************************************"
        write(gtcout,*)"lapperp_initial starts"
        write(gtcout,*)""
    endif
    call lapperp_initial

end subroutine grids_initial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine field_gradient
    use precision
    use grids_fields, only: phi, phi0,gradphi0, gradphi, gradient
    use fluid, only: suepara, suipara, spi_para, &
                     graduepara, graduipara,gradpi_para
    use particle_array, only: flowi, flowe, flowfki
    use read_in_parameters, only:tstep, ifullk,izonal
    use global_parameters, only:istep, zeta0, zeta1

    implicit none
    real(lk) :: o_ext
    
    !phi(1,:)=1.0e-4*cos(zeta1-tstep*istep*0.01)
! grad phi_electrostatic
    call gradient(phi(1,:),gradphi)
#ifndef GPU_UM
!$acc update device(gradphi)
#endif
    if(izonal==1)then
      call gradient(phi0(1,:),gradphi0)
    endif
#ifndef GPU_UM
!$acc update device(gradphi0)
#endif
    if(ifullk==0)then
    call gradient(flowi(1,:),graduipara)
    else
    call gradient(flowfki(1,:),graduipara)
    endif
!  call gradient(flowe(1,:),graduepara)

!  call gradient(suepara(1,:),graduepara)
!  call gradient(suipara(1,:),graduipara)
!  call gradient(spi_para(1,:),gradpi_para)

end subroutine field_gradient
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine gyroinitial(ngyro)
    use precision
    use spline_function, only: spline2d

!    use grids_setup, only: mgrid, mpsi_core, mpsi_sol
    use grids_setup
    use simulation_domain, only: sfcore, sfsol, couplefrc
    use mapping_splines, only: rspsol_psi_s, zspsol_psi_s, sspsol_psi_z, msp_psi_sol, aspsol_psi, &
                               lspsol, lszsol, del_psi_sol, del_z_sol, del_s_sol

    use read_in_parameters, only: Zsol_domain0, Zsol_domain1   !, ngyroi
    use global_parameters, only: rho0,mype
    use grids_eqdata, only: bmesh

    use gyro_array, only: deltap, pgyro, zetagyro, tgyro, pgyro2, tgyro2
    implicit none

! dummy parameters
    integer :: ngyro

! local parameters
    integer :: i,j,ij,ijsol,ngyromax,igyro   !,ngyromax0
    real(lk) :: bf,rhos,delr,delz,pdum0,pdum1,delpsi,sdum0,sdum1,dels

    ! for N-point gyro-averaging, N=1,4,8 for now.
    ! ngyro = 2   !ngyroi
    ngyromax = min(8,max(2,ngyro))
    allocate(deltap(mpsi),pgyro(ngyromax,mgrid),zetagyro(ngyromax,mgrid),tgyro(ngyromax,mgrid),pgyro2(ngyromax,mgrid),tgyro2(ngyromax,mgrid))   !, STAT=mstat)
!    call check_allocation(mstat, "pgyro in poisson.F90")
    pgyro = 0.0_lk
    tgyro = 0.0_lk
    zetagyro = 0.0_lk

    if(sfsol==1)then
        do i=2,mpsi_sol
            deltap(mpsi_core+i)= psimesh_sol(i)-psimesh_sol(i-1)
        enddo
    endif

! Starting from outermost point, counter-clockwise: 1,2,3,4.
! Calculate coeffecients here for rho_local=rho0*sqrt(2B_0/B), _0 represent quantities on axis.
! In subroutines poisson and locatei: rho_particle=rho_local*sqrt(m_i/m_proton)*sqrt(mu)/q_i.

    if(sfsol==1)then
        
        do i=1,mpsi_sol
            ! psi index of the grid
            pdum0=psimesh_sol(i)
            if(i==mpsi_sol)then
                pdum1=psimesh_sol(i-1)
            else
                pdum1=psimesh_sol(i+1)
            endif
            ! z index of the grid
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
                ijsol=igrid_sol(i)+j

            !!! compute distance between two radial grids along constant S
                sdum0=Ssol(ijsol)
                ! sdum=real(j-1)/real(mpara_sol(i)-1)

                ! delr=R(pdum1,sdum0)-R(pdum0,sdum0); delz=Z(pdum1,sdum0)-Z(pdum1,sdum0)
                delr=spline2d(pdum1-msp_psi_sol(1),sdum0,0,lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)   &
                         - spline2d(pdum0-msp_psi_sol(1),sdum0,0,lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)
                delz=spline2d(pdum1-msp_psi_sol(1),sdum0,0,lspsol,lszsol,del_psi_sol,del_s_sol,zspsol_psi_s,0,0)   &
                     -spline2d(pdum0-msp_psi_sol(1),sdum0,0,lspsol,lszsol,del_psi_sol,del_s_sol,zspsol_psi_s,0,0)
                ! Ssol(ij)=spline2d(psimesh_sol(i)-msp_psi_sol(1),Zsol(ij)-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)

            !!! compute change in radial flux delta_psi
                bf=bmesh(ij)
                rhos=rho0*sqrt(2.0_lk/bf) ! gyroradius of thermal proton with mu*B=T_e0

                delpsi=rhos*deltap(mpsi_core+min(mpsi_sol,i+1))/sqrt(delr**2+delz**2)

            !!! first two points perpendicular to field line on poloidal surface
                pgyro(1,ij)=delpsi*0.7071067812_lk
                pgyro(2,ij)=0.0_lk-pgyro(1,ij)
                !pgyro(1,ij)=delpsi
                !pgyro(2,ij)=0.0_lk-pgyro(1,ij)
                !pgyro(3,ij)=0.0_lk   !delpsi
                !pgyro(4,ij)=0.0_lk   !-pgyro(1,ij)
                ! non-orthorgonality between psi and theta
                tgyro(1,ij)=rhos !/sqrt(delr**2+delz**2) !0.0_lk
                tgyro(2,ij)=0.0_lk-tgyro(1,ij)
                !tgyro(3,ij)=0.0_lk   !rhos !/sqrt(delr**2+delz**2) !0.0_lk
                !tgyro(4,ij)=0.0_lk   !-tgyro(1,ij)
                ! gyromotion happens in psi-zeta plane
                zetagyro(1,ij)=rhos*0.7071067812_lk
                zetagyro(2,ij)=0.0_lk-zetagyro(1,ij)
                !zetagyro(1,ij)=0.0_lk
                !zetagyro(2,ij)=0.0_lk
                !zetagyro(3,ij)=rhos
                !zetagyro(4,ij)=-zetagyro(3,ij)
#ifndef GPU_UM
                !$acc update device(pgyro, zetagyro)
#endif

                if(mype==0 .and. ij==int(mpsi_sol/2))then
                    write(*,*)'!!!!!!!!!!!!! Gyro initial !!!!!!!!!!!!!'
                    write(*,*)'psi_diag =',pdum0,'r(psi_diag) =',spline2d(pdum0-msp_psi_sol(1),sdum0,0,lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)
                    write(*,*)'bf =',bf,'rho0 =',rho0,'rhos =',rhos
                    write(*,*)'deltap(i+1) =',deltap(mpsi_core+min(mpsi_sol,i+1)),'pgyro(delpsi) =',delpsi
                    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                endif
            enddo
        enddo

    endif

end subroutine gyroinitial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine grids_mpi_initial
use precision
use global_parameters
use simulation_domain
use equilibrium
use grids_eqdata
implicit none

integer numtmp,avr,res

ij_min_core=0
ij_min_sol=0
ij_max_core=0
ij_max_sol=0
ij_min_tor_core=0
ij_max_tor_core=0
ij_min_tor_sol=0
ij_max_tor_sol=0
if(sfcore==1)then
    numtmp = mgrid_core + nproc_partd - 1
    avr = numtmp/nproc_partd-1
    res = mgrid_core - avr * nproc_partd
    if(myrank_partd<res)then
      ij_min_core = myrank_partd*(avr+1)+1
      ij_max_core = ij_min_core + avr
    else
      ij_min_core = myrank_partd*avr + 1 + res
      ij_max_core = ij_min_core + avr - 1
    endif
    mygrid_core=ij_max_core-ij_min_core+1
    numtmp = mygrid_core + mtoroidal - 1
    avr = numtmp/mtoroidal-1
    res = mygrid_core - avr * mtoroidal
    if(myrank_toroidal<res)then
      ij_min_tor_core = ij_min_core + myrank_toroidal*(avr+1)
      ij_max_tor_core = ij_min_tor_core + avr
    else
      ij_min_tor_core = ij_min_core + myrank_toroidal*avr + res
      ij_max_tor_core = ij_min_tor_core + avr - 1
    endif
    mygrid_tor_core = ij_max_tor_core - ij_min_tor_core + 1
endif
if(sfsol==1)then
    numtmp = mgrid_sol + nproc_partd - 1
    avr = numtmp/nproc_partd-1
    res = mgrid_sol - avr * nproc_partd
    if(myrank_partd<res)then
      ij_min_sol = mgrid_core + myrank_partd*(avr+1)+1
      ij_max_sol = ij_min_sol + avr
    else
      ij_min_sol = mgrid_core + myrank_partd*avr + 1 + res
      ij_max_sol = ij_min_sol + avr -1
    endif
    mygrid_sol=ij_max_sol-ij_min_sol+1
    numtmp = mygrid_sol + mtoroidal - 1
    avr = numtmp/mtoroidal-1
    res = mygrid_sol - avr * mtoroidal
    if(myrank_toroidal<res)then
      ij_min_tor_sol = ij_min_sol + myrank_toroidal*(avr+1)
      ij_max_tor_sol = ij_min_tor_sol + avr
    else
      ij_min_tor_sol = ij_min_sol + myrank_toroidal*avr + res
      ij_max_tor_sol = ij_min_tor_sol + avr - 1
    endif
endif


end subroutine grids_mpi_initial
