diff --git a/sources/boundaries.F90 b/sources/boundaries.F90
index 921e9e3..ce4769f 100644
--- a/sources/boundaries.F90
+++ b/sources/boundaries.F90
@@ -57,7 +57,7 @@ module boundaries
 ! supported boundary types
 !
   enum, bind(c)
-    enumerator bnd_user
+    enumerator bnd_custom
     enumerator bnd_periodic
     enumerator bnd_open
     enumerator bnd_outflow
@@ -98,163 +98,98 @@ module boundaries
 ! subroutine INITIALIZE_BOUNDARIES:
 ! --------------------------------
 !
-!   Subroutine initializes the module BOUNDARIES by setting its parameters.
+!   Subroutine initializes the module BOUNDARIES.
 !
 !   Arguments:
 !
-!     status  - return flag of the procedure execution status;
+!     status  - the subroutine call status;
 !
 !===============================================================================
 !
   subroutine initialize_boundaries(status)
 
-! import external procedures and variables
-!
     use coordinates, only : periodic
+    use helpers    , only : print_message
 #ifdef MPI
     use mpitools   , only : npmax
 #endif /* MPI */
     use parameters , only : get_parameter
 
-! local variables are not implicit by default
-!
     implicit none
 
-! subroutine arguments
-!
     integer, intent(out) :: status
 
-! module parameters for the boundary update order and boundary type
-!
-    character(len = 32) :: xlbndry = "periodic"
-    character(len = 32) :: xubndry = "periodic"
-    character(len = 32) :: ylbndry = "periodic"
-    character(len = 32) :: yubndry = "periodic"
-    character(len = 32) :: zlbndry = "periodic"
-    character(len = 32) :: zubndry = "periodic"
+    character(len=64) :: str, bnd
+
+    integer :: n, s
+
+    character(len=*), parameter :: loc = 'BOUNDARIES::initialize_boundaries()'
 
-! local variables
-!
-    integer :: n
-!
 !-------------------------------------------------------------------------------
 !
     status = 0
 
-! get runtime values for the boundary types
-!
-    call get_parameter("xlbndry", xlbndry)
-    call get_parameter("xubndry", xubndry)
-    call get_parameter("ylbndry", ylbndry)
-    call get_parameter("yubndry", yubndry)
-    call get_parameter("zlbndry", zlbndry)
-    call get_parameter("zubndry", zubndry)
-
-! fill the boundary type flags
-!
-    select case(xlbndry)
-    case("open")
-      bnd_type(1,1) = bnd_open
-    case("outflow", "out")
-      bnd_type(1,1) = bnd_outflow
-    case("reflective", "reflecting", "reflect")
-      bnd_type(1,1) = bnd_reflective
-    case("hydrostatic", "gravity")
-      bnd_type(1,1) = bnd_gravity
-    case("user", "custom")
-      bnd_type(1,1) = bnd_user
-    case default
-      bnd_type(1,1) = bnd_periodic
-    end select
-
-    select case(xubndry)
-    case("open")
-      bnd_type(1,2) = bnd_open
-    case("outflow", "out")
-      bnd_type(1,2) = bnd_outflow
-    case("reflective", "reflecting", "reflect")
-      bnd_type(1,2) = bnd_reflective
-    case("hydrostatic", "gravity")
-      bnd_type(1,2) = bnd_gravity
-    case("user", "custom")
-      bnd_type(1,2) = bnd_user
-    case default
-      bnd_type(1,2) = bnd_periodic
-    end select
-
-    select case(ylbndry)
-    case("open")
-      bnd_type(2,1) = bnd_open
-    case("outflow", "out")
-      bnd_type(2,1) = bnd_outflow
-    case("reflective", "reflecting", "reflect")
-      bnd_type(2,1) = bnd_reflective
-    case("hydrostatic", "gravity")
-      bnd_type(2,1) = bnd_gravity
-    case("user", "custom")
-      bnd_type(2,1) = bnd_user
-    case default
-      bnd_type(2,1) = bnd_periodic
-    end select
-
-    select case(yubndry)
-    case("open")
-      bnd_type(2,2) = bnd_open
-    case("outflow", "out")
-      bnd_type(2,2) = bnd_outflow
-    case("reflective", "reflecting", "reflect")
-      bnd_type(2,2) = bnd_reflective
-    case("hydrostatic", "gravity")
-      bnd_type(2,2) = bnd_gravity
-    case("user", "custom")
-      bnd_type(2,2) = bnd_user
-    case default
-      bnd_type(2,2) = bnd_periodic
-    end select
-
-    select case(zlbndry)
-    case("open")
-      bnd_type(3,1) = bnd_open
-    case("outflow", "out")
-      bnd_type(3,1) = bnd_outflow
-    case("reflective", "reflecting", "reflect")
-      bnd_type(3,1) = bnd_reflective
-    case("hydrostatic", "gravity")
-      bnd_type(3,1) = bnd_gravity
-    case("user", "custom")
-      bnd_type(3,1) = bnd_user
-    case default
-      bnd_type(3,1) = bnd_periodic
-    end select
-
-    select case(zubndry)
-    case("open")
-      bnd_type(3,2) = bnd_open
-    case("outflow", "out")
-      bnd_type(3,2) = bnd_outflow
-    case("reflective", "reflecting", "reflect")
-      bnd_type(3,2) = bnd_reflective
-    case("hydrostatic", "gravity")
-      bnd_type(3,2) = bnd_gravity
-    case("user", "custom")
-      bnd_type(3,2) = bnd_user
-    case default
-      bnd_type(3,2) = bnd_periodic
-    end select
-
-! set domain periodicity
-!
     do n = 1, NDIMS
