diff --git a/src/evolution.F90 b/src/evolution.F90
index f07d0ed..ab6be92 100644
--- a/src/evolution.F90
+++ b/src/evolution.F90
@@ -31,10 +31,22 @@
 !
 module evolution
 
+#ifdef PROFILE
+! import external subroutines
+!
+  use timers, only : set_timer, start_timer, stop_timer
+#endif /* PROFILE */
+
 ! module variables are not implicit by default
 !
   implicit none
 
+#ifdef PROFILE
+! timer indices
+!
+  integer     , save :: imi, ima, imt, imu, imf, imn, imv
+#endif /* PROFILE */
+
 ! pointer to the temporal integration subroutine
 !
   procedure(evolve_euler), pointer, save :: evolve => null()
@@ -116,6 +128,22 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! set timer descriptions
+!
+    call set_timer('evolution:: initialization'  , imi)
+    call set_timer('evolution:: solution advance', ima)
+    call set_timer('evolution:: new time step'   , imt)
+    call set_timer('evolution:: solution update' , imu)
+    call set_timer('evolution:: flux update'     , imf)
+    call set_timer('evolution:: increment update', imn)
+    call set_timer('evolution:: variable update' , imv)
+
+! start accounting time for module initialization/finalization
+!
+    call start_timer(imi)
+#endif /* PROFILE */
+
 ! get the integration method and the value of the CFL coefficient
 !
     call get_parameter_string ("time_advance", integration)
@@ -185,6 +213,12 @@ module evolution
 
     end if
 
+#ifdef PROFILE
+! stop accounting time for module initialization/finalization
+!
+    call stop_timer(imi)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine initialize_evolution
@@ -213,9 +247,23 @@ module evolution
     integer, intent(inout) :: iret
 !
 !-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for module initialization/finalization
+!
+    call start_timer(imi)
+#endif /* PROFILE */
+
+! nullify pointer to integration subroutine
 !
     nullify(evolve)
 
+#ifdef PROFILE
+! stop accounting time for module initialization/finalization
+!
+    call stop_timer(imi)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine finalize_evolution
@@ -256,6 +304,12 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for solution advance
+!
+    call start_timer(ima)
+#endif /* PROFILE */
+
 ! find new time step
 !
     call new_time_step(dtnext)
@@ -296,6 +350,12 @@ module evolution
 
     end if ! toplev > 1
 
+#ifdef PROFILE
+! stop accounting time for solution advance
+!
+    call stop_timer(ima)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine advance
@@ -334,18 +394,17 @@ module evolution
 !
     implicit none
 
-! input variables
+! subroutine arguments
 !
-    real(kind=8), intent(in) :: dtnext
+    real(kind=8), intent(in)  :: dtnext
 
 ! local pointers
 !
-    type(block_data), pointer :: pblock
+    type(block_data), pointer :: pdata
 
 ! local variables
 !
-    integer                   :: iret
-    integer(kind=4)           :: lev
+    integer                   :: iret, mlev
     real(kind=8)              :: cm, dx_min
 
 ! local parameters
@@ -354,53 +413,62 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for new time step estimation
+!
+    call start_timer(imt)
+#endif /* PROFILE */
+
 ! reset the maximum speed, and the highest level
 !
     cmax   = eps
-    lev    = 1
+    mlev   = 1
+
+! assign pdata with the first block on the data block list
+!
+    pdata => list_data
 
 ! iterate over all data blocks in order to find the maximum speed among them
-! and the highest level which is required to obtain the spatial step
+! and the highest level which is required to obtain the minimum spacial step
 !
-    pblock => list_data
-    do while (associated(pblock))
+    do while (associated(pdata))
 
-! find the maximum level occupied by blocks (can be smaller than toplev)
+! update the maximum level
 !
-      lev = max(lev, pblock%meta%level)
+      mlev = max(mlev, pdata%meta%level)
 
-! obtain the maximum speed for the current block
+! get the maximum speed for the current block
 !
-      cm = maxspeed(pblock%q(:,:,:,:))
+      cm = maxspeed(pdata%q(:,:,:,:))
 
-! compare global and local maximum speeds
+! update the maximum speed
 !
       cmax = max(cmax, cm)
 
-! assiociate the pointer with the next block
+! assign pdata to the next block
 !
-      pblock => pblock%next
+      pdata => pdata%next
 
-    end do
+    end do ! over data blocks
 
 #ifdef MPI
-! find maximum speed in the system from all processors
+! reduce maximum speed and level over all processors
 !
     call reduce_maximum_real   (cmax, iret)
-    call reduce_maximum_integer(lev , iret)
+    call reduce_maximum_integer(mlev, iret)
 #endif /* MPI */
 
-! calculate squared cmax
+! calculate the squared cmax
 !
     cmax2 = cmax * cmax
 
-! find the smallest spatial step
+! find the smallest spacial step
 !
 #if NDIMS == 2
-    dx_min = min(adx(lev), ady(lev))
+    dx_min = min(adx(mlev), ady(mlev))
 #endif /* NDIMS == 2 */
 #if NDIMS == 3
-    dx_min = min(adx(lev), ady(lev), adz(lev))
+    dx_min = min(adx(mlev), ady(mlev), adz(mlev))
 #endif /* NDIMS == 3 */
 
 ! calculate the new time step
