diff --git a/sources/user_problem.F90 b/sources/user_problem.F90
index 4c493ac..dfe415c 100644
--- a/sources/user_problem.F90
+++ b/sources/user_problem.F90
@@ -45,7 +45,7 @@ module user_problem
   real(kind=8), save :: lund = 1.00d+00
 
   real(kind=8), save :: dlta = 1.00d-16
-  real(kind=8), save :: blim = 1.00d+00
+  real(kind=8), save :: tdec = 1.00d+00
 
   integer     , save :: pert = 0
   integer     , save :: nper = 10
@@ -107,7 +107,8 @@ module user_problem
 #endif /* NDIMS == 3 */
     use constants  , only : pi2
     use coordinates, only : ng => nghosts, ady, xlen, zlen, ymax
-    use equations  , only : magnetized, csnd, csnd2
+    use equations  , only : magnetized, csnd, csnd2, qpbnd
+    use equations  , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr
     use helpers    , only : print_section, print_parameter, print_message
     use mesh       , only : setup_problem
     use parameters , only : get_parameter
@@ -182,7 +183,13 @@ module user_problem
     call get_parameter("bamp", bamp)
     call get_parameter("bgui", bgui)
 
-    call get_parameter("blimit", blim)
+    call get_parameter("boundary_decay_time", tdec)
+    if (tdec <= 0.0d+00) then
+      if (verbose) &
+        call print_message(loc, "Boundary decay time must be larger than zero!")
+      status = 1
+      return
+    end if
 
 ! calculate the maximum magnetic pressure, thermal pressure from the plasma-β
 ! parameters, and the sound speed in the case of isothermal equations of state
@@ -195,7 +202,18 @@ module user_problem
     valf  = sqrt(2.0d+00 * pmag / dens)
     lund  = valf / max(tiny(eta), eta)
     dlta  = lund**(- 1.0d+00 / 3.0d+00)
-    blim  = max(blim, ng * ady(1))
+
+! fill up the boundary values
+!
+    qpbnd(idn,2,:) =   dens
+    qpbnd(ivx,2,:) = 0.0d+00
+    qpbnd(ivy,2,:) = 0.0d+00
+    qpbnd(ivz,2,:) = 0.0d+00
+    qpbnd(ibx,2,:) = [ - bamp, bamp ]
+    qpbnd(iby,2,:) = 0.0d+00
+    qpbnd(ibz,2,:) =   bgui
+    qpbnd(ibp,2,:) = 0.0d+00
+    if (ipr > 0) qpbnd(ipr,2,:) = pres
 
 ! get the geometry parameters
 !
@@ -382,7 +400,6 @@ module user_problem
       call print_parameter(verbose, '( ) S (Lundquist number)', lund)
     end if
     call print_parameter(verbose, '(*) delta (thickness)', dlta)
-    call print_parameter(verbose, '(*) blim'             , blim)
     call print_parameter(verbose, '(*) perturbation', perturbation)
     call print_parameter(verbose, '(*) bper'        , bper)
     call print_parameter(verbose, '(*) vper'        , vper)
@@ -396,6 +413,7 @@ module user_problem
     call print_parameter(verbose, '(*) ycut'        , ycut)
     call print_parameter(verbose, '(*) xdec'        , xdec)
     call print_parameter(verbose, '(*) ydec'        , ydec)
+    call print_parameter(verbose, '(*) boundary decay time', tdec)
     call print_parameter(verbose, '(*) ydistance (Vrec)', ydis)
     if (absorption) then
       call print_parameter(verbose, '(*) tabs (absorption)', tabs)
@@ -1296,8 +1314,7 @@ module user_problem
   subroutine user_boundary_y(js, il, iu, kl, ku, t, dt, x, y, z, qn)
 
     use coordinates, only : nn => bcells, nb, ne, nbl, neu
-    use equations  , only : magnetized, nv
-    use equations  , only : idn, ivy, ipr, ibx, iby, ibz, ibp
+    use equations  , only : magnetized, nv, ivy, qpbnd
 
     implicit none
 
@@ -1306,16 +1323,7 @@ module user_problem
     real(kind=8), dimension(:)      , intent(in)    :: x, y, z
     real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
 
-    integer      :: i, im1, ip1
-    integer      :: j, jm1, jp1, jm2, jp2
-    integer      :: k
-#if NDIMS == 3
-    integer      :: km1, kp1
-#endif /* NDIMS == 3 */
-    real(kind=8) :: dx, dy, dyx
-#if NDIMS == 3
-    real(kind=8) :: dz, dyz
-#endif /* NDIMS == 3 */
+    integer      :: i, j, k, jm, jp
     real(kind=8) :: fl, fr
 
 !-------------------------------------------------------------------------------
@@ -1326,130 +1334,44 @@ module user_problem
       k = 1
 #endif /* NDIMS == 2 */
 
-      dx  = x(2) - x(1)
-      dy  = y(2) - y(1)
-#if NDIMS == 3
-      dz  = z(2) - z(1)
-#endif /* NDIMS == 3 */
-      dyx = dy / dx
-#if NDIMS == 3
-      dyz = dy / dz
-#endif /* NDIMS == 3 */
+      fl = exp(- dt / tdec)
+      fr = 1.0d+00 - fl
 
