processing_mod.F90

MODULE processing
   USE prec_const,  ONLY: xp, imagu, SQRT2, SQRT3, onetwelfth, twothird
   USE grid,        ONLY: &
      local_na, local_np, local_nj, local_nky, local_nkx, local_nz, Ngz,Ngj,Ngp,nzgrid, &
      parray,pmax,ip0, iodd, ieven, total_nz, &
      CONTAINSp0,ip1,CONTAINSp1,ip2,CONTAINSp2,ip3,CONTAINSp3,&
      jarray,jmax,ij0, dmax,&
      kyarray, AA_y,&
      kxarray, AA_x,&
      zarray, zweights_SR, ieven, iodd, inv_deltaz,&
      kroneck_p0, kroneck_p1
   USE fields,           ONLY: moments, phi, psi
   USE array,            ONLY : kernel, nadiab_moments, &
      ddz_napj, ddzND_Napj, interp_napj,&
      dens, upar, uper, Tpar, Tper, temp
   USE geometry,         ONLY: Jacobian, iInt_Jacobian
   USE time_integration, ONLY: updatetlevel
   USE calculus,         ONLY: simpson_rule_z, grad_z, grad_z_5D, grad_z2, grad_z4, grad_z4_5D, interp_z
   USE model,            ONLY: EM, beta, HDz_h
   USE species,          ONLY: tau,q_tau,q_o_sqrt_tau_sigma,sqrt_tau_o_sigma
   USE parallel,         ONLY: num_procs_ky, rank_ky, comm_ky
   USE mpi
   implicit none

   REAL(xp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: pflux_x, gflux_x
   REAL(xp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: hflux_x
   INTEGER :: ierr
   PUBLIC :: init_process
   PUBLIC :: compute_nadiab_moments, compute_gradients_z, compute_interp_z
   PUBLIC :: compute_density, compute_upar, compute_uperp
   PUBLIC :: compute_Tpar, compute_Tperp, compute_fluid_moments
   PUBLIC :: compute_radial_transport
   PUBLIC :: compute_radial_heatflux
   PUBLIC :: compute_Napjz_spectrum
CONTAINS

   SUBROUTINE init_process
      USE grid,       ONLY: local_na
      IMPLICIT NONE
      ALLOCATE( pflux_x(local_na))
      ALLOCATE( gflux_x(local_na))
      ALLOCATE( hflux_x(local_na))
   END SUBROUTINE init_process
   
!------------------------------ HIGH FREQUENCY ROUTINES -------------
! The following routines (nadiab computing, interp and grads) must be
! as fast as possible since they are called every RK step
   ! evaluate the non-adiabatique ion moments
   !
   ! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0
   !
   SUBROUTINE compute_nadiab_moments
      IMPLICIT NONE
      INTEGER :: ia,ij
      !non adiab moments
      ! default : same as moments
      nadiab_moments(:,:,:,:,:,:) = moments(:,:,:,:,:,:,updatetlevel)
      ! add phi and psi contribution
      IF (CONTAINSp0) THEN
         DO ij=1,local_nj+ngj
         DO ia = 1,local_na
            nadiab_moments(ia,ip0,ij,:,:,:) = moments(ia,ip0,ij,:,:,:,updatetlevel) &
            +q_tau(ia)*kernel(ia,ij,:,:,:,ieven)*phi(:,:,:)
         ENDDO 
         ENDDO
      ENDIF
      IF (CONTAINSp1 .AND. (beta .GT. 0)) THEN
         DO ij=1,local_nj+ngj
         DO ia = 1,local_na
            nadiab_moments(ia,ip1,ij,:,:,:) = moments(ia,ip1,ij,:,:,:,updatetlevel) &
               -q_o_sqrt_tau_sigma(ia)*kernel(ia,ij,:,:,:,iodd)*psi(:,:,:)
         ENDDO 
         ENDDO
      ENDIF
   END SUBROUTINE compute_nadiab_moments

   ! ! z grid gradients
   SUBROUTINE compute_gradients_z
      IMPLICIT NONE
      INTEGER :: eo, p_int, ia,ip,ij,iky,ikx
      COMPLEX(xp), DIMENSION(local_nz+ngz) :: f_in
      COMPLEX(xp), DIMENSION(local_nz)     :: f_out
   IF(total_nz .GT. 4) THEN
      DO ikx = 1,local_nkx
      DO iky = 1,local_nky
      DO ij = 1,local_nj+ngj
      DO ip = 1,local_np+ngp
      DO ia = 1,local_na
         IF(nzgrid .GT. 1) THEN
            p_int = parray(ip+ngp/2)
            eo    = MODULO(p_int,2)+1 ! Indicates if we are on even or odd z grid
         ELSE
            eo = 0
         ENDIF
         ! Compute z first derivative
         f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
         CALL   grad_z(eo,local_nz,ngz,inv_deltaz,f_in,f_out)
         ddz_napj(ia,ip,ij,iky,ikx,:) = f_out
         ! Parallel numerical diffusion
         IF (.NOT. HDz_h) &
            f_in = moments(ia,ip,ij,iky,ikx,:,updatetlevel)
         CALL  grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
         ! get output
         ddzND_Napj(ia,ip,ij,iky,ikx,:) = f_out
      ENDDO
      ENDDO
      ENDDO
      ENDDO
      ENDDO
   ELSE ! 2D Z-pinch
      ddz_napj   = 0._xp
      ddzND_Napj = 0._xp
   ENDIF
   END SUBROUTINE compute_gradients_z

      ! z grid interpolation
   SUBROUTINE compute_interp_z
      IMPLICIT NONE
      INTEGER :: eo, ia,ip,ij,iky,ikx,iz
      COMPLEX(xp), DIMENSION(local_nz+ngz) :: f_in
      COMPLEX(xp), DIMENSION(local_nz)     :: f_out
      IF(nzgrid .GT. 1) THEN
         DO ikx = 1,local_nkx
         DO iky = 1,local_nky
         DO ij = 1,local_nj+ngj
         DO ip = 1,local_np+ngp
         DO ia = 1,local_na
            ! Compute even odd grids interpolation
            f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
            CALL interp_z(eo,local_nz,ngz,f_in,f_out)
            DO iz = 1,local_nz
               interp_napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
            ENDDO
         ENDDO
         ENDDO
         ENDDO
         ENDDO
         ENDDO
      ELSE
         interp_napj(:,:,:,:,:,1:local_nz) = nadiab_moments(:,:,:,:,:,(1+ngz/2):(local_nz+ngz/2))
      ENDIF
   END SUBROUTINE compute_interp_z

   !--------------------- LOW FREQUENCY PROCESSING ROUTINES -------------------!
   ! The following routines are called by the diagnose routine every nsave3D steps 
   ! (does not need to be optimized)
   ! 1D diagnostic to compute the average radial particle transport <n_a v_ExB_x>_xyz
   SUBROUTINE compute_radial_transport
      IMPLICIT NONE
      COMPLEX(xp) :: pflux_local, gflux_local, integral
      REAL(xp)    :: buffer(2)
      INTEGER     :: i_, root, iky, ikx, ia, iz, in, izi, ini
      COMPLEX(xp), DIMENSION(local_nz) :: integrant
      DO ia = 1,local_na
         pflux_local = 0._xp ! particle flux
         gflux_local = 0._xp ! gyrocenter flux
         integrant   = 0._xp ! auxiliary variable for z integration
         !!---------- Gyro center flux (drift kinetic) ------------
         ! Electrostatic part
         IF(CONTAINSp0) THEN
            DO iz = 1,local_nz
            izi = iz + ngz/2 !interior index for ghosted arrays
            DO ikx = 1,local_nkx
            DO iky = 1,local_nky
               integrant(iz) = integrant(iz) &
               +Jacobian(izi,ieven)*moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
               *imagu*kyarray(iky)*CONJG(phi(iky,ikx,izi))
            ENDDO
            ENDDO
            ENDDO
         ENDIF
         ! Electromagnetic part
         IF( EM .AND. CONTAINSp1 ) THEN
            DO iz = 1,local_nz ! we take interior points only
            izi = iz + ngz/2 !interior index for ghosted arrays
            DO ikx = 1,local_nkx
            DO iky = 1,local_nky
               integrant(iz) = integrant(iz)&
                  -Jacobian(izi,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,ij0,iky,ikx,izi,updatetlevel)&
                  *imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))
            ENDDO
            ENDDO
            ENDDO
         ENDIF
         ! Integrate over z
         call simpson_rule_z(local_nz,zweights_SR,integrant,integral)
         ! Get process local gyrocenter flux with a factor two to account for the negative ky modes
         gflux_local = 2._xp*integral*iInt_Jacobian
         !
         !!---------- Particle flux (gyrokinetic) ------------
         integrant   = 0._xp ! reset auxiliary variable
         ! Electrostatic part
         IF(CONTAINSp0) THEN
            DO iz = 1,local_nz ! we take interior points only
            izi = iz + ngz/2 !interior index for ghosted arrays
            DO ikx = 1,local_nkx
            DO iky = 1,local_nky
            DO in = 1, local_nj
            ini = in + ngj/2 !interior index for ghosted arrays
               integrant(iz) = integrant(iz)+ &
                  Jacobian(izi,ieven)*moments(ia,ip0,ini,iky,ikx,izi,updatetlevel)&
                  *imagu*kyarray(iky)*kernel(ia,ini,iky,ikx,izi,ieven)*CONJG(phi(iky,ikx,izi))
            ENDDO
            ENDDO
            ENDDO
            ENDDO
         ENDIF
         ! Electromagnetic part
         IF( EM .AND. CONTAINSp1 ) THEN
            DO iz = 1,local_nz ! we take interior points only
            izi = iz + ngz/2 !interior index for ghosted arrays
            DO ikx = 1,local_nkx
            DO iky = 1,local_nky
            DO in = 1, local_nj
            ini = in + ngj/2 !interior index for ghosted arrays
                  integrant(iz) = integrant(iz) - &
                  Jacobian(izi,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,ini,iky,ikx,izi,updatetlevel)&
                  *imagu*kyarray(iky)*kernel(ia,ini,iky,ikx,izi,iodd)*CONJG(psi(iky,ikx,izi))
            ENDDO
            ENDDO
            ENDDO
            ENDDO
         ENDIF
         ! Integrate over z
         call simpson_rule_z(local_nz,zweights_SR,integrant,integral)
         ! Get process local particle flux with a factor two to account for the negative ky modes
         pflux_local = 2._xp*integral*iInt_Jacobian
         !!!!---------- Sum over all processes ----------
         buffer(1) = REAL(gflux_local,xp)
         buffer(2) = REAL(pflux_local,xp)
         root = 0
         !Gather manually among the rank_p=0 processes and perform the sum
         gflux_x(ia) = 0
         pflux_x(ia) = 0
         IF (num_procs_ky .GT. 1) THEN
            !! Everyone sends its local_sum to root = 0
            IF (rank_ky .NE. root) THEN
               CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
            ELSE
               ! Recieve from all the other processes
               DO i_ = 0,num_procs_ky-1
                  IF (i_ .NE. rank_ky) &
                     CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
                  gflux_x(ia) = gflux_x(ia) + buffer(1)
                  pflux_x(ia) = pflux_x(ia) + buffer(2)
               ENDDO
            ENDIF
         ELSE
            gflux_x(ia) = REAL(gflux_local,xp)
            pflux_x(ia) = REAL(pflux_local,xp)
         ENDIF
      ENDDO
   END SUBROUTINE compute_radial_transport