@@ -408,13 +476,15 @@ module evolution
     dtn = cfl * dx_min / max(cmax                                              &
                         + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps)
 
-! substitute the new time step
-!
-    dt  = dtn
-
 ! round the time
 !
-    if (dtnext > 0.0d+00) dt = min(dt, dtnext)
+    if (dtnext > 0.0d+00) dt = min(dtn, dtnext)
+
+#ifdef PROFILE
+! stop accounting time for new time step estimation
+!
+    call stop_timer(imt)
+#endif /* PROFILE */
 
 !-------------------------------------------------------------------------------
 !
@@ -434,6 +504,12 @@ module evolution
 !   Subroutine advances the solution by one time step using the 1st order
 !   Euler integration method.
 !
+!   References:
+!
+!     [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,
+!         "Numerical Recipes in Fortran",
+!         Cambridge University Press, Cambridge, 1992
+!
 !===============================================================================
 !
   subroutine evolve_euler()
@@ -455,7 +531,7 @@ module evolution
 
 ! local pointers
 !
-    type(block_data), pointer    :: pblock
+    type(block_data), pointer            :: pdata
 
 ! local arrays
 !
@@ -463,40 +539,51 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
-! update fluxes for the first step of the RK2 integration
+#ifdef PROFILE
+! start accounting time for one step update
+!
+    call start_timer(imu)
+#endif /* PROFILE */
+
+! update fluxes
 !
     call update_fluxes()
 
-! update the solution using numerical fluxes stored in the data blocks
+! assign pdata with the first block on the data block list
 !
-    pblock => list_data
-    do while (associated(pblock))
+    pdata => list_data
 
-! calculate variable increment for the current block
+! iterate over all data blocks
 !
-      call update_increment(pblock, du(:,:,:,:))
+    do while (associated(pdata))
 
-! add source terms
+! calculate the variable increment
 !
-      call update_sources(pblock, du(:,:,:,:))
+      call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
 
-! update the solution for the fluid variables
+! add the source terms
 !
-      pblock%u0(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + dt * du(1:nv,:,:,:)
+      call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
+
+! update the solution
+!
+      pdata%u0(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km)            &
+                                     + dt * du(1:nv,1:im,1:jm,1:km)
 
 ! update the conservative variable pointer
 !
-      pblock%u => pblock%u0
+      pdata%u => pdata%u0
 
-! update ψ by its source term
+! update ψ with its source term
 !
-      if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:)
+      if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) =                               &
+                                           decay * pdata%u(ibp,1:im,1:jm,1:km)
 
-! assign pointer to the next block
+! assign pdata to the next block
 !
-      pblock => pblock%next
+      pdata => pdata%next
 
-    end do
+    end do ! over data blocks
 
 ! update primitive variables
 !
@@ -506,6 +593,12 @@ module evolution
 !
     call boundary_variables()
 
+#ifdef PROFILE
+! stop accounting time for one step update
+!
+    call stop_timer(imu)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine evolve_euler
@@ -545,7 +638,7 @@ module evolution
 
 ! local pointers
 !
-    type(block_data), pointer    :: pblock
+    type(block_data), pointer            :: pdata
 
 ! local arrays
 !
@@ -553,36 +646,46 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
-! update fluxes for the first step of the RK2 integration
+#ifdef PROFILE
+! start accounting time for one step update
+!
+    call start_timer(imu)
+#endif /* PROFILE */
+
+! update fluxes
 !
     call update_fluxes()
 
-! update the solution using numerical fluxes stored in the data blocks
+! assign pdata with the first block on the data block list
 !
-    pblock => list_data
-    do while (associated(pblock))
+    pdata => list_data
 
-! calculate variable increment for the current block
+! iterate over all data blocks
 !
-      call update_increment(pblock, du(:,:,:,:))
+    do while (associated(pdata))
 
-! add source terms
+! calculate the variable increment
 !
-      call update_sources(pblock, du(:,:,:,:))
+      call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
 
-! update the solution for the fluid variables
+! add the source terms
 !
-      pblock%u1(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + dt * du(1:nv,:,:,:)
+      call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
+
+! update the intermediate solution
+!
+      pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km)            &
+                                     + dt * du(1:nv,1:im,1:jm,1:km)
 
 ! update the conservative variable pointer
 !
-      pblock%u => pblock%u1
+      pdata%u => pdata%u1
 
-! assign pointer to the next block
+! assign pdata to the next block
 !
-      pblock => pblock%next
+      pdata => pdata%next
 
-    end do
+    end do ! over data blocks
 
 ! update primitive variables
 !
@@ -592,41 +695,46 @@ module evolution
 !
     call boundary_variables()
 
-! update fluxes for the second step of the RK2 integration
+! update fluxes from the intermediate stage
 !
     call update_fluxes()
 
-! update the solution using numerical fluxes stored in the data blocks
+! assign pdata with the first block on the data block list
 !
-    pblock => list_data
-    do while (associated(pblock))
+    pdata => list_data
 
-! calculate variable increment for the current block
+! iterate over all data blocks
 !
-      call update_increment(pblock, du(:,:,:,:))
+    do while (associated(pdata))
 
-! add source terms
+! calculate the variable increment
 !
-      call update_sources(pblock, du(:,:,:,:))
+      call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
 
