diff --git a/build/makefile b/build/makefile
index b1ed814..4793f3d 100644
--- a/build/makefile
+++ b/build/makefile
@@ -49,16 +49,6 @@ $(info )
 $(shell sleep 15)
 endif
 
-#-------------------------------------------------------------------------------
-#
-# generate object file dependencies
-#
-ifneq ($(wildcard makedeps),makedeps)
-$(info Generating dependencies.)
-$(shell ./mkdeps.sh $(SRCSDIR) $(OBJSDIR) > makedeps)
-$(info )
-endif
-
 #-------------------------------------------------------------------------------
 #
 # check flag conditions
@@ -80,6 +70,14 @@ FFLAGS  += ${CPPPREFIX}-DNDIMS=${NDIMS}
 #
 FFLAGS  += ${CPPPREFIX}-D${OUTPUT}
 
+# add module path to compiler options
+#
+ifeq ($(COMPILER),GNU)
+FFLAGS  += -J $(OBJSDIR)
+else
+FFLAGS  += -module $(OBJSDIR)
+endif
+
 #-------------------------------------------------------------------------------
 
 .SUFFIXES:
@@ -89,21 +87,22 @@ FFLAGS  += ${CPPPREFIX}-D${OUTPUT}
 
 name    = amun
 
-modules = algebra blocks boundaries constants coordinates domains driver       \
+files   = algebra blocks boundaries constants coordinates domains driver       \
           equations error evolution gravity integrals interpolations io mesh   \
           mpitools operators parameters problems random refinement schemes     \
           shapes sources timers user_problem utils
 
-sources := $(addprefix $(SRCSDIR)/,$(addsuffix .F90, $(modules)))
-objects := $(addprefix $(OBJSDIR)/,$(addsuffix .o,   $(modules)))
+sources := $(addprefix $(SRCSDIR)/,$(addsuffix .F90, $(files)))
+objects := $(addprefix $(OBJSDIR)/,$(addsuffix .o,   $(files)))
+modules := $(addprefix $(OBJSDIR)/,$(addsuffix .mod, $(files)))
 
 all: $(name).x
 
-$(name).x: $(deps) $(objects) | $(DESTDIR)
+$(name).x: $(objects) | $(DESTDIR)
 	$(LD) $(LDFLAGS) $(objects) $(LIBS) -o $(DESTDIR)/$(name).x
 
 $(OBJSDIR)/%.o : $(SRCSDIR)/%.F90
-	$(FC) -c $(FFLAGS) -J $(OBJSDIR) $< -o $@
+	$(FC) -c $(FFLAGS) $< -o $@
 
 $(objects): | $(OBJSDIR)
 
@@ -113,8 +112,14 @@ $(OBJSDIR):
 $(DESTDIR):
 	mkdir -p $(DESTDIR)
 
+makedeps: $(sources)
+	./mkdeps.sh $(SRCSDIR) $(OBJSDIR) > makedeps
+
+.PHONY: clean
+
 clean:
-	rm -rf $(OBJSDIR) $(DESTDIR)/$(name).x
+	rm -f $(objects) $(modules) $(DESTDIR)/$(name).x makedeps
+	if [ -d $(OBJSDIR) ]; then rmdir $(OBJSDIR); fi
 
 #-------------------------------------------------------------------------------
 #
diff --git a/python/amun.py b/python/amun.py
index 863e064..76261e9 100644
--- a/python/amun.py
+++ b/python/amun.py
@@ -418,6 +418,103 @@ def rebin(a, newshape):
       a = np.repeat(a, newshape[n] / a.shape[n], axis = n)
     return(a)
 
