diff --git a/src/gravity.F90 b/src/gravity.F90
index 200fb64..82824ec 100644
--- a/src/gravity.F90
+++ b/src/gravity.F90
@@ -216,11 +216,11 @@ module gravity
 !
 !     t, dt   - time and the time increment;
 !     x, y, z - rectangular coordinates;
-!     gacc    - vector of the gravitational acceleration;
+!     acc     - vector of the gravitational acceleration;
 !
 !===============================================================================
 !
-  subroutine gacc_none(t, dt, x, y, z, gacc)
+  subroutine gacc_none(t, dt, x, y, z, acc)
 
 ! local variables are not implicit by default
 !
@@ -230,7 +230,7 @@ module gravity
 !
     real(kind=8)              , intent(in)  :: t, dt
     real(kind=8)              , intent(in)  :: x, y, z
-    real(kind=8), dimension(3), intent(out) :: gacc
+    real(kind=8), dimension(3), intent(out) :: acc
 !
 !-------------------------------------------------------------------------------
 !
@@ -240,6 +240,10 @@ module gravity
     call start_timer(imc)
 #endif /* PROFILE */
 
+! reset gravitational acceleration
+!
+    acc(:) = 0.0d+00
+
 #ifdef PROFILE
 ! stop accounting time for the gravitational acceleration calculation
 !
diff --git a/src/interpolations.F90 b/src/interpolations.F90
index c04521d..80fa8c1 100644
--- a/src/interpolations.F90
+++ b/src/interpolations.F90
@@ -86,6 +86,7 @@ module interpolations
 
 ! flags for reconstruction corrections
 !
+  logical     , save :: mlp        = .false.
   logical     , save :: positivity = .false.
   logical     , save :: clip       = .false.
 
@@ -137,6 +138,7 @@ module interpolations
     character(len=255) :: tlimiter        = "mm"
     character(len=255) :: plimiter        = "mm"
     character(len=255) :: climiter        = "mm"
+    character(len=255) :: mlp_limiting    = "off"
     character(len=255) :: positivity_fix  = "off"
     character(len=255) :: clip_extrema    = "off"
     character(len=255) :: name_rec        = ""
@@ -169,6 +171,7 @@ module interpolations
     call get_parameter_string ("clip_extrema"        , clip_extrema   )
     call get_parameter_string ("extrema_limiter"     , climiter       )
     call get_parameter_string ("prolongation_limiter", plimiter       )
+    call get_parameter_string ("mlp_limiting"        , mlp_limiting   )
     call get_parameter_integer("nghosts"             , ng             )
     call get_parameter_integer("ngp"                 , ngp            )
     call get_parameter_real   ("sgp"                 , sgp            )
@@ -390,6 +393,12 @@ module interpolations
 
 ! check additional reconstruction limiting
 !
+    select case(trim(mlp_limiting))
+    case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
+      mlp = .true.
+    case default
+      mlp = .false.
+    end select
     select case(trim(positivity_fix))
     case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
       positivity = .true.
@@ -410,6 +419,7 @@ module interpolations
       write (*,"(4x,a14, 9x,'=',1x,a)") "reconstruction"      , trim(name_rec)
       write (*,"(4x,a11,12x,'=',1x,a)") "TVD limiter"         , trim(name_tlim)
       write (*,"(4x,a20, 3x,'=',1x,a)") "prolongation limiter", trim(name_plim)
+      write (*,"(4x,a12,11x,'=',1x,a)") "MLP limiting"        , trim(mlp_limiting)
       write (*,"(4x,a14, 9x,'=',1x,a)") "fix positivity"      , trim(positivity_fix)
       write (*,"(4x,a12,11x,'=',1x,a)") "clip extrema"        , trim(clip_extrema)
       if (clip) then
@@ -673,9 +683,12 @@ module interpolations
 
 ! local variables
 !
-    integer                        :: i, im1, ip1
-    integer                        :: j, jm1, jp1
-    integer                        :: k, km1, kp1
+    integer :: i, im1, ip1
+    integer :: j, jm1, jp1
+    integer :: k, km1, kp1
+
+! local vectors
+!
     real(kind=8), dimension(NDIMS) :: dql, dqr, dq
 !
 !-------------------------------------------------------------------------------
@@ -751,6 +764,10 @@ module interpolations
       end do ! j = jbl, jeu
     end do ! k = kbl, keu
 
+! apply Multi-dimensional Limiting Process
+!
+    if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:))
+
 !-------------------------------------------------------------------------------
 !
   end subroutine interfaces_tvd
@@ -858,6 +875,10 @@ module interpolations
 
     end if
 
+! apply Multi-dimensional Limiting Process
+!
+    if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:))
+
 !-------------------------------------------------------------------------------
 !
   end subroutine interfaces_dir
@@ -868,7 +889,7 @@ module interpolations
 ! -------------------------
 !
 !   Subroutine reconstructs both side interfaces of variable using
-!   multidimentional Gaussian Process method.
+!   multidimensional Gaussian Process method.
 !
 !   Arguments:
 !