-! update the solution for the fluid variables
+! add the source terms
 !
-      pblock%u0(1:nv,:,:,:) = 0.5d+00 * (pblock%u0(1:nv,:,:,:)                 &
-                                + pblock%u1(1:nv,:,:,:) + dt * du(1:nv,:,:,:))
+      call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
+
+! update the final solution
+!
+      pdata%u0(1:nv,1:im,1:jm,1:km) = 0.5d+00 * (pdata%u0(1:nv,1:im,1:jm,1:km) &
+                                               + pdata%u1(1:nv,1:im,1:jm,1:km) &
+                                               +  dt * du(1:nv,1:im,1:jm,1:km))
 
 ! update the conservative variable pointer
 !
-      pblock%u => pblock%u0
+      pdata%u => pdata%u0
 
-! update ψ by its source term
+! update ψ with its source term
 !
-      if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:)
+      if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) =                               &
+                                           decay * pdata%u(ibp,1:im,1:jm,1:km)
 
-! assign pointer to the next block
+! assign pdata to the next block
 !
-      pblock => pblock%next
+      pdata => pdata%next
 
-    end do
+    end do ! over data blocks
 
 ! update primitive variables
 !
@@ -636,6 +744,12 @@ module evolution
 !
     call boundary_variables()
 
+#ifdef PROFILE
+! stop accounting time for one step update
+!
+    call stop_timer(imu)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine evolve_rk2
@@ -696,6 +810,12 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for one step update
+!
+    call start_timer(imu)
+#endif /* PROFILE */
+
 ! prepare things which don't change later
 !
     if (first) then
@@ -819,7 +939,7 @@ module evolution
 !
       pdata%u => pdata%u0
 
-! update ψ by its source term
+! update ψ with its source term
 !
       if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) =                               &
                                            decay * pdata%u(ibp,1:im,1:jm,1:km)
@@ -838,6 +958,12 @@ module evolution
 !
     call boundary_variables()
 
+#ifdef PROFILE
+! stop accounting time for one step update
+!
+    call stop_timer(imu)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine evolve_ssprk2
@@ -878,7 +1004,7 @@ module evolution
 
 ! local pointers
 !
-    type(block_data), pointer    :: pblock
+    type(block_data), pointer            :: pdata
 
 ! local variables
 !
@@ -895,42 +1021,52 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for one step update
+!
+    call start_timer(imu)
+#endif /* PROFILE */
+
 !! 1st substep of integration
 !!
-! prepare fractional time step
+! prepare the fractional time step
 !
     ds = dt
 
-! update fluxes for the first step of the RK2 integration
+! update fluxes
 !
     call update_fluxes()
 
-! update the solution using numerical fluxes stored in the data blocks
+! assign pdata with the first block on the data block list
 !
-    pblock => list_data
-    do while (associated(pblock))
+    pdata => list_data
 
-! calculate variable increment for the current block
+! iterate over all data blocks
 !
-      call update_increment(pblock, du(:,:,:,:))
+    do while (associated(pdata))
 
-! add source terms
+! calculate the variable increment
 !
-      call update_sources(pblock, du(:,:,:,:))
+      call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
 
-! update the solution for the fluid variables
+! add the source terms
 !
-      pblock%u1(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
+      call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
+
+! update the first intermediate solution
+!
+      pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km)            &
+                                     + ds * du(1:nv,1:im,1:jm,1:km)
 
 ! update the conservative variable pointer
 !
-      pblock%u => pblock%u1
+      pdata%u => pdata%u1
 
-! assign pointer to the next block
+! assign pdata to the next block
 !
-      pblock => pblock%next
+      pdata => pdata%next
 
-    end do
+    end do ! over data blocks
 
 ! update primitive variables
 !
@@ -942,37 +1078,41 @@ module evolution
 
 !! 2nd substep of integration
 !!
-! prepare fractional time step
+! prepare the fractional time step
 !
     ds = f22 * dt
 
-! update fluxes for the first step of the RK2 integration
+! update fluxes
 !
     call update_fluxes()
 
-! update the solution using numerical fluxes stored in the data blocks
+! assign pdata with the first block on the data block list
 !
-    pblock => list_data
-    do while (associated(pblock))
+    pdata => list_data
 
-! calculate variable increment for the current block
+! iterate over all data blocks
 !
-      call update_increment(pblock, du(:,:,:,:))
+    do while (associated(pdata))
 
-! add source terms
+! calculate the variable increment
 !
-      call update_sources(pblock, du(:,:,:,:))
+      call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
 
-! update the solution for the fluid variables
+! add the source terms
 !
