diff --git a/src/equations.F90 b/src/equations.F90
index e75506c..f843f1b 100644
--- a/src/equations.F90
+++ b/src/equations.F90
@@ -1165,8 +1165,8 @@ module equations
 !
       write(sid,'(i15)') id
       write(snc,'(i15)') nc
-      write(msg,'("Unphysical cells in block ID:",a," (",a," counted)")')      &
-                                                     trim(sid), trim(snc)
+      write(msg,'("Unphysical cells in block ID:",a," (",a," counted).")')     &
+                                        trim(adjustl(sid)), trim(adjustl(snc))
       call print_warning(loc, trim(msg))
 
 ! allocate temporary vectors for unphysical states
@@ -3522,7 +3522,7 @@ module equations
       en  = u(ien,i) + u(idn,i)
       dn  = u(idn,i)
 
-! find the exact W using an Newton-Ralphson interative method
+! find the exact W using the Newton-Ralphson interative method
 !
       call nr_iterate(mm, bb, mb, en, dn, w, vv, info)
 
@@ -3873,124 +3873,135 @@ module equations
     call start_timer(imp)
 #endif /* PROFILE */
 
+! check if the state is physical; this can save some time if unphysical state
+! is considered
+!
+    wl = sqrt(mm + dn * dn)
+    if (en > wl) then
+
 ! prepare the initial brackets
 !
-    wl   = sqrt(mm + dn * dn) + gammaxi * pmin
-    wu   = en + pmin
+      wl = wl + gammaxi * pmin
+      wu = en + pmin
 
 ! calculate the value of function for the lower bracket
 !
-    call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
+      call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
 
 ! the lower bracket gives negative function, so there is chance it bounds
 ! the root
 !
-    if (fl < 0.0d+00) then
+      if (fl < 0.0d+00) then
 
 ! make sure that the upper bracket is larger than the lower one and
 ! the function has positive value
 !
-      if (wu <= wl) wu = 2.0d+00 * wl
+        if (wu <= wl) wu = 2.0d+00 * wl
 
 ! check if the brackets bound the root region, if not proceed until
 ! opposite function signs are found for the brackets
 !
-      call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
-      it   = nrmax
-      keep = fl * fu > 0.0d+00
-      do while (keep)
-        it = it - 1
-        wl = wu
-        fl = fu
-        wu = 2.0d+00 * wu
         call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
-        keep = (fl * fu > 0.0d+00) .and. it > 0
-      end do
+        it   = nrmax
+        keep = fl * fu > 0.0d+00
+        do while (keep)
+          it = it - 1
+          wl = wu
+          fl = fu
+          wu = 2.0d+00 * wu
+          call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
+          keep = (fl * fu > 0.0d+00) .and. it > 0
+        end do
 
 ! the upper bracket was found, so proceed with determining the root
 !
-      if (it > 0 .and. fu >= 0.0d+00) then
+        if (it > 0 .and. fu >= 0.0d+00) then
 
 ! estimate the value of enthalpy close to the root and corresponding v²
 !
-        w = wl - fl * (wu - wl) / (fu - fl)
+          w = wl - fl * (wu - wl) / (fu - fl)
 
 ! initialize iteration parameters
 !
-        info = .true.
-        keep = .true.
-        it   = nrmax
-        cn   = nrext
+          info = .true.
+          keep = .true.
+          it   = nrmax
+          cn   = nrext
 
 ! iterate using the Newton-Raphson method in order to find a root w of the
 ! function
 !
-        do while(keep)
+          do while(keep)
 
 ! calculate F(W) and dF(W)/dW
 !
-          call nr_function_srhd_adi_1d(mm, en, dn, w, f, df)
+            call nr_function_srhd_adi_1d(mm, en, dn, w, f, df)
 
 ! update brackets
 !
-          if (f > fl .and. f < 0.0d+00) then
-            wl = w
-            fl = f
-          end if
-          if (f < fu .and. f > 0.0d+00) then
-            wu = w
-            fu = f
-          end if
+            if (f > fl .and. f < 0.0d+00) then
+              wl = w
+              fl = f
+            end if
+            if (f < fu .and. f > 0.0d+00) then
+              wu = w
+              fu = f
+            end if
 
 ! calculate the increment dW, update the solution, and estimate the error
 !
