diff --git a/python/amunpy/src/amunpy/amun.py b/python/amunpy/src/amunpy/amun.py
index c9d00f1..12f1447 100644
--- a/python/amunpy/src/amunpy/amun.py
+++ b/python/amunpy/src/amunpy/amun.py
@@ -215,11 +215,12 @@ class Amun:
         if all(v in self.variables for v in ['dens', 'pres']):
             self.variables['temperature']                   = 'temp'
         if self.attributes['eqsys'] in [ 'hd', 'mhd' ] \
-                        and all(v in self.variables for v in ['dens','velo']):
+                        and all(v in self.variables for v in ['dens','velx','vely','velz']):
             self.variables['kinetic energy']                = 'ekin'
         if self.attributes['eqsys'] in [ 'mhd', 'srmhd' ] \
                         and all(v in self.variables for v in ['magx','magy','magz']):
             self.variables['magnetic energy']               = 'emag'
+            self.variables['magnetic pressure']             = 'emag'
         if all(v in self.variables for v in ['velx','vely','velz', 'magx','magy','magz']) in self.variables:
             self.variables['electric field']                = 'evec'
             self.variables['electric field magnitude']      = 'elec'
diff --git a/sources/equations.F90 b/sources/equations.F90
index 9f187ed..d2e120f 100644
--- a/sources/equations.F90
+++ b/sources/equations.F90
@@ -75,6 +75,10 @@ module equations
       real(kind=8), dimension(:,:,:,:), intent(in) :: qq
       real(kind=8) :: maxspeed
     end function
+    subroutine get_maximum_speeds_iface(qq, um, cm)
+      real(kind=8), dimension(:,:,:,:), intent(in) :: qq
+      real(kind=8)                    , intent(out):: um, cm
+    end subroutine
     subroutine esystem_roe_iface(x, y, q, c, r, l)
       real(kind=8)                , intent(in)    :: x, y
       real(kind=8), dimension(:)  , intent(in)    :: q
@@ -101,6 +105,9 @@ module equations
 !
   procedure(maxspeed_iface)   , pointer, save :: maxspeed        => null()
 
+  procedure(get_maximum_speeds_iface), pointer, save :: &
+                                                 get_maximum_speeds => null()
+
 ! pointer to the Roe eigensystem procedure
 !
   procedure(esystem_roe_iface), pointer, save :: eigensystem_roe => null()
@@ -164,9 +171,10 @@ module equations
 !
   real(kind=8), save :: csnd    = 1.0d+00, csnd2   = 1.0d+00
 
-! maximum speed in the system
+! maximum speeds in the system
 !
   real(kind=8), save :: cmax    = 0.0d+00, cmax2   = 0.0d+00
+  real(kind=8), save :: umax    = 0.0d+00, cglm    = 0.0d+00
 
 ! the lower limits for density and pressure to be treated as physical
 !
@@ -222,13 +230,13 @@ module equations
   public :: prim2cons_mhd_adi, fluxspeed_mhd_adi
   public :: prim2cons_srhd_adi, fluxspeed_srhd_adi
   public :: prim2cons_srmhd_adi, fluxspeed_srmhd_adi
-  public :: maxspeed, reset_maxspeed, get_maxspeed
+  public :: maxspeed, reset_maxspeed, get_maximum_speeds
   public :: eigensystem_roe
   public :: update_primitive_variables
   public :: fix_unphysical_cells, correct_unphysical_states
   public :: adiabatic_index, relativistic, magnetized
   public :: csnd, csnd2
-  public :: cmax, cmax2
+  public :: cmax, cmax2, umax, cglm
   public :: nv, nf, ns
   public :: inx, iny, inz
   public :: idn, ivx, ivy, ivz, imx, imy, imz
@@ -361,11 +369,12 @@ module equations
 
 ! set pointers to subroutines
 !
-        prim2cons       => prim2cons_hd_iso
-        cons2prim       => cons2prim_hd_iso
-        fluxspeed       => fluxspeed_hd_iso
-        maxspeed        => maxspeed_hd_iso
-        eigensystem_roe => esystem_roe_hd_iso
+        prim2cons          => prim2cons_hd_iso
+        cons2prim          => cons2prim_hd_iso
+        fluxspeed          => fluxspeed_hd_iso
+        maxspeed           => maxspeed_hd_iso
+        eigensystem_roe    => esystem_roe_hd_iso
+        get_maximum_speeds => get_maximum_speeds_hd_iso
 
       case("adi", "ADI", "adiabatic", "ADIABATIC")
 
@@ -401,11 +410,12 @@ module equations
 
 ! set pointers to subroutines
 !
-        prim2cons       => prim2cons_hd_adi
-        cons2prim       => cons2prim_hd_adi
-        fluxspeed       => fluxspeed_hd_adi
-        maxspeed        => maxspeed_hd_adi
-        eigensystem_roe => esystem_roe_hd_adi
+        prim2cons          => prim2cons_hd_adi
+        cons2prim          => cons2prim_hd_adi
+        fluxspeed          => fluxspeed_hd_adi
+        maxspeed           => maxspeed_hd_adi
+        eigensystem_roe    => esystem_roe_hd_adi
+        get_maximum_speeds => get_maximum_speeds_hd_adi
 
 ! warn about the unimplemented equation of state
 !
@@ -527,11 +537,12 @@ module equations
 
 ! set pointers to the subroutines
 !
-        prim2cons       => prim2cons_mhd_iso
-        cons2prim       => cons2prim_mhd_iso
-        fluxspeed       => fluxspeed_mhd_iso
-        maxspeed        => maxspeed_mhd_iso
-        eigensystem_roe => esystem_roe_mhd_iso
+        prim2cons          => prim2cons_mhd_iso
+        cons2prim          => cons2prim_mhd_iso
+        fluxspeed          => fluxspeed_mhd_iso
+        maxspeed           => maxspeed_mhd_iso
+        eigensystem_roe    => esystem_roe_mhd_iso
+        get_maximum_speeds => get_maximum_speeds_mhd_iso
 
       case("adi", "ADI", "adiabatic", "ADIABATIC")
 
@@ -574,11 +585,12 @@ module equations
 
 ! set pointers to subroutines
 !
-        prim2cons       => prim2cons_mhd_adi
-        cons2prim       => cons2prim_mhd_adi
-        fluxspeed       => fluxspeed_mhd_adi
-        maxspeed        => maxspeed_mhd_adi
-        eigensystem_roe => esystem_roe_mhd_adi
+        prim2cons          => prim2cons_mhd_adi
+        cons2prim          => cons2prim_mhd_adi
+        fluxspeed          => fluxspeed_mhd_adi
+        maxspeed           => maxspeed_mhd_adi
+        eigensystem_roe    => esystem_roe_mhd_adi
+        get_maximum_speeds => get_maximum_speeds_mhd_adi
 
       case default
 
@@ -703,10 +715,11 @@ module equations
 
 ! set pointers to subroutines
 !
-        prim2cons       => prim2cons_srhd_adi
-        cons2prim       => cons2prim_srhd_adi
-        fluxspeed       => fluxspeed_srhd_adi
-        maxspeed        => maxspeed_srhd_adi
+        prim2cons          => prim2cons_srhd_adi
+        cons2prim          => cons2prim_srhd_adi
+        fluxspeed          => fluxspeed_srhd_adi
+        maxspeed           => maxspeed_srhd_adi
+        get_maximum_speeds => get_maximum_speeds_srhd_adi
 
 ! warn about the unimplemented equation of state
 !
@@ -884,10 +897,11 @@ module equations
 
 ! set pointers to subroutines
 !
-        prim2cons       => prim2cons_srmhd_adi
-        cons2prim       => cons2prim_srmhd_adi
-        fluxspeed       => fluxspeed_srmhd_adi
-        maxspeed        => maxspeed_srmhd_adi
+        prim2cons          => prim2cons_srmhd_adi
+        cons2prim          => cons2prim_srmhd_adi
+        fluxspeed          => fluxspeed_srmhd_adi
+        maxspeed           => maxspeed_srmhd_adi
+        get_maximum_speeds => get_maximum_speeds_srmhd_adi
 
 ! warn about the unimplemented equation of state
 !