-      pblock%u1(1:nv,:,:,:) = f21 * pblock%u0(1:nv,:,:,:)                      &
-                            + f22 * pblock%u1(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
+      call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
 
-! assign pointer to the next block
+! update the second intermediate solution
 !
-      pblock => pblock%next
+      pdata%u1(1:nv,1:im,1:jm,1:km) = f21 * pdata%u0(1:nv,1:im,1:jm,1:km)      &
+                                    + f22 * pdata%u1(1:nv,1:im,1:jm,1:km)      &
+                                           + ds * du(1:nv,1:im,1:jm,1:km)
 
-    end do
+! assign pdata to the next block
+!
+      pdata => pdata%next
+
+    end do ! over data blocks
 
 ! update primitive variables
 !
@@ -984,45 +1124,50 @@ module evolution
 
 !! 3rd substep of integration
 !!
-! prepare fractional time step
+! prepare the fractional time step
 !
     ds = f32 * dt
 
-! update fluxes for the second step of the RK2 integration
+! update fluxes
 !
     call update_fluxes()
 
-! update the solution using numerical fluxes stored in the data blocks
+! assign pdata with the first block on the data block list
 !
-    pblock => list_data
-    do while (associated(pblock))
+    pdata => list_data
 
-! calculate variable increment for the current block
+! iterate over all data blocks
 !
-      call update_increment(pblock, du(:,:,:,:))
+    do while (associated(pdata))
 
-! add source terms
+! calculate the variable increment
 !
-      call update_sources(pblock, du(:,:,:,:))
+      call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
 
-! update the solution for the fluid variables
+! add the source terms
 !
-      pblock%u0(1:nv,:,:,:) = f31 * pblock%u0(1:nv,:,:,:)                      &
-                            + f32 * pblock%u1(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
+      call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
+
+! update the final solution
+!
+      pdata%u0(1:nv,1:im,1:jm,1:km) = f31 * pdata%u0(1:nv,1:im,1:jm,1:km)      &
+                                    + f32 * pdata%u1(1:nv,1:im,1:jm,1:km)      &
+                                           + ds * du(1:nv,1:im,1:jm,1:km)
 
 ! update the conservative variable pointer
 !
-      pblock%u => pblock%u0
+      pdata%u => pdata%u0
 
-! update ψ by its source term
+! update ψ with its source term
 !
-      if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:)
+      if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) =                               &
+                                           decay * pdata%u(ibp,1:im,1:jm,1:km)
 
-! assign pointer to the next block
+! assign pdata to the next block
 !
-      pblock => pblock%next
+      pdata => pdata%next
 
-    end do
+    end do ! over data blocks
 
 ! update primitive variables
 !
@@ -1032,6 +1177,12 @@ module evolution
 !
     call boundary_variables()
 
+#ifdef PROFILE
+! stop accounting time for one step update
+!
+    call stop_timer(imu)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine evolve_rk3
@@ -1090,6 +1241,12 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for one step update
+!
+    call start_timer(imu)
+#endif /* PROFILE */
+
 != 1st step: U(1) = U(n) + 1/2 dt F[U(n)]
 !
 ! calculate the fractional time step
@@ -1261,7 +1418,7 @@ module evolution
 !
       pdata%u => pdata%u0
 
-! update ψ by its source term
+! update ψ with its source term
 !
       if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) =                               &
                                            decay * pdata%u(ibp,1:im,1:jm,1:km)
@@ -1280,6 +1437,12 @@ module evolution
 !
     call boundary_variables()
 
+#ifdef PROFILE
+! stop accounting time for one step update
+!
+    call stop_timer(imu)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine evolve_ssprk34
@@ -1346,6 +1509,12 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for one step update
+!
+    call start_timer(imu)
+#endif /* PROFILE */
+
 != 1st step: U(1) = U(n) + b1 dt F[U(n)]
 !
 ! calculate the fractional time step
@@ -1564,7 +1733,7 @@ module evolution
                                     + a55 * pdata%u0(1:nv,1:im,1:jm,1:km)      &
                                            + ds * du(1:nv,1:im,1:jm,1:km)
 
-! update ψ by its source term
+! update ψ with its source term
 !
       if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) =                               &
                                            decay * pdata%u(ibp,1:im,1:jm,1:km)
@@ -1583,6 +1752,12 @@ module evolution
 !
     call boundary_variables()
 
+#ifdef PROFILE
+! stop accounting time for one step update
+!
+    call stop_timer(imu)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine evolve_ssprk35
@@ -1610,6 +1785,8 @@ module evolution
 !
     use blocks        , only : block_data, list_data
     use coordinates   , only : adx, ady, adz
+    use coordinates   , only : im, jm, km
+    use equations     , only : nv
 
 ! local variables are not implicit by default
 !
@@ -1617,7 +1794,7 @@ module evolution
 
 ! local pointers
 !
-    type(block_data), pointer  :: pblock
+    type(block_data), pointer  :: pdata
 
 ! local vectors
 !
@@ -1629,34 +1806,50 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for fluxe update
+!
+    call start_timer(imf)
+#endif /* PROFILE */
+
+! assign pdata with the first block on the data block list
+!
+    pdata => list_data
+
 ! iterate over all data blocks
 !
-    pblock => list_data
-    do while (associated(pblock))
+    do while (associated(pdata))
 
 ! obtain dx, dy, and dz for the current block
 !
-      dx(1) = adx(pblock%meta%level)
-      dx(2) = ady(pblock%meta%level)
-      dx(3) = adz(pblock%meta%level)
+      dx(1) = adx(pdata%meta%level)
+      dx(2) = ady(pdata%meta%level)
+      dx(3) = adz(pdata%meta%level)
 
-! update the flux for the current block
+! update fluxes for the current block
 !
       do n = 1, NDIMS