+
+def amun_integrals(field, filename, pathlist):
+    '''
+        get_integral: iterate over pathlist and read and merge field values from filename files in the provided paths
+    '''
+    # Initiate the return values with empty array and file number.
+    #
+    vals = array([])
+    num  = 1
+
+    # Iterate over all paths provided in the list 'pathlist'.
+    #
+    for path in pathlist:
+
+      # Iterate over all files in the current path.
+      #
+      while True:
+
+        # Generate file name.
+        #
+        dfile = path + '/' + filename + '_' + str(num).zfill(2) + '.dat'
+
+        # Check if the file exists.
+        #
+        if isfile(dfile):
+
+          # Read values from the current integrals file.
+          #
+          lvals = read_integrals(dfile, field)
+
+          # Append to the return array.
+          #
+          vals = append(vals, lvals)
+
+          # Increase the number file.
+          #
+          num = num + 1
+
+        else:
+
+          # File does not exists, so go to the next path.
+          #
+          break
+
+    # Return appended values.
+    #
+    return vals
+
+def read_integrals(filename, column):
+    '''
+        read_integrals: reads a given column from an integral file.
+    '''
+    # Open the given file and check if it is text file.
+    #
+    f = open(filename, 'r')
+
+    # Read fist line and store it in h, since it will be used to obtain the
+    # column headers.
+    #
+    l = f.readline()
+    h = l
+
+    # Read first line which is not comment in order to determine the number of
+    # columns and store the number of columns in nc.  Calculate the column width
+    # and store it in wc.
+    #
+    while l.startswith('#'):
+        l = f.readline()
+    nc = len(l.rsplit())
+    wc = int((len(l) - 9) / (nc - 1))
+
+    # Split header line into a list.
+    #
+    lh = [h[1:9].strip()]
+    for i in range(nc - 1):
+      ib = i * wc + 10
+      ie = ib + wc - 1
+      lh.append(h[ib:ie].strip())
+
+
+    ic = lh.index(column)
+
+    # Read given column.
+    #
+    if (ic > -1):
+      lc = [float(l.split()[ic])]
+      for l in f:
+        lc.append(float(l.split()[ic]))
+
+    # Close the file.
+    #
+    f.close()
+
+    # Return values.
+    #
+    return(array(lc))
+
 if __name__ == "__main__":
   fname = './p000030_00000.h5'
 
diff --git a/src/blocks.F90 b/src/blocks.F90
index 181fd64..f529178 100644
--- a/src/blocks.F90
+++ b/src/blocks.F90
@@ -343,7 +343,7 @@ module blocks
   public :: set_last_id, get_last_id, get_mblocks, get_dblocks, get_nleafs
   public :: set_blocks_update
   public :: change_blocks_process
-  public :: set_neighbors_refine
+  public :: set_neighbors_refine, set_neighbors_update
   public :: metablock_set_id, metablock_set_process, metablock_set_level
   public :: metablock_set_configuration, metablock_set_refinement
   public :: metablock_set_position, metablock_set_coordinates
diff --git a/src/boundaries.F90 b/src/boundaries.F90
index 32986d2..fc49aca 100644
--- a/src/boundaries.F90
+++ b/src/boundaries.F90
@@ -5765,6 +5765,7 @@ module boundaries
     use coordinates    , only : im , jm , km
     use coordinates    , only : faces_dp
     use equations      , only : nv, idn, ipr
+    use error          , only : print_warning
     use interpolations , only : limiter_prol
 
 ! local variables are not implicit by default
@@ -5842,9 +5843,15 @@ module boundaries
             dq(3) = limiter_prol(0.5d+00, dql, dqr)
 
             if (p == idn .or. p == ipr) then
-              do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
-                dq(:) = 0.5d+00 * dq(:)
-              end do
+              if (qn(p,i,j,k) > 0.0d+00) then
+                do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
+                  dq(:) = 0.5d+00 * dq(:)
+                end do
+              else
+                call print_warning("boundaries::block_face_prolong"            &
+                                       , "Positive variable is not positive!")
+                dq(:) = 0.0d+00
+              end if
             end if
 
             dq(:) = 0.5d+00 * dq(:)
@@ -6063,6 +6070,7 @@ module boundaries
     use coordinates    , only : im , jm , km
     use coordinates    , only : edges_dp
     use equations      , only : nv, idn, ipr
+    use error          , only : print_warning
     use interpolations , only : limiter_prol
 
 ! local variables are not implicit by default
@@ -6155,9 +6163,15 @@ module boundaries
 #endif /* NDIMS == 3 */
 
             if (p == idn .or. p == ipr) then
-              do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
-                dq(:) = 0.5d+00 * dq(:)
-              end do
+              if (qn(p,i,j,k) > 0.0d+00) then
+                do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
+                  dq(:) = 0.5d+00 * dq(:)
+                end do
+              else
+                call print_warning("boundaries::block_edge_prolong"            &
+                                       , "Positive variable is not positive!")
+                dq(:) = 0.0d+00
+              end if
             end if
 
             dq(:) = 0.5d+00 * dq(:)
@@ -6330,6 +6344,7 @@ module boundaries
     use coordinates    , only : im , jm , km
     use coordinates    , only : corners_dp
     use equations      , only : nv, idn, ipr