@@ -1881,6 +1895,58 @@ module equations
 !
 !===============================================================================
 !
+! subroutine GET_MAXIMUM_SPEEDS_HD_ISO:
+! ------------------------------------
+!
+!   Subroutine determines the maximum characteristic speed and eigenvalue
+!   in the input array of primitive variables.
+!
+!   Arguments:
+!
+!     qq - the input array of primitive variables;
+!     vm - the maximum physical speed;
+!     cm - the maximum eigenvalue;
+!
+!===============================================================================
+!
+  subroutine get_maximum_speeds_hd_iso(qq, vm, cm)
+
+    use coordinates, only : nb, ne
+
+    implicit none
+
+    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
+    real(kind=8)                    , intent(out) :: vm, cm
+
+    integer      :: i, j, k = 1
+    real(kind=8) :: vl, vu
+
+!-------------------------------------------------------------------------------
+!
+    vm = 0.0d+00
+    cm = 0.0d+00
+
+#if NDIMS == 3
+    do k = nb, ne
+#endif /* NDIMS == 3 */
+      do j = nb, ne
+        do i = nb, ne
+          vl = minval(qq(ivx:ivz,i,j,k))
+          vu = maxval(qq(ivx:ivz,i,j,k))
+          vm = max(vm, abs(vl), abs(vu))
+          cm = max(cm, abs(vl - csnd), abs(vu + csnd))
+        end do
+      end do
+#if NDIMS == 3
+    end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine get_maximum_speeds_hd_iso
+!
+!===============================================================================
+!
 ! subroutine ESYSTEM_ROE_HD_ISO:
 ! -----------------------------
 !
@@ -2255,6 +2321,60 @@ module equations
 !
 !===============================================================================
 !
+! subroutine GET_MAXIMUM_SPEEDS_HD_ADI:
+! ------------------------------------
+!
+!   Subroutine determines the maximum characteristic speed and eigenvalue
+!   in the input array of primitive variables.
+!
+!   Arguments:
+!
+!     qq - the input array of primitive variables;
+!     vm - the maximum physical speed;
+!     cm - the maximum eigenvalue;
+!
+!===============================================================================
+!
+  subroutine get_maximum_speeds_hd_adi(qq, vm, cm)
+
+    use coordinates, only : nb, ne
+
+    implicit none
+
+    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
+    real(kind=8)                    , intent(out) :: vm, cm
+
+    integer      :: i, j, k = 1
+    real(kind=8) :: vl, vu, cc
+
+!-------------------------------------------------------------------------------
+!
+    vm = 0.0d+00
+    cm = 0.0d+00
+
+#if NDIMS == 3
+    do k = nb, ne
+#endif /* NDIMS == 3 */
+      do j = nb, ne
+        do i = nb, ne
+          vl = minval(qq(ivx:ivz,i,j,k))
+          vu = maxval(qq(ivx:ivz,i,j,k))
+          vm = max(vm, abs(vl), abs(vu))
+
+          cc = sqrt(adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k))
+          cm = max(cm, abs(vl - cc), abs(vu + cc))
+        end do
+      end do
+#if NDIMS == 3
+    end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine get_maximum_speeds_hd_adi
+!
+!===============================================================================
+!
 ! subroutine ESYSTEM_ROE_HD_ADI:
 ! -----------------------------
 !
@@ -2295,7 +2415,7 @@ module equations
     real(kind=8), dimension(5,5), save :: lvec, rvec
 !$omp threadprivate(first, lvec, rvec)
 
-    real(kind=8)  :: vv, vh, c2, cc, vc, nd, fc, fh, f1, f2, fv, fx
+    real(kind=8)  :: vv, vh, c2, cc, vc, fc, fh, f1, f2, fv, fx
 
 !-------------------------------------------------------------------------------
 !
@@ -2682,6 +2802,65 @@ module equations
 !
 !===============================================================================
 !
+! subroutine GET_MAXIMUM_SPEEDS_MHD_ISO:
+! -------------------------------------
+!
+!   Subroutine determines the maximum characteristic speed and eigenvalue
+!   in the input array of primitive variables.
+!
+!   Arguments:
+!
+!     qq - the input array of primitive variables;
+!     vm - the maximum physical speed;
+!     cm - the maximum eigenvalue;
+!
+!===============================================================================
+!
+  subroutine get_maximum_speeds_mhd_iso(qq, vm, cm)
+
+    use coordinates, only : nb, ne
+
+    implicit none
+
+    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
+    real(kind=8)                    , intent(out) :: vm, cm
+
+    integer      :: i, j, k = 1
+    real(kind=8) :: vl, vu, cc, xx
+
+    real(kind=8), dimension(3) :: bb
+
+!-------------------------------------------------------------------------------
+!
+    vm = 0.0d+00
+    cm = 0.0d+00
+
+#if NDIMS == 3
+    do k = nb, ne
+#endif /* NDIMS == 3 */
+      do j = nb, ne
+        do i = nb, ne
+          vl = minval(qq(ivx:ivz,i,j,k))
+          vu = maxval(qq(ivx:ivz,i,j,k))
+          vm = max(vm, abs(vl), abs(vu))
+
+          bb = qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k) / qq(idn,i,j,k)
+          xx = csnd2 + bb(1) + bb(2) + bb(3)
+          cc = max(0.0d+00, xx**2 - 4.0d+00 * csnd2 * bb(1))
+          cc = sqrt(5.0d-01 * (xx + sqrt(cc)))
+          cm = max(cm, abs(vl - cc), abs(vu + cc))
+        end do
+      end do
+#if NDIMS == 3
+    end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine get_maximum_speeds_mhd_iso
+!
+!===============================================================================
+!
 ! subroutine ESYSTEM_ROE_MHD_ISO:
 ! ------------------------------
 !
@@ -2960,24 +3139,16 @@ module equations
 !
   subroutine prim2cons_mhd_adi(q, u, s)
 
-! local variables are not implicit by default
-!
     implicit none
 
-! input/output arguments
-!
     real(kind=8), dimension(:,:), intent(in)  :: q
     real(kind=8), dimension(:,:), intent(out) :: u
     logical     , optional      , intent(in)  :: s
 
-! local variables
-!
     integer      :: i, p
-    real(kind=8) :: ei, ek, em
-!
+    real(kind=8) :: ei, ek, em, ep
+
 !-------------------------------------------------------------------------------
-!
-! iterate over all positions
 !
     do i = 1, size(q,2)
 
@@ -2990,15 +3161,13 @@ module equations
       u(ibz,i) = q(ibz,i)
       u(ibp,i) = q(ibp,i)
       ei       = gammam1i * q(ipr,i)
-      ek       = 0.5d+00 * (u(imx,i) * q(ivx,i) + u(imy,i) * q(ivy,i)          &
-                                                + u(imz,i) * q(ivz,i))
-      em       = 0.5d+00 * sum(q(ibx:ibz,i) * q(ibx:ibz,i))
-      u(ien,i) = ei + ek + em
+      ek       = sum(u(imx:imz,i) * q(ivx:ivz,i))
+      em       = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
+      ep       =     q(ibp,i)     * q(ibp,i)
+      u(ien,i) = ei + 5.0d-01 * (ek + em + ep)
 
     end do
 
-! update primitive passive scalars
-!
     if (ns > 0 .and. present(s)) then
       if (s) then
         do p = isl, isu
@@ -3036,7 +3205,7 @@ module equations
     integer                     , intent(out) :: s
 
     integer      :: i, p
-    real(kind=8) :: ei, ek, em
+    real(kind=8) :: ei, ek, em, ep
 
 !-------------------------------------------------------------------------------
 !
@@ -3053,9 +3222,10 @@ module equations
         q(iby,i) = u(iby,i)
         q(ibz,i) = u(ibz,i)
         q(ibp,i) = u(ibp,i)