-        call update_flux(n, dx(n), pblock%q(:,:,:,:), pblock%f(n,:,:,:,:))
+        call update_flux(n, dx(n), pdata%q(1:nv,1:im,1:jm,1:km)                &
+                                             , pdata%f(n,1:nv,1:im,1:jm,1:km))
       end do
 
-! assign pointer to the next block
+! assign pdata to the next block
 !
-      pblock => pblock%next
+      pdata => pdata%next
 
-    end do
+    end do ! over data blocks
 
 ! correct the numerical fluxes of the blocks which have neighbours at higher
-! level
+! levels
 !
     call boundary_fluxes()
 
+#ifdef PROFILE
+! stop accounting time for flux update
+!
+    call stop_timer(imf)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine update_fluxes
@@ -1701,6 +1894,12 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for increment update
+!
+    call start_timer(imn)
+#endif /* PROFILE */
+
 ! reset the increment array du
 !
     du(:,:,:,:) = 0.0d+00
@@ -1736,6 +1935,12 @@ module evolution
     end do
 #endif /* NDIMS == 3 */
 
+#ifdef PROFILE
+! stop accounting time for increment update
+!
+    call stop_timer(imn)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine update_increment
@@ -1774,6 +1979,12 @@ module evolution
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for variable update
+!
+    call start_timer(imv)
+#endif /* PROFILE */
+
 ! associate the pointer with the first block on the data block list
 !
     pdata => list_data
@@ -1800,6 +2011,12 @@ module evolution
 
     end do
 
+#ifdef PROFILE
+! stop accounting time for variable update
+!
+    call stop_timer(imv)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine update_variables
diff --git a/src/makefile b/src/makefile
index d81ae3b..5bf97dd 100644
--- a/src/makefile
+++ b/src/makefile
@@ -148,7 +148,7 @@ mesh.o           : mesh.F90 blocks.o coordinates.o domains.o equations.o       \
                    error.o interpolations.o mpitools.o problems.o refinement.o \
                    timers.o
 mpitools.o       : mpitools.F90 timers.o
-operators.o      : operators.F90
+operators.o      : operators.F90 timers.o
 parameters.o     : parameters.F90 mpitools.o
 problems.o       : problems.F90 blocks.o constants.o coordinates.o equations.o \
                    error.o operators.o parameters.o random.o timers.o
diff --git a/src/mesh.F90 b/src/mesh.F90
index 2c3524f..6eb6758 100644
--- a/src/mesh.F90
+++ b/src/mesh.F90
@@ -878,7 +878,7 @@ module mesh
 
 ! local buffer for data block exchange
 !
-    real(kind=8)   , dimension(2,nv,im,jm,km) :: rbuf
+    real(kind=8)   , dimension(3,nv,im,jm,km) :: rbuf
 #endif /* MPI */
 
 !-------------------------------------------------------------------------------
@@ -931,8 +931,9 @@ module mesh
 
 ! copy data to buffer
 !
-              rbuf(1,:,:,:,:) = pmeta%data%u(:,:,:,:)
-              rbuf(2,:,:,:,:) = pmeta%data%q(:,:,:,:)
+              rbuf(1,:,:,:,:) = pmeta%data%q (:,:,:,:)
+              rbuf(2,:,:,:,:) = pmeta%data%u0(:,:,:,:)
+              rbuf(3,:,:,:,:) = pmeta%data%u1(:,:,:,:)
 
 ! send data
 !
@@ -961,8 +962,9 @@ module mesh
 
 ! coppy the buffer to data block
 !
-              pmeta%data%u(:,:,:,:) = rbuf(1,:,:,:,:)
-              pmeta%data%q(:,:,:,:) = rbuf(2,:,:,:,:)
+              pmeta%data%q (:,:,:,:) = rbuf(1,:,:,:,:)
+              pmeta%data%u0(:,:,:,:) = rbuf(2,:,:,:,:)
+              pmeta%data%u1(:,:,:,:) = rbuf(3,:,:,:,:)
 
             end if ! nproc == n
 
@@ -1031,11 +1033,11 @@ module mesh
 !
 !   Arguments:
 !
-!     pblock - the input meta block;
+!     pmeta - the input meta block;
 !
 !===============================================================================
 !
-  subroutine prolong_block(pblock)
+  subroutine prolong_block(pmeta)
 
 ! import external procedures and variables
 !
@@ -1051,7 +1053,7 @@ module mesh
 
 ! input arguments
 !
-    type(block_meta), pointer, intent(inout) :: pblock
+    type(block_meta), pointer, intent(inout) :: pmeta
 
 ! local variables
 !
@@ -1083,7 +1085,7 @@ module mesh
 
 ! assign the pdata pointer
 !
-    pdata => pblock%data
+    pdata => pmeta%data
 
 ! prepare dimensions
 !
@@ -1175,7 +1177,7 @@ module mesh
 
 ! assign pointer to the current child
 !
-      pchild => pblock%child(p)%ptr
+      pchild => pmeta%child(p)%ptr
 
 ! obtain the position of child in the parent block
 !
@@ -1236,11 +1238,11 @@ module mesh
 !
 !   Arguments:
 !
-!     pblock - the input meta block;
+!     pmeta - the input meta block;
 !
 !===============================================================================
 !
-  subroutine restrict_block(pblock)
+  subroutine restrict_block(pmeta)
 
 ! import external procedures and variables
 !
