Performance of tokamak fusion plasmas is heavily linked to the radial heat and particle transport, which is known to be mainly produced by turbulence driven by micro-instabilities. Understanding such processes is thus of key importance for the design and operation of future fusion reactors such as ITER. One of the most commonly used tool to study turbulent transport in the very hot core of tokamaks is gyrokinetic theory, which aims at reducing the physical complexity by averaging out the fast gyromotion, while keeping its averaged effect. The \orb{} program is a numerical code solving the gyrokinetic Maxwell-Boltzmann system of equations to study various phenomena such as turbulent transport. It is a global particle-in-cell code using a finite element representation with B-spline basis function for the fields and various noise control schemes. The main objective of this work is to improve the performance and physical model of \orb{} in order to be able to carry out simulations relevant for the Tokamak à Configuration Variable (TCV) at the Swiss Plasma Center. This is done in three major steps. First, the code undergoes a complete refactoring and optimization process with the aim of increasing its maintability and performance. Using a simple test bed being a simplification of the PIC method while retaining the key elements, various algorithm and parallelization schemes are developed. Furthermore, it is shown that one of the most critical kernels of the PIC algorithm, the charge deposition, can be accelerated by a factor 2.4 in \orb{}. Second, the hybrid electron model, which is used to simulate certain electron modes at a lower numerical cost than the fully kinetic electron model, needs to be corrected. Indeed, in this model, trapped electrons are kinetic while passing electrons have an adiabatic response. This is problematic because it does not satisfy the ambipolarity condition and produces spurious sources of particles and momentum. A corrected model in which the flux-surface averaged contribution of the passing electrons is treated as kinetic is developed and implemented in \orb{}. As a first step, an approximation of the flux-surface average consisting of a Fourier filter retaining only the $n=m=0$ mode is tested. Linearly, both the proper flux-surface average and its $n=m=0$ approximation lead to an accurate description of the geodesic acoustic mode dynamics. Nonlinearly, the $n=m=0$ model exhibit interesting physical features such as sharp radial structures at the mode rational surface $q=1$ and a transport regime transition. Finally, a moment-conserving heat source is implemented in the \orb{} code to allow for flux-driven simulations. The heat source conservation properties are verified and it is shown that conservation of density, parallel momentum, and zonal flows is satisfied up to machine precision. A ``mode-switching'' approach consisting in using an effective heat source computed from a gradient-driven simulation as an input for the flux-driven run is also tested and consistent results are obtained between gradient-driven, flux-driven, and ``mode-switched'' simulations. Finally, the heat source is characterized, along with variations, using a 2D phase-space diagnostics. It is shown that in some cases, the use of this fixed heat source can lead to a sampling problem causing an artificial source/sink of moments.