Merge branch 'reconnection' of /home/.snapshots/current/kowal/Research/software/codes/amun/ into reconnection

src/boundaries.F90

@@ -55,6 +55,7 @@ module boundaries
   integer, parameter            :: bnd_open         = 1
   integer, parameter            :: bnd_outflow      = 2
   integer, parameter            :: bnd_reflective   = 3
+  integer, parameter            :: bnd_reconnection = 4
 
 ! variable to store boundary type flags
 !
@@ -172,6 +173,8 @@ module boundaries
       bnd_type(1,1) = bnd_outflow
     case("reflective", "reflecting", "reflect")
       bnd_type(1,1) = bnd_reflective
+    case("reconnection", "recon", "rec")
+      bnd_type(1,1) = bnd_reconnection
     case default
       bnd_type(1,1) = bnd_periodic
     end select
@@ -183,6 +186,8 @@ module boundaries
       bnd_type(1,2) = bnd_outflow
     case("reflective", "reflecting", "reflect")
       bnd_type(1,2) = bnd_reflective
+    case("reconnection", "recon", "rec")
+      bnd_type(1,2) = bnd_reconnection
     case default
       bnd_type(1,2) = bnd_periodic
     end select
@@ -194,6 +199,8 @@ module boundaries
       bnd_type(2,1) = bnd_outflow
     case("reflective", "reflecting", "reflect")
       bnd_type(2,1) = bnd_reflective
+    case("reconnection", "recon", "rec")
+      bnd_type(2,1) = bnd_reconnection
     case default
       bnd_type(2,1) = bnd_periodic
     end select
@@ -205,6 +212,8 @@ module boundaries
       bnd_type(2,2) = bnd_outflow
     case("reflective", "reflecting", "reflect")
       bnd_type(2,2) = bnd_reflective
+    case("reconnection", "recon", "rec")
+      bnd_type(2,2) = bnd_reconnection
     case default
       bnd_type(2,2) = bnd_periodic
     end select
@@ -1118,6 +1127,7 @@ module boundaries
 !
                   if (.not. associated(pmeta%edges(i,j,m)%ptr))                &
                             call block_boundary_specific(i, j, k, n            &
+                                          , pmeta%level                        &
                                           , pmeta%data%q(1:nv,1:im,1:jm,1:km))
 
                 end do ! i = 1, sides
@@ -1146,6 +1156,7 @@ module boundaries
 !
                     if (.not. associated(pmeta%faces(i,j,k,n)%ptr))            &
                             call block_boundary_specific(i, j, k, n            &
+                                          , pmeta%level                        &
                                           , pmeta%data%q(1:nv,1:im,1:jm,1:km))
 
                   end do ! i = 1, sides
@@ -4944,20 +4955,23 @@ module boundaries
 !
 !     nc         - the edge direction;
 !     ic, jc, kc - the corner position;
+!     lv         - the block level;
 !     qn         - the variable array;
 !
 !===============================================================================
 !
-  subroutine block_boundary_specific(ic, jc, kc, nc, qn)
+  subroutine block_boundary_specific(ic, jc, kc, nc, lv, qn)
 
 ! import external procedures and variables
 !
     use coordinates    , only : im , jm , km , ng
     use coordinates    , only : ib , jb , kb , ie , je , ke
     use coordinates    , only : ibl, jbl, kbl, ieu, jeu, keu
+    use coordinates    , only : adx, ady, adxi, adyi, adzi
     use equations      , only : nv
-    use equations      , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp
+    use equations      , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp
     use error          , only : print_error, print_warning
+    use parameters     , only : get_parameter_real
 
 ! local variables are not implicit by default
 !
@@ -4966,9 +4980,21 @@ module boundaries
 ! subroutine arguments
 !
     integer                                     , intent(in)    :: ic, jc, kc
