Numerical codes for electromagnetic wave propagation in fusion plasmas are mainly based on frequency-domain asymptotic methods, which provide a fast solution and are thus valuable for experiment design and control applications. However, in several cases of practical interest (like O-X-B mode conversion, mm-diagnostics) these tools run close to their limits of validity and should be compared to analytic or numerical full-wave solutions. The Fortran code RFFW (successor of FWTOR) solves Maxwell’s equations for the evolution of electromagnetic waves in plasmas. The propagation is computed with the FDTD method, whereas the medium response is described in terms of the generated plasma current. The plasma equilibrium can be analytic or experimental, including non-axisymmetric perturbations and density fluctuations, and the wave electric field spectrum can be arbitrary. In fusion-relevant applications, the code may conduct investigations of RF propagation and absorption in plasma geometry relevant to H&CD, reflectometry and MHD control. The results may provide a more detailed picture of the physics involved, and also benchmark the validity of the currently established tools in this field, which are based on simpler models. The code has been parallelized under EUROfusion HLST support, which has allowed the exploitation of the much larger CPU power and memory of supercomputers. The first step was to obtain a hybrid parallel (MPI + OpenMP) version while retaining the original code structure, with an overall speedup factor of 2500 as compared to the serial code version. Afterwards, the work was focused on MPI communication improvements like e.g., piggy-backing from multiple array exchanges in one and finer grain exchanges, further optimization of the memory usage from loop structures, and introduction of parallel I/O and restart functionality. All these improvements, apart from the simulation of cases where a small grid scale is required (e.g., EC absorption in high-density plasma), allow to upgrade the physics model of the code with the inclusion of processes that require themselves a very large computational burden, like the convolutional computation of the plasma current response in time domain based on the Fourier transform of the frequency-domain dielectric tensor. The latter may open the way for the exploration of EC physics problems for which asymptotic methods cannot be applied. In this work, first we present the main aspects of the code physics and technical setup, and also refer to details of the parallelization scheme. Then, we show results that exhibit the strong scaling of the code, as well as examine cases relevant to ECRH/ECCD application in medium-sized tokamaks (AUG, TCV).