-        ek       = 0.5d+00 * sum(u(imx:imz,i) * q(ivx:ivz,i))
-        em       = 0.5d+00 * sum(q(ibx:ibz,i) * q(ibx:ibz,i))
-        ei       = u(ien,i) - (ek + em)
+        ek       = sum(u(imx:imz,i) * q(ivx:ivz,i))
+        em       = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
+        ep       =     q(ibp,i)     * q(ibp,i)
+        ei       = u(ien,i) - 5.0d-01 * (ek + em + ep)
         if (ei > 0.0d+00) then
           q(ipr,i) = gammam1 * ei
         else
@@ -3245,6 +3415,66 @@ module equations
 !
 !===============================================================================
 !
+! subroutine GET_MAXIMUM_SPEEDS_MHD_ADI:
+! -------------------------------------
+!
+!   Subroutine determines the maximum characteristic speed and eigenvalue
+!   in the input array of primitive variables.
+!
+!   Arguments:
+!
+!     qq - the input array of primitive variables;
+!     vm - the maximum physical speed;
+!     cm - the maximum eigenvalue;
+!
+!===============================================================================
+!
+  subroutine get_maximum_speeds_mhd_adi(qq, vm, cm)
+
+    use coordinates, only : nb, ne
+
+    implicit none
+
+    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
+    real(kind=8)                    , intent(out) :: vm, cm
+
+    integer      :: i, j, k = 1
+    real(kind=8) :: vl, vu, aa, cc, xx
+
+    real(kind=8), dimension(3) :: bb
+
+!-------------------------------------------------------------------------------
+!
+    vm = 0.0d+00
+    cm = 0.0d+00
+
+#if NDIMS == 3
+    do k = nb, ne
+#endif /* NDIMS == 3 */
+      do j = nb, ne
+        do i = nb, ne
+          vl = minval(qq(ivx:ivz,i,j,k))
+          vu = maxval(qq(ivx:ivz,i,j,k))
+          vm = max(vm, abs(vl), abs(vu))
+
+          aa = adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k)
+          bb = qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k) / qq(idn,i,j,k)
+          xx = aa + bb(1) + bb(2) + bb(3)
+          cc = max(0.0d+00, xx**2 - 4.0d+00 * aa * bb(1))
+          cc = sqrt(5.0d-01 * (xx + sqrt(cc)))
+          cm = max(cm, abs(vl - cc), abs(vu + cc))
+        end do
+      end do
+#if NDIMS == 3
+    end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine get_maximum_speeds_mhd_adi
+!
+!===============================================================================
+!
 ! subroutine ESYSTEM_ROE_MHD_ADI:
 ! ------------------------------
 !
@@ -3869,6 +4099,65 @@ module equations
 !
 !===============================================================================
 !
+! subroutine GET_MAXIMUM_SPEEDS_SRHD_ADI:
+! --------------------------------------
+!
+!   Subroutine determines the maximum characteristic speed and eigenvalue
+!   in the input array of primitive variables.
+!
+!   Arguments:
+!
+!     qq - the input array of primitive variables;
+!     vm - the maximum physical speed;
+!     cm - the maximum eigenvalue;
+!
+!===============================================================================
+!
+  subroutine get_maximum_speeds_srhd_adi(qq, vm, cm)
+
+    use coordinates, only : nb, ne
+
+    implicit none
+
+    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
+    real(kind=8)                    , intent(out) :: vm, cm
+
+    integer      :: i, j, k = 1
+    real(kind=8) :: vl, vu, vv, ww, aa, cc, ss, fc
+
+!-------------------------------------------------------------------------------
+!
+    vm = 0.0d+00
+    cm = 0.0d+00
+
+#if NDIMS == 3
+    do k = nb, ne
+#endif /* NDIMS == 3 */
+      do j = nb, ne
+        do i = nb, ne
+          vl = minval(qq(ivx:ivz,i,j,k))
+          vu = maxval(qq(ivx:ivz,i,j,k))
+          vm = max(vm, abs(vl), abs(vu))
+          vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k))
+
+          ww = qq(idn,i,j,k) + qq(ipr,i,j,k) / gammaxi
+          aa = adiabatic_index * qq(ipr,i,j,k) / ww
+          ss = aa * (1.0d+00 - vv) / (1.0d+00 - aa)
+          fc = 1.0d+00 + ss
+          cc = sqrt(ss * (fc - vv))
+          cm = max(cm, abs(vl - cc), abs(vu + cc))
+        end do
+      end do
+#if NDIMS == 3
+    end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine get_maximum_speeds_srhd_adi
+!
+!===============================================================================
+!
 ! subroutine NR_FUNCTION_SRHD_ADI_1D:
 ! ----------------------------------
 !
@@ -5073,6 +5362,57 @@ module equations
 !
 !===============================================================================
 !
+! subroutine GET_MAXIMUM_SPEEDS_SRMHD_ADI:
+! ---------------------------------------
+!
+!   Subroutine determines the maximum characteristic speed and eigenvalue
+!   in the input array of primitive variables.
+!
+!   Arguments:
+!
+!     qq - the input array of primitive variables;
+!     vm - the maximum physical speed;
+!     cm - the maximum eigenvalue;
+!
+!===============================================================================
+!
+  subroutine get_maximum_speeds_srmhd_adi(qq, vm, cm)
+
+    use coordinates, only : nb, ne
+
+    implicit none
+
+    real(kind=8), dimension(:,:,:,:), intent(in)  :: qq
+    real(kind=8)                    , intent(out) :: vm, cm
+
+    integer      :: i, j, k = 1
+    real(kind=8) :: vl, vu
+
+!-------------------------------------------------------------------------------
+!
+    vm = 0.0d+00
+    cm = 1.0d+00
+
+#if NDIMS == 3
+    do k = nb, ne
+#endif /* NDIMS == 3 */
+      do j = nb, ne
+        do i = nb, ne
+          vl = minval(qq(ivx:ivz,i,j,k))
+          vu = maxval(qq(ivx:ivz,i,j,k))
+          vm = max(vm, abs(vl), abs(vu))
+        end do
+      end do
+#if NDIMS == 3
+    end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine get_maximum_speeds_srmhd_adi
+!
+!===============================================================================
+!
 ! subroutine NR_VELOCITY_SRMHD_ADI_1D:
 ! -----------------------------------
 !
diff --git a/sources/evolution.F90 b/sources/evolution.F90
index 613b0a3..4a265c4 100644
--- a/sources/evolution.F90
+++ b/sources/evolution.F90
@@ -1170,7 +1170,7 @@ module evolution
 #if NDIMS == 3
     use coordinates, only : adz
 #endif /* NDIMS == 3 */
-    use equations  , only : nf, maxspeed, cmax, cmax2
+    use equations  , only : nf, get_maximum_speeds, cmax2, cmax, umax, cglm
     use helpers    , only : print_message
 #ifdef MPI
     use mpitools   , only : reduce_maximum, reduce_sum
@@ -1183,7 +1183,7 @@ module evolution
 
     logical                    :: flag
     integer                    :: l, m, n, status
-    real(kind=8)               :: cm, dx_min, fnorm, h0, h1
+    real(kind=8)               :: um, cm, dx_min, fnorm, h0, h1
     real(kind=8), dimension(3) :: d
 
     real(kind=8), dimension(:,:), allocatable :: s
@@ -1200,6 +1200,7 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+    umax = 0.0d+00
     cmax = eps
     m    = 1
     n    = get_dblocks()
@@ -1208,7 +1209,9 @@ module evolution
     do l = 1, n
       pdata => data_blocks(l)%ptr
 
-      cm   = maxspeed(pdata%q)
+      call get_maximum_speeds(pdata%q, um, cm)
+!$omp atomic update
+      umax = max(umax, um)
 !$omp atomic update
       cmax = max(cmax, cm)
 !$omp atomic update
@@ -1217,10 +1220,12 @@ module evolution
 !$omp end parallel do
 
 #ifdef MPI
+    call reduce_maximum(umax)
     call reduce_maximum(cmax)
     call reduce_maximum(m)
 #endif /* MPI */
 
