diff --git a/sources/integrals.F90 b/sources/integrals.F90
index 47db729..7bf6bd6 100644
--- a/sources/integrals.F90
+++ b/sources/integrals.F90
@@ -331,13 +331,14 @@ module integrals
! write the integral file header
!
- write(runit,'("#",a8,17a25)') &
+ write(runit,'("#",a8,20a25)') &
'step', 'time', '|Bx| int', '|Bx| inf' &
, '|Bx| y-adv', '|Bx| y-shr', '|Bx| y-dif' &
, '|Bx| z-adv', '|Bx| z-shr', '|Bx| z-dif' &
, 'Vin lower' , 'Vin upper' &
, 'Emag', 'Emag(inf)' &
- , 'Emag(x-adv)', 'Emag(y-adv)', 'Emag(z-adv)'
+ , 'Emag(x-adv)', 'Emag(y-adv)', 'Emag(z-adv)' &
+ , 'Emag(x-dif)', 'Emag(y-dif)', 'Emag(z-dif)'
write(runit,"('#')")
end if ! store
@@ -426,10 +427,10 @@ module integrals
use coordinates , only : nb, nbl, nbu, ne, nel, neu
use coordinates , only : adx, ady, adz, advol, voli
use coordinates , only : periodic
- use coordinates , only : xmin, xmax, xarea
+ use coordinates , only : xmin, xmax
use coordinates , only : ymin, ymax, yarea
#if NDIMS == 3
- use coordinates , only : zmin, zmax, zarea
+ use coordinates , only : zmin, zmax
#endif /* NDIMS == 3 */
use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp
use equations , only : ien, imx, imy, imz
@@ -523,11 +524,11 @@ module integrals
!
dvol = advol(pdata%meta%level)
dvolh = 0.5d+00 * dvol
- dyz = ady(pdata%meta%level) * adz(pdata%meta%level) / xarea
+ dyz = ady(pdata%meta%level) * adz(pdata%meta%level)
#if NDIMS == 3
- dxy = adx(pdata%meta%level) * ady(pdata%meta%level) / zarea
+ dxy = adx(pdata%meta%level) * ady(pdata%meta%level)
#endif /* NDIMS == 3 */
- dxz = adx(pdata%meta%level) * adz(pdata%meta%level) / yarea
+ dxz = adx(pdata%meta%level) * adz(pdata%meta%level)
dh(1) = adx(pdata%meta%level)
dh(2) = ady(pdata%meta%level)
dh(3) = adz(pdata%meta%level)
@@ -693,6 +694,158 @@ module integrals
!
inarr(21) = inarr(7)
+! fluxes across the X boundaries
+!
+ if (.not. periodic(1)) then
+
+! lower X-boundary
+!
+ if (pdata%meta%xmin <= xmin + dh(1)) then
+
+! advection of magnetic energy through the lower X boundary
+!
+#if NDIMS == 3
+ inarr(23) = inarr(23) + sum((pdata%q(iby,nbl:nb,nb:ne,nb:ne)**2 &
+ + pdata%q(ibz,nbl:nb,nb:ne,nb:ne)**2) &
+ * pdata%q(ivx,nbl:nb,nb:ne,nb:ne) &
+ - (pdata%q(ivy,nbl:nb,nb:ne,nb:ne) &
+ * pdata%q(iby,nbl:nb,nb:ne,nb:ne) &
+ + pdata%q(ivz,nbl:nb,nb:ne,nb:ne) &
+ * pdata%q(ibz,nbl:nb,nb:ne,nb:ne)) &
+ * pdata%q(ibx,nbl:nb,nb:ne,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+ inarr(23) = inarr(23) + sum((pdata%q(iby,nbl:nb,nb:ne, : )**2 &
+ + pdata%q(ibz,nbl:nb,nb:ne, : )**2) &
+ * pdata%q(ivx,nbl:nb,nb:ne, : ) &
+ - (pdata%q(ivy,nbl:nb,nb:ne, : ) &
+ * pdata%q(iby,nbl:nb,nb:ne, : ) &
+ + pdata%q(ivz,nbl:nb,nb:ne, : ) &
+ * pdata%q(ibz,nbl:nb,nb:ne, : )) &
+ * pdata%q(ibx,nbl:nb,nb:ne, : )) * dyz
+#endif /* NDIMS == 3 */
+
+ if (resistivity > 0.0d+00) then
+
+#if NDIMS == 3
+! the Y component of current density
+!
+ tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,nb ,nb:ne,nbu:neu) &
+ - pdata%q(ibx,nb ,nb:ne,nbl:nel)) / dh(3) &
+ - (pdata%q(ibz,nb ,nb:ne,nb :ne ) &
+ - pdata%q(ibz,nbl,nb:ne,nb :ne )) / dh(1)
+
+! the Z component of current density
+!
+ tmp(2,:,:) = (pdata%q(iby,nb ,nb :ne ,nb:ne) &
+ - pdata%q(iby,nbl,nb :ne ,nb:ne)) / dh(1) &
+ - 0.5d+00 * (pdata%q(ibx,nb ,nbu:neu,nb:ne) &
+ - pdata%q(ibx,nb ,nbl:nel,nb:ne)) / dh(2)
+
+! diffusion of magnetic energy through the lower X boundary
+!
+ inarr(26) = inarr(26) &
+ + sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne,nb:ne) &
+ - tmp(2,:,:) * pdata%q(iby,nb,nb:ne,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+! the Y component of current density
+!
+ tmp(1,:,:) = (pdata%q(ibz,nbl,nb:ne, : ) &
+ - pdata%q(ibz,nb ,nb:ne, : )) / dh(1)
+
+! the Z component of current density
+!
+ tmp(2,:,:) = (pdata%q(iby,nb ,nb :ne , : ) &
+ - pdata%q(iby,nbl,nb :ne , : )) / dh(1) &
+ - 0.5d+00 * (pdata%q(ibx,nb ,nbu:neu, : ) &
+ - pdata%q(ibx,nb ,nbl:nel, : )) / dh(2)
+
+! diffusion of magnetic energy through the lower X boundary
+!
+ inarr(26) = inarr(26) &
+ + sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne, : ) &
+ - tmp(2,:,:) * pdata%q(iby,nb,nb:ne, : )) * dyz
+#endif /* NDIMS == 3 */
+
+ end if ! resistivity
+
+ end if ! xmin
+
+! upper X-boundary
+!
+ if (pdata%meta%xmax >= xmax - dh(1)) then
+
+! advection of magnetic energy through the upper X boundary
+!
+#if NDIMS == 3
+ inarr(23) = inarr(23) + sum((pdata%q(iby,ne:neu,nb:ne,nb:ne)**2 &
+ + pdata%q(ibz,ne:neu,nb:ne,nb:ne)**2) &
+ * pdata%q(ivx,ne:neu,nb:ne,nb:ne) &
+ - (pdata%q(ivy,ne:neu,nb:ne,nb:ne) &
+ * pdata%q(iby,ne:neu,nb:ne,nb:ne) &
+ + pdata%q(ivz,ne:neu,nb:ne,nb:ne) &
+ * pdata%q(ibz,ne:neu,nb:ne,nb:ne)) &
+ * pdata%q(ibx,ne:neu,nb:ne,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+ inarr(23) = inarr(23) + sum((pdata%q(iby,ne:neu,nb:ne, : )**2 &
+ + pdata%q(ibz,ne:neu,nb:ne, : )**2) &
+ * pdata%q(ivx,ne:neu,nb:ne, : ) &
+ - (pdata%q(ivy,ne:neu,nb:ne, : ) &
+ * pdata%q(iby,ne:neu,nb:ne, : ) &
+ + pdata%q(ivz,ne:neu,nb:ne, : ) &
+ * pdata%q(ibz,ne:neu,nb:ne, : )) &
+ * pdata%q(ibx,ne:neu,nb:ne, : )) * dyz
+#endif /* NDIMS == 3 */
+
+ if (resistivity > 0.0d+00) then
+
+#if NDIMS == 3
+! the Y component of current density
+!
+ tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,ne ,nb:ne,nbu:neu) &
+ - pdata%q(ibx,ne ,nb:ne,nbl:nel)) / dh(3) &
+ - (pdata%q(ibz,neu,nb:ne,nb :ne ) &
+ - pdata%q(ibz,ne ,nb:ne,nb :ne )) / dh(1)
+
+! the Z component of current density
+!
+ tmp(2,:,:) = (pdata%q(iby,neu,nb :ne ,nb:ne) &
+ - pdata%q(iby,ne ,nb :ne ,nb:ne)) / dh(1) &
+ - 0.5d+00 * (pdata%q(ibx,ne ,nbu:neu,nb:ne) &
+ - pdata%q(ibx,ne ,nbl:nel,nb:ne)) / dh(2)
+
+! diffusion of magnetic energy through the upper X boundary
+!
+ inarr(26) = inarr(26) &
+ - sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne,nb:ne) &
+ - tmp(2,:,:) * pdata%q(iby,ne,nb:ne,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+! the Y component of current density
+!
+ tmp(1,:,:) = (pdata%q(ibz,neu,nb:ne, : ) &
+ - pdata%q(ibz,ne ,nb:ne, : )) / dh(1)
+
+! the Z component of current density
+!
+ tmp(2,:,:) = (pdata%q(iby,neu,nb :ne , : ) &
+ - pdata%q(iby,ne ,nb :ne , : )) / dh(1) &
+ - 0.5d+00 * (pdata%q(ibx,ne ,nbu:neu, : ) &
+ - pdata%q(ibx,ne ,nbl:nel, : )) / dh(2)
+
+! diffusion of magnetic energy through the upper X boundary
+!
+ inarr(26) = inarr(26) &
+ - sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne, : ) &
+ - tmp(2,:,:) * pdata%q(iby,ne,nb:ne, : )) * dyz
+#endif /* NDIMS == 3 */
+
+ end if ! resistivity
+
+ end if ! xmax
+
+ end if ! not periodic in X
+
+! fluxes across the Y boundaries
+!
! lower Y-boundary
!
if (pdata%meta%ymin <= ymin + dh(2)) then
@@ -733,60 +886,86 @@ module integrals
pdata%q(ibx,nb:ne,nbl:nb, : ))) * dxz
#endif /* NDIMS == 3 */
-! diffusion along Y
+! mean magnetic energy at the lower Y boundary
!
+#if NDIMS == 3
+ inarr(22) = inarr(22) + sum(pdata%q(ibx:ibz,nb:ne,nb,nb:ne)**2) * dxz
+#else /* NDIMS == 3 */
+ inarr(22) = inarr(22) + sum(pdata%q(ibx:ibz,nb:ne,nb, : )**2) * dxz
+#endif /* NDIMS == 3 */
+
+! advection of magnetic energy through the lower Y boundary
+!
+#if NDIMS == 3
+ inarr(24) = inarr(24) + sum((pdata%q(ibx,nb:ne,nbl:nb,nb:ne)**2 &
+ + pdata%q(ibz,nb:ne,nbl:nb,nb:ne)**2) &
+ * pdata%q(ivy,nb:ne,nbl:nb,nb:ne) &
+ - (pdata%q(ivx,nb:ne,nbl:nb,nb:ne) &
+ * pdata%q(ibx,nb:ne,nbl:nb,nb:ne) &
+ + pdata%q(ivz,nb:ne,nbl:nb,nb:ne) &
+ * pdata%q(ibz,nb:ne,nbl:nb,nb:ne)) &
+ * pdata%q(iby,nb:ne,nbl:nb,nb:ne)) * dxz
+#else /* NDIMS == 3 */
+ inarr(24) = inarr(24) + sum((pdata%q(ibx,nb:ne,nbl:nb, : )**2 &
+ + pdata%q(ibz,nb:ne,nbl:nb, : )**2) &
+ * pdata%q(ivy,nb:ne,nbl:nb, : ) &
+ - (pdata%q(ivx,nb:ne,nbl:nb, : ) &
+ * pdata%q(ibx,nb:ne,nbl:nb, : ) &
+ + pdata%q(ivz,nb:ne,nbl:nb, : ) &
+ * pdata%q(ibz,nb:ne,nbl:nb, : )) &
+ * pdata%q(iby,nb:ne,nbl:nb, : )) * dxz
+#endif /* NDIMS == 3 */
+
if (resistivity > 0.0d+00) then
#if NDIMS == 3
-! current density Z component
+! the Z component of current density
!
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb ,nb:ne) &
- pdata%q(iby,nbl:nel,nb ,nb:ne)) / dh(1) &
- (pdata%q(ibx,nb :ne ,nb ,nb:ne) &
- pdata%q(ibx,nb :ne ,nbl,nb:ne)) / dh(2)
+! the X component of current density
+!
+ tmp(:,2,:) = (pdata%q(ibz,nb:ne,nb ,nb :ne ) &
+ - pdata%q(ibz,nb:ne,nbl,nb :ne )) / dh(2) &
+ - 0.5d+00 * (pdata%q(iby,nb:ne,nb ,nbu:neu) &
+ - pdata%q(iby,nb:ne,nb ,nbl:nel)) / dh(3)
+
inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), &
- pdata%q(ibx,nb:ne,nb ,nb:ne))) * dxz
+ pdata%q(ibx,nb:ne,nb,nb:ne))) * dxz
+
+! diffusion of magnetic energy through the lower Y boundary
+!
+ inarr(27) = inarr(27) &
+ + sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb,nb:ne) &
+ - tmp(:,2,:) * pdata%q(ibz,nb:ne,nb,nb:ne)) * dxz
#else /* NDIMS == 3 */
-! current density Z component
+! the Z component of current density
!
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb , : ) &
- pdata%q(iby,nbl:nel,nb , : )) / dh(1) &
- (pdata%q(ibx,nb :ne ,nb , : ) &
- pdata%q(ibx,nb :ne ,nbl, : )) / dh(2)
+! the X component of current density
+!
+ tmp(:,2,:) = (pdata%q(ibz,nb:ne,nb , : ) &
+ - pdata%q(ibz,nb:ne,nbl, : )) / dh(2)
+
inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), &
pdata%q(ibx,nb:ne,nb , : ))) * dxz
+
+! diffusion of magnetic energy through the lower Y boundary
+!
+ inarr(27) = inarr(27) &
+ + sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb, : ) &
+ - tmp(:,2,:) * pdata%q(ibz,nb:ne,nb, : )) * dxz
#endif /* NDIMS == 3 */
end if ! resistivity
-! mean magnetic energy at boundary
-!
-#if NDIMS == 3
- inarr(22) = inarr(22) &
- + 2.5d-01 * sum((pdata%q(ibx,nb:ne,nb,nb:ne)**2 &
- + pdata%q(ibz,nb:ne,nb,nb:ne)**2)) * dxz
-#else /* NDIMS == 3 */
- inarr(22) = inarr(22) &
- + 2.5d-01 * sum((pdata%q(ibx,nb:ne,nb, : )**2 &
- + pdata%q(ibz,nb:ne,nb, : )**2)) * dxz
-#endif /* NDIMS == 3 */
-
-! advection of magnetic field along Y
-!
-#if NDIMS == 3
- inarr(24) = inarr(24) &
- + 5.0d-01 * sum((pdata%q(ibx,nb:ne,nbl:nb,nb:ne)**2 &
- + pdata%q(ibz,nb:ne,nbl:nb,nb:ne)**2) &
- * pdata%q(ivy,nb:ne,nbl:nb,nb:ne)) * dxz
-#else /* NDIMS == 3 */
- inarr(24) = inarr(24) &
- + 5.0d-01 * sum((pdata%q(ibx,nb:ne,nbl:nb, : )**2 &
- + pdata%q(ibz,nb:ne,nbl:nb, : )**2) &
- * pdata%q(ivy,nb:ne,nbl:nb, : )) * dxz
-#endif /* NDIMS == 3 */
-
end if ! ymin
! upper Y-boundary
@@ -829,107 +1008,91 @@ module integrals
pdata%q(ibx,nb:ne,ne:neu, : ))) * dxz
#endif /* NDIMS == 3 */
-! diffusion along Y
+! mean magnetic energy at the upper Y boundary
!
+#if NDIMS == 3
+ inarr(22) = inarr(22) + sum(pdata%q(ibx:ibz,nb:ne,ne,nb:ne)**2) * dxz
+#else /* NDIMS == 3 */
+ inarr(22) = inarr(22) + sum(pdata%q(ibx:ibz,nb:ne,ne, : )**2) * dxz
+#endif /* NDIMS == 3 */
+
+! advection of magnetic energy through the upper Y boundary
+!
+#if NDIMS == 3
+ inarr(24) = inarr(24) - sum((pdata%q(ibx,nb:ne,ne:neu,nb:ne)**2 &
+ + pdata%q(ibz,nb:ne,ne:neu,nb:ne)**2) &
+ * pdata%q(ivy,nb:ne,ne:neu,nb:ne) &
+ - (pdata%q(ivx,nb:ne,ne:neu,nb:ne) &
+ * pdata%q(ibx,nb:ne,ne:neu,nb:ne) &
+ + pdata%q(ivz,nb:ne,ne:neu,nb:ne) &
+ * pdata%q(ibz,nb:ne,ne:neu,nb:ne)) &
+ * pdata%q(iby,nb:ne,ne:neu,nb:ne)) * dxz
+#else /* NDIMS == 3 */
+ inarr(24) = inarr(24) - sum((pdata%q(ibx,nb:ne,ne:neu, : )**2 &
+ + pdata%q(ibz,nb:ne,ne:neu, : )**2) &
+ * pdata%q(ivy,nb:ne,ne:neu, : ) &
+ - (pdata%q(ivx,nb:ne,ne:neu, : ) &
+ * pdata%q(ibx,nb:ne,ne:neu, : ) &
+ + pdata%q(ivz,nb:ne,ne:neu, : ) &
+ * pdata%q(ibz,nb:ne,ne:neu, : )) &
+ * pdata%q(iby,nb:ne,ne:neu, : )) * dxz
+#endif /* NDIMS == 3 */
+
if (resistivity > 0.0d+00) then
#if NDIMS == 3
-! current density Z component
+! the Z component of current density
!
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,ne ,nb:ne) &
- pdata%q(iby,nbl:nel,ne ,nb:ne)) / dh(1) &
- (pdata%q(ibx,nb :ne ,neu,nb:ne) &
- pdata%q(ibx,nb :ne ,ne ,nb:ne)) / dh(2)
+! the X component of current density
+
+ tmp(:,2,:) = 0.5d+00 * (pdata%q(ibx,nb :ne ,ne,nbu:neu) &
+ - pdata%q(ibx,nb :ne ,ne,nbl:nel)) / dh(3) &
+ - 0.5d+00 * (pdata%q(ibz,nbu:neu,ne,nb :ne ) &
+ - pdata%q(ibz,nbl:nel,ne,nb :ne )) / dh(1)
+
inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), &
pdata%q(ibx,nb:ne,ne,nb:ne))) * dxz
+
+! diffusion of magnetic energy through the upper Y boundary
+!
+ inarr(27) = inarr(27) &
+ - sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne,nb:ne) &
+ - tmp(:,2,:) * pdata%q(ibz,nb:ne,ne,nb:ne)) * dxz
#else /* NDIMS == 3 */
-! current density Z component
+! the Z component of current density
!
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,ne , : ) &
- pdata%q(iby,nbl:nel,ne , : )) / dh(1) &
- (pdata%q(ibx,nb :ne ,neu, : ) &
- pdata%q(ibx,nb :ne ,ne , : )) / dh(2)
+! the X component of current density
+!
+ tmp(:,2,:) = 0.5d+00 * (pdata%q(ibz,nbl:nel,ne, : ) &
+ - pdata%q(ibz,nbu:neu,ne, : )) / dh(1)
+
inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), &
pdata%q(ibx,nb:ne,ne, : ))) * dxz
+
+! diffusion of magnetic energy through the upper Y boundary
+!
+ inarr(27) = inarr(27) &
+ - sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne, : ) &
+ - tmp(:,2,:) * pdata%q(ibz,nb:ne,ne, : )) * dxz
#endif /* NDIMS == 3 */
end if ! resistivity
-! mean magnetic energy at boundary
-!
-#if NDIMS == 3
- inarr(22) = inarr(22) &
- + 2.5d-01 * sum((pdata%q(ibx,nb:ne,ne,nb:ne)**2 &
- + pdata%q(ibz,nb:ne,ne,nb:ne)**2)) * dxz
-#else /* NDIMS == 3 */
- inarr(22) = inarr(22) &
- + 2.5d-01 * sum((pdata%q(ibx,nb:ne,ne, : )**2 &
- + pdata%q(ibz,nb:ne,ne, : )**2)) * dxz
-#endif /* NDIMS == 3 */
-
-! advection of magnetic field along Y
-!
-#if NDIMS == 3
- inarr(24) = inarr(24) &
- - 5.0d-01 * sum((pdata%q(ibx,nb:ne,ne:neu,nb:ne)**2 &
- + pdata%q(ibz,nb:ne,ne:neu,nb:ne)**2) &
- * pdata%q(ivy,nb:ne,ne:neu,nb:ne)) * dxz
-#else /* NDIMS == 3 */
- inarr(24) = inarr(24) &
- - 5.0d-01 * sum((pdata%q(ibx,nb:ne,ne:neu, : )**2 &
- + pdata%q(ibz,nb:ne,ne:neu, : )**2) &
- * pdata%q(ivy,nb:ne,ne:neu, : )) * dxz
-#endif /* NDIMS == 3 */
-
end if ! ymax
- if (.not. periodic(1)) then
-
-! lower X-boundary
-!
- if (pdata%meta%xmin <= xmin + dh(1)) then
-
-! advection of magnetic field along X
-!
#if NDIMS == 3
- inarr(23) = inarr(23) &
- + 5.0d-01 * sum((pdata%q(iby,nbl:nb,nb:ne,nb:ne)**2 &
- + pdata%q(ibz,nbl:nb,nb:ne,nb:ne)**2) &
- * pdata%q(ivx,nbl:nb,nb:ne,nb:ne)) * dyz
-#else /* NDIMS == 3 */
- inarr(23) = inarr(23) &
- + 5.0d-01 * sum((pdata%q(iby,nbl:nb,nb:ne, : )**2 &
- + pdata%q(ibz,nbl:nb,nb:ne, : )**2) &
- * pdata%q(ivx,nbl:nb,nb:ne, : )) * dyz
-#endif /* NDIMS == 3 */
-
- end if ! xmin
-
-! upper X-boundary
+! fluxes across the Z boundaries
!
- if (pdata%meta%xmax >= xmax - dh(1)) then
-
-! advection of magnetic field along X
-!
-#if NDIMS == 3
- inarr(23) = inarr(23) &
- - 5.0d-01 * sum((pdata%q(iby,ne:neu,nb:ne,nb:ne)**2 &
- + pdata%q(ibz,ne:neu,nb:ne,nb:ne)**2) &
- * pdata%q(ivx,ne:neu,nb:ne,nb:ne)) * dyz
-#else /* NDIMS == 3 */
- inarr(23) = inarr(23) &
- - 5.0d-01 * sum((pdata%q(iby,ne:neu,nb:ne, : )**2 &
- + pdata%q(ibz,ne:neu,nb:ne, : )**2) &
- * pdata%q(ivx,ne:neu,nb:ne, : )) * dyz
-#endif /* NDIMS == 3 */
-
- end if ! xmax
-
- end if ! not periodic in X
-
-#if NDIMS == 3
if (.not. periodic(3)) then
! lower Z-boundary
@@ -949,28 +1112,43 @@ module integrals
* pdata%q(ivx,nb:ne,nb:ne,nbl:nb), &
pdata%q(ibx,nb:ne,nb:ne,nbl:nb))) * dxy
-! diffusion along Z
+! advection of magnetic energy through the lower Z boundary
!
+ inarr(25) = inarr(25) + sum((pdata%q(ibx,nb:ne,nb:ne,nbl:nb)**2 &
+ + pdata%q(iby,nb:ne,nb:ne,nbl:nb)**2) &
+ * pdata%q(ivz,nb:ne,nb:ne,nbl:nb) &
+ - (pdata%q(ivx,nb:ne,nb:ne,nbl:nb) &
+ * pdata%q(ibx,nb:ne,nb:ne,nbl:nb) &
+ + pdata%q(ivy,nb:ne,nb:ne,nbl:nb) &
+ * pdata%q(iby,nb:ne,nb:ne,nbl:nb)) &
+ * pdata%q(ibz,nb:ne,nb:ne,nbl:nb)) * dxy
+
if (resistivity > 0.0d+00) then
-! current density Y component
+! the X component of current density
!
- tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,nb ) &
- - pdata%q(ibz,nbl:nel,nb:ne,nb )) / dh(1) &
- - (pdata%q(ibx,nb :ne ,nb:ne,nb ) &
- - pdata%q(ibx,nb :ne ,nb:ne,nbl)) / dh(3)
+ tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbu:neu,nb ) &
+ - pdata%q(ibz,nb:ne,nbl:nel,nb )) / dh(2) &
+ - (pdata%q(iby,nb:ne,nb :ne ,nb ) &
+ - pdata%q(iby,nb:ne,nb :ne ,nbl)) / dh(3)
- inarr(18) = inarr(18) + resistivity * sum(sign(tmp(:,:,1), &
+! the Y component of current density
+!
+ tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,nb ) &
+ - pdata%q(ibx,nb :ne ,nb:ne,nbl)) / dh(3) &
+ - 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,nb ) &
+ - pdata%q(ibz,nbl:nel,nb:ne,nb )) / dh(1)
+
+ inarr(18) = inarr(18) - resistivity * sum(sign(tmp(:,:,2), &
pdata%q(ibx,nb:ne,nb:ne,nb))) * dxy
- end if ! resistivity
-
-! advection of magnetic field along Z
+! diffusion of magnetic energy through the lower Z boundary
!
- inarr(25) = inarr(25) &
- + 0.5d+00 * sum((pdata%q(ibx,nb:ne,nb:ne,nbl:nb)**2 &
- + pdata%q(iby,nb:ne,nb:ne,nbl:nb)**2) &
- * pdata%q(ivz,nb:ne,nb:ne,nbl:nb)) * dxy
+ inarr(28) = inarr(28) &
+ + sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,nb) &
+ - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,nb)) * dxy
+
+ end if ! resistivity
end if ! zmin
@@ -991,28 +1169,43 @@ module integrals
* pdata%q(ivx,nb:ne,nb:ne,ne:neu), &
pdata%q(ibx,nb:ne,nb:ne,ne:neu))) * dxy
-! diffusion along Z
+! advection of magnetic energy through the upper Z boundary
!
+ inarr(25) = inarr(25) + sum((pdata%q(ibx,nb:ne,nb:ne,ne:neu)**2 &
+ + pdata%q(iby,nb:ne,nb:ne,ne:neu)**2) &
+ * pdata%q(ivz,nb:ne,nb:ne,ne:neu) &
+ - (pdata%q(ivx,nb:ne,nb:ne,ne:neu) &
+ * pdata%q(ibx,nb:ne,nb:ne,ne:neu) &
+ + pdata%q(ivy,nb:ne,nb:ne,ne:neu) &
+ * pdata%q(iby,nb:ne,nb:ne,ne:neu)) &
+ * pdata%q(ibz,nb:ne,nb:ne,ne:neu)) * dxy
+
if (resistivity > 0.0d+00) then
-! current density Y component
+! the X component of current density
!
- tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,ne ) &
- - pdata%q(ibz,nbl:nel,nb:ne,ne )) / dh(1) &
- - (pdata%q(ibx,nb :ne ,nb:ne,neu) &
- - pdata%q(ibx,nb :ne ,nb:ne,ne )) / dh(3)
+ tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbu:neu,ne ) &
+ - pdata%q(ibz,nb:ne,nbl:nel,ne )) / dh(2) &
+ - (pdata%q(iby,nb:ne,nb :ne ,neu ) &
+ - pdata%q(iby,nb:ne,nb :ne ,ne )) / dh(3)
- inarr(18) = inarr(18) - resistivity * sum(sign(tmp(:,:,1), &
+! the Y component of current density
+!
+ tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,neu) &
+ - pdata%q(ibx,nb :ne ,nb:ne,ne )) / dh(3) &
+ - 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,ne ) &
+ - pdata%q(ibz,nbl:nel,nb:ne,ne )) / dh(1)
+
+ inarr(18) = inarr(18) + resistivity * sum(sign(tmp(:,:,1), &
pdata%q(ibx,nb:ne,nb:ne,ne))) * dxy
- end if ! resistivity
-
-! advection of magnetic field along Z
+! diffusion of magnetic energy through the upper Z boundary
!
- inarr(25) = inarr(25) &
- - 0.5d+00 * sum((pdata%q(ibx,nb:ne,nb:ne,ne:neu)**2 &
- + pdata%q(iby,nb:ne,nb:ne,ne:neu)**2) &
- * pdata%q(ivz,nb:ne,nb:ne,ne:neu)) * dxy
+ inarr(28) = inarr(28) &
+ - sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,ne) &
+ - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,ne)) * dxy
+
+ end if ! resistivity
end if ! zmax
@@ -1064,6 +1257,13 @@ module integrals
!
avarr(:) = avarr(:) * voli
+! apply factors to the reconnection rate terms
+!
+ inarr(12) = inarr(12) / yarea
+ inarr(22) = 2.5d-01 * inarr(22) / yarea
+ inarr(23:25) = 5.0d-01 * inarr(23:25)
+ inarr(26:28) = resistivity * inarr(26:28)
+
! write down the integrals and statistics to appropriate files
!
if (stored) then
@@ -1077,7 +1277,7 @@ module integrals
, avarr(6), mnarr(6), mxarr(6) &
, avarr(7), mnarr(7), mxarr(7)
write(eunit,efmt) step, time, dtn, dte, maxval(errors(:)), errors(:)
- write(runit,"(i9, 16es25.15e3)") step, time, inarr(11:25)
+ write(runit,"(i9, 19es25.15e3)") step, time, inarr(11:28)
call flush_and_sync(funit)
call flush_and_sync(sunit)