diff --git a/sources/boundaries.F90 b/sources/boundaries.F90
index a78abac..9a352d1 100644
--- a/sources/boundaries.F90
+++ b/sources/boundaries.F90
@@ -5333,8 +5333,8 @@ module boundaries
     use equations      , only : csnd2
     use gravity        , only : gravitational_acceleration
     use iso_fortran_env, only : error_unit
-    use user_problem   , only : boundary_user_x, boundary_user_y               &
-                              , boundary_user_z
+    use user_problem   , only : user_boundary_x, user_boundary_y               &
+                              , user_boundary_z
 
 ! local variables are not implicit by default
 !
@@ -5535,7 +5535,7 @@ module boundaries
 !
       case(bnd_user)
 
-        call boundary_user_x(side(1), jl, ju, kl, ku                           &
+        call user_boundary_x(side(1), jl, ju, kl, ku                           &
                            , t, dt, x(:), y(:), z(:), qn(:,:,:,:))
 
 ! wrong boundary conditions
@@ -5713,7 +5713,7 @@ module boundaries
 !
       case(bnd_user)
 
-        call boundary_user_y(side(2), il, iu, kl, ku                           &
+        call user_boundary_y(side(2), il, iu, kl, ku                           &
                            , t, dt, x(:), y(:), z(:), qn(:,:,:,:))
 
 ! wrong boundary conditions
@@ -5887,7 +5887,7 @@ module boundaries
 !
       case(bnd_user)
 
-        call boundary_user_z(side(3), il, iu, jl, ju                           &
+        call user_boundary_z(side(3), il, iu, jl, ju                           &
                            , t, dt, x(:), y(:), z(:), qn(:,:,:,:))
 
 ! wrong boundary conditions
diff --git a/sources/driver.F90 b/sources/driver.F90
index 8be0daa..974b68e 100644
--- a/sources/driver.F90
+++ b/sources/driver.F90
@@ -86,6 +86,7 @@ program amun
   use timers         , only : start_timer, stop_timer, set_timer, get_timer
   use timers         , only : get_timer_total, timer_enabled, timer_description
   use timers         , only : get_count, ntimers
+  use user_problem   , only : user_time_statistics
 
 ! module variables are not implicit by default
 !
@@ -564,6 +565,7 @@ program amun
   if (.not. restart_from_snapshot()) then
 
     call store_integrals()
+    call user_time_statistics()
     call write_snapshot()
 
   end if
@@ -655,9 +657,10 @@ program amun
 !
     call store_mesh_stats(step, time)
 
-! store integrals
+! store time statistics
 !
     call store_integrals()
+    call user_time_statistics()
 
 ! write down the restart snapshot
 !
diff --git a/sources/evolution.F90 b/sources/evolution.F90
index 86abe1e..93e2f8c 100644
--- a/sources/evolution.F90
+++ b/sources/evolution.F90
@@ -350,12 +350,13 @@ module evolution
       cfl       = 2.5d+00 * cfl
       name_int  = "Optimized 3ʳᵈ-order 5-step embedded RK3(2)5"
 
-      allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
+      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
       c (:)   = [ 0.000000000000000000000000000000000000d+00,                  &
                   2.300285062878154351930669430512780706d-01,                  &
                   4.050049049262914975700372321130661410d-01,                  &
                   8.947823877926760224705450466361360720d-01,                  &
-                  7.235108137218888081489570284485201518d-01 ]
+                  7.235108137218888081489570284485201518d-01,                  &
+                  1.000000000000000000000000000000000000d+00 ]
       dl(:)   = [ 1.000000000000000000000000000000000000d+00,                  &
                   3.407687209321455242558804921815861422d-01,                  &
                   3.414399280584625023244387687873774697d-01,                  &
@@ -415,12 +416,13 @@ module evolution
       cfl       = 2.5d+00 * cfl
       name_int  = "Optimized 3ʳᵈ-order 5-step embedded RK3(2)5 FSAL"
 
-      allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
+      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
       c (:)   = [ 0.000000000000000000000000000000000000d+00,                  &
                   2.300298624518076223899418286314123354d-01,                  &
                   4.050046072094990912268498160116125481d-01,                  &
                   8.947822893693433545220710894560512805d-01,                  &
-                  7.235136928826589010272834603680114769d-01 ]
+                  7.235136928826589010272834603680114769d-01,                  &
+                  1.000000000000000000000000000000000000d+00 ]
       dl(:)   = [ 1.000000000000000000000000000000000000d+00,                  &
                   3.407655879334525365094815965895763636d-01,                  &
                   3.414382655003386206551709871126405331d-01,                  &
@@ -479,7 +481,7 @@ module evolution
       cfl       = 3.5d+00 * cfl
       name_int  = "Optimized 4ᵗʰ-order 9-step embedded RK4(3)9"
 
-      allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
+      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
       c (:)   = [  0.000000000000000000000000000000000000d+00,                 &
                    2.836343531977826022543660465926414772d-01,                 &
                    5.484073767552486705240014599676811834d-01,                 &
@@ -488,7 +490,8 @@ module evolution
                    3.518526451892056368706593492732753284d-01,                 &
                    1.665941920204672094647868254892387293d+00,                 &
                    9.715276989307335935187466054546761665d-01,                 &
-                   9.051569554420045339601721625247585643d-01 ]
+                   9.051569554420045339601721625247585643d-01,                 &
+                   1.000000000000000000000000000000000000d+00 ]
       dl(:)   = [  1.000000000000000000000000000000000000d+00,                 &
                    1.262923854387806460989545005598562667d+00,                 &
                    7.574967177560872438940839460448329992d-01,                 &
@@ -572,7 +575,7 @@ module evolution
       cfl       = 3.5d+00 * cfl
       name_int  = "Optimized 4ᵗʰ-order 9-step embedded RK4(3)9 FSAL"
 
-      allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
+      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
       c (:)   = [  0.000000000000000000000000000000000000d+00,                 &
                    2.836343005184365275160654678626695428d-01,                 &
                    5.484076570002894365286665352032296535d-01,                 &
@@ -581,7 +584,8 @@ module evolution
                    3.518526124230705801739919476290327750d-01,                 &
                    1.665941994879593315477304663913129942d+00,                 &
                    9.715279295934715835299192116436237065d-01,                 &
-                   9.051569840159589594903399929316959062d-01 ]
+                   9.051569840159589594903399929316959062d-01,                 &
+                   1.000000000000000000000000000000000000d+00 ]
       dl(:)   = [  1.000000000000000000000000000000000000d+00,                 &
                    1.262923876648114432874834923838556100d+00,                 &
                    7.574967189685911558308119415539596711d-01,                 &
@@ -664,7 +668,7 @@ module evolution
       cfl       = 4.2d+00 * cfl
       name_int  = "Optimized 5ᵗʰ-order 10-step embedded RK5(4)10"
 
-      allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
+      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
       c (:)   = [  0.000000000000000000000000000000000000d+00,                 &
                    2.597883575710995826783320802193635406d-01,                 &
                    9.904573115730917688557891428202061598d-02,                 &
@@ -674,7 +678,8 @@ module evolution
                    5.449986973408778242805929551952000165d-01,                 &
                    7.615224662599497796472095353126697300d-01,                 &
                    8.427062083059167761623893618875787414d-01,                 &
-                   9.152209807185253394871325258038753352d-01 ]
+                   9.152209807185253394871325258038753352d-01,                 &
+                   1.000000000000000000000000000000000000d+00 ]
       dl(:)   = [  1.000000000000000000000000000000000000d+00,                 &
                   -1.331778409133849616712007380176762548d-01,                 &
                    8.260422785246030254485064732649153253d-01,                 &
@@ -764,7 +769,7 @@ module evolution
       cfl       = 4.2d+00 * cfl
       name_int  = "Optimized 5ᵗʰ-order 10-step embedded RK5(4)10 FSAL"
 
-      allocate(c(stages), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
+      allocate(c(stages+1), dl(stages), bt(stages), bh(stages+1), gm(3,stages))
       c (:)   = [  0.000000000000000000000000000000000000d+00,                 &
                    2.597883554788674084039539165398464630d-01,                 &
                    9.904573247592460887087003212056568980d-02,                 &
@@ -774,7 +779,8 @@ module evolution
                    5.449986978853637084972622392134732553d-01,                 &
                    7.615224694532590139829150720490417596d-01,                 &
                    8.427062083267360939805493320684741215d-01,                 &
-                   9.152209805057669959657927210873423883d-01 ]
+                   9.152209805057669959657927210873423883d-01,                 &
+                   1.000000000000000000000000000000000000d+00 ]
       dl(:)   = [  1.000000000000000000000000000000000000d+00,                 &
                   -1.331778419508803397033287009506932673d-01,                 &
                    8.260422814750207498262063505871077303d-01,                 &
@@ -3217,14 +3223,15 @@ module evolution
         pdata => pdata%next
       end do
 
-      call update_variables(tm, dtm)
+      dtm = c(2) * dt
+      call update_variables(time, dtm)
 
 != remaining stages
 !
       do i = 2, stages
 
-          tm  = time + c(i) * dt
-          dtm = (c(i) - c(i-1)) * dt
+          dtm = c(i) * dt
+          tm  = time + dtm
 
           pdata => list_data
           do while (associated(pdata))
@@ -3252,6 +3259,7 @@ module evolution
           pdata => pdata%next
         end do
 
+        dtm = c(i+1) * dt
         call update_variables(tm, dtm)
 
       end do ! i = 2, stages
@@ -3348,9 +3356,6 @@ module evolution
 
 != final step: U(n+1) = U(1)
 !
-    dtm = dt
-    tm  = time + dt
-
     pdata => list_data
     do while (associated(pdata))
 
@@ -3361,7 +3366,8 @@ module evolution
       pdata => pdata%next
     end do
 
-    call update_variables(tm, dtm)
+    tm = time + dt
+    call update_variables(tm, dt)
 
 #ifdef PROFILE
     call stop_timer(imu)
diff --git a/sources/gravity.F90 b/sources/gravity.F90
index 205b648..5ad75b6 100644
--- a/sources/gravity.F90
+++ b/sources/gravity.F90
@@ -103,8 +103,8 @@ module gravity
 ! include external procedures and variables
 !
     use parameters     , only : get_parameter
-    use user_problem   , only : gravitational_acceleration_user                &
-                              , gravity_enabled_user
+    use user_problem   , only : user_gravitational_acceleration                &
+                              , user_gravity_enabled
 
 ! local variables are not implicit by default
 !
@@ -152,8 +152,8 @@ module gravity
 
 ! in case of other problems, gravity is calculated by user
 !
-      gravitational_acceleration => gravitational_acceleration_user
-      gravity_enabled            =  gravity_enabled_user
+      gravitational_acceleration => user_gravitational_acceleration
+      gravity_enabled            =  user_gravity_enabled
 
     end select
 
diff --git a/sources/problems.F90 b/sources/problems.F90
index 8a4cb3d..be80efa 100644
--- a/sources/problems.F90
+++ b/sources/problems.F90
@@ -102,7 +102,7 @@ module problems
   subroutine initialize_problems(problem, verbose, status)
 
     use parameters  , only : get_parameter
-    use user_problem, only : initialize_user_problem, setup_problem_user
+    use user_problem, only : initialize_user_problem, setup_user_problem
 
     implicit none
 
@@ -150,7 +150,7 @@ module problems
       setup_problem => setup_problem_tearing
 
     case default
-      setup_problem => setup_problem_user
+      setup_problem => setup_user_problem
 
       call initialize_user_problem(problem, verbose, status)
 
diff --git a/sources/shapes.F90 b/sources/shapes.F90
index 72372a0..65bba40 100644
--- a/sources/shapes.F90
+++ b/sources/shapes.F90
@@ -102,7 +102,7 @@ module shapes
 ! include external procedures and variables
 !
     use parameters  , only : get_parameter
-    use user_problem, only : update_shapes_user
+    use user_problem, only : update_user_shapes
 
 ! local variables are not implicit by default
 !
@@ -160,7 +160,7 @@ module shapes
 ! no shape update
 !
       case default
-        update_shapes => update_shapes_user
+        update_shapes => update_user_shapes
 
       end select
 
diff --git a/sources/sources.F90 b/sources/sources.F90
index 523d9e9..55a22cc 100644
--- a/sources/sources.F90
+++ b/sources/sources.F90
@@ -330,7 +330,7 @@ module sources
     use equations      , only : ibx, iby, ibz, ibp
     use gravity        , only : gravity_enabled, gravitational_acceleration
     use operators      , only : divergence, gradient, laplace, curl
-    use user_problem   , only : update_sources_user
+    use user_problem   , only : update_user_sources
 
 ! local variables are not implicit by default
 !
@@ -750,7 +750,7 @@ module sources
 
 ! add user defined source terms
 !
-    call update_sources_user(pdata, t, dt, du(:,:,:,:))
+    call update_user_sources(pdata, t, dt, du(:,:,:,:))
 
 #ifdef PROFILE
 ! stop accounting time for source terms
diff --git a/sources/user_problem.F90 b/sources/user_problem.F90
index 8dc3c3d..fe15cf2 100644
--- a/sources/user_problem.F90
+++ b/sources/user_problem.F90
@@ -76,12 +76,21 @@ module user_problem
 
 ! flag indicating if the gravitational source term is enabled
 !
-  logical, save :: gravity_enabled_user = .false.
+  logical, save :: user_gravity_enabled = .false.
 
 ! allocatable arrays for velocity perturbation
 !
   real(kind=8), dimension(:), allocatable :: kx, ky, kz, ux, uy, uz, ph
 
+! export subroutines
+!
+  private
+  public :: initialize_user_problem, finalize_user_problem
+  public :: setup_user_problem, update_user_sources, update_user_shapes
+  public :: user_gravitational_acceleration, user_gravity_enabled
+  public :: user_boundary_x, user_boundary_y, user_boundary_z
+  public :: user_time_statistics
+
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 !
   contains
@@ -134,7 +143,7 @@ module user_problem
 #if NDIMS == 3
     real(kind=8)      :: thv, tx, ty, tz, tt
 #endif /* NDIMS == 3 */
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
@@ -411,7 +420,7 @@ module user_problem
 ! subroutine arguments
 !
     integer, intent(out) :: status
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
@@ -442,7 +451,7 @@ module user_problem
 !
 !===============================================================================
 !
-! subroutine SETUP_PROBLEM_USER:
+! subroutine SETUP_USER_PROBLEM:
 ! -----------------------------
 !
 !   Subroutine sets the initial conditions for the user specific problem.
@@ -454,7 +463,7 @@ module user_problem
 !
 !===============================================================================
 !
-  subroutine setup_problem_user(pdata)
+  subroutine setup_user_problem(pdata)
 
 ! include external procedures and variables
 !
@@ -512,7 +521,7 @@ module user_problem
     real(kind=8), dimension(nn)          :: zc
 #endif /* NDIMS == 3 */
     real(kind=8), dimension(3)           :: dh
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
@@ -878,104 +887,11 @@ module user_problem
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine setup_problem_user
+  end subroutine setup_user_problem
 !
 !===============================================================================
 !
-! subroutine LOG_COSH:
-! -------------------
-!
-!   Function calculates the logarithm of the hyperbolic cosine, which is
-!   the result of the integration of tanh(x). Direct calculation using
-!   Fortran intrinsic subroutines fails for large values of x, therefore
-!   the logarithm of cosh is approximated as |x| + log(1/2) for
-!   |x| > threshold.
-!
-!   Arguments:
-!
-!     x - function argument;
-!
-!===============================================================================
-!
-  function log_cosh(x) result(y)
-
-! local variables are not implicit by default
-!
-    implicit none
-
-! function arguments
-!
-    real(kind=8), intent(in) :: x
-    real(kind=8)             :: y
-
-! local parameters
-!
-    real(kind=8), parameter :: th = acosh(huge(x)), lh = log(0.5d+00)
-!
-!-------------------------------------------------------------------------------
-!
-    if (abs(x) < th) then
-      y = log(cosh(x))
-    else
-      y = abs(x) + lh
-    end if
-
-!-------------------------------------------------------------------------------
-!
-  end function log_cosh
-!
-!===============================================================================
-!
-! subroutine UPDATE_SHAPES_USER:
-! -----------------------------
-!
-!   Subroutine defines the regions updated by user.
-!
-!   Arguments:
-!
-!     pdata - pointer to the data block structure of the currently initialized
-!             block;
-!     time  - time at the moment of update;
-!     dt    - time step since the last update;
-!
-!===============================================================================
-!
-  subroutine update_shapes_user(pdata, time, dt)
-
-! include external procedures and variables
-!
-    use blocks, only : block_data
-
-! local variables are not implicit by default
-!
-    implicit none
-
-! subroutine arguments
-!
-    type(block_data), pointer, intent(inout) :: pdata
-    real(kind=8)             , intent(in)    :: time, dt
-!
-!-------------------------------------------------------------------------------
-!
-#ifdef PROFILE
-! start accounting time for the shape update
-!
-    call start_timer(ims)
-#endif /* PROFILE */
-
-#ifdef PROFILE
-! stop accounting time for the shape update
-!
-    call stop_timer(ims)
-#endif /* PROFILE */
-
-!-------------------------------------------------------------------------------
-!
-  end subroutine update_shapes_user
-!
-!===============================================================================
-!
-! subroutine UPDATE_SOURCES_USER:
+! subroutine UPDATE_USER_SOURCES:
 ! ------------------------------
 !
 !   Subroutine adds the user defined source terms.
@@ -988,7 +904,7 @@ module user_problem
 !
 !===============================================================================
 !
-  subroutine update_sources_user(pdata, t, dt, du)
+  subroutine update_user_sources(pdata, t, dt, du)
 
 ! include external variables
 !
@@ -1003,7 +919,7 @@ module user_problem
     type(block_data), pointer       , intent(inout) :: pdata
     real(kind=8)                    , intent(in)    :: t, dt
     real(kind=8), dimension(:,:,:,:), intent(inout) :: du
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
@@ -1020,11 +936,60 @@ module user_problem
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine update_sources_user
+  end subroutine update_user_sources
 !
 !===============================================================================
 !
-! subroutine GRAVITATIONAL_ACCELERATION_USER:
+! subroutine UPDATE_USER_SHAPES:
+! -----------------------------
+!
+!   Subroutine defines the regions updated by user.
+!
+!   Arguments:
+!
+!     pdata - pointer to the data block structure of the currently initialized
+!             block;
+!     time  - time at the moment of update;
+!     dt    - time step since the last update;
+!
+!===============================================================================
+!
+  subroutine update_user_shapes(pdata, time, dt)
+
+! include external procedures and variables
+!
+    use blocks, only : block_data
+
+! local variables are not implicit by default
+!
+    implicit none
+
+! subroutine arguments
+!
+    type(block_data), pointer, intent(inout) :: pdata
+    real(kind=8)             , intent(in)    :: time, dt
+
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for the shape update
+!
+    call start_timer(ims)
+#endif /* PROFILE */
+
+#ifdef PROFILE
+! stop accounting time for the shape update
+!
+    call stop_timer(ims)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine update_user_shapes
+!
+!===============================================================================
+!
+! subroutine USER_GRAVITATIONAL_ACCELERATION:
 ! ------------------------------------------
 !
 !   Subroutine returns the user defined gravitational acceleration.
@@ -1037,7 +1002,7 @@ module user_problem
 !
 !===============================================================================
 !
-  subroutine gravitational_acceleration_user(t, dt, x, y, z, acc)
+  subroutine user_gravitational_acceleration(t, dt, x, y, z, acc)
 
 ! include external procedures and variables
 !
@@ -1052,7 +1017,7 @@ module user_problem
     real(kind=8)              , intent(in)  :: t, dt
     real(kind=8)              , intent(in)  :: x, y, z
     real(kind=8), dimension(3), intent(out) :: acc
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
@@ -1073,11 +1038,11 @@ module user_problem
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine gravitational_acceleration_user
+  end subroutine user_gravitational_acceleration
 !
 !===============================================================================
 !
-! subroutine BOUNDARY_USER_X:
+! subroutine USER_BOUNDARY_X:
 ! --------------------------
 !
 !   Subroutine updates ghost zones within the specific region along
@@ -1094,7 +1059,7 @@ module user_problem
 !
 !===============================================================================
 !
-  subroutine boundary_user_x(ic, jl, ju, kl, ku, t, dt, x, y, z, qn)
+  subroutine user_boundary_x(ic, jl, ju, kl, ku, t, dt, x, y, z, qn)
 
 ! import external procedures and variables
 !
@@ -1131,7 +1096,7 @@ module user_problem
 #if NDIMS == 3
     real(kind=8) :: dz, dxz
 #endif /* NDIMS == 3 */
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
@@ -1288,11 +1253,11 @@ module user_problem
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine boundary_user_x
+  end subroutine user_boundary_x
 !
 !===============================================================================
 !
-! subroutine BOUNDARY_USER_Y:
+! subroutine USER_BOUNDARY_Y:
 ! --------------------------
 !
 !   Subroutine updates ghost zones within the specific region along
@@ -1309,7 +1274,7 @@ module user_problem
 !
 !===============================================================================
 !
-  subroutine boundary_user_y(jc, il, iu, kl, ku, t, dt, x, y, z, qn)
+  subroutine user_boundary_y(jc, il, iu, kl, ku, t, dt, x, y, z, qn)
 
 ! import external procedures and variables
 !
@@ -1345,7 +1310,7 @@ module user_problem
     real(kind=8) :: dz, dyz
 #endif /* NDIMS == 3 */
     real(kind=8) :: fl, fr
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
@@ -1515,11 +1480,11 @@ module user_problem
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine boundary_user_y
+  end subroutine user_boundary_y
 !
 !===============================================================================
 !
-! subroutine BOUNDARY_USER_Z:
+! subroutine USER_BOUNDARY_Z:
 ! --------------------------
 !
 !   Subroutine updates ghost zones within the specific region along
@@ -1536,7 +1501,7 @@ module user_problem
 !
 !===============================================================================
 !
-  subroutine boundary_user_z(kc, il, iu, jl, ju, t, dt, x, y, z, qn)
+  subroutine user_boundary_z(kc, il, iu, jl, ju, t, dt, x, y, z, qn)
 
 ! local variables are not implicit by default
 !
@@ -1552,7 +1517,7 @@ module user_problem
     real(kind=8), dimension(:)      , intent(in)    :: y
     real(kind=8), dimension(:)      , intent(in)    :: z
     real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
@@ -1569,7 +1534,80 @@ module user_problem
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine boundary_user_z
+  end subroutine user_boundary_z
+!
+!===============================================================================
+!
+! subroutine USER_TIME_STATISTICS:
+! -------------------------------
+!
+!   Subroutine can be use to store user defined time statistics. The file to
+!   store these statistics should be properly created in subroutine
+!   initialize_user_problem() and closed in finalize_user_problem().
+!
+!===============================================================================
+!
+  subroutine user_time_statistics()
+
+    implicit none
+
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+    call start_timer(imt)
+#endif /* PROFILE */
+
+#ifdef PROFILE
+    call stop_timer(imt)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine user_time_statistics
+!
+!===============================================================================
+!
+! subroutine LOG_COSH:
+! -------------------
+!
+!   Function calculates the logarithm of the hyperbolic cosine, which is
+!   the result of the integration of tanh(x). Direct calculation using
+!   Fortran intrinsic subroutines fails for large values of x, therefore
+!   the logarithm of cosh is approximated as |x| + log(1/2) for
+!   |x| > threshold.
+!
+!   Arguments:
+!
+!     x - function argument;
+!
+!===============================================================================
+!
+  function log_cosh(x) result(y)
+
+! local variables are not implicit by default
+!
+    implicit none
+
+! function arguments
+!
+    real(kind=8), intent(in) :: x
+    real(kind=8)             :: y
+
+! local parameters
+!
+    real(kind=8), parameter :: th = acosh(huge(x)), lh = log(0.5d+00)
+
+!-------------------------------------------------------------------------------
+!
+    if (abs(x) < th) then
+      y = log(cosh(x))
+    else
+      y = abs(x) + lh
+    end if
+
+!-------------------------------------------------------------------------------
+!
+  end function log_cosh
 
 !===============================================================================
 !