-! process left and right side boundary separatelly
-!
       if (js == 1) then
 
-! iterate over left-side ghost layers
-!
         do j = nbl, 1, -1
-
-! calculate neighbor cell indices
-!
-          jp1 = min(nn, j + 1)
-          jp2 = min(nn, j + 2)
-
-! calculate variable decay coefficients
-!
-          fr  = (dy * (nb - j - 5.0d-01)) / blim
-          fl  = 1.0d+00 - fr
-
-! iterate over boundary layer
-!
+          jp = j + 1
 #if NDIMS == 3
           do k = kl, ku
-            km1 = max( 1, k - 1)
-            kp1 = min(nn, k + 1)
 #endif /* NDIMS == 3 */
             do i = il, iu
-              im1 = max( 1, i - 1)
-              ip1 = min(nn, i + 1)
-
-! make normal derivatives zero
-!
-              qn(1:nv,i,j,k) = qn(1:nv,i,nb,k)
-
-! decay density and pressure to their limits
-!
-              qn(idn,i,j,k) = fl * qn(idn,i,nb,k) + fr * dens
-              if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,nb,k) + fr * pres
-
-! decay magnetic field to its limit
-!
-              qn(ibx,i,j,k) = fl * qn(ibx,i,nb,k) - fr * bamp
-              qn(ibz,i,j,k) = fl * qn(ibz,i,nb,k) + fr * bgui
-
-! update By from div(B)=0
-!
-              qn(iby,i,j,k) = qn(iby,i,jp2,k)                              &
-                           + (qn(ibx,ip1,jp1,k) - qn(ibx,im1,jp1,k)) * dyx
-#if NDIMS == 3
-              qn(iby,i,j,k) = qn(iby,i,j  ,k)                              &
-                           + (qn(ibz,i,jp1,kp1) - qn(ibz,i,jp1,km1)) * dyz
-#endif /* NDIMS == 3 */
-              qn(ibp,i,j,k) = 0.0d+00
+              qn( : ,i,j,k) = fl * qn( : ,i,jp,k) + fr * qpbnd(:,2,1)
+              qn(ivy,i,j,k) =      qn(ivy,i,jp,k)
             end do ! i = il, iu
 #if NDIMS == 3
           end do ! k = kl, ku
 #endif /* NDIMS == 3 */
         end do ! j = nbl, 1, -1
+
       else ! js = 1
 
-! iterate over right-side ghost layers
-!
         do j = neu, nn
-
-! calculate neighbor cell indices
-!
-          jm1 = max( 1, j - 1)
-          jm2 = max( 1, j - 2)
-
-! calculate variable decay coefficients
-!
-          fr  = (dy * (j - ne - 5.0d-01)) / blim
-          fl  = 1.0d+00 - fr
-
-! iterate over boundary layer
-!
+          jm = j - 1
 #if NDIMS == 3
           do k = kl, ku
-            km1 = max( 1, k - 1)
-            kp1 = min(nn, k + 1)
 #endif /* NDIMS == 3 */
             do i = il, iu
-              im1 = max( 1, i - 1)
-              ip1 = min(nn, i + 1)
-
-! make normal derivatives zero
-!
-              qn(1:nv,i,j,k) = qn(1:nv,i,ne,k)
-
-! decay density and pressure to their limits
-!
-              qn(idn,i,j,k) = fl * qn(idn,i,ne,k) + fr * dens
-              if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,ne,k) + fr * pres
-
-! decay magnetic field to its limit
-!
-              qn(ibx,i,j,k) = fl * qn(ibx,i,ne,k) + fr * bamp
-              qn(ibz,i,j,k) = fl * qn(ibz,i,ne,k) + fr * bgui
-
-! update By from div(B)=0
-!
-              qn(iby,i,j,k) = qn(iby,i,jm2,k)                              &
-                           + (qn(ibx,im1,jm1,k) - qn(ibx,ip1,jm1,k)) * dyx
-#if NDIMS == 3
-              qn(iby,i,j,k) = qn(iby,i,j  ,k)                              &
-                           + (qn(ibz,i,jm1,km1) - qn(ibz,i,jm1,kp1)) * dyz
-#endif /* NDIMS == 3 */
-              qn(ibp,i,j,k) = 0.0d+00
+              qn( : ,i,j,k) = fl * qn( : ,i,jm,k) + fr * qpbnd(:,2,2)
+              qn(ivy,i,j,k) =      qn(ivy,i,jm,k)
             end do ! i = il, iu
 #if NDIMS == 3
           end do ! k = kl, ku
 #endif /* NDIMS == 3 */
         end do ! j = neu, nn
       end if ! js = 1
-    else ! ibx > 0
+
+    else ! magnetized
+
       if (js == 1) then
         do j = nbl, 1, -1
 #if NDIMS == 3
@@ -1471,7 +1393,8 @@ module user_problem
 #endif /* NDIMS == 3 */
         end do ! j = neu, nn
       end if
-    end if ! ibx > 0
+
+    end if ! magnetized
 
 !-------------------------------------------------------------------------------
 !