@@ -1028,12 +1049,194 @@ module interpolations
       end do ! j = jbl, jeu
     end do ! k = kbl, keu
 
+! apply Multi-dimensional Limiting Process
+!
+    if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:))
+
 !-------------------------------------------------------------------------------
 !
   end subroutine interfaces_mgp
 !
 !===============================================================================
 !
+! subroutine MLP_LIMITING:
+! -----------------------
+!
+!   Subroutine applies the multi-dimensional limiting process to
+!   the reconstructed states in order to control oscillations due to
+!   the Gibbs phenomena near discontinuities.
+!
+!   Arguments:
+!
+!     q        - the variable array;
+!     qi       - the array of reconstructed interfaces (2 in each direction);
+!
+!   References:
+!
+!     [1] Gerlinger, P.,
+!         "Multi-dimensional limiting for high-order schemes including
+!          turbulence and combustion",
+!         Journal of Computational Physics,
+!         2012, vol. 231, pp. 2199-2228,
+!         http://dx.doi.org/10.1016/j.jcp.2011.10.024
+!
+!===============================================================================
+!
+  subroutine mlp_limiting(q, qi)
+
+! include external procedures
+!
+    use coordinates    , only : im , jm , km
+    use coordinates    , only : ibl, jbl, kbl, ieu, jeu, keu
+
+! local variables are not implicit by default
+!
+    implicit none
+
+! subroutine arguments
+!
+    real(kind=8), dimension(im,jm,km)        , intent(in)    :: q
+    real(kind=8), dimension(im,jm,km,2,NDIMS), intent(inout) :: qi
+
+! local variables
+!
+    integer      :: i, im1, ip1
+    integer      :: j, jm1, jp1
+    integer      :: k, km1, kp1
+    integer      :: m
+#if NDIMS == 3
+    integer      :: n, np1, np2
+#endif /* NDIMS == 3 */
+    real(kind=8) :: qmn, qmx, dqc
+    real(kind=8) :: fc, fl
+
+! local vectors
+!
+    real(kind=8), dimension(NDIMS) :: dql, dqr, dq
+    real(kind=8), dimension(NDIMS) :: dqm, ap
+#if NDIMS == 3
+    real(kind=8), dimension(NDIMS) :: hh, uu
+#endif /* NDIMS == 3 */
+!
+!-------------------------------------------------------------------------------
+!
+    do k = kbl, keu
+#if NDIMS == 3
+      km1 = k - 1
+      kp1 = k + 1
+#endif /* NDIMS == 3 */
+      do j = jbl, jeu
+        jm1 = j - 1
+        jp1 = j + 1
+        do i = ibl, ieu
+          im1 = i - 1
+          ip1 = i + 1
+
+! calculate the minmod TVD derivatives
+!
+          dql(1) = q(i  ,j,k) - q(im1,j,k)
+          dqr(1) = q(ip1,j,k) - q(i  ,j,k)
+          dq (1) = limiter_minmod(0.5d+00, dql(1), dqr(1))
+
+          dql(2) = q(i,j  ,k) - q(i,jm1,k)
+          dqr(2) = q(i,jp1,k) - q(i,j  ,k)
+          dq (2) = limiter_minmod(0.5d+00, dql(2), dqr(2))
+
+#if NDIMS == 3
+          dql(3) = q(i,j,k  ) - q(i,j,km1)
+          dqr(3) = q(i,j,kp1) - q(i,j,k  )
+          dq (3) = limiter_minmod(0.5d+00, dql(3), dqr(3))
+#endif /* NDIMS == 3 */
+
+! calculate dqc
+!
+#if NDIMS == 3
+          qmn = minval(q(im1:ip1,jm1:jp1,km1:kp1))
+          qmx = maxval(q(im1:ip1,jm1:jp1,km1:kp1))
+#else /* NDIMS == 3 */
+          qmn = minval(q(im1:ip1,jm1:jp1,k))
+          qmx = maxval(q(im1:ip1,jm1:jp1,k))
+#endif /* NDIMS == 3 */
+          dqc = min(qmx - q(i,j,k), q(i,j,k) - qmn)
+
+! if needed, apply the multi-dimensional limiting process
+!
+          if (sum(abs(dq(1:NDIMS))) > dqc) then
+
+            dqm(1) = abs(q(ip1,j,k) - q(im1,j,k))
+            dqm(2) = abs(q(i,jp1,k) - q(i,jm1,k))
+#if NDIMS == 3
+            dqm(3) = abs(q(i,j,kp1) - q(i,j,km1))
+#endif /* NDIMS == 3 */
+
+            fc = dqc / max(1.0d-16, sum(abs(dqm(1:NDIMS))))
+            do m = 1, NDIMS
+              ap(m) = fc * abs(dqm(m))
+            end do
+
+#if NDIMS == 3
+            do n = 1, NDIMS
+              hh(n)   = max(ap(n) - abs(dq(n)), 0.0d+00)
+              np1     = mod(n    , NDIMS) + 1
+              np2     = mod(n + 1, NDIMS) + 1
+              uu(n)   = ap(n)   -           hh(n)
+              uu(np1) = ap(np1) + 0.5d+00 * hh(n)
+              uu(np2) = ap(np2) + 0.5d+00 * hh(n)
+              fc      = hh(n) / (hh(n) + 1.0d-16)
+              fl      = fc * (max(ap(np1) - abs(dq(np1)), 0.0d+00)             &
+                            - max(ap(np2) - abs(dq(np2)), 0.0d+00))
+              ap(n  ) = uu(n)
+              ap(np1) = uu(np1) - fl
+              ap(np2) = uu(np2) + fl
+            end do
+#else /* NDIMS == 3 */
+            fl     = max(ap(1) - abs(dq(1)), 0.0d+00)                          &
+                   - max(ap(2) - abs(dq(2)), 0.0d+00)
+            ap(1)  = ap(1) - fl
+            ap(2)  = ap(2) + fl
+#endif /* NDIMS == 3 */
+
+            do m = 1, NDIMS
+              dq(m) = sign(ap(m), dq(m))
+            end do
+
+! update the interpolated states
+!
+            dql(1) = qi(i  ,j,k,1,1) - q(i,j,k)
+            dqr(1) = qi(im1,j,k,2,1) - q(i,j,k)
+            if (max(abs(dql(1)), abs(dqr(1))) > abs(dq(1))) then
+              qi(i  ,j,k,1,1) = q(i,j,k) + dq(1)
+              qi(im1,j,k,2,1) = q(i,j,k) - dq(1)
+            end if
+
+            dql(2) = qi(i,j  ,k,1,2) - q(i,j,k)
+            dqr(2) = qi(i,jm1,k,2,2) - q(i,j,k)
+            if (max(abs(dql(2)), abs(dqr(2))) > abs(dq(2))) then
+              qi(i,j  ,k,1,2) = q(i,j,k) + dq(2)
+              qi(i,jm1,k,2,2) = q(i,j,k) - dq(2)
+            end if
+
+#if NDIMS == 3
+            dql(3) = qi(i,j,k  ,1,3) - q(i,j,k))
+            dqr(3) = qi(i,j,km1,2,3) - q(i,j,k))
+            if (max(abs(dql(3)), abs(dqr(3))) > abs(dq(3))) then
+              qi(i,j,k  ,1,3) = q(i,j,k) + dq(3)
+              qi(i,j,km1,2,3) = q(i,j,k) - dq(3)
+            end if
+#endif /* NDIMS == 3 */
+
+          end if
+
+        end do ! i = ibl, ieu
+      end do ! j = jbl, jeu
+    end do ! k = kbl, keu
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine mlp_limiting
+!
+!===============================================================================
+!
 ! subroutine RECONSTRUCT:
 ! ----------------------
 !
