diff --git a/sources/user_problem.F90 b/sources/user_problem.F90
index bfa6e14..d71b7da 100644
--- a/sources/user_problem.F90
+++ b/sources/user_problem.F90
@@ -31,6 +31,8 @@ module user_problem
 
   implicit none
 
+  integer            :: profile = 1
+
   real(kind=8), save :: beta = 2.00d+00
   real(kind=8), save :: zeta = 0.00d+00
   real(kind=8), save :: dens = 1.00d+00
@@ -88,13 +90,21 @@ module user_problem
     logical          , intent(in)  :: verbose
     integer          , intent(out) :: status
 
-    character(len=64) :: append = "off", fname
+    character(len=64) :: tmp, append = "off", fname
     logical           :: flag
 
     character(len=*), parameter :: &
                                loc = 'USER_PROBLEM::initialize_user_problem()'
 !-------------------------------------------------------------------------------
 !
+    call get_parameter("profile", tmp)
+    select case(tmp)
+      case('bhattacharjee', 'sincos', 'trigonometric')
+        profile = 2
+      case default
+        profile = 1
+    end select
+
     call get_parameter("dens", dens)
     call get_parameter("bamp", bamp)
     call get_parameter("bgui", bgui)
@@ -192,6 +202,14 @@ module user_problem
     write(runit,"('#')")
 
     call print_section(verbose, "Parameters (* - set, otherwise calculated)")
+    select case(profile)
+    case(2)
+      call print_parameter(verbose, '(*) flux tubes profile', &
+                                                  'Bhattacharjee et al. (2009)')
+    case default
+      call print_parameter(verbose, '(*) flux tubes profile', &
+                                                  'Kowal et al. (2023)')
+    end select
     call print_parameter(verbose, '(*) beta (plasma-beta)'   , beta)
     call print_parameter(verbose, '(*) zeta'                 , zeta)
     call print_parameter(verbose, '(*) dens'                 , dens)
@@ -260,7 +278,7 @@ module user_problem
 
     use blocks     , only : block_data
     use constants  , only : pi, pi2
-    use coordinates, only : nn => bcells
+    use coordinates, only : nn => bcells, xmin, xmax
     use coordinates, only : ax, ay
     use equations  , only : prim2cons, csnd2
     use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
@@ -270,34 +288,68 @@ module user_problem
     type(block_data), pointer, intent(inout) :: pdata
 
     integer                     :: i, j, k
-    real(kind=8)                :: bx, by, bz
+    real(kind=8)                :: r1, r2, bx, by, bz
     real(kind=8), dimension(nn) :: x, y
     real(kind=8), dimension(nn) :: snx, csx, sny, csy, thy, ch2
+    real(kind=8), dimension(nn,nn) :: az
 
 !-------------------------------------------------------------------------------
 !
     x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
     y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
 
-    snx(:) = sin(pi * x(:))
-    csx(:) = cos(pi * x(:))
-    sny(:) = sin(pi2 * y(:))
-    csy(:) = cos(pi2 * y(:))
-    thy(:) = tanh(y(:) / dlta)
-    ch2(:) = cosh(y(:) / dlta)**2
+    select case(profile)
 
-    do j = 1, nn
-      do i = 1, nn
-        bx = bamp * csx(i) * (csy(j) * thy(j) + sny(j) / (pi2 * dlta * ch2(j)))
-        by = bamp * 5.0d-01 * snx(i) * sny(j) * thy(j)
-        bz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - bx**2 - by**2)))
+    case(2)
+      snx(:) = sin(pi * x(:))
+      csx(:) = cos(pi * x(:))
+      sny(:) = sin(pi2 * y(:))
+      csy(:) = cos(pi2 * y(:))
+      thy(:) = tanh(y(:) / dlta)
+      ch2(:) = cosh(y(:) / dlta)**2
 
-        pdata%q(ibx,i,j,:) = bx
-        pdata%q(iby,i,j,:) = by
-        pdata%q(ibz,i,j,:) = bz
+      do j = 1, nn
+        do i = 1, nn
+          if (abs(x(i)) <= 5.0d-01) then
+            bx = bamp * csx(i) * (csy(j) * thy(j) + sny(j) / (pi2 * dlta * ch2(j)))
+            by = bamp * 5.0d-01 * snx(i) * sny(j) * thy(j)
+          else
+            bx = 0.0d+00
+            by = 0.0d+00
+          end if
+          bz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - bx**2 - by**2)))
+
+          pdata%q(ibx,i,j,:) = bx
+          pdata%q(iby,i,j,:) = by
+          pdata%q(ibz,i,j,:) = bz
+        end do
       end do
-    end do
-    pdata%q(ibp,:,:,:) = 0.0d+00
+      pdata%q(ibp,:,:,:) = 0.0d+00
+
+    case default
+      do j = 1, nn
+        do i = 1, nn
+          r1 = 5.0d-01 - sqrt(x(i)**2 + ((2.0d+00 * y(j) - 5.0d-01))**2)
+          r2 = 5.0d-01 - sqrt(x(i)**2 + ((2.0d+00 * y(j) + 5.0d-01))**2)
+          az(i,j) = 2.5d-01 * (r1 * (tanh(r1 / dlta) + 1.0d+00) &
+                            +  r2 * (tanh(r2 / dlta) + 1.0d+00))
+        end do
+      end do
+
+      do j = 2, nn - 1
+        do i = 2, nn - 1
+          bx = (az(i,j+1) - az(i,j-1)) / (y(j+1) - y(j-1))
+          by = (az(i-1,j) - az(i+1,j)) / (x(i+1) - x(i-1))
+          bz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - bx**2 - by**2)))
+
+          pdata%q(ibx,i,j,:) = bx
+          pdata%q(iby,i,j,:) = by
+          pdata%q(ibz,i,j,:) = bz
+        end do
+      end do
+      pdata%q(ibp,:,:,:) = 0.0d+00
+
+    end select
 
     if (ipr > 0) then
       pdata%q(idn,:,:,:) = dens
@@ -306,6 +358,7 @@ module user_problem
       pdata%q(idn,:,:,:) = (ptot - 5.0d-01 * sum(pdata%q(ibx:ibz,:,:,:)**2,1)) &
                          / csnd2
     end if
+
     pdata%q(ivx,:,:,:) = 0.0d+00
     pdata%q(ivy,:,:,:) = 0.0d+00
     pdata%q(ivz,:,:,:) = 0.0d+00