-    integer                                     , intent(in)    :: nc
+    integer                                     , intent(in)    :: nc, lv
     real(kind=8), dimension(1:nv,1:im,1:jm,1:km), intent(inout) :: qn
 
+! default parameter values
+!
+    real(kind=8), save :: dens = 1.00d+00
+    real(kind=8), save :: pres = 1.00d+00
+    real(kind=8), save :: bamp = 1.00d+00
+    real(kind=8), save :: bgui = 0.00d+00
+    real(kind=8), save :: blim = 1.00d+00
+
+! local saved parameters
+!
+    logical     , save :: first = .true.
+
 ! local variables
 !
     integer :: i , j , k
@@ -4976,9 +5002,38 @@ module boundaries
     integer :: iu, ju, ku
     integer :: is, js, ks
     integer :: it, jt, kt
+    integer :: im2, im1, ip1, ip2
+    integer :: jm2, jm1, jp1, jp2
+#if NDIMS == 3
+    integer :: km2, km1, kp1, kp2
+#endif /* NDIMS == 3 */
+    real(kind=8) :: dxy, dxz, dyx, dyz
+    real(kind=8) :: fl, fr, ds
 !
 !-------------------------------------------------------------------------------
 !
+! prepare problem constants during the first subroutine call
+!
+    if (first) then
+
+! get problem parameters
+!
+      call get_parameter_real("dens"  , dens)
+      call get_parameter_real("pres"  , pres)
+      call get_parameter_real("bamp"  , bamp)
+      call get_parameter_real("bgui"  , bgui)
+      call get_parameter_real("blimit", blim)
+
+! upper limit for blim
+!
+      blim = max(blim, ng * ady(1))
+
+! reset the first execution flag
+!
+      first = .false.
+
+    end if ! first call
+
 ! apply specific boundaries depending on the direction
 !
     select case(nc)
@@ -5040,6 +5095,180 @@ module boundaries
           end do ! i = ieu, im
         end if
 