diff --git a/src/io.F90 b/src/io.F90
index 1381bb2..44e9df0 100644
--- a/src/io.F90
+++ b/src/io.F90
@@ -33,9 +33,7 @@ module io
 ! import external subroutines
 !
   use blocks, only : pointer_meta
-#ifdef PROFILE
   use timers, only : set_timer, start_timer, stop_timer
-#endif /* PROFILE */
 
 ! module variables are not implicit by default
 !
@@ -89,9 +87,10 @@ module io
 #endif /* HDF5 */
   end interface
 
-#ifdef PROFILE
 ! timer indices
 !
+  integer            , save :: iio
+#ifdef PROFILE
   integer            , save :: ioi, iow, ios
 #endif /* PROFILE */
 
@@ -204,9 +203,10 @@ module io
 !
 !-------------------------------------------------------------------------------
 !
-#ifdef PROFILE
 ! set timer descriptions
 !
+    call set_timer('SNAPSHOTS I/O'        , iio)
+#ifdef PROFILE
     call set_timer('io:: initialization'  , ioi)
     call set_timer('io:: snapshot writing', iow)
     call set_timer('io:: snapshot reading', ios)
@@ -360,12 +360,24 @@ module io
     call start_timer(ios)
 #endif /* PROFILE */
 
+! reset the return flag
+!
+    iret = 0
+
+! start accounting time for I/O
+!
+    call start_timer(iio)
+
 #ifdef HDF5
 ! read HDF5 restart file and rebuild the meta and data block structures
 !
     call read_restart_snapshot_h5(iret)
 #endif /* HDF5 */
 
+! stop accounting time for I/O
+!
+    call stop_timer(iio)
+
 ! calculate the shift of the snapshot counter, and the next snapshot time
 !
     ishift = int(time / hsnap) - isnap + 1
@@ -411,6 +423,10 @@ module io
 !
 !-------------------------------------------------------------------------------
 !
+! reset the return flag
+!
+    iret = 0
+
 ! check if conditions for storing the restart snapshot have been met
 !
     if (hrest < 5.0d-02 .or. thrs < irest * hrest) return
@@ -421,12 +437,20 @@ module io
     call start_timer(iow)
 #endif /* PROFILE */
 