+    use error          , only : print_warning
     use interpolations , only : limiter_prol
 
 ! local variables are not implicit by default
@@ -6427,9 +6442,15 @@ module boundaries
 #endif /* NDIMS == 3 */
 
             if (p == idn .or. p == ipr) then
-              do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
-                dq(:) = 0.5d+00 * dq(:)
-              end do
+              if (qn(p,i,j,k) > 0.0d+00) then
+                do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
+                  dq(:) = 0.5d+00 * dq(:)
+                end do
+              else
+                call print_warning("boundaries::block_corner_prolong"          &
+                                       , "Positive variable is not positive!")
+                dq(:) = 0.0d+00
+              end if
             end if
 
             dq(:) = 0.5d+00 * dq(:)
diff --git a/src/equations.F90 b/src/equations.F90
index 6c8c3c4..380677e 100644
--- a/src/equations.F90
+++ b/src/equations.F90
@@ -139,7 +139,7 @@ module equations
 
 ! the lower limits for density and pressure to be treated as physical
 !
-  real(kind=8)     , save :: dmin    = 1.0d-08, pmin    = 1.0d-08
+  real(kind=8)     , save :: dmin    = 1.0d-16, pmin    = 1.0d-16
 
 ! the upper limits for the Lorentz factor and corresponding |v|²
 !
@@ -158,9 +158,12 @@ module equations
   integer          , save :: nrmax   = 100
   integer          , save :: nrext   = 2
 
-! flags for corrections
+! flag for unphysical cells correction, the maximum distance of neighbors for
+! averaging region, and the minimum number of cells for averaging
 !
   logical          , save :: fix_unphysical_cells = .false.
+  integer          , save :: ngavg   = 2
+  integer          , save :: npavg   = 4
 
 ! by default everything is private
 !
@@ -868,6 +871,16 @@ module equations
       fix_unphysical_cells = .false.
     end select
 
+! get parameters for unphysical cells correction
+!
+    call get_parameter_integer("ngavg", ngavg)
+    call get_parameter_integer("npavg", npavg)
+
+! correct the above parameters to reasonable values
+!
+    ngavg = max(1, ngavg)
+    npavg = max(2, npavg)
+
 ! print information about the equation module
 !
     if (verbose) then
@@ -879,6 +892,10 @@ module equations
       end if
       write (*,"(4x,a20, 3x,'=',1x,a)") "fix unphysical cells"                 &
                                                         , trim(unphysical_fix)
+      if (fix_unphysical_cells) then
+        write (*,"(4x,a20, 3x,'=',1x,i4)") "ngavg               ", ngavg
+        write (*,"(4x,a20, 3x,'=',1x,i4)") "npavg               ", npavg
+      end if
 
     end if
 
@@ -1193,7 +1210,7 @@ module equations
 !
               np = 0
               p  = 1
-              do while (np <= 2 .and. p <= 4)
+              do while (np <= npavg .and. p <= ngavg)
                 il = max( 1, i - p)
                 iu = min(im, i + p)
                 jl = max( 1, j - p)
@@ -1208,7 +1225,7 @@ module equations
 
 ! average primitive variables
 !