+! "reconnection" boundary conditions
+!
+      case(bnd_reconnection)
+
+! process case with magnetic field, otherwise revert to standard outflow
+!
+        if (ibx > 0) then
+
+! get the cell size ratios
+!
+          dxy = adx(lv) * adyi(lv)
+          dxz = adx(lv) * adzi(lv)
+
+! process left and right side boundary separatelly
+!
+          if (ic == 1) then
+
+! iterate over left-side ghost layers
+!
+            do i = ibl, 1, -1
+
+! calculate neighbor cell indices
+!
+              ip1 = min(im, i + 1)
+              ip2 = min(im, i + 2)
+
+! iterate over boundary layer
+!
+              do k = kl, ku
+#if NDIMS == 3
+                km2 = max( 1, k - 2)
+                km1 = max( 1, k - 1)
+                kp1 = min(km, k + 1)
+                kp2 = min(km, k + 2)
+#endif /* NDIMS == 3 */
+                do j = jl, ju
+                  jm2 = max( 1, j - 2)
+                  jm1 = max( 1, j - 1)
+                  jp1 = min(jm, j + 1)
+                  jp2 = min(jm, j + 2)
+
+! make normal derivatives zero
+!
+                  qn(1:nv,i,j,k) = qn(1:nv,ib,j,k)
+
+! prevent the inflow
+!
+                  qn(ivx,i,j,k) = min(0.0d+00, qn(ivx,ib,j,k))
+
+! prevent from creating strong reconnection at the boundary (only near
+! the plane of Bx polarity change)
+!
+#if NDIMS == 3
+! detect the current sheet
+!
+                  ds = min(qn(ibx,ib,jm2,k) * qn(ibx,ib,jp2,k)                 &
+                         , qn(ibx,ib,j,km2) * qn(ibx,ib,j,kp2))
+
+! apply curl-free condition in the vicinity of current sheet
+!
+                  if (ds < 0.0d+00) then
+                    qn(iby,i,j,k) = qn(iby,ip2,j,k)                            &
+                                 + (qn(ibx,ip1,jm1,k) - qn(ibx,ip1,jp1,k)) * dxy
+                    qn(ibz,i,j,k) = qn(ibz,ip2,j,k)                            &
+                                 + (qn(ibx,ip1,j,km1) - qn(ibx,ip1,j,kp1)) * dxz
+                  end if
+#else /* NDIMS == 3 */
+! apply curl-free condition in the vicinity of current sheet
+!
+                  if (qn(ibx,ib,jm2,k) * qn(ibx,ib,jp2,k) < 0.0d+00) then
+                    qn(iby,i,j,k) = qn(iby,ip2,j,k)                            &
+                                 + (qn(ibx,ip1,jm1,k) - qn(ibx,ip1,jp1,k)) * dxy
+                  end if
+#endif /* NDIMS == 3 */
+
+! update Bx from div(B)=0
+!
+                  qn(ibx,i,j,k) = qn(ibx,ip2,j,k)                              &
+                               + (qn(iby,ip1,jp1,k) - qn(iby,ip1,jm1,k)) * dxy
+#if NDIMS == 3
+                  qn(ibx,i,j,k) = qn(ibx,i  ,j,k)                              &
+                               + (qn(ibz,ip1,j,kp1) - qn(ibz,ip1,j,km1)) * dxz
+#endif /* NDIMS == 3 */
+                  qn(ibp,i,j,k) = 0.0d+00
+                end do ! j = jl, ju
+              end do ! k = kl, ku
+            end do ! i = ibl, 1, -1
+          else ! ic == 1
+
+! iterate over right-side ghost layers
+!
+            do i = ieu, im
+
+! calculate neighbor cell indices
+!
+              im1 = max( 1, i - 1)
+              im2 = max( 1, i - 2)
+
+! iterate over boundary layer
+!
+              do k = kl, ku
+#if NDIMS == 3
+                km1 = max( 1, k - 1)
+                kp1 = min(km, k + 1)
+                km2 = max( 1, k - 2)
+                kp2 = min(km, k + 2)
+#endif /* NDIMS == 3 */
+                do j = jl, ju
+                  jm1 = max( 1, j - 1)
+                  jp1 = min(jm, j + 1)
+                  jm2 = max( 1, j - 2)
+                  jp2 = min(jm, j + 2)
+
+! make normal derivatives zero
+!
+                  qn(1:nv,i,j,k) = qn(1:nv,ie,j,k)
+
+! prevent the inflow
+!
+                  qn(ivx,i,j,k) = max(0.0d+00, qn(ivx,ie,j,k))
+
+! prevent from creating strong reconnection at the boundary (only near
+! the plane of Bx polarity change)
+!
+#if NDIMS == 3
+! detect the current sheet
+!
+                  ds = min(qn(ibx,ie,jm2,k) * qn(ibx,ie,jp2,k)                 &
+                         , qn(ibx,ie,j,km2) * qn(ibx,ie,j,kp2))
+
+! apply curl-free condition in the vicinity of current sheet
+!
+                  if (ds < 0.0d+00) then
+                    qn(iby,i,j,k) = qn(iby,im2,j,k)                            &
+                                 + (qn(ibx,im1,jp1,k) - qn(ibx,im1,jm1,k)) * dxy
+                    qn(ibz,i,j,k) = qn(ibz,im2,j,k)                            &
+                                 + (qn(ibx,im1,j,kp1) - qn(ibx,im1,j,km1)) * dxz
+                  end if
+#else /* NDIMS == 3 */
+! apply curl-free condition in the vicinity of current sheet
+!
+                  if (qn(ibx,ie,jm2,k) * qn(ibx,ie,jp2,k) < 0.0d+00) then
+                    qn(iby,i,j,k) = qn(iby,im2,j,k)                            &
+                                 + (qn(ibx,im1,jp1,k) - qn(ibx,im1,jm1,k)) * dxy
+                  end if
+#endif /* NDIMS == 3 */
+
+! update Bx from div(B)=0
+!
+                  qn(ibx,i,j,k) = qn(ibx,im2,j,k)                              &
+                               + (qn(iby,im1,jm1,k) - qn(iby,im1,jp1,k)) * dxy
+#if NDIMS == 3
+                  qn(ibx,i,j,k) = qn(ibx,i  ,j,k)                              &
+                               + (qn(ibz,im1,j,km1) - qn(ibz,im1,j,kp1)) * dxz
+#endif /* NDIMS == 3 */
+                  qn(ibp,i,j,k) = 0.0d+00
+                end do ! j = jl, ju
+              end do ! k = kl, ku
+            end do ! i = ieu, im
+          end if ! ic == 1
+        else ! ibx > 0
+          if (ic == 1) then
+            do i = ibl, 1, -1
+              qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,ib,jl:ju,kl:ku)
+              qn(ivx ,i,jl:ju,kl:ku) = min(0.0d+00, qn(ivx,ib,jl:ju,kl:ku))
+            end do ! i = ibl, 1, -1
+          else
+            do i = ieu, im
+              qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,ie,jl:ju,kl:ku)
+              qn(ivx ,i,jl:ju,kl:ku) = max(0.0d+00, qn(ivx,ie,jl:ju,kl:ku))
+            end do ! i = ieu, im
+          end if
+        end if ! ibx > 0
+
 ! "reflective" boundary conditions
 !
       case(bnd_reflective)
