bin_variables Subroutine

  SUBROUTINE bin_variables(params,spp)
    IMPLICIT NONE
    TYPE(KORC_PARAMS), INTENT(IN) :: params
    TYPE(SPECIES), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: spp
    CHARACTER(MAX_STRING_LENGTH) :: var_name
    REAL(rp), DIMENSION(3) :: X
    REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: eta
    REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: g
    REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: N
    REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: array3D
    REAL(rp), DIMENSION(:,:), ALLOCATABLE :: array2D
    REAL(rp), DIMENSION(:), ALLOCATABLE :: array1D
    REAL(rp) :: R,Z,phi
    REAL(rp) :: Dtor
    REAL(rp) :: q,m
    INTEGER :: ii,jj,kk,ss,pp
    REAL(rp), DIMENSION(:), ALLOCATABLE :: send_buffer, receive_buffer
    INTEGER :: numel, mpierr
    REAL(rp) :: units

    ALLOCATE(eta(binning_params%num_bins(1),binning_params%num_bins(2),binning_params%ntor_sections,params%num_species))
    ALLOCATE(g(binning_params%num_bins(1),binning_params%num_bins(2),binning_params%ntor_sections,params%num_species))
    ALLOCATE(N(binning_params%num_bins(1),binning_params%num_bins(2),binning_params%ntor_sections,params%num_species))

    eta = 0.0_rp
    g = 0.0_rp
    N = 0.0_rp

    Dtor = 2.0_rp*C_PI/REAL(binning_params%ntor_sections,rp)

    do ss=1_idef,params%num_species
       q = ABS(spp(ss)%q)*params%cpp%charge
       m = spp(ss)%m*params%cpp%mass
       !$OMP PARALLEL DO FIRSTPRIVATE(q,m,Dtor) PRIVATE(X,ii,jj,kk,pp)&
       !$OMP& SHARED(params,spp,ss,eta,g,N)
       do pp=1_idef,spp(ss)%ppp
          if ( spp(ss)%vars%flag(pp) .EQ. 1_idef ) then
             X = spp(ss)%vars%X(pp,:)*params%cpp%length

             R = SQRT(SUM(X(1:2)**2))
             Z = X(3)

             ii = FLOOR((R - binning_params%rmin)/binning_params%dr) + 1_idef
             jj = FLOOR((Z + ABS(binning_params%zmin))/binning_params%dz) + 1_idef

             phi = ATAN2(X(2),X(1))
             if (phi.LT.0.0_rp) phi = 2.0_rp*C_PI + phi
             kk = floor(phi/Dtor) + 1_idef

             eta(ii,jj,kk,ss) = eta(ii,jj,kk,ss) + spp(ss)%vars%eta(pp)
             g(ii,jj,kk,ss) = g(ii,jj,kk,ss) + spp(ss)%vars%g(pp)
             N(ii,jj,kk,ss) = N(ii,jj,kk,ss) + 1.0_rp
          end if ! if confined
       end do ! particles
       !$OMP END PARALLEL DO
    end do ! species

    if (params%mpi_params%rank.EQ.0_idef) then
       if (.NOT.binning_params%toroidal_sections) then
          ALLOCATE(array3D(binning_params%num_bins(1),binning_params%num_bins(2),params%num_species))
       end if
    end if

    if (params%mpi_params%nmpi.GT.1_idef) then
       numel = binning_params%num_bins(1)*binning_params%num_bins(2)*binning_params%ntor_sections*params%num_species

       ALLOCATE(send_buffer(numel))
       ALLOCATE(receive_buffer(numel))

       send_buffer = RESHAPE(eta,(/numel/))
       receive_buffer = 0.0_rp
       CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
       if (params%mpi_params%rank.EQ.0_idef) then
          eta = RESHAPE(receive_buffer,(/binning_params%num_bins(1),binning_params%num_bins(2),&
               binning_params%ntor_sections,params%num_species/))

          var_name = 'eta'
          if (binning_params%toroidal_sections) then
             call save_snapshot_var(params,eta,var_name)
          else
             array3D = SUM(eta,3)
             call save_snapshot_var(params,array3D,var_name)
          end if
       end if

       DEALLOCATE(send_buffer)
       DEALLOCATE(receive_buffer)

       ALLOCATE(send_buffer(numel))
       ALLOCATE(receive_buffer(numel))

       send_buffer = RESHAPE(g,(/numel/))
       receive_buffer = 0.0_rp
       CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
       if (params%mpi_params%rank.EQ.0_idef) then
          g = RESHAPE(receive_buffer,(/binning_params%num_bins(1),binning_params%num_bins(2),&
               binning_params%ntor_sections,params%num_species/))

          var_name = 'g'
          if (binning_params%toroidal_sections) then
             call save_snapshot_var(params,g,var_name)
          else
             array3D = SUM(g,3)
             call save_snapshot_var(params,array3D,var_name)
          end if
       end if

       DEALLOCATE(send_buffer)
       DEALLOCATE(receive_buffer)

       ALLOCATE(send_buffer(numel))
       ALLOCATE(receive_buffer(numel))

       send_buffer = RESHAPE(N,(/numel/))
       receive_buffer = 0.0_rp
       CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
       if (params%mpi_params%rank.EQ.0_idef) then
          N = RESHAPE(receive_buffer,(/binning_params%num_bins(1),binning_params%num_bins(2),&
               binning_params%ntor_sections,params%num_species/))

          var_name = 'N'
          if (binning_params%toroidal_sections) then
             call save_snapshot_var(params,N,var_name)
          else
             array3D = SUM(N,3)
             call save_snapshot_var(params,array3D,var_name)
          end if
       end if

       DEALLOCATE(send_buffer)
       DEALLOCATE(receive_buffer)
    else
       var_name = 'eta'
       if (binning_params%toroidal_sections) then
          call save_snapshot_var(params,eta,var_name)
       else
          array3D = SUM(eta,3)
          call save_snapshot_var(params,array3D,var_name)
       end if

       var_name = 'g'
       if (binning_params%toroidal_sections) then
          call save_snapshot_var(params,g,var_name)
       else
          array3D = SUM(g,3)
          call save_snapshot_var(params,array3D,var_name)
       end if

       var_name = 'N'
       if (binning_params%toroidal_sections) then
          call save_snapshot_var(params,N,var_name)
       else
          array3D = SUM(N,3)
          call save_snapshot_var(params,array3D,var_name)
       end if
    end if

    DEALLOCATE(eta)
    DEALLOCATE(g)
    DEALLOCATE(N)

    if (ALLOCATED(array3D)) DEALLOCATE(array3D)
    if (ALLOCATED(array2D)) DEALLOCATE(array2D)
    if (ALLOCATED(array1D)) DEALLOCATE(array1D)
  END SUBROUTINE bin_variables