module fluid
use precision
implicit none

public fluid_initial, pushfluid,&
       sdeltni,deltni,deltni0,suipara,uipara,uipara0, &
       sdeltne,deltne,deltne0,suepara,uepara,uepara0, &
       spi_para,pi_para,pi_para0,&

       graduepara, graduipara,gradpi_para
private

real(lk), dimension(:,:), allocatable :: sdeltni, deltni, deltni0, &
                                         suipara, uipara, uipara0, &
                                         spi_para, pi_para, pi_para0, &
                                         sdeltne, deltne, deltne0, &
                                         suepara, uepara, uepara0
                                           
real(lk), dimension(:,:,:), allocatable :: graduepara, graduipara, gradpi_para
save
contains
!===================================================================
subroutine fluid_initial
    use grids_setup, only: mgrid
    use global_parameters, only: mype, gtcout
    if(mype==0) then
        write(gtcout, *) "============================================="
        write(gtcout, *) "Fluid_initial Starts:"
        write(gtcout, *) "============================================="
        write(gtcout, *) "" !this is a new line
        flush(gtcout)
    endif
    allocate(sdeltni(0:1,mgrid),deltni0(0:1,mgrid),deltni(0:1,mgrid),&
             suipara(0:1,mgrid),uipara0(0:1,mgrid),uipara(0:1,mgrid),&
             sdeltne(0:1,mgrid),deltne0(0:1,mgrid),deltne(0:1,mgrid),&
             suepara(0:1,mgrid),uepara0(0:1,mgrid),uepara(0:1,mgrid),&
             spi_para(0:1,mgrid),pi_para(0:1,mgrid),pi_para0(0:1,mgrid),&
             gradpi_para(3,0:1,mgrid),graduepara(3,0:1,mgrid),graduipara(3,0:1,mgrid))

    sdeltni=0.0_lk
    deltni0=0.0_lk
    deltni=0.0_lk

    suipara=0.0_lk
    uipara0=0.0_lk
    uipara=0.0_lk

    sdeltne=0.0_lk
    deltne0=0.0_lk
    deltne=0.0_lk

    suepara=0.0_lk
    uepara0=0.0_lk
    uepara=0.0_lk

    spi_para=0.0_lk
    pi_para=0.0_lk
    pi_para0=0.0_lk
 
    graduepara=0.0_lk
    graduipara=0.0_lk
    gradpi_para=0.0_lk


end subroutine fluid_initial

!===================================================================
subroutine pushfluid
    use precision
    use global_parameters
    use particle_array
    use grids_fields
    use simulation_domain, only: sfsol
    use grids_setup
    use grids_eqdata
    use read_in_parameters, only: qion, aion, qelectron, aelectron, tstep, min_r
    implicit none
  
    integer :: i,ij,j,ij0
    real(lk) :: dtime,q,perturb,cmratioi,cmratioe,wpara,wstar,wresistivity,wpressure,coll_freq,&
                mass_ion_inv,&
                ni, b_inv, br, bz, dtdr, dtdz, r_inv
  
    perturb=0.0!real(nonlinear)
    cmratioi=qion/aion
    cmratioe=qelectron/aelectron
    mass_ion_inv=1.0_lk/aion
  
    coll_freq=0.01_lk
  
    if(irk==1)then
! 1st step of Runge-Kutta method
        dtime=0.5*tstep
!$omp parallel do private(i,j,ij)
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
                deltni0(1,ij)=deltni(1,ij)
                uipara0(1,ij)=uipara(1,ij)
                pi_para0(1,ij)=pi_para(1,ij)
                deltne0(1,ij)=deltne(1,ij)
                uepara0(1,ij)=uepara(1,ij)
            enddo
        enddo
! 2nd step of Runge-Kutta method
    else
        dtime=tstep
    endif
  
    if(sfsol==1)then
!$omp parallel do private(i,j,ij,ij0,wpara,wstar,wresistivity)
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
                ij0=igrid_sol(i)+j
! parallel compression
                wpara=meshni(ij)*graduipara(3,1,ij)
! wstar
                wstar=meshni(ij)*(-bmeshz_unit(ij)*kapanir(ij)+bmeshr_unit(ij)*kapaniz(ij))*gradphi(2,1,ij)/bmesh(ij)/max(min_r,Rsol(ij0))
                wresistivity=-coll_freq*sdeltni(1,ij)
                deltni(1,ij)=deltni0(1,ij)-dtime*(1.0_lk*wpara + 1.0_lk*wstar + 0.0_lk*wresistivity)
            enddo
        enddo
    endif
  
    sdeltni=deltni
    call toroidal_filter(sdeltni)
  
  
    if(sfsol==1)then
!$omp parallel do private(i,j,ij,wpara,wpressure)
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
! parallel electric field
                wpara=cmratioi*gradphi(3,1,ij)
! parallel compression 
                wpressure=gradpi_para(3,1,ij)*mass_ion_inv/meshni(ij)
                uipara(1,ij)=uipara0(1,ij)-dtime*(1.0_lk*wpara + 1.0_lk*wpressure)
            enddo
        enddo
    endif
    suipara=uipara
    call toroidal_filter(suipara)
  
  
    if(sfsol==1)then
!$omp parallel do private(i,j,ij,ni,b_inv,br,bz,dtdr,dtdz,r_inv)
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
                ij0=igrid_sol(i)+j
                ni=meshni(ij)
                b_inv=1.0_lk/bmesh(ij)
                br=bmeshr_unit(ij)
                bz=bmeshz_unit(ij)
                dtdr=meshti(ij)*kapatir(ij)
                dtdz=meshti(ij)*kapatiz(ij)
                r_inv=1.0_lk/max(min_r,Rsol(ij0))
! omega_temperature_star term
                wstar=ni*b_inv*(-bz*dtdr+br*dtdz)*gradphi(2,1,ij)*r_inv
                pi_para(1,ij) = pi_para0(1,ij)-dtime*(wstar)
          !     pi_para(1,ij)=sdeltni(1,ij)*meshti(ij)
            enddo
        enddo
    endif
    spi_para=pi_para
    call toroidal_filter(spi_para)
  
!==================================================================
    if(sfsol==1)then
!$omp parallel do private(i,j,ij,wpara,wstar)
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
! parallel compression
                wpara=meshne(ij)*graduepara(3,1,ij)
! omega_star
                wstar=meshne(ij)*(-bmeshz_unit(ij)*kapaner(ij)+bmeshr_unit(ij)*kapanez(ij))*gradphi(2,1,ij)/bmesh(ij)/max(min_r,Rsol(ij0))
                deltne(1,ij)=deltne0(1,ij)-dtime*(wpara + wstar)
            enddo
        enddo
    endif  
    sdeltne=deltne
    call toroidal_filter(sdeltne)

    if(sfsol==1)then  
!$omp parallel do private(i,j,ij,wpara)
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
! parallel compression
                wpara=cmratioe*gradphi(3,1,ij)
      
                uepara(1,ij)=uepara0(1,ij)-dtime*wpara
            enddo
        enddo
    endif   
    suepara=uepara
    call toroidal_filter(suepara)
  
end subroutine pushfluid
end module fluid