+    cglm  = cmax - umax
     cmax2 = cmax * cmax
 
 #if NDIMS == 2
@@ -1398,7 +1403,7 @@ module evolution
 #if NDIMS == 3
     use coordinates   , only : adz
 #endif /* NDIMS == 3 */
-    use equations     , only : maxspeed, cmax, cmax2
+    use equations     , only : get_maximum_speeds, cmax, cmax2, umax, cglm
 #ifdef MPI
     use mpitools      , only : reduce_maximum
 #endif /* MPI */
@@ -1409,12 +1414,13 @@ module evolution
     type(block_data), pointer :: pdata
 
     integer      :: l, n, m
-    real(kind=8) :: cm, dx_min
+    real(kind=8) :: um, cm, dx_min
 
     real(kind=8), parameter :: eps = tiny(cmax)
 
 !-------------------------------------------------------------------------------
 !
+    umax = 0.0d+00
     cmax = eps
     m    = 1
     n    = get_dblocks()
@@ -1423,7 +1429,9 @@ module evolution
     do l = 1, n
       pdata => data_blocks(l)%ptr
 
-      cm   = maxspeed(pdata%q)
+      call get_maximum_speeds(pdata%q, um, cm)
+!$omp atomic update
+      umax = max(umax, um)
 !$omp atomic update
       cmax = max(cmax, cm)
 !$omp atomic update
@@ -1432,10 +1440,12 @@ module evolution
 !$omp end parallel do
 
 #ifdef MPI
+    call reduce_maximum(umax)
     call reduce_maximum(cmax)
     call reduce_maximum(m)
 #endif /* MPI */
 
+    cglm  = cmax - umax
     cmax2 = cmax * cmax
 
 #if NDIMS == 2
diff --git a/sources/schemes.F90 b/sources/schemes.F90
index c916316..628cdad 100644
--- a/sources/schemes.F90
+++ b/sources/schemes.F90
@@ -2060,11 +2060,12 @@ module schemes
               wcl(:) = (sml - sl) * uil(:) + wl(:)
               wcr(:) = (smr - sr) * uir(:) + wr(:)
 
-              f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
+              dvl = smr - sml
+              f(:,i) = sml * wcr(:) / dvl - smr * wcl(:) / dvl
 
             else ! Bx = 0 -> Sₘ = 0
 
-              f(:,i) = (sl * wr(:) - sr *wl(:)) / srml
+              f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
 
             end if
 
@@ -2085,11 +2086,12 @@ module schemes
 
             wcl(:) = (sml - sl) * uil(:) + wl(:)
 
-            f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml)
+            dvl = sr - sml
+            f(:,i) = sml * wr(:) / dvl - sr * wcl(:) / dvl
 
           else ! Bx = 0 -> Sₘ = 0
 
-            f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
+            f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
 
           end if ! sl* < 0
 
@@ -2108,17 +2110,18 @@ module schemes
 
             wcr(:) = (smr - sr) * uir(:) + wr(:)
 
-            f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl)
+            dvr = smr - sl
+            f(:,i) = sl * wcr(:) / dvr - smr * wl(:) / dvr
 
           else ! Bx = 0 -> Sₘ = 0
 
-            f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
+            f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
 
           end if ! sr* > 0
 
         else ! sl* < sl & sr* > sr
 
-          f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
+          f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
 
         end if
 
@@ -3244,8 +3247,8 @@ module schemes
 !
   subroutine riemann_hd_iso_kepes(ql, qr, ul, ur, fl, fr, cl, cr, f)
 
-    use equations, only : idn, ivx, ivy, ivz
-    use equations, only : csnd, csnd2
+    use equations, only : idn, ivx, ivy, ivz, imx, imy, imz
+    use equations, only : nf, csnd, csnd2
 
     implicit none
 
@@ -3256,18 +3259,18 @@ module schemes
     real(kind=8), dimension(4,4), save :: rm
 !$omp threadprivate(first, rm)
 
-    integer      :: i
-    real(kind=8) :: dnl, vxa, vya, vza
+    integer      :: i, p
+    real(kind=8) :: dnl, dna, vxa, vya, vza
 
-    real(kind=8), dimension(4) :: v, lm, tm
+    real(kind=8), dimension(nf) :: v, lm, tm, lv
 
 !-------------------------------------------------------------------------------
 !
     if (first) then