@@ -5134,6 +5363,142 @@ module boundaries
           end do ! j = jeu, jm
         end if
 
+! "reconnection" boundary conditions
+!
+      case(bnd_reconnection)
+
+! process case with magnetic field, otherwise revert to standard outflow
+!
+        if (ibx > 0) then
+
+! get the cell size ratios
+!
+          dyx = ady(lv) * adxi(lv)
+          dyz = ady(lv) * adzi(lv)
+
+! process left and right side boundary separatelly
+!
+          if (jc == 1) then
+
+! iterate over left-side ghost layers
+!
+            do j = jbl, 1, -1
+
+! calculate neighbor cell indices
+!
+              jp1 = min(jm, j + 1)
+              jp2 = min(jm, j + 2)
+
+! calculate variable decay coefficients
+!
+              fr  = (ady(lv) * (jb - j - 5.0d-01)) / blim
+              fl  = 1.0d+00 - fr
+
+! iterate over boundary layer
+!
+              do k = kl, ku
+#if NDIMS == 3
+                km1 = max( 1, k - 1)
+                kp1 = min(km, k + 1)
+#endif /* NDIMS == 3 */
+                do i = il, iu
+                  im1 = max( 1, i - 1)
+                  ip1 = min(im, i + 1)
+
+! make normal derivatives zero
+!
+                  qn(1:nv,i,j,k) = qn(1:nv,i,jb,k)
+
+! decay density and pressure to their limits
+!
+                  qn(idn,i,j,k) = fl * qn(idn,i,jb,k) + fr * dens
+                  if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,jb,k) + fr * pres
+
+! decay magnetic field to its limit
+!
+                  qn(ibx,i,j,k) = fl * qn(ibx,i,jb,k) - fr * bamp
+                  qn(ibz,i,j,k) = fl * qn(ibz,i,jb,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
+                end do ! i = il, iu
+              end do ! k = kl, ku
+            end do ! j = jbl, 1, -1
+          else ! jc = 1
+
+! iterate over right-side ghost layers
+!
+            do j = jeu, jm
+
+! calculate neighbor cell indices
+!
+              jm1 = max( 1, j - 1)
+              jm2 = max( 1, j - 2)
+
+! calculate variable decay coefficients
+!
+              fr  = (ady(lv) * (j - je - 5.0d-01)) / blim
+              fl  = 1.0d+00 - fr
+
+! iterate over boundary layer
+!
+              do k = kl, ku
+#if NDIMS == 3
+                km1 = max( 1, k - 1)
+                kp1 = min(km, k + 1)
+#endif /* NDIMS == 3 */
+                do i = il, iu
+                  im1 = max( 1, i - 1)
+                  ip1 = min(im, i + 1)
+
+! make normal derivatives zero
+!
+                  qn(1:nv,i,j,k) = qn(1:nv,i,je,k)
+
+! decay density and pressure to their limits
+!
+                  qn(idn,i,j,k) = fl * qn(idn,i,je,k) + fr * dens
+                  if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,je,k) + fr * pres
+
+! decay magnetic field to its limit
+!
+                  qn(ibx,i,j,k) = fl * qn(ibx,i,je,k) + fr * bamp
+                  qn(ibz,i,j,k) = fl * qn(ibz,i,je,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
+                end do ! i = il, iu
+              end do ! k = kl, ku
+            end do ! j = jeu, jm
+          end if ! jc = 1
+        else ! ibx > 0
+          if (jc == 1) then
+            do j = jbl, 1, -1
+              qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,jb,kl:ku)
+              qn(ivy ,il:iu,j,kl:ku) = min(0.0d+00, qn(ivy,il:iu,jb,kl:ku))
+            end do ! j = jbl, 1, -1
+          else
+            do j = jeu, jm
+              qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,je,kl:ku)
+              qn(ivy ,il:iu,j,kl:ku) = max(0.0d+00, qn(ivy,il:iu,je,kl:ku))
+            end do ! j = jeu, jm
+          end if
+        end if ! ibx > 0
+
 ! "reflective" boundary conditions
 !
       case(bnd_reflective)