-          dw  = f / df
-          w   = w - dw
-          err = abs(dw / w)
+            dw  = f / df
+            w   = w - dw
+            err = abs(dw / w)
 
 ! check the convergence, if the convergence is not reached, iterate until
 ! the maximum number of iteration is reached
 !
-          if (err < tol) then
-            keep = cn > 0
-            cn   = cn - 1
-          else
-            keep = it > 0
-          end if
+            if (err < tol) then
+              keep = cn > 0
+              cn   = cn - 1
+            else
+              keep = it > 0
+            end if
 
 ! if new W leaves the brackets, use the bisection method to estimate
 ! the new guess
 !
-          if (w < wl .or. w > wu) then
-            w = 0.5d+00 * (wl + wu)
-          end if
+            if (w < wl .or. w > wu) then
+              w = 0.5d+00 * (wl + wu)
+            end if
 
-          it = it - 1
-        end do ! NR iterations
+            it = it - 1
+          end do ! NR iterations
 
 ! print information about failed convergence or unphysical variables
 !
-        if (err >= tol) then
-          call print_warning(loc, "Convergence not reached!")
-          write(*,"(a,1x,1e24.16e3)") "Error: ", err
-        end if
+          if (err >= tol) then
+            call print_warning(loc, "Convergence not reached!")
+            write(*,"(a,1x,1e24.16e3)") "Error: ", err
+          end if
 
 ! calculate |V|² from W
 !
-        vv  = mm / (w * w)
+          vv  = mm / (w * w)
 
-      else ! the upper brack not found
-        call print_warning(loc, "Could not find the upper bracket!")
+        else ! the upper brack not found
+          call print_warning(loc, "Could not find the upper bracket!")
+          info = .false.
+
+        end if
+
+      else if (fl == 0.0d+00) then ! the lower bracket is a root, so return it
+        w    = wl
+        info = .true.
+
+      else ! the root cannot be found, since it is below the lower bracket
+        call print_warning(loc, "Positive function for lower bracket!")
         info = .false.
-
       end if
 
-    else if (fl == 0.0d+00) then ! the lower bracket is a root, so return it
-      w    = wl
-      info = .true.
-
-    else ! the root cannot be found, since it is below the lower bracket
-      call print_warning(loc, "Positive function for lower bracket!")
+    else ! the state is unphysical
+      call print_warning(loc, "The state is not physical!")
       info = .false.
     end if
 
@@ -5363,7 +5374,8 @@ module equations
 
 ! include external procedures
 !
-    use algebra        , only : cubic_normalized, quartic
+    use algebra, only : cubic_normalized, quartic
+    use error  , only : print_warning
 
 ! local variables are not implicit by default
 !
@@ -5387,6 +5399,11 @@ module equations
 !
     real(kind=8), dimension(5) :: a
     real(kind=8), dimension(4) :: x
+
+! local parameters
+!
+    character(len=*), parameter :: loc =                                       &
+                                  'EQUATIONS::nr_initial_brackets_srmhd_adi()'
 !
 !-------------------------------------------------------------------------------
 !
@@ -5450,29 +5467,19 @@ module equations
       keep = .true.
       it   = nrmax
       wl   = wu
-
       do while(keep)
-
         call nr_pressure_srmhd_adi_1d(mm, bb, mb, ec, dn, wl, f, df)
-
         dw   = f / df
         wl   = wl - dw
-
         err  = abs(dw / wl)
         it   = it - 1
         keep = (err > tol) .and. it > 0
-
       end do
-
       if (it <= 0) then
-        write(*,*)
-        write(*,"(a,1x,a)") "ERROR in"                                         &
-                          , "EQUATIONS::nr_initial_brackets_srmhd_adi()"
-        write(*,"(a)"     ) "Could not find the lower limit for the enthalpy!"
+        call print_warning(loc,                                                &
+                         "Iterative solver failed to find the lower bracket!")
         write(*,"(a,5(1x,1e24.16e3))") " D, |m|², m.B, |B|², E = "             &
                                                           , dn, mm, mb, bb, en
-        info = .false.
-        return
       end if
 
     end if ! nr > 0
@@ -5481,71 +5488,53 @@ module equations
 !
     call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wl, fl)
 