-      rm(:,:) = 0.0d+00
       rm(1,:) = [ 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
-      rm(3,2) = 1.0d+00
-      rm(4,3) = 1.0d+00
+      rm(2,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
+      rm(3,:) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
+      rm(4,:) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
 
       first = .false.
     end if
@@ -3277,16 +3280,25 @@ module schemes
 ! the intermediate state averages
 !
       dnl = lmean(ql(idn,i), qr(idn,i))
+      dna = amean(ql(idn,i), qr(idn,i))
       vxa = amean(ql(ivx,i), qr(ivx,i))
       vya = amean(ql(ivy,i), qr(ivy,i))
       vza = amean(ql(ivz,i), qr(ivz,i))
 
+! the difference in entropy variables between states
+!
+      v(ivx) = qr(ivx,i) - ql(ivx,i)
+      v(ivy) = qr(ivy,i) - ql(ivy,i)
+      v(ivz) = qr(ivz,i) - ql(ivz,i)
+      v(idn) = csnd2 * (log(qr(idn,i)) - log(ql(idn,i))) &
+             - (vxa * v(ivx) + vya * v(ivy) + vza * v(ivz))
+
 ! the eigenvalue diagonal matrix
 !
-      lm(1) = vxa - csnd
+      lm(1) = vxa + csnd
       lm(2) = vxa
       lm(3) = vxa
-      lm(4) = vxa + csnd
+      lm(4) = vxa - csnd
 
 ! the scaling diagonal matrix
 !
@@ -3297,32 +3309,32 @@ module schemes
 
 ! the average of the right eigenvector matrix
 !
-      rm(2,1) = vxa - csnd
-      rm(2,4) = vxa + csnd
+      rm(2,1) = vxa + csnd
       rm(3,1) = vya
-      rm(3,4) = vya
       rm(4,1) = vza
+      rm(2,4) = vxa - csnd
+      rm(3,4) = vya
       rm(4,4) = vza
 
-! the difference in entropy variables between states
-!
-      v(1) =   csnd2 * diff(log(ql(idn,i))       , log(qr(idn,i))       ) &
-           - 5.0d-01 * diff(sum(ql(ivx:ivz,i)**2), sum(qr(ivx:ivz,i)**2))
-      v(2) = diff(ql(ivx,i), qr(ivx,i))
-      v(3) = diff(ql(ivy,i), qr(ivy,i))
-      v(4) = diff(ql(ivz,i), qr(ivz,i))
-
 ! KEPEC flux
 !
       f(idn,i) = dnl * vxa
-      f(ivx,i) = f(idn,i) * vxa + csnd2 * dnl
-      f(ivy,i) = f(idn,i) * vya
-      f(ivz,i) = f(idn,i) * vza
+      f(imx,i) = f(idn,i) * vxa + csnd2 * dna
+      f(imy,i) = f(idn,i) * vya
+      f(imz,i) = f(idn,i) * vza
 
 ! KEPES correction
 !
-      f(:,i) = f(:,i) &
-             - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(transpose(rm), v))
+      lv = 5.0d-01 * (abs(lm) * tm) * matmul(v, rm)
+      if (vxa >= 0.0d+00) then
+        do p = 1, nf
+          f(:,i) = f(:,i) - lv(p) * rm(:,p)
+        end do
+      else
+        do p = nf, 1, -1
+          f(:,i) = f(:,i) - lv(p) * rm(:,p)
+        end do
+      end if
 
     end do
 
@@ -3357,8 +3369,8 @@ module schemes
 !
   subroutine riemann_hd_adi_kepes(ql, qr, ul, ur, fl, fr, cl, cr, f)
 
-    use equations, only : idn, ivx, ivy, ivz, ipr
-    use equations, only : adiabatic_index
+    use equations, only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien
+    use equations, only : nf, adiabatic_index
 
     implicit none
 
@@ -3370,11 +3382,11 @@ module schemes
     real(kind=8), dimension(5,5), save :: rm
 !$omp threadprivate(first, adi_m1, adi_m1x, rm)
 
-    integer      :: i
-    real(kind=8) :: btl, btr, bta, dnl, prl, pra, vxa, vya, vza, vva
-    real(kind=8) :: csn, ent, sl, sr, tl, tr
+    integer      :: i, p
+    real(kind=8) :: bl, br, btl, dnl, prl, bta, dna, vxa, vya, vza, vva, pra
+    real(kind=8) :: csn, ent
 
-    real(kind=8), dimension(5) :: v, lm, tm
+    real(kind=8), dimension(nf) :: v, lm, tm, lv
 
 !-------------------------------------------------------------------------------
 !
@@ -3382,39 +3394,58 @@ module schemes
       adi_m1  = adiabatic_index - 1.0d+00
       adi_m1x = adi_m1 / adiabatic_index
 
-      rm(:,:) = 0.0d+00
       rm(1,:) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
-      rm(3,3) = 1.0d+00
-      rm(4,4) = 1.0d+00
+      rm(2,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
+      rm(3,:) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
+      rm(4,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
+      rm(5,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
 
       first = .false.
     end if
 
     do i = 1, size(ql, 2)
 
+! various parameters
+!
+      bl = ql(idn,i) / ql(ipr,i)
+      br = qr(idn,i) / qr(ipr,i)
+
 ! the intermediate state averages
 !
-      btl = ql(idn,i) / ql(ipr,i)
-      btr = qr(idn,i) / qr(ipr,i)
-      bta = lmean(btl, btr)
+      btl = lmean(bl, br)
+      bta = amean(bl, br)
       dnl = lmean(ql(idn,i), qr(idn,i))
+      dna = amean(ql(idn,i), qr(idn,i))
       vxa = amean(ql(ivx,i), qr(ivx,i))
       vya = amean(ql(ivy,i), qr(ivy,i))
       vza = amean(ql(ivz,i), qr(ivz,i))
-      pra = dnl / amean(btl, btr)
-      prl = dnl / bta
+      pra = dna / bta
+      prl = dnl / btl
       vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
+      ent = 1.0d+00 / adi_m1x / btl + vva
 
-! the eigenvalue diagonal matrix
+! the difference in entropy variables between states
+!
+      v(idn) = (log(qr(idn,i)) - log(ql(idn,i))) &
+             + (log(br) - log(bl)) / adi_m1 &
+             - 5.0d-01 * (br * sum(qr(ivx:ivz,i)**2) &
+                        - bl * sum(ql(ivx:ivz,i)**2))
+      v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
+      v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
+      v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
+      v(ipr) = bl - br
+
+! the eigenvalues
 !
       csn = sqrt(adiabatic_index * pra / dnl)
-      ent = 1.0d+00 / adi_m1x / bta + vva
 
-      lm(1) = vxa - csn
+! the diagonal matrix of eigenvalues
+!
+      lm(1) = vxa + csn
       lm(2) = vxa
       lm(3) = vxa
       lm(4) = vxa
-      lm(5) = vxa + csn
+      lm(5) = vxa - csn
 
 ! the scaling diagonal matrix
 !
@@ -3426,43 +3457,37 @@ module schemes
 
 ! the average of the right eigenvector matrix
 !
-      rm(2,1) = vxa - csn
-      rm(2,2) = vxa
-      rm(2,5) = vxa + csn
+      rm(2,1) = vxa + csn
       rm(3,1) = vya
-      rm(3,2) = vya
-      rm(3,5) = vya
       rm(4,1) = vza
+      rm(2,2) = vxa
+      rm(3,2) = vya
       rm(4,2) = vza
+      rm(2,5) = vxa - csn
+      rm(3,5) = vya
       rm(4,5) = vza
-      rm(5,:) = [ ent - vxa * csn, vva, vya, vza, ent + vxa * csn ]
-
-! the difference in entropy variables between states
-!
-      sl = - adi_m1 * log(ql(idn,i)) - log(btl)
-      sr = - adi_m1 * log(qr(idn,i)) - log(btr)
-      tl = (adiabatic_index - sl) / adi_m1 &
-                            - 5.0d-01 * btl * sum(ql(ivx:ivz,i)**2)
-      tr = (adiabatic_index - sr) / adi_m1 &
-                            - 5.0d-01 * btr * sum(qr(ivx:ivz,i)**2)
-      v(1) = diff(tl, tr)
-      v(2) = diff(btl * ql(ivx,i), btr * qr(ivx,i))
-      v(3) = diff(btl * ql(ivy,i), btr * qr(ivy,i))
-      v(4) = diff(btl * ql(ivz,i), btr * qr(ivz,i))
-      v(5) = - diff(btl, btr)
+      rm(5,:) = [ ent + vxa * csn, vva, vya, vza, ent - vxa * csn ]
 
 ! KEPEC flux
 !
       f(idn,i) = dnl * vxa
-      f(ivx,i) = f(idn,i) * vxa + pra
-      f(ivy,i) = f(idn,i) * vya
-      f(ivz,i) = f(idn,i) * vza
-      f(ipr,i) = (prl / adi_m1 + pra + dnl * vva) * vxa
+      f(imx,i) = f(idn,i) * vxa + pra
+      f(imy,i) = f(idn,i) * vya
+      f(imz,i) = f(idn,i) * vza
+      f(ien,i) = (prl / adi_m1 + pra + dnl * vva) * vxa
 
 ! KEPES correction
 !
-      f(:,i) = f(:,i) &
-             - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(transpose(rm), v))
+      lv = 5.0d-01 * (abs(lm) * tm) * matmul(v, rm)
+      if (vxa >= 0.0d+00) then
+        do p = 1, nf
+          f(:,i) = f(:,i) - lv(p) * rm(:,p)
+        end do
+      else
+        do p = nf, 1, -1
+          f(:,i) = f(:,i) - lv(p) * rm(:,p)
+        end do
+      end if
 
     end do
 
@@ -3500,7 +3525,7 @@ module schemes
   subroutine riemann_mhd_iso_kepes(ql, qr, ul, ur, fl, fr, cl, cr, f)
 
     use equations, only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp
-    use equations, only : csnd, csnd2, cmax
+    use equations, only : nf, csnd, csnd2, ch => cglm
 
     implicit none
 
@@ -3511,15 +3536,15 @@ module schemes
     real(kind=8), dimension(8,8), save :: rm
 !$omp threadprivate(first, rm)
 
-    integer      :: i
+    integer      :: i, p
     real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa
-    real(kind=8) :: dnl, bp
+    real(kind=8) :: dnl, sqd, bp
     real(kind=8) :: bb2, bp2, b1, b2, b3, bb
     real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3
     real(kind=8) :: alf2, als2, alf, als
     real(kind=8) :: f1, f2, f3, sgnb1, v2l, v2r, b2l, b2r
 
-    real(kind=8), dimension(8) :: v, lm, tm
+    real(kind=8), dimension(nf) :: v, lm, tm, lv
 
 !-------------------------------------------------------------------------------
 !
@@ -3535,24 +3560,23 @@ module schemes
 
     do i = 1, size(ql, 2)
 
-! the difference in entropy variables between the states
+! various parameters
 !
       v2l = sum(ql(ivx:ivz,i)**2)
       v2r = sum(qr(ivx:ivz,i)**2)
       b2l = sum(ql(ibx:ibz,i)**2)
       b2r = sum(qr(ibx:ibz,i)**2)
 
-      f1 = csnd2 * (log(ql(idn,i)) + 1.0d+00) - 5.0d-01 * v2l
-      f2 = csnd2 * (log(qr(idn,i)) + 1.0d+00) - 5.0d-01 * v2r
-
-      v(idn) = diff(f1, f2)
-      v(ivx) = diff(ql(ivx,i), qr(ivx,i))
-      v(ivy) = diff(ql(ivy,i), qr(ivy,i))
-      v(ivz) = diff(ql(ivz,i), qr(ivz,i))
-      v(ibx) = diff(ql(ibx,i), qr(ibx,i))
-      v(iby) = diff(ql(iby,i), qr(iby,i))
-      v(ibz) = diff(ql(ibz,i), qr(ibz,i))
-      v(ibp) = diff(ql(ibp,i), qr(ibp,i))
+! the difference in entropy variables between the states
+!
+      v(idn) = csnd2 * (log(qr(idn,i)) - log(ql(idn,i))) - 5.0d-01 * (v2r - v2l)
+      v(ivx) = qr(ivx,i) - ql(ivx,i)
+      v(ivy) = qr(ivy,i) - ql(ivy,i)
+      v(ivz) = qr(ivz,i) - ql(ivz,i)
+      v(ibx) = qr(ibx,i) - ql(ibx,i)
+      v(iby) = qr(iby,i) - ql(iby,i)
+      v(ibz) = qr(ibz,i) - ql(ibz,i)
+      v(ibp) = qr(ibp,i) - ql(ibp,i)
 
 ! the intermediate state averages
 !
@@ -3569,9 +3593,10 @@ module schemes
 
 ! the eigenvalues and related parameters
 !
-      b1  = bxa / sqrt(dnl)
-      b2  = bya / sqrt(dnl)
-      b3  = bza / sqrt(dnl)
+      sqd = sqrt(dnl)
+      b1  = bxa / sqd
+      b2  = bya / sqd
+      b3  = bza / sqd
       bp2 = b2**2 + b3**2
       bp  = sqrt(bp2)
       bb2 = b1**2 + bp2
@@ -3624,7 +3649,7 @@ module schemes
 !
       f1 = dnl * alf
       f2 = dnl * als * cs * sgnb1
-      f3 = als * csnd * sqrt(dnl)
+      f3 = als * csnd * sqd
       rm(1,1) = f1
       rm(2,1) = f1 * (vxa + cf)
       rm(3,1) = f1 * vya - f2 * x2
@@ -3656,7 +3681,7 @@ module schemes
 !
       f1 = dnl * als
       f2 = dnl * alf * cf * sgnb1
-      f3 = - alf * csnd * sqrt(dnl)
+      f3 = - alf * csnd * sqd
       rm(1,3) = f1
       rm(2,3) = f1 * (vxa + cs)
       rm(3,3) = f1 * vya + f2 * x2
@@ -3676,8 +3701,8 @@ module schemes
       lm(1) = vxa + cf
       lm(2) = vxa + ca
       lm(3) = vxa + cs
-      lm(4) = vxa + cmax
-      lm(5) = vxa - cmax
+      lm(4) = vxa + ch
+      lm(5) = vxa - ch
       lm(6) = vxa - cs
       lm(7) = vxa - ca
       lm(8) = vxa - cf
@@ -3688,15 +3713,23 @@ module schemes
       f(ivx,i) = f(idn,i) * vxa - bxa * bxa + pta
       f(ivy,i) = f(idn,i) * vya - bxa * bya
       f(ivz,i) = f(idn,i) * vza - bxa * bza
-      f(ibx,i) = cmax * bpa
+      f(ibx,i) = ch  * bpa
       f(iby,i) = vxa * bya - bxa * vya
       f(ibz,i) = vxa * bza - bxa * vza
-      f(ibp,i) = cmax * bxa
+      f(ibp,i) = ch  * bxa
 
 ! KEPES correction
 !
-      f(:,i) = f(:,i) &
-             - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(transpose(rm), v))
+      lv = 5.0d-01 * (abs(lm) * tm) * matmul(v, rm)
+      if (vxa >= 0.0d+00) then
+        do p = 1, nf
+          f(:,i) = f(:,i) - lv(p) * rm(:,p)
+        end do
+      else
+        do p = nf, 1, -1
+          f(:,i) = f(:,i) - lv(p) * rm(:,p)
+        end do
+      end if
 
     end do
 
@@ -3707,7 +3740,7 @@ module schemes
 !===============================================================================
 !
 ! subroutine RIEMANN_MHD_ADI_KEPES:
-! -------------------------------
+! --------------------------------
 !
 !   Subroutine solves one dimensional adiabatic magnetohydrodynamic
 !   Riemann problem using the entropy stable KEPES method.
@@ -3732,7 +3765,8 @@ module schemes
   subroutine riemann_mhd_adi_kepes(ql, qr, ul, ur, fl, fr, cl, cr, f)
 
     use equations, only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
-    use equations, only : adiabatic_index, cmax
+    use equations, only : imx, imy, imz, ien
+    use equations, only : nf, adiabatic_index, ch => cglm
 
     implicit none
 
@@ -3744,7 +3778,7 @@ module schemes
     real(kind=8), dimension(9,9), save :: rm
 !$omp threadprivate(first, adi_m1, adi_m1x, rm)
 
-    integer      :: i
+    integer      :: i, p
     real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa
     real(kind=8) :: dnl, prl, pta, vva, br, bl, bp
     real(kind=8) :: btl, bta, eka, ema, ub2, uba
@@ -3753,7 +3787,7 @@ module schemes
     real(kind=8) :: alf2, als2, alf, als, wps, wms, wpf, wmf
     real(kind=8) :: f1, f2, f3, f4, sgnb1, v2l, v2r, b2l, b2r
 
-    real(kind=8), dimension(9) :: v, lm, tm
+    real(kind=8), dimension(nf) :: v, lm, tm
 
 !-------------------------------------------------------------------------------
 !
@@ -3775,27 +3809,25 @@ module schemes
 
 ! the difference in entropy variables between the states
 !
-      v2l = sum(ql(ivx:ivz,i)**2)
-      v2r = sum(qr(ivx:ivz,i)**2)
-      b2l = sum(ql(ibx:ibz,i)**2)
-      b2r = sum(qr(ibx:ibz,i)**2)
+      v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
+      v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
+      b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
+      b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
 
       bl  = ql(idn,i) / ql(ipr,i)
       br  = qr(idn,i) / qr(ipr,i)
-      f1  = - adi_m1 * log(ql(idn,i)) - log(bl)
-      f2  = - adi_m1 * log(qr(idn,i)) - log(br)
-      f1  = (adiabatic_index - f1) / adi_m1 - 5.0d-01 * bl * v2l
-      f2  = (adiabatic_index - f2) / adi_m1 - 5.0d-01 * br * v2r
 
-      v(idn) = diff(f1, f2)
-      v(ivx) = diff(bl * ql(ivx,i), br * qr(ivx,i))
-      v(ivy) = diff(bl * ql(ivy,i), br * qr(ivy,i))
-      v(ivz) = diff(bl * ql(ivz,i), br * qr(ivz,i))
-      v(ipr) = -diff(bl, br)
-      v(ibx) = diff(bl * ql(ibx,i), br * qr(ibx,i))
-      v(iby) = diff(bl * ql(iby,i), br * qr(iby,i))
-      v(ibz) = diff(bl * ql(ibz,i), br * qr(ibz,i))
-      v(ibp) = diff(bl * ql(ibp,i), br * qr(ibp,i))
+      v(idn) = (log(qr(idn,i)) - log(ql(idn,i))) &
+             + (log(br) - log(bl)) / adi_m1 &
+             - 5.0d-01 * (br * v2r - bl * v2l)
+      v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
+      v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
+      v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
+      v(ipr) = bl - br
+      v(ibx) = br * qr(ibx,i) - bl * ql(ibx,i)
+      v(iby) = br * qr(iby,i) - bl * ql(iby,i)
+      v(ibz) = br * qr(ibz,i) - bl * ql(ibz,i)
+      v(ibp) = br * qr(ibp,i) - bl * ql(ibp,i)
 
 ! the intermediate state averages
 !
@@ -3907,8 +3939,8 @@ module schemes
       rm(7,1) = f3 * x2
       rm(8,1) = f3 * x3
 
-      rm(1,9) = rm(1,1)
-      rm(2,9) = rm(1,9) * (vxa - cf)
+      rm(1,9) = f1
+      rm(2,9) = f1 * (vxa - cf)
       rm(3,9) = f1 * vya + f2 * x2
       rm(4,9) = f1 * vza + f2 * x3
       rm(5,9) = wmf
@@ -3968,9 +4000,9 @@ module schemes
       lm(1) = vxa + cf
       lm(2) = vxa + ca
       lm(3) = vxa + cs
-      lm(4) = vxa + cmax
+      lm(4) = vxa + ch
       lm(5) = vxa
-      lm(6) = vxa - cmax
+      lm(6) = vxa - ch
       lm(7) = vxa - cs
       lm(8) = vxa - ca
       lm(9) = vxa - cf
@@ -3978,23 +4010,22 @@ module schemes
 ! KEPEC flux
 !
       f(idn,i) = dnl * vxa
-      f(ivx,i) = f(idn,i) * vxa - bxa * bxa + pta
-      f(ivy,i) = f(idn,i) * vya - bxa * bya
-      f(ivz,i) = f(idn,i) * vza - bxa * bza
-      f(ibx,i) = cmax * bpa
+      f(imx,i) = f(idn,i) * vxa - bxa * bxa + pta
+      f(imy,i) = f(idn,i) * vya - bxa * bya
+      f(imz,i) = f(idn,i) * vza - bxa * bza
+      f(ibx,i) = ch  * bpa
       f(iby,i) = vxa * bya - bxa * vya
       f(ibz,i) = vxa * bza - bxa * vza
-      f(ibp,i) = cmax * bxa
-      f(ipr,i) = f(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) &
-               + f(ivx,i) * vxa + f(ivy,i) * vya + f(ivz,i) * vza &
-               + f(ibx,i) * bxa + f(iby,i) * bya + f(ibz,i) * bza &
-               + f(ibp,i) * bpa - ub2 + bxa * uba &
-               - cmax * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i))
+      f(ibp,i) = ch  * bxa
+      f(ien,i) = f(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) &
+               + (f(imx,i) * vxa + f(imy,i) * vya + f(imz,i) * vza) &
+               + (f(ibx,i) * bxa + f(iby,i) * bya + f(ibz,i) * bza) &
+               +  f(ibp,i) * bpa - ub2 + bxa * uba &
+               - ch * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i))
 
 ! KEPES correction
 !