src/coordinates.F90

@@ -94,6 +94,10 @@ module coordinates
   real(kind=8), save :: yarea = 1.0d+00
   real(kind=8), save :: zarea = 1.0d+00
 
+! the characteristic decay length
+!
+  real(kind=8), save :: ldec = 1.0d-03
+
 ! the block coordinates for all levels of refinement
 !
   real(kind=8), dimension(:,:), allocatable, save :: ax  , ay  , az
@@ -131,6 +135,10 @@ module coordinates
   type(rectangular), dimension(2,2,2)      , save :: corners_dr, corners_gr
 #endif /* NDIMS == 3 */
 
+! the decay factors for all levels
+!
+  real(kind=8), dimension(:  ), allocatable, save :: adec
+
 ! by default everything is private
 !
   public
@@ -304,6 +312,7 @@ module coordinates
     allocate(adyi (toplev))
     allocate(adzi (toplev))
     allocate(advol(toplev))
+    allocate(adec (toplev))
 
 ! reset all coordinate variables to initial values
 !
@@ -714,6 +723,16 @@ module coordinates
     end do ! k = 1, 2
 #endif /* NDIMS == 3 */
 
+! get the characteristic decay scale
+!
+    call get_parameter_real("ldecay", ldec)
+
+! calculate the decay factors
+!
+    do l = 1, toplev
+      adec(l) = exp(- adx(l) / ldec)
+    end do ! l = 1, toplev
+
 ! print general information about the level resolutions
 !
     if (verbose) then
@@ -794,6 +813,7 @@ module coordinates
     if (allocated(adyi) ) deallocate(adyi)
     if (allocated(adzi) ) deallocate(adzi)
     if (allocated(advol)) deallocate(advol)
+    if (allocated(adec) ) deallocate(adec)
 
 !-------------------------------------------------------------------------------
 !

src/integrals.F90

@@ -55,6 +55,7 @@ module integrals
 !
   integer(kind=4), save :: funit = 7
   integer(kind=4), save :: sunit = 8
+  integer(kind=4), save :: runit = 9
   integer(kind=4), save :: iintd = 1
 
 ! by default everything is private
@@ -193,6 +194,30 @@ module integrals
                                  , 'mean(Malf)', 'min(Malf)', 'max(Malf)'
       write(sunit,"('#')")
 