@@ -1256,7 +1258,7 @@ module mesh
 
 ! subroutine arguments
 !
-    type(block_meta), pointer, intent(inout) :: pblock
+    type(block_meta), pointer, intent(inout) :: pmeta
 
 ! local variables
 !
@@ -1282,7 +1284,7 @@ module mesh
 
 ! assign the parent data pointer
 !
-    pparent => pblock%data
+    pparent => pmeta%data
 
 ! iterate over all children
 !
@@ -1290,7 +1292,7 @@ module mesh
 
 ! assign a pointer to the current child
 !
-      pchild  => pblock%child(p)%ptr%data
+      pchild  => pmeta%child(p)%ptr%data
 
 ! obtain the child position in the parent block
 !
diff --git a/src/operators.F90 b/src/operators.F90
index b94b4ef..917e0a0 100644
--- a/src/operators.F90
+++ b/src/operators.F90
@@ -30,10 +30,22 @@
 !
 module operators
 
+#ifdef PROFILE
+! import external subroutines
+!
+  use timers, only : set_timer, start_timer, stop_timer
+#endif /* PROFILE */
+
 ! module variables are not implicit by default
 !
   implicit none
 
+#ifdef PROFILE
+! timer indices
+!
+  integer     , save :: imi, imd, img, imc, iml, im1, im2
+#endif /* PROFILE */
+
 ! by default everything is public
 !
   private
@@ -41,7 +53,7 @@ module operators
 ! declare public subroutines
 !
   public :: initialize_operators, finalize_operators
-  public :: divergence, curl, laplace
+  public :: divergence, gradient, curl, laplace
 
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 !
@@ -80,6 +92,27 @@ module operators
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! set timer descriptions
+!
+    call set_timer('operators:: initialization', imi)
+    call set_timer('operators:: divergence'    , imd)
+    call set_timer('operators:: gradient'      , img)
+    call set_timer('operators:: curl'          , imc)
+    call set_timer('operators:: laplace'       , iml)
+    call set_timer('operators:: 1st derivative', im1)
+    call set_timer('operators:: 2nd derivative', im2)
+
+! start accounting time for the module initialization/finalization
+!
+    call start_timer(imi)
+#endif /* PROFILE */
+
+#ifdef PROFILE
+! stop accounting time for the module initialization/finalization
+!
+    call stop_timer(imi)
+#endif /* PROFILE */
 
 !-------------------------------------------------------------------------------
 !
@@ -110,6 +143,17 @@ module operators
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for the module initialization/finalization
+!
+    call start_timer(imi)
+#endif /* PROFILE */
+
+#ifdef PROFILE
+! stop accounting time for the module initialization/finalization
+!
+    call stop_timer(imi)
+#endif /* PROFILE */
 
 !-------------------------------------------------------------------------------
 !
@@ -155,6 +199,12 @@ module operators
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for the divergence operator calculation
+!
+    call start_timer(imd)
+#endif /* PROFILE */
+
 ! allocate temporary array
 !
     allocate(w(size(u,2), size(u,3), size(u,4)))
@@ -167,7 +217,7 @@ module operators
 !
     do dir = 1, NDIMS
 
-! calculate contribution from the Y derivative of the Y component
+! calculate derivative along the current direction
 !
       call derivative_1st(dir, dh(dir), u(dir,:,:,:), w(:,:,:))
 
@@ -181,12 +231,83 @@ module operators
 !
     deallocate(w)
 
+#ifdef PROFILE
+! stop accounting time for the divergence operator calculation
+!
+    call stop_timer(imd)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine divergence
 !
 !===============================================================================
 !
+! subroutine GRADIENT:
+! -------------------
+!
+!   Subroutine calculates the cell centered gradient of the input scalar field.
+!
+!      grad(U) = ∇ U = [ ∂x U, ∂y U, ∂z U ]
+!
+!   Arguments:
+!
+!     dh  - the spacial intervals in all direction;
+!     u   - the input scalar field;
+!     v   - the output gradient vector field;
+!
+!===============================================================================
+!
+  subroutine gradient(dh, u, v)
+
+! local variables are not implicit by default
+!
+    implicit none
+
+! input and output variables
+!
+    real(kind=8), dimension(3)      , intent(in)  :: dh
+    real(kind=8), dimension(:,:,:)  , intent(in)  :: u
+    real(kind=8), dimension(:,:,:,:), intent(out) :: v
+
+! local variables
+!
+    integer :: dir
+!
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for the gradient operator calculation
+!
+    call start_timer(img)
+#endif /* PROFILE */
+
+! reset the output array
+!
+    v(:,:,:,:) = 0.0d+00
+
+! iterate over directions and calculate gradient components
+!
+    do dir = 1, NDIMS
+
+! calculate derivative along the current direction
+!
+      call derivative_1st(dir, dh(dir), u(:,:,:), v(dir,:,:,:))
+
+    end do ! directions
+
+#ifdef PROFILE
+! stop accounting time for the gradient operator calculation
+!
+    call stop_timer(img)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine gradient
+!
+!===============================================================================
+!
 ! subroutine CURL:
 ! ---------------
 !
@@ -221,6 +342,12 @@ module operators
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for the rotation operator calculation
+!
+    call start_timer(imc)
+#endif /* PROFILE */
+
 ! allocate temporary array
 !
     allocate(w(size(u,2), size(u,3), size(u,4)))