-      f(:,i) = f(:,i) &
-             - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(transpose(rm), v))
+      f(:,i) = f(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
 
     end do
 
@@ -4107,19 +4138,25 @@ module schemes
 
     real(kind=8), intent(in) :: l, r
 
-    real(kind=8) :: u, rr
+    real(kind=8) :: u, d, s
 
     real(kind=8), parameter :: c1 = 2.0d+00, c2 = 2.0d+00 / 3.0d+00, &
                                c3 = 4.0d-01, c4 = 2.0d+00 / 7.0d+00
 
 !-------------------------------------------------------------------------------
 !
-    rr = r * r
-    u = (l * (l - 2.0d+00 * r) + rr) / (l * (l + 2.0d+00 * r) + rr)
+    d = r - l
+    s = r + l
+    u = (d / s)**2
     if (u < 1.0d-04) then
-      lmean = (l + r) / (c1 + u * (c2 + u * (c3 + u * c4)))
+      lmean = s / (c1 + u * (c2 + u * (c3 + u * c4)))
     else
-      lmean = (r - l) / log(r / l)
+      if (d >= 0.0d+00) then
+        s = log(r / l)
+      else
+        s = log(l / r)
+      end if
+      lmean = abs(d) / s
     end if
 
 !-------------------------------------------------------------------------------
@@ -4152,33 +4189,6 @@ module schemes
 !-------------------------------------------------------------------------------
 !
   end function amean
-!
-!===============================================================================
-!
-! function DIFF:
-! -------------
-!
-!   Function calculate the difference.
-!
-!   Arguments:
-!
-!     l, r - the values to calculate the difference of;
-!
-!===============================================================================
-!
-  real(kind=8) function diff(l, r)
-
-    implicit none
-
-    real(kind=8), intent(in) :: l, r
-
-!-------------------------------------------------------------------------------
-!
-    diff = r - l
-
-!-------------------------------------------------------------------------------
-!
-  end function diff
 
 !===============================================================================
 !
diff --git a/sources/sources.F90 b/sources/sources.F90
index 8976909..9fa7399 100644
--- a/sources/sources.F90
+++ b/sources/sources.F90
@@ -121,6 +121,9 @@ module sources
     case("heglm", "HEGLM")
       glm_type = 2
       glm_name = "HEGLM"
+    case("kepes", "KEPES")
+      glm_type = 3
+      glm_name = "KEPES"
     case default
       glm_type = 0
       glm_name = "none"
@@ -556,7 +559,7 @@ module sources
       dh(2) = ady(pdata%meta%level)
       dh(3) = adz(pdata%meta%level)
 
-! add the EGLM-MHD source terms
+! add the EGLM-MHD and KEPES source terms
 !
       if (glm_type > 0) then
 
@@ -571,23 +574,7 @@ module sources
         du(imy,:,:,:) = du(imy,:,:,:) - db(:,:,:) * pdata%q(iby,:,:,:)
         du(imz,:,:,:) = du(imz,:,:,:) - db(:,:,:) * pdata%q(ibz,:,:,:)
 
-! update the energy equation
-!
-        if (ien > 0 .and. ibp > 0) then
-
-! calculate the gradient of divergence potential
-!
-          call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))
-
-! add the divergence potential source term to the energy equation, i.e.
-!     d/dt E + ∇.F = - B.(∇ψ)
-!
-          du(ien,:,:,:) = du(ien,:,:,:)                                        &
-                     - sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
-
-        end if ! ien > 0
-
-! add the HEGLM-MHD source terms
+! add the HEGLM-MHD and KEPES source terms
 !
         if (glm_type > 1) then
 