+! generate the reconnection file name
+!
+      write(fname, "('reconnection_',i2.2,'.dat')") irun
+
+! create a new statistics file
+!
+#ifdef GNU
+      open (newunit = runit, file = fname, form = 'formatted'                  &
+                                       , status = 'replace')
+#endif /* GNU */
+#ifdef INTEL
+      open (newunit = runit, file = fname, form = 'formatted'                  &
+                                       , status = 'replace', buffered = 'yes')
+#endif /* INTEL */
+#ifdef IBM
+      open (unit = runit, file = fname, form = 'formatted'                     &
+                                       , status = 'replace')
+#endif /* IBM */
+
+! write the integral file header
+!
+      write(runit,"('#',a8,3(1x,a18))") 'step', 'time', 'Vrec_l', 'Vrec_u'
+      write(runit,"('#')")
+
     end if ! master
 
 #ifdef PROFILE
@@ -238,6 +263,7 @@ module integrals
     if (master) then
       close(funit)
       close(sunit)
+      close(runit)
     end if
 
 #ifdef PROFILE
@@ -267,7 +293,8 @@ module integrals
 !
     use blocks         , only : block_meta, block_data, list_data
     use coordinates    , only : in, jn, kn, ib, jb, kb, ie, je, ke
-    use coordinates    , only : advol, voli
+    use coordinates    , only : adx, adz, advol, voli
+    use coordinates    , only : ymin, ymax, xlen, zlen, yarea
     use equations      , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp
     use equations      , only : ien, imx, imy, imz
     use equations      , only : gamma, csnd
@@ -286,7 +313,7 @@ module integrals
 ! local variables
 !
     integer                       :: iret
-    real(kind=8)                  :: dvol, dvolh
+    real(kind=8)                  :: dvol, dvolh, dxz
 
 ! local pointers
 !
@@ -352,6 +379,7 @@ module integrals
 !
       dvol  = advol(pdata%meta%level)
       dvolh = 0.5d+00 * dvol
+      dxz   = adx(pdata%meta%level) * adz(pdata%meta%level) / yarea
 
 ! sum up density and momenta components
 !
@@ -444,6 +472,15 @@ module integrals
         mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:)))
       end if
 
+! get the inflow speed
+!
+      if (pdata%meta%ymin == ymin) then
+        inarr(10) = inarr(10) + sum(pdata%q(ivy,ib:ie,jb,kb:ke)) * dxz
+      end if
+      if (pdata%meta%ymax == ymax) then
+        inarr(11) = inarr(11) + sum(pdata%q(ivy,ib:ie,je,kb:ke)) * dxz
+      end if
+
 ! associate the pointer with the next block on the list
 !
       pdata => pdata%next
@@ -482,6 +519,7 @@ module integrals
                                             , avarr(5), mnarr(5), mxarr(5)     &
                                             , avarr(6), mnarr(6), mxarr(6)     &
                                             , avarr(7), mnarr(7), mxarr(7)
+      write(runit,"(i9, 3(1x,1e18.8e3))") step, time, inarr(10:11)
     end if
 
 #ifdef PROFILE

src/makefile

@@ -149,7 +149,7 @@ mpitools.o       : mpitools.F90 error.o timers.o
 operators.o      : operators.F90 timers.o
 parameters.o     : parameters.F90 mpitools.o
 problems.o       : problems.F90 blocks.o constants.o coordinates.o equations.o \
-                   error.o parameters.o random.o timers.o
+                   error.o operators.o parameters.o random.o timers.o
 random.o         : random.F90 mpitools.o parameters.o
 refinement.o     : refinement.F90 blocks.o coordinates.o equations.o           \
                    operators.o parameters.o timers.o

src/problems.F90

@@ -143,6 +143,9 @@ module problems
     case("current_sheet")
       setup_problem => setup_problem_current_sheet
 