@@ -295,6 +422,12 @@ module operators
 !
     deallocate(w)
 
+#ifdef PROFILE
+! stop accounting time for the rotation operator calculation
+!
+    call stop_timer(imc)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine curl
@@ -334,6 +467,12 @@ module operators
 !
 !-------------------------------------------------------------------------------
 !
+#ifdef PROFILE
+! start accounting time for the laplace operator calculation
+!
+    call start_timer(iml)
+#endif /* PROFILE */
+
 ! allocate temporary array
 !
     allocate(w(size(u,1), size(u,2), size(u,3)))
@@ -368,6 +507,12 @@ module operators
 !
     deallocate(w)
 
+#ifdef PROFILE
+! stop accounting time for the laplace operator calculation
+!
+    call stop_timer(iml)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine laplace
@@ -420,6 +565,12 @@ module operators
     if (dir < 1 .or. dir > NDIMS .or. dh == 0.0d+00) return
 #endif /* DEBUG */
 
+#ifdef PROFILE
+! start accounting time for the 1st derivative calculation
+!
+    call start_timer(im1)
+#endif /* PROFILE */
+
 ! prepare index limits
 !
     m0 = size(u, dir)
@@ -458,6 +609,12 @@ module operators
 
     end select
 
+#ifdef PROFILE
+! stop accounting time for the 1st derivative calculation
+!
+    call stop_timer(im1)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine derivative_1st
@@ -504,6 +661,12 @@ module operators
     if (dir < 1 .or. dir > NDIMS .or. dh == 0.0d+00) return
 #endif /* DEBUG */
 
+#ifdef PROFILE
+! start accounting time for the 2nd derivative calculation
+!
+    call start_timer(im2)
+#endif /* PROFILE */
+
 ! prepare index limits
 !
     m0 = size(u, dir)
@@ -542,6 +705,12 @@ module operators
 
     end select
 
+#ifdef PROFILE
+! stop accounting time for the 2nd derivative calculation
+!
+    call stop_timer(im2)
+#endif /* PROFILE */
+
 !-------------------------------------------------------------------------------
 !
   end subroutine derivative_2nd
diff --git a/src/sources.F90 b/src/sources.F90
index 753888f..bcd9b7e 100644
--- a/src/sources.F90
+++ b/src/sources.F90
@@ -45,6 +45,10 @@ module sources
   integer, save :: imi, imu
 #endif /* PROFILE */
 
+! GLM-MHD source terms type (1 - EGLM, 2 - HEGLM)
+!
+  integer     , save :: glm_type    = 0
+
 ! gravitational acceleration coefficient
 !
   real(kind=8), save :: gpoint      = 0.0d+00
@@ -95,7 +99,7 @@ module sources
 
 ! include external procedures and variables
 !
-    use parameters , only : get_parameter_real
+    use parameters , only : get_parameter_string, get_parameter_real
 
 ! local variables are not implicit by default
 !
@@ -105,6 +109,10 @@ module sources
 !
     logical, intent(in)    :: verbose
     integer, intent(inout) :: iret
+
+! local variables
+!
+    character(len=8)       :: tglm = "none"
 !
 !-------------------------------------------------------------------------------
 !
@@ -119,6 +127,23 @@ module sources
     call start_timer(imi)
 #endif /* PROFILE */
 
+! get the type of the GLM source terms
+!
+    call get_parameter_string("glm_source_terms", tglm)
+
+! set the glm_type variable to correct value
+!
+    select case(trim(tglm))
+    case("eglm", "EGLM")
+      glm_type = 1
+
+    case("heglm", "HEGLM")
+      glm_type = 2
+
+    case default
+      glm_type = 0
+    end select
+
 ! get acceleration coefficient
 !
       call get_parameter_real("gpoint"     , gpoint   )
@@ -135,6 +160,7 @@ module sources
 !
     if (verbose) then
 
+      write (*,"(4x,a,1x,a)    ") "glm source terms       =", trim(tglm)
       write (*,"(4x,a,1x,1e9.2)") "point mass constant    =", gpoint
       write (*,"(4x,a,1x,1e9.2)") "viscosity              =", viscosity
       write (*,"(4x,a,1x,1e9.2)") "resistivity            =", resistivity
@@ -215,8 +241,8 @@ module sources
     use coordinates    , only : ax, ay, az, adx, ady, adz, adxi, adyi, adzi
     use equations      , only : nv, inx, iny, inz
     use equations      , only : idn, ivx, ivy, ivz, imx, imy, imz, ien
-    use equations      , only : ibx, iby, ibz
-    use operators      , only : laplace, curl
+    use equations      , only : ibx, iby, ibz, ibp
+    use operators      , only : divergence, gradient, laplace, curl
 
 ! local variables are not implicit by default
 !
@@ -241,6 +267,7 @@ module sources
     real(kind=8), dimension(im) :: x
     real(kind=8), dimension(jm) :: y
     real(kind=8), dimension(km) :: z
+    real(kind=8), dimension(im,jm,km)   :: db
     real(kind=8), dimension(3,im,jm,km) :: jc
 !
 !-------------------------------------------------------------------------------
@@ -440,9 +467,9 @@ module sources
 
     end if ! viscosity is not zero
 