-      periodic(n) = (bnd_type(n,1) == bnd_periodic) .and.                      &
+      do s = 1, 2
+
+        bnd = "periodic"
+
+        if (s == 1) then
+          str = char(119+n) // "lbndry"
+        else
+          str = char(119+n) // "ubndry"
+        end if
+        call get_parameter(str, bnd)
+
+        if (s == 1) then
+          str = char(119+n) // "_lower_boundary"
+        else
+          str = char(119+n) // "_upper_boundary"
+        end if
+        call get_parameter(str, bnd)
+
+        select case(bnd)
+        case("open")
+          bnd_type(n,s) = bnd_open
+        case("outflow", "out")
+          bnd_type(n,s) = bnd_outflow
+        case("reflective", "reflecting", "reflect")
+          bnd_type(n,s) = bnd_reflective
+        case("hydrostatic", "gravity")
+          bnd_type(n,s) = bnd_gravity
+        case("user", "custom")
+          bnd_type(n,s) = bnd_custom
+        case default
+          bnd_type(n,s) = bnd_periodic
+        end select
+
+      end do
+
+      if ((bnd_type(n,1) == bnd_periodic .and. &
+           bnd_type(n,2) /= bnd_periodic) .or. &
+          (bnd_type(n,2) == bnd_periodic .and. &
+           bnd_type(n,1) /= bnd_periodic)) then
+
+        call print_message(loc, char(87+n) // &
+                                "-boundary cannot be periodic on one " // &
+                                "side and non-periodic on another!")
+        status = 1
+        return
+      end if
+
+      periodic(n) = (bnd_type(n,1) == bnd_periodic) .and. &
                     (bnd_type(n,2) == bnd_periodic)
+
     end do
 
-#ifdef MPI
-! allocate the exchange arrays
-!
-    allocate(barray(0:npmax,0:npmax), bcount(0:npmax,0:npmax), stat = status)
 
-! prepare the exchange arrays
+#ifdef MPI
+! allocate the arrays for the list of data blocks at different processes which
+! need to exchange boundaries, and for the count of the blocks in each list
+!
+    allocate(barray(0:npmax,0:npmax), bcount(0:npmax,0:npmax), stat=status)
+
+! generate the list of blocks on different processes
 !
     if (status == 0) call prepare_exchange_array()
 #endif /* MPI */
@@ -5377,7 +5312,7 @@ module boundaries
 
 ! user specific boundary conditions
 !
-      case(bnd_user)
+      case(bnd_custom)
 
         if (associated(custom_boundary_x)) then
           call custom_boundary_x(side(1), jl, ju, kl, ku, &
@@ -5560,7 +5495,7 @@ module boundaries
 
 ! user specific boundary conditions
 !
-      case(bnd_user)
+      case(bnd_custom)
 
         if (associated(custom_boundary_y)) then
           call custom_boundary_y(side(2), il, iu, kl, ku, &
@@ -5739,7 +5674,7 @@ module boundaries
 
 ! user specific boundary conditions
 !
-      case(bnd_user)
+      case(bnd_custom)
 
         if (associated(custom_boundary_z)) then
           call custom_boundary_z(side(3), il, iu, jl, ju, &