+! start accounting time for I/O
+!
+    call start_timer(iio)
+
 #ifdef HDF5
 ! store restart file
 !
     call write_restart_snapshot_h5(nrun, iret)
 #endif /* HDF5 */
 
+! stop accounting time for I/O
+!
+    call stop_timer(iio)
+
 ! increase the restart snapshot counter
 !
     irest = irest + 1
@@ -476,6 +500,10 @@ module io
     call start_timer(iow)
 #endif /* PROFILE */
 
+! start accounting time for I/O
+!
+    call start_timer(iio)
+
 #ifdef HDF5
 ! store variable snapshot file
 !
@@ -486,6 +514,10 @@ module io
     end if
 #endif /* HDF5 */
 
+! stop accounting time for I/O
+!
+    call stop_timer(iio)
+
 ! increase the snapshot counter and calculate the next snapshot time
 !
     isnap = isnap + 1
diff --git a/src/schemes.F90 b/src/schemes.F90
index d07c90d..95f4611 100644
--- a/src/schemes.F90
+++ b/src/schemes.F90
@@ -2091,10 +2091,9 @@ module schemes
         dnl = wl(idn) / slmm
         dnr = wr(idn) / srmm
 
-! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to
-! the HLL one
+! apply full HLLD solver only if the Alfvén wave is strong enough
 !
-        if (b2 > 0.0d+00) then ! Bₓ² > 0
+        if (b2 > 1.0d-08 * max(dnl * sl**2, dnr * sr**2)) then ! Bₓ² > ε
 
 ! left and right Alfvén speeds
 !
@@ -2123,7 +2122,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
               ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
-              ui(ibp) = ql(ibp,i)
+              ui(ibp) = ul(ibp,i)
 
 ! the outer left intermediate flux
 !
@@ -2140,7 +2139,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
               ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
-              ui(ibp) = qr(ibp,i)
+              ui(ibp) = ur(ibp,i)
 
 ! the outer right intermediate flux
 !
@@ -2157,7 +2156,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
               ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
-              ui(ibp) = ql(ibp,i)
+              ui(ibp) = ul(ibp,i)
 
 ! vector of the left-going Alfvén wave
 !
@@ -2172,7 +2171,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
               ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
-              ui(ibp) = qr(ibp,i)
+              ui(ibp) = ur(ibp,i)
 
 ! vector of the right-going Alfvén wave
 !
@@ -2199,7 +2198,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
               ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
-              ui(ibp) = ql(ibp,i)
+              ui(ibp) = ul(ibp,i)
 
 ! choose the correct state depending on the speed signs
 !
@@ -2232,7 +2231,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
               ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
-              ui(ibp) = qr(ibp,i)
+              ui(ibp) = ur(ibp,i)
 
 ! choose the correct state depending on the speed signs
 !
@@ -2264,12 +2263,16 @@ module schemes
 
           end if ! one degeneracy
 
-        else ! Bₓ² = 0
+        else if (b2 > 0.0d+00) then ! Bₓ² > 0
 
-! in the vase of vanishing Bₓ there is no Alfvén wave, density is constant, and
-! still we can have a discontinuity in the perpendicular components
+! the Alfvén wave is too weak, so ignore it in order to not introduce any
+! numerical instabilities; since Bₓ is not zero, the perpendicular components
+! of velocity and magnetic field are continuous; we have only one intermediate
+! state, therefore simply apply the HLL solver
 !
-          dn   = dn / srml
+          f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
+
+        else ! Bₓ² = 0
 
 ! take the right flux depending on the sign of the advection speed
 !
@@ -2277,14 +2280,14 @@ module schemes
 
 ! conservative variables for the left intermediate state
 !
-            ui(idn) = dn
-            ui(imx) = dn *  sm
+            ui(idn) = dnl
+            ui(imx) = dnl * sm
             ui(imy) = wl(imy) / slmm
             ui(imz) = wl(imz) / slmm
             ui(ibx) = 0.0d+00
             ui(iby) = wl(iby) / slmm
             ui(ibz) = wl(ibz) / slmm
-            ui(ibp) = ql(ibp,i)
+            ui(ibp) = ul(ibp,i)
 
 ! the left intermediate flux
 !
@@ -2294,14 +2297,14 @@ module schemes
 
 ! conservative variables for the right intermediate state
 !
-            ui(idn) = dn
-            ui(imx) = dn *  sm
+            ui(idn) = dnr
+            ui(imx) = dnr * sm
             ui(imy) = wr(imy) / srmm
             ui(imz) = wr(imz) / srmm
             ui(ibx) = 0.0d+00
             ui(iby) = wr(iby) / srmm
             ui(ibz) = wr(ibz) / srmm
-            ui(ibp) = qr(ibp,i)
+            ui(ibp) = ur(ibp,i)
 
 ! the right intermediate flux
 !
@@ -2310,16 +2313,9 @@ module schemes
           else ! sm = 0
 
 ! the intermediate flux; since the advection speed is zero, perpendicular
-! components do not change
+! components do not change, so revert it to HLL
 !