@@ -616,7 +603,41 @@ module sources
 
         end if ! glm_type > 1
 
-      end if ! glmtype > 0
+        if  (ibp > 0 .and. (ien > 0 .or. glm_type == 3)) then
+
+! calculate the gradient of divergence correcting field
+!
+          call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))
+
+          db(:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
+
+! update the divergence correcting field
+!     d/dt ψ+ ∇.F = - (v.∇)ψ
+!
+          du(ibp,:,:,:) = du(ibp,:,:,:) - db(:,:,:)
+
+! update the energy equation
+!
+          if (ien > 0) then
+            if  (glm_type == 3) then
+
+! add the divergence potential source term to the energy equation, i.e.
+!     d/dt E + ∇.F = - ψ(v.∇)ψ
+!
+              du(ien,:,:,:) = du(ien,:,:,:) - pdata%q(ibp,:,:,:) * db(:,:,:)
+
+            else
+
+! add the divergence potential source term to the energy equation, i.e.
+!     d/dt E + ∇.F = - B.(∇ψ)
+!
+              du(ien,:,:,:) = du(ien,:,:,:) &
+                       - sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
+
+            end if
+          end if ! ien > 0
+        end if ! glm_type == 3
+      end if ! glm_type > 0
 
 ! if anomalous resistivity is enabled
 !
diff --git a/sources/statistics.F90 b/sources/statistics.F90
index 0f29495..d593949 100644
--- a/sources/statistics.F90
+++ b/sources/statistics.F90
@@ -38,12 +38,14 @@ module statistics
 !   munit - a file handler for the mesh statistics file;
 !   cunit - a file handler for the conserved variable integrals file;
 !   sunit - a file handler for the variable statistics file;
+!   funit - a file handler for the forcing statistics file;
 !   eunit - a file handler for the integration errors file;
 !
   integer(kind=4), save :: munit = 10
   integer(kind=4), save :: cunit = 11
   integer(kind=4), save :: sunit = 12
-  integer(kind=4), save :: eunit = 13
+  integer(kind=4), save :: funit = 13
+  integer(kind=4), save :: eunit = 14
 
 !   iintd - the number of steps between subsequent intervals;
 !
@@ -96,6 +98,7 @@ module statistics
     use coordinates, only : toplev, ncells, bcells, nghosts, domain_base_dims
     use equations  , only : pvars, nf
     use evolution  , only : error_control
+    use forcing    , only : forcing_enabled
     use mpitools   , only : master
 #ifdef MPI
     use mpitools   , only : nprocs
@@ -234,10 +237,9 @@ module statistics
 #endif /* __INTEL_COMPILER */
       end if
 
-      write(cunit,"('#',a8,13(1x,a18))") 'step', 'time', 'dt'                  &
+      write(cunit,"('#',a8,11(1x,a18))") 'step', 'time', 'dt'                  &
                                        , 'mass', 'momx', 'momy', 'momz'        &
-                                       , 'ener', 'ekin', 'emag', 'eint'        &
-                                       , 'einj', 'erat', 'arms'
+                                       , 'ener', 'ekin', 'emag', 'eint', 'entr'
       write(cunit,"('#')")
 
 ! prepare the variable statistics file
@@ -279,6 +281,40 @@ module statistics
                                  , 'mean(Malf)', 'min(Malf)', 'max(Malf)'
       write(sunit,"('#')")
 
+! prepare the forcing statistics file
+!
+      if (forcing_enabled) then
+        if (append) then
+          write(fname, "('forcing-statistics.dat')")
+          inquire(file=fname, exist=flag)
+          flag = flag .and. nrun > 1
+        else
+          write(fname, "('forcing-statistics_',i2.2,'.dat')") nrun
+          flag = .false.
+        end if
+
+        if (flag) then
+#ifdef __INTEL_COMPILER
+          open(newunit=funit, file=fname, form='formatted', status='old', &
+                                          position='append', buffered='yes')
+#else /* __INTEL_COMPILER */
+          open(newunit=funit, file=fname, form='formatted', status='old', &
+                                          position='append')
+#endif /* __INTEL_COMPILER */
+          write(cunit,"('#')")
+        else
+#ifdef __INTEL_COMPILER
+          open(newunit=funit, file=fname, form='formatted', &
+                                          status='replace', buffered='yes')
+#else /* __INTEL_COMPILER */
+          open(newunit=funit, file=fname, form='formatted', status='replace')
+#endif /* __INTEL_COMPILER */
+        end if
+
+        write(funit,"('#',a18,4(1x,a18))") 'time', 'einj', 'erat', 'arms'
+        write(funit,"('#')")
+      end if
+
 ! prepare the integration errors file
 !
       if (error_control) then
@@ -344,6 +380,7 @@ module statistics
   subroutine finalize_statistics(status)
 
     use evolution, only : error_control
+    use forcing  , only : forcing_enabled
     use mpitools , only : master
 
     implicit none
@@ -358,6 +395,7 @@ module statistics
       close(munit)
       close(cunit)
       close(sunit)
+      if (forcing_enabled) close(funit)
       if (error_control) close(eunit)
     end if
 
@@ -386,11 +424,11 @@ module statistics
     use coordinates, only : toplev
     use equations  , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp
     use equations  , only : ien, imx, imy, imz
-    use equations  , only : magnetized, adiabatic_index, csnd
+    use equations  , only : magnetized, adiabatic_index, csnd, csnd2
     use equations  , only : errors
     use evolution  , only : error_control
     use evolution  , only : step, time, dt, dth, dte
-    use forcing    , only : einj, rinj, arms
+    use forcing    , only : forcing_enabled, einj, rinj, arms
     use helpers    , only : print_message, flush_and_sync
     use mesh       , only : changed
     use mpitools   , only : master
@@ -566,7 +604,7 @@ module statistics
       inarr(4) = inarr(4) + sum(pdata%u(imz,nb:ne,nb:ne,  :  )) * dvol
 #endif /* NDIMS == 3 */
 
-! sum up total energy
+! sum up total energy and entropy
 !
       if (ien > 0) then
 #if NDIMS == 3
@@ -574,7 +612,26 @@ module statistics
 #else /* NDIMS == 3 */
         inarr(5) = inarr(5) + sum(pdata%u(ien,nb:ne,nb:ne,  :  )) * dvol
 #endif /* NDIMS == 3 */
+
+#if NDIMS == 3
+        tmp(:,:,:) = - log(pdata%q(ipr,nb:ne,nb:ne,nb:ne) &
+                   / pdata%q(idn,nb:ne,nb:ne,nb:ne)**adiabatic_index) &
+                   * pdata%q(idn,nb:ne,nb:ne,nb:ne) / (adiabatic_index - 1.0d+00)
+#else /* NDIMS == 3 */
+        tmp(:,:,:) = - log(pdata%q(ipr,nb:ne,nb:ne,  :  ) &
+                   / pdata%q(idn,nb:ne,nb:ne,  :  )**adiabatic_index) &
+                   * pdata%q(idn,nb:ne,nb:ne,  :  ) / (adiabatic_index - 1.0d+00)
+#endif /* NDIMS == 3 */
+      else
+#if NDIMS == 3
+        tmp(:,:,:) = csnd2 * log(pdata%q(idn,nb:ne,nb:ne,nb:ne)) &
+                           *     pdata%q(idn,nb:ne,nb:ne,nb:ne)
+#else /* NDIMS == 3 */
+        tmp(:,:,:) = csnd2 * log(pdata%q(idn,nb:ne,nb:ne,  :  )) &
+                           *     pdata%q(idn,nb:ne,nb:ne,  :  )
+#endif /* NDIMS == 3 */
       end if
+      inarr(9) = inarr(9) + sum(tmp(:,:,:)) * dvol
 
 ! sum up kinetic energy
 !
@@ -602,12 +659,21 @@ module statistics
                                 + pdata%u(iby,nb:ne,nb:ne,  :  )**2            &
                                 + pdata%u(ibz,nb:ne,nb:ne,  :  )**2) * dvolh
 #endif /* NDIMS == 3 */