-              if (np > 2) then
+              if (np >= npavg) then
 
                 do p = 1, nv
                   q(p,n) = sum(qq(p,il:iu,jl:ju,kl:ku),                        &
@@ -1218,13 +1235,17 @@ module equations
 
               else
 
-! print error, since no physical cells found for averaging
+! limit density or pressure to minimum value, since the averaging over
+! neighbours failed
 !
                 write(msg,'(a,1x,a)')                                          &
-                              "Cannot correct the unphysical cell."            &
-                            , "Not sufficient number of physical neighbors!"
-                call print_error(loc, trim(msg))
-                stop
+                              "Not sufficient number of physical neighbors!"   &
+                            , "Applying lower bounds for positive variables."
+                call print_warning(loc, trim(msg))
+
+                q(1:nv,n) = qq(1:nv,i,j,k)
+                q(idn ,n) = max(dmin, qq(idn,i,j,k))
+                if (ipr > 0) q(ipr,n) = max(pmin, qq(ipr,i,j,k))
 
               end if ! not sufficient number of physical cells for averaging
 
diff --git a/src/evolution.F90 b/src/evolution.F90
index c38b3f4..dc86c5c 100644
--- a/src/evolution.F90
+++ b/src/evolution.F90
@@ -1988,6 +1988,7 @@ module evolution
 
 ! include external procedures
 !
+    use blocks        , only : set_neighbors_update
     use boundaries    , only : boundary_variables
     use equations     , only : update_primitive_variables
     use equations     , only : fix_unphysical_cells, correct_unphysical_states
@@ -2037,12 +2038,28 @@ module evolution
 ! correct unphysical states if detected
 !
     if (fix_unphysical_cells) then
+
+! if an unphysical cell appeared in a block while updating its primitive
+! variables it could be propagated to its neighbors through boundary update;
+! mark all neighbors of such a block to be verified and corrected for
+! unphysical cells too
+!
+      pmeta => list_meta
+      do while (associated(pmeta))
+
+        if (pmeta%update) call set_neighbors_update(pmeta)
+
+        pmeta => pmeta%next
+      end do
+
+! verify and correct, if necessary, unphysical cells in recently updated blocks
+!
       pdata => list_data
       do while (associated(pdata))
         pmeta => pdata%meta
 
-        if (pmeta%update) call correct_unphysical_states(pmeta%id              &
-                                                           , pdata%q, pdata%u)
+        if (pmeta%update)                                                      &
+                    call correct_unphysical_states(pmeta%id, pdata%q, pdata%u)
 
         pdata => pdata%next
       end do
diff --git a/src/interpolations.F90 b/src/interpolations.F90
index 617d486..40af9e2 100644
--- a/src/interpolations.F90
+++ b/src/interpolations.F90
@@ -695,6 +695,7 @@ module interpolations
     use coordinates    , only : im , jm , km
     use coordinates    , only : ib , jb , kb , ie , je , ke
     use coordinates    , only : ibl, jbl, kbl, ieu, jeu, keu
+    use error          , only : print_warning
 
 ! local variables are not implicit by default
 !
@@ -768,9 +769,15 @@ module interpolations
 ! variables
 !
           if (positive) then
-            do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
-              dq(:) = 0.5d+00 * dq(:)
-            end do
+            if (q(i,j,k) > 0.0d+00) then
+              do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
+                dq(:) = 0.5d+00 * dq(:)
+              end do
+            else
+              call print_warning("interpolations::interfaces_tvd"              &
+                                       , "Positive variable is not positive!")
+              dq(:) = 0.0d+00
+            end if
           end if
 
 ! interpolate states
@@ -933,6 +940,7 @@ module interpolations
     use coordinates    , only : im , jm , km
     use coordinates    , only : ib , jb , kb , ie , je , ke
     use coordinates    , only : ibl, jbl, kbl, ieu, jeu, keu
+    use error          , only : print_warning
 
 ! local variables are not implicit by default
 !
@@ -1051,9 +1059,15 @@ module interpolations
 ! variables
 !
             if (positive) then
-              do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
-                dq(:) = 0.5d+00 * dq(:)
-              end do
+              if (q(i,j,k) > 0.0d+00) then
+                do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
+                  dq(:) = 0.5d+00 * dq(:)
+                end do
+              else
+                call print_warning("interpolations::interfaces_mgp"            &
+                                       , "Positive variable is not positive!")
+                dq(:) = 0.0d+00
+              end if
             end if
 
 ! interpolate states
diff --git a/src/mesh.F90 b/src/mesh.F90
index 5327e65..f11343e 100644
--- a/src/mesh.F90
+++ b/src/mesh.F90
@@ -1037,6 +1037,7 @@ module mesh
     use coordinates    , only : ng, nh, in, jn, kn, im, jm, km
     use coordinates    , only : ib, ie, jb, je, kb, ke
     use equations      , only : nv, idn, ien
+    use error          , only : print_warning
     use interpolations , only : limiter_prol
 
 ! local variables are not implicit by default
@@ -1137,9 +1138,15 @@ module mesh
 #endif /* NDIMS == 3 */
 
             if (p == idn .or. p == ien) then
-              do while (pdata%u(p,i,j,k) <= sum(abs(du(1:NDIMS))))
-                du(:) = 0.5d+00 * du(:)
-              end do
+              if (pdata%u(p,i,j,k) > 0.0d+00) then
+                do while (pdata%u(p,i,j,k) <= sum(abs(du(1:NDIMS))))
+                  du(:) = 0.5d+00 * du(:)
+                end do
+              else
+                call print_warning("mesh::prolong_block"                       &
+                                       , "Positive variable is not positive!")
+                du(:) = 0.0d+00
+              end if
             end if
 
             du(:) = 0.5d+00 * du(:)