-            f(idn,i) = (sl * wr(idn) - sr * wl(idn)) / srml
-            f(imx,i) = (sl * wr(imx) - sr * wl(imx)) / srml
-            f(imy,i) = 0.0d+00
-            f(imz,i) = 0.0d+00
-            f(ibx,i) = (sl * wr(ibx) - sr * wl(ibx)) / srml
-            f(iby,i) = 0.0d+00
-            f(ibz,i) = 0.0d+00
-            f(ibp,i) = (sl * wr(ibp) - sr * wl(ibp)) / srml
+            f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
 
           end if ! sm = 0
 
@@ -2384,7 +2380,7 @@ module schemes
 !
     integer                       :: i
     real(kind=8)                  :: sl, sr, sm, sml, smr, srml, slmm, srmm
-    real(kind=8)                  :: bx, bp, b2, dn, dnl, dnr, dvl, dvr, ca
+    real(kind=8)                  :: bx, bp, b2, dn, dnl, dnr, dvl, dvr, ca, ca2
 
 ! local arrays to store the states
 !
@@ -2475,14 +2471,17 @@ module schemes
 !
         dn = dn / srml
 
-! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to
-! the HLL one
+! get the square of the Alfvén speed
 !
-        if (b2 > 0.0d+00) then ! Bₓ² > 0
+        ca2 = b2 / dn
+
+! apply full HLLD solver only if the Alfvén wave is strong enough
+!
+        if (ca2 > 1.0d-08 * max(sl**2, sr**2)) then ! Bₓ² > ε
 
 ! left and right Alfvén speeds
 !
-          ca   = sqrt(b2 / dn)
+          ca   = sqrt(ca2)
           sml  = sm - ca
           smr  = sm + ca
 
@@ -2508,7 +2507,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl
               ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl
-              ui(ibp) = ql(ibp,i)
+              ui(ibp) = ul(ibp,i)
 
 ! the outer left intermediate flux
 !
@@ -2525,7 +2524,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr
               ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr
-              ui(ibp) = qr(ibp,i)
+              ui(ibp) = ur(ibp,i)
 
 ! the outer right intermediate flux
 !
@@ -2542,7 +2541,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl
               ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl
-              ui(ibp) = ql(ibp,i)
+              ui(ibp) = ul(ibp,i)
 
 ! vector of the left-going Alfvén wave
 !
@@ -2557,7 +2556,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr
               ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr
-              ui(ibp) = qr(ibp,i)
+              ui(ibp) = ur(ibp,i)
 
 ! vector of the right-going Alfvén wave
 !
@@ -2584,7 +2583,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl
               ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl
-              ui(ibp) = ql(ibp,i)
+              ui(ibp) = ul(ibp,i)
 
 ! choose the correct state depending on the speed signs
 !
@@ -2617,7 +2616,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr
               ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr
-              ui(ibp) = qr(ibp,i)
+              ui(ibp) = ur(ibp,i)
 
 ! choose the correct state depending on the speed signs
 !
@@ -2649,6 +2648,15 @@ module schemes
 
           end if ! one degeneracy
 
+        else if (b2 > 0.0d+00) then ! Bₓ² > 0
+
+! the Alfvén wave is very weak, so ignore it in order to not introduce any
+! numerical instabilities; since Bₓ is not zero, the perpendicular components
+! of velocity and magnetic field are continuous; we have only one intermediate
+! state, therefore simply apply the HLL solver
+!
+          f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
+
         else ! Bₓ² = 0
 
 ! in the vase of vanishing Bₓ there is no Alfvén wave, density is constant, and
@@ -2667,7 +2675,7 @@ module schemes
             ui(ibx) = 0.0d+00
             ui(iby) = wl(iby) / slmm
             ui(ibz) = wl(ibz) / slmm
-            ui(ibp) = ql(ibp,i)
+            ui(ibp) = ul(ibp,i)
 
 ! the left intermediate flux
 !
@@ -2684,7 +2692,7 @@ module schemes
             ui(ibx) = 0.0d+00
             ui(iby) = wr(iby) / srmm
             ui(ibz) = wr(ibz) / srmm
-            ui(ibp) = qr(ibp,i)
+            ui(ibp) = ur(ibp,i)
 
 ! the right intermediate flux
 !
@@ -2692,17 +2700,9 @@ module schemes
 
           else ! sm = 0
 
-! the intermediate flux; since the advection speed is zero, perpendicular
-! components do not change
+! both states are equal so revert to the HLL flux
 !
-            f(idn,i) = (sl * wr(idn) - sr * wl(idn)) / srml
-            f(imx,i) = (sl * wr(imx) - sr * wl(imx)) / srml
-            f(imy,i) = 0.0d+00
-            f(imz,i) = 0.0d+00
-            f(ibx,i) = (sl * wr(ibx) - sr * wl(ibx)) / srml
-            f(iby,i) = 0.0d+00
-            f(ibz,i) = 0.0d+00
-            f(ibp,i) = (sl * wr(ibp) - sr * wl(ibp)) / srml
+            f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
 
           end if ! sm = 0
 