! 1D diagnostic to compute the average radial particle heatflux <T_i v_ExB_x>_xyz
   SUBROUTINE compute_radial_heatflux
      IMPLICIT NONE
      COMPLEX(xp) :: hflux_local, integral
      REAL(xp)    :: buffer(2), n_xp
      INTEGER     :: i_, root, in, ia, iky, ikx, iz, izi, ini
      COMPLEX(xp), DIMENSION(local_nz) :: integrant        ! charge density q_a n_a
      DO ia = 1,local_na
         hflux_local = 0._xp ! heat flux
         integrant   = 0._xp ! z integration auxiliary variable
         !!----------------ELECTROSTATIC CONTRIBUTION---------------------------
         IF(CONTAINSp0 .AND. CONTAINSp2) THEN
            ! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Na00 * phi
            DO iz = 1,local_nz ! we take interior points only
            izi = iz + ngz/2 !interior index for ghosted arrays
            DO ikx = 1,local_nkx
            DO iky = 1,local_nky
            DO in = 1, local_nj
            ini  = in+ngj/2 !interior index for ghosted arrays
            n_xp = jarray(ini)
               integrant(iz) = integrant(iz) &
                  -Jacobian(izi,ieven)*tau(ia)*imagu*kyarray(iky)*phi(iky,ikx,izi)&
                  *kernel(ia,ini,iky,ikx,izi,ieven)*CONJG(&
                              0.5_xp*SQRT2*moments(ia,ip2,ini  ,iky,ikx,izi,updatetlevel)&
                  +(2._xp*n_xp + 1.5_xp)*moments(ia,ip0,ini  ,iky,ikx,izi,updatetlevel)&
                           -(n_xp+1._xp)*moments(ia,ip0,ini+1,iky,ikx,izi,updatetlevel)&
                                    -n_xp*moments(ia,ip0,ini-1,iky,ikx,izi,updatetlevel))
            ENDDO
            ENDDO
            ENDDO
            ENDDO
         ELSEIF(CONTAINSp0) THEN
            ERROR STOP "Degrees p=0 and p=2 should be owned by the same process"
         ENDIF
         IF(EM .AND. CONTAINSp1 .AND. CONTAINSp3) THEN
            !!----------------ELECTROMAGNETIC CONTRIBUTION--------------------
            DO iz  = 1,local_nz
            izi = iz + ngz/2 !interior index for ghosted arrays
            DO ikx = 1,local_nkx
            DO iky = 1,local_nky
            DO in = 1, local_nj
            ini = in + ngj/2 !interior index for ghosted arrays
            n_xp = jarray(ini)
               integrant(iz) = integrant(iz) &
                     +Jacobian(izi,iodd)*tau(ia)*sqrt_tau_o_sigma(ia)*imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))&
                  *kernel(ia,ini,iky,ikx,izi,iodd)*(&
                  0.5_xp*SQRT2*SQRT3*moments(ia,ip3,ini  ,iky,ikx,izi,updatetlevel)&
                              +1.5_xp*moments(ia,ip1,ini  ,iky,ikx,izi,updatetlevel)&
                  +(2._xp*n_xp+1._xp)*moments(ia,ip1,ini  ,iky,ikx,izi,updatetlevel)&
                        -(n_xp+1._xp)*moments(ia,ip1,ini+1,iky,ikx,izi,updatetlevel)&
                                 -n_xp*moments(ia,ip1,ini-1,iky,ikx,izi,updatetlevel))
            ENDDO
            ENDDO
            ENDDO
            ENDDO
         ENDIF
         ! Add polarisation contribution
         ! integrant(iz) = integrant(iz) + tau_i*imagu*ky_&
         ! *CONJG(phi(iky,ikx,iz))*phi(iky,ikx,iz) * HF_phi_correction_operator(iky,ikx,iz)
         ! Integrate over z
         call simpson_rule_z(local_nz,zweights_SR,integrant,integral)
         ! Double it for taking into account the other half plane
         hflux_local = 2._xp*integral*iInt_Jacobian
         buffer(2)   = REAL(hflux_local,xp)
         root = 0
         !Gather manually among the rank_p=0 processes and perform the sum
         hflux_x(ia) = 0
         IF (num_procs_ky .GT. 1) THEN
            !! Everyone sends its local_sum to root = 0
            IF (rank_ky .NE. root) THEN
               CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
            ELSE
               ! Recieve from all the other processes
               DO i_ = 0,num_procs_ky-1
                  IF (i_ .NE. rank_ky) &
                     CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
                  hflux_x(ia) = hflux_x(ia) + buffer(2)
               ENDDO
            ENDIF
         ELSE
            hflux_x(ia) = REAL(hflux_local,xp)
         ENDIF
      ENDDO
   END SUBROUTINE compute_radial_heatflux

   SUBROUTINE compute_Napjz_spectrum
      USE fields, ONLY : moments
      USE array,  ONLY : Napjz
      USE time_integration, ONLY : updatetlevel
      IMPLICIT NONE
      REAL(xp), DIMENSION(local_np,local_nj,local_nz) :: local_sum, global_sum, buffer
      INTEGER  :: i_, root, count, ia, ip, ij, iky, ikx, iz
      root = 0
      DO ia=1,local_na
         ! z-moment spectrum
         ! build local sum
         local_sum = 0._xp
         DO iz  = 1,local_nz
         DO ikx = 1,local_nkx
         DO iky = 1,local_nky
         DO ij  = 1,local_nj
         DO ip  = 1,local_np
            local_sum(ip,ij,iz)  = local_sum(ip,ij,iz)  + &
               REAL(moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel) &
               * CONJG(moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel)),xp)
         ENDDO
         ENDDO
         ENDDO
         ENDDO
         ENDDO
         ! sum reduction
         buffer     = local_sum
         global_sum = 0._xp
         count = local_np*local_nj*local_nz
         IF (num_procs_ky .GT. 1) THEN
            !! Everyone sends its local_sum to root = 0
            IF (rank_ky .NE. root) THEN
               CALL MPI_SEND(buffer, count , MPI_DOUBLE_PRECISION, root, 5678, comm_ky, ierr)
            ELSE
               ! Recieve from all the other processes
               DO i_ = 0,num_procs_ky-1
                  IF (i_ .NE. rank_ky) &
                     CALL MPI_RECV(buffer, count , MPI_DOUBLE_PRECISION, i_, 5678, comm_ky, MPI_STATUS_IGNORE, ierr)
                  global_sum = global_sum + buffer
               ENDDO
            ENDIF
         ELSE
            global_sum = local_sum
         ENDIF
         Napjz(ia,:,:,:) = global_sum
      ENDDO
   END SUBROUTINE compute_Napjz_spectrum