-! proceed only if the resistivity coefficient is not zero
+!=== add magnetic field related source terms ===
 !
-    if (resistivity > 0.0d+00 .and. ibx > 0) then
+    if (ibx > 0) then
 
 ! prepare coordinate increments
 !
@@ -450,43 +477,111 @@ module sources
       dh(2) = ady(pdata%meta%level)
       dh(3) = adz(pdata%meta%level)
 
+! add the EGLM-MHD source terms
+!
+      if (glm_type > 0) then
+
+! calculate the magnetic field divergence
+!
+        call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:))
+
+! update the momentum component increments, i.e.
+!     d/dt (ρv) + ∇.F = - (∇.B)B
+!
+        du(imx,:,:,:) = du(imx,:,:,:) - db(:,:,:) * pdata%q(ibx,:,:,:)
+        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,:,:,:), jc(inx:inz,:,:,:))
+
+! 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,:,:,:) * jc(inx:inz,:,:,:), 1)
+
+        end if ! ien > 0
+
+! add the HEGLM-MHD source terms
+!
+        if (glm_type > 1) then
+
+! update magnetic field component increments, i.e.
+!     d/dt B + ∇.F = - (∇.B)v
+!
+          du(ibx,:,:,:) = du(ibx,:,:,:) - db(:,:,:) * pdata%q(ivx,:,:,:)
+          du(iby,:,:,:) = du(iby,:,:,:) - db(:,:,:) * pdata%q(ivy,:,:,:)
+          du(ibz,:,:,:) = du(ibz,:,:,:) - db(:,:,:) * pdata%q(ivz,:,:,:)
+
+! update the energy equation
+!
+          if (ien > 0) then
+
+! calculate scalar product of velocity and magnetic field
+!
+            jc(inx,:,:,:) = sum(pdata%q(ivx:ivz,:,:,:)                         &
+                                                  * pdata%q(ibx:ibz,:,:,:), 1)
+
+! add the divergence potential source term to the energy equation, i.e.
+!     d/dt E + ∇.F = - (∇.B) (v.B)
+!
+            du(ien,:,:,:) = du(ien,:,:,:) - db(:,:,:) * jc(inx,:,:,:)
+
+          end if ! ien > 0
+
+        end if ! glm_type > 1
+
+      end if ! glmtype > 0
+
+! proceed only if the resistivity coefficient is not zero
+!
+      if (resistivity > 0.0d+00) then
+
 ! calculate the Laplace operator of B, i.e. Δ(B)
 !
-      call laplace(dh(:), pdata%q(ibx,:,:,:), jc(inx,:,:,:))
-      call laplace(dh(:), pdata%q(iby,:,:,:), jc(iny,:,:,:))
-      call laplace(dh(:), pdata%q(ibz,:,:,:), jc(inz,:,:,:))
+        call laplace(dh(:), pdata%q(ibx,:,:,:), jc(inx,:,:,:))
+        call laplace(dh(:), pdata%q(iby,:,:,:), jc(iny,:,:,:))
+        call laplace(dh(:), pdata%q(ibz,:,:,:), jc(inz,:,:,:))
 
 ! update magnetic field component increments
 !
-      du(ibx,:,:,:) = du(ibx,:,:,:) + resistivity * jc(inx,:,:,:)
-      du(iby,:,:,:) = du(iby,:,:,:) + resistivity * jc(iny,:,:,:)
-      du(ibz,:,:,:) = du(ibz,:,:,:) + resistivity * jc(inz,:,:,:)
+        du(ibx,:,:,:) = du(ibx,:,:,:) + resistivity * jc(inx,:,:,:)
+        du(iby,:,:,:) = du(iby,:,:,:) + resistivity * jc(iny,:,:,:)
+        du(ibz,:,:,:) = du(ibz,:,:,:) + resistivity * jc(inz,:,:,:)
 
 ! update energy equation
 !
-      if (ien > 0) then
+        if (ien > 0) then
 
 ! add the first resistive source term to the energy equation, i.e.
 !     d/dt E + ∇.F = η B.[Δ(B)]
 !
-        du(ien,:,:,:) = du(ien,:,:,:)                                          &
+          du(ien,:,:,:) = du(ien,:,:,:)                                        &
                           + resistivity * (pdata%q(ibx,:,:,:) * jc(inx,:,:,:)  &
                                          + pdata%q(iby,:,:,:) * jc(iny,:,:,:)  &
                                          + pdata%q(ibz,:,:,:) * jc(inz,:,:,:))
 
 ! calculate current density J = ∇xB
 !
-        call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(inx:inz,:,:,:))
+          call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(inx:inz,:,:,:))
 
 ! add the second resistive source term to the energy equation, i.e.
 !     d/dt E + ∇.F = η J²
 !
-        du(ien,:,:,:) = du(ien,:,:,:)                                          &
+          du(ien,:,:,:) = du(ien,:,:,:)                                        &
                              + resistivity * sum(jc(:,:,:,:) * jc(:,:,:,:), 1)
 
-      end if ! energy equation present
+        end if ! energy equation present
 
-    end if ! resistivity is not zero
+      end if ! resistivity is not zero
+
+    end if ! ibx > 0
 
 #ifdef PROFILE
 ! stop accounting time for source terms