-    if (fl > 0.0d+00) then
-      write(*,*)
-      write(*,"(a,1x,a)") "ERROR in"                                           &
-                        , "EQUATIONS::nr_initial_brackets_srmhd_adi()"
-      write(*,"(a)"     ) "Lower limit positive!"
-      write(*,"(a,6(1x,1e24.16e3))") " D, |m|², m.B, |B|², E, W = "            &
-                                                      , dn, mm, mb, bb, en, wl
-      info = .false.
-      return
-    end if
+    if (fl <= 0.0d+00) then
 
 ! make sure that the upper limit is larger than the lower one
 !
-    keep = wl >= wu
-    it   = nrmax
-    do while(keep)
-      wu   = 2.0d+00 * wu
-      it   = it - 1
-      keep = (wl >= wu) .and. it > 0
-    end do
-    if (it <= 0) then
-      write(*,*)
-      write(*,"(a,1x,a)") "ERROR in"                                           &
-                        , "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
-      write(*,"(a)"     ) "Could not find the upper limit for enthalpy!"
-      write(*,"(a,6(1x,1e24.16e3))") " D, |m|², m.B, |B|², E, W = "            &
-                                                      , dn, mm, mb, bb, en, wl
-      info = .false.
-      return
-    end if
+      if (wu <= wl) wu = 2.0d+00 * wl
 
 ! check if the brackets bound the root region, if not proceed until
 ! opposite function signs are found for the brackets
 !
-    call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wu, fu)
-
-    keep = (fl * fu > 0.0d+00)
-    it   = nrmax
-
-    do while (keep)
-
-      wl = wu
-      fl = fu
-      wu = 2.0d+00 * wu
-
       call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wu, fu)
+      it   = nrmax
+      keep = fl * fu > 0.0d+00
+      do while (keep)
+        it = it - 1
+        wl = wu
+        fl = fu
+        wu = 2.0d+00 * wu
+        call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wu, fu)
+        keep = (fl * fu > 0.0d+00) .and. it > 0
+      end do
 
-      it   = it - 1
-      keep = (fl * fu > 0.0d+00) .and. it > 0
-
-    end do
-    if (it <= 0) then
-      write(*,*)
-      write(*,"(a,1x,a)") "ERROR in"                                           &
-                        , "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
-      write(*,"(a)"     ) "No initial brackets found!"
-      write(*,"(a,5(1x,1e24.16e3))") " D, |m|², m.B, |B|², E = "               &
-                                                          , dn, mm, mb, bb, en
-      info = .false.
-      return
-    end if
+! the upper bracket was found, so proceed with determining the root
+!
+      if (it > 0 .and. fu >= 0.0d+00) then
 
 ! estimate the enthalpy value close to the root
 !
-    wc   = wl - fl * (wu - wl) / (fu - fl)
+        wc = wl - fl * (wu - wl) / (fu - fl)
+
+! we have good brackets and guess, so good to go
+!
+        info = .true.
+
+      else ! the upper brack not found
+        call print_warning(loc, "Could not find the upper bracket!")
+        write(*,"(a,5(1x,1e24.16e3))") " D, |m|², m.B, |B|², E = "             &
+                                                          , dn, mm, mb, bb, en
+        info = .false.
+
+      end if
+
+    else ! the root cannot be found, since it is below the lower bracket
+      call print_warning(loc, "Positive function for lower bracket!")
+      write(*,"(a,6(1x,1e24.16e3))") " D, |m|², m.B, |B|², E, W = "            &
+                                                      , dn, mm, mb, bb, en, wl
+      info = .false.
+    end if
 
 #ifdef PROFILE
 ! stop accounting time for the initial brackets
@@ -5587,6 +5576,10 @@ module equations
 !
   subroutine nr_iterate_srmhd_adi_1dw(mm, bb, mb, en, dn, w, vv, info)
 
+! include external procedures
+!
+    use error, only : print_warning
+
 ! local variables are not implicit by default
 !
     implicit none
@@ -5604,6 +5597,10 @@ module equations
     real(kind=8) :: wl, wu, fl, fu
     real(kind=8) :: f , df, dw
     real(kind=8) :: err
+
+! local parameters
+!
+    character(len=*), parameter :: loc = 'EQUATIONS::nr_iterate_srmhd_adi_1dw()'
 !
 !-------------------------------------------------------------------------------
 !