+    case("reconnection")
+      setup_problem => setup_problem_reconnection
+
     case("jet")
       setup_problem => setup_problem_jet
 
@@ -1367,6 +1370,245 @@ module problems
 !
 !===============================================================================
 !
+! subroutine SETUP_PROBLEM_RECONNECTION:
+! -------------------------------------
+!
+!   Subroutine sets the initial conditions for the reconnection problem.
+!
+!   Arguments:
+!
+!     pdata - pointer to the datablock structure of the currently initialized
+!             block;
+!
+!===============================================================================
+!
+  subroutine setup_problem_reconnection(pdata)
+
+! include external procedures and variables
+!
+    use blocks     , only : block_data
+    use constants  , only : pi2
+    use coordinates, only : im, jm, km
+    use coordinates, only : ax, ay, adx, ady, adz
+    use equations  , only : prim2cons
+    use equations  , only : nv
+    use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
+    use equations  , only : csnd2
+    use operators  , only : curl
+    use parameters , only : get_parameter_real
+    use random     , only : randomn
+
+! local variables are not implicit by default
+!
+    implicit none
+
+! input arguments
+!
+    type(block_data), pointer, intent(inout) :: pdata
+
+! default parameter values
+!
+    real(kind=8), save :: dens = 1.00d+00
+    real(kind=8), save :: pres = 1.00d+00
+    real(kind=8), save :: bamp = 1.00d+00
+    real(kind=8), save :: bper = 0.00d+00
+    real(kind=8), save :: bgui = 0.00d+00
+    real(kind=8), save :: vper = 1.00d-02
+    real(kind=8), save :: xcut = 4.00d-01
+    real(kind=8), save :: ycut = 1.00d-02
+    real(kind=8), save :: yth  = 1.00d-16
+    real(kind=8), save :: pth  = 1.00d-02
+    real(kind=8), save :: pmag = 5.00d-01
+
+! local saved parameters
+!
+    logical     , save :: first = .true.
+
+! local variables
+!
+    integer            :: i, j, k
+    real(kind=8)       :: yt, yp
+
+! local arrays
+!
+    real(kind=8), dimension(nv,jm) :: q, u
+    real(kind=8), dimension(im)    :: x
+    real(kind=8), dimension(jm)    :: y, pm
+    real(kind=8), dimension(3)     :: dh
+!
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for the problem setup
+!
+    call start_timer(imu)
+#endif /* PROFILE */
+
+! prepare problem constants during the first subroutine call
+!
+    if (first) then
+
+! get problem parameters
+!
+      call get_parameter_real("dens", dens)
+      call get_parameter_real("pres", pres)
+      call get_parameter_real("bamp", bamp)
+      call get_parameter_real("bper", bper)
+      call get_parameter_real("bgui", bgui)
+      call get_parameter_real("vper", vper)
+      call get_parameter_real("xcut", xcut)
+      call get_parameter_real("ycut", ycut)
+      call get_parameter_real("yth" , yth )
+      call get_parameter_real("pth" , pth )
+
+! calculate the maximum magnetic pressure
+!
+      pmag = 0.5d+00 * (bamp**2 + bgui**2)
+
+! reset the first execution flag
+!
+      first = .false.
+
+    end if ! first call
+
+! prepare block coordinates
+!
+    x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im)
+    y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm)
+
+! prepare cell sizes
+!
+    dh(1) = adx(pdata%meta%level)
+    dh(2) = ady(pdata%meta%level)
+    dh(3) = adz(pdata%meta%level)
+
+! calculate the perturbation of magnetic field
+!
+    if (bper /= 0.0d+00) then
+
+! initiate the vector potential (we use velocity components to store vector
+! potential temporarily, and we store magnetic field perturbation in U)
+!
+      do k = 1, km
+        do j = 1, jm
+          do i = 1, im
+            yt = abs(y(j) / yth)
+            yt = log(exp(yt) + exp(- yt))
+            yp = y(j) / pth
+
+            pdata%q(ivx,i,j,k) = 0.0d+00
+            pdata%q(ivy,i,j,k) = 0.0d+00
+            pdata%q(ivz,i,j,k) = bper * cos(pi2 * x(i)) * exp(- yp * yp) / pi2
+          end do ! i = 1, im
+        end do ! j = 1, jm
+      end do ! k = 1, km
+
+! calculate magnetic field components from vector potential
+!
+      call curl(dh(1:3), pdata%q(ivx:ivz,1:im,1:jm,1:km)                       &
+                                            , pdata%q(ibx:ibz,1:im,1:jm,1:km))
+
+    else
+
+! reset magnetic field components
+!
+      pdata%q(ibx:ibz,:,:,:) = 0.0d+00
+
+    end if ! bper /= 0.0
+
+! iterate over all positions in the XZ plane
+!
+    do k = 1, km
+      do i = 1, im
+
+! if magnetic field is present, set its initial configuration
+!
+        if (ibx > 0) then
+
+! set antiparallel magnetic field component
+!
+          do j = 1, jm
+            q(ibx,j) = bamp * tanh(y(j) / yth)
+          end do ! j = 1, jm
+
+! set tangential magnetic field components
+!
+          q(iby,1:jm) = 0.0d+00
+          q(ibz,1:jm) = bgui
+          q(ibp,1:jm) = 0.0d+00
+
+! calculate local magnetic pressure
+!
+          pm(1:jm)    = 0.5d+00 * sum(q(ibx:ibz,:) * q(ibx:ibz,:), 1)
+
+! add magnetic field perturbation
+!
+          if (bper /= 0.0d+00) then
+            q(ibx,1:jm) = q(ibx,1:jm) + pdata%q(ibx,i,1:jm,k)
+            q(iby,1:jm) = q(iby,1:jm) + pdata%q(iby,i,1:jm,k)
+            q(ibz,1:jm) = q(ibz,1:jm) + pdata%q(ibz,i,1:jm,k)
+          end if ! bper /= 0.0
+
+        end if ! ibx > 0
+
+! set the uniform density and pressure
+!
+        if (ipr > 0) then
+          q(idn,1:jm) = dens
+          q(ipr,1:jm) = pres + (pmag - pm(1:jm))
+        else
+          q(idn,1:jm) = dens + (pmag - pm(1:jm)) / csnd2
+        end if
+
+! reset velocity components
+!
+        q(ivx,1:jm) = 0.0d+00
+        q(ivy,1:jm) = 0.0d+00
+        q(ivz,1:jm) = 0.0d+00
+
+! set the random velocity field in a layer near current sheet
+!
+        if (abs(x(i)) <= xcut) then
+          do j = 1, jm
+            if (abs(y(j)) <= ycut) then
+
+              q(ivx,j) = vper * randomn()
+              q(ivy,j) = vper * randomn()
+#if NDIMS == 3
+              q(ivz,j) = vper * randomn()
+#endif /* NDIMS == 3 */
+
+            end if ! |y| < ycut
+          end do ! j = 1, jm
+        end if ! |x| < xcut
+
+! convert the primitive variables to conservative ones
+!
+        call prim2cons(jm, q(1:nv,1:jm), u(1:nv,1:jm))
+
+! copy the conserved variables to the current block
+!
+        pdata%u(1:nv,i,1:jm,k) = u(1:nv,1:jm)
+
+! copy the primitive variables to the current block
+!
+        pdata%q(1:nv,i,1:jm,k) = q(1:nv,1:jm)
+
+      end do ! i = 1, im
+    end do ! k = 1, km
+
+#ifdef PROFILE
+! stop accounting time for the problems setup
+!
+    call stop_timer(imu)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine setup_problem_reconnection
+!
+!===============================================================================
+!
 ! subroutine SETUP_PROBLEM_JET:
 ! ----------------------------
 !