@@ -3668,35 +3668,35 @@ module schemes
         dn =  wr(idn) - wl(idn)
         sm = (wr(imx) - wl(imx)) / dn
 
+! speed differences
+!
+        slmm = sl - sm
+        srmm = sr - sm
+
+! left and right state densities
+!
+        dnl = wl(idn) / slmm
+        dnr = wr(idn) / srmm
+
 ! square of Bₓ, i.e. Bₓ²
 !
         bx = ql(ibx,i)
         b2 = ql(ibx,i) * qr(ibx,i)
 
-! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to
-! the HLLC solver
-!
-        if (b2 > 0.0d+00) then ! Bₓ² > 0
-
 ! the total pressure, constant across the contact discontinuity and Alfvén waves
 !
-          pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
+        pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
 
-! speed differences
+! if Alfvén wave is strong enough, apply the full HLLD solver, otherwise
+! revert to the HLLC one
 !
-          slmm = sl - sm
-          srmm = sr - sm
-
-! left and right state densities
-!
-          dnl = wl(idn) / slmm
-          dnr = wr(idn) / srmm
+        if (b2 > 1.0d-08 * max(dnl * sl**2, dnr * sr**2)) then ! Bₓ² > ε
 
 ! left and right Alfvén speeds
 !
-          cal = - sqrt(b2 / dnl)
-          car =   sqrt(b2 / dnr)
-          sml = sm + cal
+          cal = sqrt(b2 / dnl)
+          car = sqrt(b2 / dnr)
+          sml = sm - cal
           smr = sm + car
 
 ! calculate division factors (also used to determine degeneracies)
@@ -3729,7 +3729,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = by
               ui(ibz) = bz
-              ui(ibp) = ql(ibp,i)
+              ui(ibp) = ul(ibp,i)
               ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
 
 ! the outer left intermediate flux
@@ -3755,7 +3755,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = by
               ui(ibz) = bz
-              ui(ibp) = qr(ibp,i)
+              ui(ibp) = ur(ibp,i)
               ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
 
 ! the outer right intermediate flux
@@ -3781,7 +3781,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = by
               ui(ibz) = bz
-              ui(ibp) = ql(ibp,i)
+              ui(ibp) = ul(ibp,i)
               ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
 
 ! vector of the left-going Alfvén wave
@@ -3805,7 +3805,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = by
               ui(ibz) = bz
-              ui(ibp) = qr(ibp,i)
+              ui(ibp) = ur(ibp,i)
               ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
 
 ! vector of the right-going Alfvén wave
@@ -3814,10 +3814,10 @@ module schemes
 
 ! prepare constant primitive variables of the intermediate states
 !
-              dv      = car * dnr - cal * dnl
+              dv      = car * dnr + cal * dnl
               vy      = (wcr(imy) - wcl(imy)) / dv
               vz      = (wcr(imz) - wcl(imz)) / dv
-              dv      = car       - cal
+              dv      = car       + cal
               by      = (wcr(iby) - wcl(iby)) / dv
               bz      = (wcr(ibz) - wcl(ibz)) / dv
               vb      = sm * bx + vy * by + vz * bz
@@ -3836,8 +3836,8 @@ module schemes
                 ui(ibx) = bx
                 ui(iby) = by
                 ui(ibz) = bz
-                ui(ibp) = ql(ibp,i)
-                ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal
+                ui(ibp) = ul(ibp,i)
+                ui(ien) = (wcl(ien) + sm * pt - bx * vb) / (sml - sm)
 
 ! the inmost left intermediate flux
 !
@@ -3854,8 +3854,8 @@ module schemes
                 ui(ibx) = bx
                 ui(iby) = by
                 ui(ibz) = bz
-                ui(ibp) = qr(ibp,i)
-                ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
+                ui(ibp) = ur(ibp,i)
+                ui(ien) = (wcr(ien) + sm * pt - bx * vb) / (smr - sm)
 
 ! the inmost right intermediate flux
 !
@@ -3895,7 +3895,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = by
               ui(ibz) = bz
-              ui(ibp) = ql(ibp,i)
+              ui(ibp) = ul(ibp,i)
               ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
 
 ! choose the correct state depending on the speed signs
@@ -3914,7 +3914,7 @@ module schemes
 
 ! primitive variables for the inner left intermediate state
 !
-                dv      = srmm * dnr - cal * dnl
+                dv      = srmm * dnr + cal * dnl
                 vy      = (wr(imy) - wcl(imy)) / dv
                 vz      = (wr(imz) - wcl(imz)) / dv
                 dv      = sr - sml
@@ -3922,23 +3922,23 @@ module schemes
                 bz      = (wr(ibz) - wcl(ibz)) / dv
                 vb      = sm * bx + vy * by + vz * bz
 
+! conservative variables for the inner left intermediate state
+!
+                ui(idn) = dnl
+                ui(imx) = dnl * sm
+                ui(imy) = dnl * vy
+                ui(imz) = dnl * vz
+                ui(ibx) = bx
+                ui(iby) = by
+                ui(ibz) = bz
+                ui(ibp) = ul(ibp,i)
+                ui(ien) = (wcl(ien) + sm * pt - bx * vb) / (sml - sm)
+
 ! choose the correct state depending on the sign of contact discontinuity
 ! advection speed
 !
                 if (sm >= 0.0d+00) then ! sm ≥ 0
 