+        if (ien <= 0) then
+#if NDIMS == 3
+          inarr(9) = inarr(9) + sum(pdata%u(ibp,nb:ne,nb:ne,nb:ne)**2) * dvolh
+#else /* NDIMS == 3 */
+          inarr(9) = inarr(9) + sum(pdata%u(ibp,nb:ne,nb:ne,  :  )**2) * dvolh
+#endif /* NDIMS == 3 */
+        end if
       end if
 
 ! sum up the injected energy and injection rate
 !
-      inarr( 9) = einj
-      inarr(10) = rinj
+      if (forcing_enabled) then
+        inarr(11) = einj
+        inarr(12) = rinj
+      end if
 
 ! get average, minimum and maximum values of density
 !
@@ -720,9 +786,13 @@ module statistics
     call reduce_maximum(mxarr(1:narr))
 #endif /* MPI */
 
-! calculate the internal energy
+! calculate the internal energy, or update the entropy
 !
-    if (ien > 0) inarr(8) = inarr(5) - inarr(6) - inarr(7)
+    if (ien > 0) then
+      inarr(8) = inarr(5) - inarr(6) - inarr(7)
+    else
+      inarr(9) = inarr(9) + inarr(6) + inarr(7)
+    end if
 
 ! normalize the averages by the volume of domain
 !
@@ -731,7 +801,7 @@ module statistics
 ! write down the integrals and statistics to appropriate files
 !
     if (master) then
-      write(cunit,"(i9,13(1x,1es18.8e3))") step, time, dt, inarr(1:10), arms
+      write(cunit,"(i9,11(1x,1es18.8e3))") step, time, dt, inarr(1:9)
       write(sunit,"(i9,23(1x,1es18.8e3))") step, time                          &
                                             , avarr(1), mnarr(1), mxarr(1)     &
                                             , avarr(2), mnarr(2), mxarr(2)     &
@@ -744,6 +814,10 @@ module statistics
       call flush_and_sync(cunit)
       call flush_and_sync(sunit)
 
+      if (forcing_enabled) then
+        write(funit,"(4(1x,1es18.8e3))") time, inarr(11:12), arms
+        call flush_and_sync(funit)
+      end if
       if (error_control) then
         write(eunit,efmt) step, time, dth, dte, maxval(errors(:)), errors(:)
         call flush_and_sync(eunit)