@@ -5616,104 +5613,80 @@ module equations
 ! find the initial brackets and estimate the initial enthalpy
 !
     call nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn                      &
-                                                     , wl, wu, w, fl, fu, info)
+                                                 , wl, wu, w, fl, fu, info)
 
-! if the brackets could not be found, return the lower bracket as the solution
+! continue if brackets found
 !
-    if (.not. info) then
-      write(*,*)
-      write(*,"(a,1x,a)"        ) "WARNING in"                                 &
-                                , "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
-      write(*,"(a,1x)"          ) "The solution lays in unphysical regime."
-      write(*,"(a,1x,1e24.16e3)") "Using the lower bracket as solution: ", wl
+    if (info) then
 
-! use the lower bracket, since it guarantees the positive pressure
+! initialize iteration parameters
 !
-      w = wl
+      info = .true.
+      keep = .true.
+      it   = nrmax
+      cn   = nrext
+
+! iterate using the Newton-Raphson method in order to find a root w of the
+! function
+!
+      do while(keep)
+
+! calculate F(W) and dF(W)/dW
+!
+        call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, w, f, df)
+
+! update brackets
+!
+        if (f > fl .and. f < 0.0d+00) then
+          wl = w
+          fl = f
+        end if
+        if (f < fu .and. f > 0.0d+00) then
+          wu = w
+          fu = f
+        end if
+
+! calculate the increment dW, update the solution, and estimate the error
+!
+        dw  = f / df
+        w   = w - dw
+        err = abs(dw / w)
+
+! check the convergence, if the convergence is not reached, iterate until
+! the maximum number of iteration is reached
+!
+        if (err < tol) then
+          keep = cn > 0
+          cn   = cn - 1
+        else
+          keep = it > 0
+        end if
+
+! if new W lays out of the brackets, use the bisection method to estimate
+! the new guess
+!
+        if (w < wl .or. w > wu) then
+          w = 0.5d+00 * (wl + wu)
+        end if
+
+! decrease the number of remaining iterations
+!
+        it = it - 1
+
+      end do ! NR iterations
+
+! let know the user if the convergence failed
+!
+      if (err >= tol) then
+        call print_warning(loc, "Convergence not reached!")
+        write(*,"(a,1x,1e24.16e3)") "Error: ", err
+      end if
 
 ! calculate |V|² from W
 !
       call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
 
-      info = .true.
-      return
-
-    end if
-
-! initialize iteration parameters
-!
-    info = .true.
-    keep = .true.
-    it   = nrmax
-    cn   = nrext
-
-! iterate using the Newton-Raphson method in order to find a root w of the
-! function
-!
-    do while(keep)
-
-! calculate F(W) and dF(W)/dW
-!
-      call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, w, f, df)
-
-! update brackets
-!
-      if (f > fl .and. f < 0.0d+00) then
-        wl = w
-        fl = f
-      end if
-      if (f < fu .and. f > 0.0d+00) then
-        wu = w
-        fu = f
-      end if
-
-! calculate the increment dW
-!
-      dw  = f / df
-
-! update the solution
-!
-      w   = w - dw
-
-! calculate the error
-!
-      err = abs(dw / w)
-
-! check the convergence, if the convergence is not reached, iterate until
-! the maximum number of iteration is reached
-!
-      if (err < tol) then
-        keep = cn > 0
-        cn   = cn - 1
-      else
-        keep = it > 0
-      end if
-
-! if new W lays out of the brackets, use the bisection method to estimate
-! the new guess
-!
-      if (w < wl .or. w > wu) then
-        w = 0.5d+00 * (wl + wu)
-      end if
-
-! decrease the number of remaining iterations
-!
-      it = it - 1
-
-    end do ! NR iterations
-
-! calculate |V|² from W
-!
-    call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
-
-! let know the user if the convergence failed
-!
-    if (err >= tol) then
-      write(*,*)
-      write(*,"(a,1x,a)"        ) "WARNING in"                                 &
-                                , "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
-      write(*,"(a,1x,1e24.16e3)") "Convergence not reached: ", err
-    end if
+    end if ! correct brackets
 
 #ifdef PROFILE
 ! stop accounting time for variable solver