-! conservative variables for the inner left intermediate state
-!
-                  ui(idn) = dnl
-                  ui(imx) = dnl * sm
-                  ui(imy) = dnl * vy
-                  ui(imz) = dnl * vz
-                  ui(ibx) = bx
-                  ui(iby) = by
-                  ui(ibz) = bz
-                  ui(ibp) = ql(ibp,i)
-                  ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal
-
 ! the inner left intermediate flux
 !
                   f(:,i)  = sml * ui(:) - wcl(:)
@@ -3947,11 +3947,11 @@ module schemes
 
 ! vector of the left-going Alfvén wave
 !
-                  wcl(:)  = (sm - sml) * ui(:) + wcl(:)
+                  wcr(:)  = (sm - sml) * ui(:) + wcl(:)
 
 ! calculate the average flux over the right inner intermediate state
 !
-                  f(:,i)  = (sm * wr(:) - sr * wcl(:)) / srmm
+                  f(:,i)  = (sm * wr(:) - sr * wcr(:)) / srmm
 
                 end if ! sm < 0
 
@@ -3976,7 +3976,7 @@ module schemes
               ui(ibx) = bx
               ui(iby) = by
               ui(ibz) = bz
-              ui(ibp) = qr(ibp,i)
+              ui(ibp) = ur(ibp,i)
               ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
 
 ! choose the correct state depending on the speed signs
@@ -4003,23 +4003,23 @@ module schemes
                 bz      = (wl(ibz) - wcr(ibz)) / dv
                 vb      = sm * bx + vy * by + vz * bz
 
+! conservative variables for the inner left intermediate state
+!
+                ui(idn) = dnr
+                ui(imx) = dnr * sm
+                ui(imy) = dnr * vy
+                ui(imz) = dnr * vz
+                ui(ibx) = bx
+                ui(iby) = by
+                ui(ibz) = bz
+                ui(ibp) = ur(ibp,i)
+                ui(ien) = (wcr(ien) + sm * pt - bx * vb) / (smr - sm)
+
 ! choose the correct state depending on the sign of contact discontinuity
 ! advection speed
 !
                 if (sm <= 0.0d+00) then ! sm ≤ 0
 
-! conservative variables for the inner left intermediate state
-!
-                  ui(idn) = dnr
-                  ui(imx) = dnr * sm
-                  ui(imy) = dnr * vy
-                  ui(imz) = dnr * vz
-                  ui(ibx) = bx
-                  ui(iby) = by
-                  ui(ibz) = bz
-                  ui(ibp) = qr(ibp,i)
-                  ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
-
 ! the inner right intermediate flux
 !
                   f(:,i)  = smr * ui(:) - wcr(:)
@@ -4028,11 +4028,11 @@ module schemes
 
 ! vector of the right-going Alfvén wave
 !
-                  wcr(:)  = (sm - smr) * ui(:) + wcr(:)
+                  wcl(:)  = (sm - smr) * ui(:) + wcr(:)
 
 ! calculate the average flux over the left inner intermediate state
 !
-                  f(:,i)  = (sm * wl(:) - sl * wcr(:)) / slmm
+                  f(:,i)  = (sm * wl(:) - sl * wcl(:)) / slmm
 
                 end if ! sm > 0
 
@@ -4048,30 +4048,83 @@ module schemes
 
           end if ! one degeneracy
 
-        else ! Bₓ² = 0
+        else if (b2 > 0.0d+00) then ! Bₓ² > 0
 
-! the total pressure, constant across the contact discontinuity
+! constant intermediate state tangential components of velocity and
+! magnetic field
 !
-          pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn
+          vy = (wr(imy) - wl(imy)) / dn
+          vz = (wr(imz) - wl(imz)) / dn
+          by = (wr(iby) - wl(iby)) / srml
+          bz = (wr(ibz) - wl(ibz)) / srml
+
+! scalar product of velocity and magnetic field for the intermediate states
+!
+          vb = sm * bx + vy * by + vz * bz
 
 ! separate intermediate states depending on the sign of the advection speed
 !
           if (sm > 0.0d+00) then ! sm > 0
 
-! the left speed difference
+! conservative variables for the left intermediate state
 !