!_____________________________________________________________________________!
!!!!! FLUID MOMENTS COMPUTATIONS !!!!!
! Compute the 2D particle density for electron and ions (sum over Laguerre)
   SUBROUTINE compute_density
   IMPLICIT NONE
   COMPLEX(xp) :: dens_
   INTEGER :: ia, iz, iky, ikx, ij
   DO ia=1,local_na
   IF ( CONTAINSp0 ) THEN
   ! Loop to compute dens_i = sum_j kernel_j Ni0j at each k
   DO iz = 1,local_nz
   DO iky = 1,local_nky
   DO ikx = 1,local_nkx
      dens_ = 0._xp
      DO ij = 1, local_nj
         dens_ = dens_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) * moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
      ENDDO
      dens(ia,iky,ikx,iz) = dens_
   ENDDO
   ENDDO
   ENDDO
   ENDIF
   ENDDO
   END SUBROUTINE compute_density

! Compute the 2D particle fluid perp velocity for electron and ions (sum over Laguerre)
   SUBROUTINE compute_uperp
      IMPLICIT NONE
      COMPLEX(xp) :: uperp_
      INTEGER :: ia, iz, iky, ikx, ij
      DO ia=1,local_na
      IF ( CONTAINSp0 ) THEN
      DO iz = 1,local_nz
      DO iky = 1,local_nky
      DO ikx = 1,local_nkx
      uperp_ = 0._xp
      DO ij = 1, local_nj
         uperp_ = uperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) *&
            0.5_xp*(moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
                     -moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
      ENDDO
      uper(ia,iky,ikx,iz) = uperp_
      ENDDO
      ENDDO
      ENDDO
      ENDIF
      ENDDO
   END SUBROUTINE compute_uperp

! Compute the 2D particle fluid par velocity for electron and ions (sum over Laguerre)
   SUBROUTINE compute_upar
      IMPLICIT NONE
      INTEGER :: ia, iz, iky, ikx, ij
      COMPLEX(xp) :: upar_
      DO ia=1,local_na
      IF ( CONTAINSp1 ) THEN
      DO iz = 1,local_nz
      DO iky = 1,local_nky
      DO ikx = 1,local_nkx
         upar_ = 0._xp
         DO ij = 1, local_nj
            upar_ = upar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,iodd)*moments(ia,ip1,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
         ENDDO
         upar(ia,iky,ikx,iz) = upar_
      ENDDO
      ENDDO
      ENDDO
      ENDIF
      ENDDO
   END SUBROUTINE compute_upar

! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
   SUBROUTINE compute_tperp
      USE time_integration, ONLY : updatetlevel
      IMPLICIT NONE
      REAL(xp)    :: j_xp
      COMPLEX(xp) :: Tperp_
      INTEGER     :: ia, iz, iky, ikx, ij
      DO ia=1,local_na
      IF ( CONTAINSp0 .AND. CONTAINSp2 ) THEN
      ! Loop to compute T = 1/3*(Tpar + 2Tperp)
      DO iz = 1,local_nz
      DO iky = 1,local_nky
      DO ikx = 1,local_nkx
      Tperp_ = 0._xp
      DO ij = 1, local_nj
         j_xp = jarray(ij+ngj/2)
         Tperp_ = Tperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
            ((2_xp*j_xp+1)*moments(ia,ip0,ij  +ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
                     -j_xp*moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
                  -(j_xp+1)*moments(ia,ip0,ij+1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
      ENDDO
      Tper(ia,iky,ikx,iz) = Tperp_
      ENDDO
      ENDDO
      ENDDO
      ENDIF
      ENDDO
   END SUBROUTINE compute_Tperp

! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
   SUBROUTINE compute_Tpar
      USE time_integration, ONLY : updatetlevel
      IMPLICIT NONE
      REAL(xp)    :: j_xp
      COMPLEX(xp) :: Tpar_
      INTEGER     :: ia, iz, iky, ikx, ij
      DO ia=1,local_na
      IF ( CONTAINSp0 .AND. CONTAINSp0 ) THEN
      ! Loop to compute T = 1/3*(Tpar + 2Tperp)
      DO iz = 1,local_nz
      DO iky = 1,local_nky
      DO ikx = 1,local_nkx
      Tpar_ = 0._xp
      DO ij = 1, local_nj
         j_xp = REAL(ij-1,xp)
         Tpar_  = Tpar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
            (SQRT2 * moments(ia,ip2,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel) &
                     + moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
      ENDDO
      Tpar(ia,iky,ikx,iz) = Tpar_
      ENDDO
      ENDDO
      ENDDO
      ENDIF
      ENDDO
   END SUBROUTINE compute_Tpar

! Compute the 2D particle fluid moments for electron and ions (sum over Laguerre)
   SUBROUTINE compute_fluid_moments
      IMPLICIT NONE
      CALL compute_density
      CALL compute_upar
      CALL compute_uperp
      CALL compute_Tpar
      CALL compute_Tperp
      ! Temperature
      temp = (Tpar - 2._xp * Tper)/3._xp - dens
   END SUBROUTINE compute_fluid_moments

END MODULE processing