  Jet-plasmas emanating from the vicinity of relativistic objects and in
gamma-ray bursts have been observed to propagate with Lorentz factors laying in
the range between one and several hundreds. On the other hand, the numerical
studies of such flows have been focussed so far mainly on the lowest possible
range of Lorentz factors $\Gamma,$ specifically, on the regime $1\leq \Gamma
\leq 5.$ Therefore, relativistic flows with high $\Gamma-$factors have poorly
studied, as most numerical methods are found to encounter severe numerical
difficulties or even become numerically unstable for $\Gamma \gg 1. In this
paper we present an implicit numerical advection scheme for modeling the
propagation of relativistic plasmas with shocks, discuss its consistency with
respect to both the internal and total energy formulation in general
relativity. Using the total energy formulation, the scheme is found to be
viable for modeling moving shocks with moderate Lorentz factors, though with
relatively small Courant numbers. In the limit of high Lorentz factors, the
internal energy formulation in combination with a fine-tuned artificial
viscosity is much more robust and efficient. We confirm our conclusions by
performing test calculations and compare the results with analytical solutions
of the relativistic shock tube problem. The aim of the present modification is
to enhance the robustness of the general relativistic implicit radiative MHD
solver: GR-I-RMHD
(http://www1.iwr.uni-heidelberg.de/groups/compastro/home/gr-i-mhd-solver) and
extend its range of applications into the high $\Gamma-$regime.}