-            slmm    =  sl - sm
+            ui(idn) =  dnl
+            ui(imx) =  dnl * sm
+            ui(imy) =  dnl * vy
+            ui(imz) =  dnl * vz
+            ui(ibx) =  bx
+            ui(iby) =  by
+            ui(ibz) =  bz
+            ui(ibp) =  ul(ibp,i)
+            ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
+
+! the left intermediate flux
+!
+            f(:,i)  = sl * ui(:) - wl(:)
+
+          else if (sm < 0.0d+00) then ! sm < 0
+
+! conservative variables for the right intermediate state
+!
+            ui(idn) =  dnr
+            ui(imx) =  dnr * sm
+            ui(imy) =  dnr * vy
+            ui(imz) =  dnr * vz
+            ui(ibx) =  bx
+            ui(iby) =  by
+            ui(ibz) =  bz
+            ui(ibp) =  ur(ibp,i)
+            ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
+
+! the right intermediate flux
+!
+            f(:,i)  = sr * ui(:) - wr(:)
+
+          else ! sm = 0
+
+! when Sₘ = 0 all variables are continuous, therefore the flux reduces
+! to the HLL one
+!
+            f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
+
+          end if ! sm = 0
+
+        else ! Bₓ² = 0
+
+! separate intermediate states depending on the sign of the advection speed
+!
+          if (sm > 0.0d+00) then ! sm > 0
 
 ! conservative variables for the left intermediate state
 !
-            ui(idn) =  wl(idn) / slmm
-            ui(imx) =  ui(idn) * sm
-            ui(imy) =  ui(idn) * ql(ivy,i)
-            ui(imz) =  ui(idn) * ql(ivz,i)
+            ui(idn) =  dnl
+            ui(imx) =  dnl * sm
+            ui(imy) =  wl(imy) / slmm
+            ui(imz) =  wl(imz) / slmm
             ui(ibx) =  0.0d+00
             ui(iby) =  wl(iby) / slmm
             ui(ibz) =  wl(ibz) / slmm
-            ui(ibp) =  ql(ibp,i)
+            ui(ibp) =  ul(ibp,i)
             ui(ien) = (wl(ien) + sm * pt) / slmm
 
 ! the left intermediate flux
@@ -4080,20 +4133,16 @@ module schemes
 
           else if (sm < 0.0d+00) then ! sm < 0
 
-! the right speed difference
-!
-            srmm    =  sr - sm
-
 ! conservative variables for the right intermediate state
 !
-            ui(idn) =  wr(idn) / srmm
-            ui(imx) =  ui(idn) * sm
-            ui(imy) =  ui(idn) * qr(ivy,i)
-            ui(imz) =  ui(idn) * qr(ivz,i)
+            ui(idn) =  dnr
+            ui(imx) =  dnr * sm
+            ui(imy) =  wr(imy) / srmm
+            ui(imz) =  wr(imz) / srmm
             ui(ibx) =  0.0d+00
             ui(iby) =  wr(iby) / srmm
             ui(ibz) =  wr(ibz) / srmm
-            ui(ibp) =  qr(ibp,i)
+            ui(ibp) =  ur(ibp,i)
             ui(ien) = (wr(ien) + sm * pt) / srmm
 
 ! the right intermediate flux
@@ -4102,18 +4151,10 @@ module schemes
 
           else ! sm = 0
 
-! intermediate flux is constant across the contact discontinuity and all except
-! the parallel momentum flux are zero
+! when Sₘ = 0 all variables are continuous, therefore the flux reduces
+! to the HLL one
 !
-            f(idn,i) =   0.0d+00
-            f(imx,i) = - 0.5d+00 * (wl(imx) + wr(imx))
-            f(imy,i) =   0.0d+00
-            f(imz,i) =   0.0d+00
-            f(ibx,i) = fl(ibx,i)
-            f(iby,i) =   0.0d+00
-            f(ibz,i) =   0.0d+00
-            f(ibp,i) = fl(ibp,i)
-            f(ien,i) =   0.0d+00
+            f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
 
           end if ! sm = 0
 
diff --git a/src/sources.F90 b/src/sources.F90
index 88d198f..1c13354 100644
--- a/src/sources.F90
+++ b/src/sources.F90
@@ -582,11 +582,17 @@ module sources
           call curl(dh(:), pdata%q(ibx:ibz,1:im,1:jm,1:km)                     &
                                             , tmp(inz,inx:inz,1:im,1:jm,1:km))
 
+! calculate J²
+!
+          db(1:im,1:jm,1:km) = tmp(inz,inx,1:im,1:jm,1:km)**2                  &
+                             + tmp(inz,iny,1:im,1:jm,1:km)**2                  &
+                             + tmp(inz,inz,1:im,1:jm,1:km)**2
+
 ! add the second resistive source term to the energy equation, i.e.
 !     d/dt E + ∇.F = η J²
 !
           du(ien,1:im,1:jm,1:km) = du(ien,1:im,1:jm,1:km)                      &
-                    + resistivity * sum(tmp(inz,inx:inz,1:im,1:jm,1:km)**2, 2)
+                                 + resistivity * db(1:im,1:jm,1:km)
 
         end if ! energy equation present
 
diff --git a/src/user_problem.F90 b/src/user_problem.F90
index 8a2332e..dab5585 100644
--- a/src/user_problem.F90
+++ b/src/user_problem.F90
@@ -831,6 +831,10 @@ module user_problem
     call start_timer(img)
 #endif /* PROFILE */
 
+! reset gravitational acceleration
+!
+    acc(:) = 0.0d+00
+
 #ifdef PROFILE
 ! stop accounting time for the gravitational acceleration calculation
 !