gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

INTEGRALS: Add remaining terms of the magnetic energy evolution.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2020-12-29 17:47:57 -03:00
parent d49b03e7e7
commit 4af2db266c
1 changed files with 275 additions and 195 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

464
sources/integrals.F90
View File

@ -331,13 +331,15 @@ module integrals
! write the integral file header ! write the integral file header
! !
write(runit,'("#",a8,20a25)') & write(runit,'("#",a8,25a25)') &
'step', 'time', '|Bx| int', '|Bx| inf' & 'step', 'time', '|Bx| int', '|Bx| inf' &
, '|Bx| y-adv', '|Bx| z-adv', '|Bx| y-shr', '|Bx| z-shr' & , '|Bx| y-adv', '|Bx| z-adv', '|Bx| y-shr', '|Bx| z-shr' &
, '|Bx| y-dif', '|Bx| z-dif', 'Vin lower' , 'Vin upper' & , '|Bx| y-dif', '|Bx| z-dif', 'Vin lower' , 'Vin upper' &
, 'Emag', 'Emag inf' & , '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' , 'Emag x-v.b', 'Emag y-v.b', 'Emag z-v.b' &
, 'Emag x-dif', 'Emag y-dif', 'Emag z-dif' &
, 'Emag-Ekin', 'Emag-Eint'
write(runit,"('#')") write(runit,"('#')")
end if ! store end if ! store
@ -423,7 +425,7 @@ module integrals
! !
use blocks , only : block_meta, block_data, list_data use blocks , only : block_meta, block_data, list_data
use coordinates , only : ni => ncells use coordinates , only : ni => ncells
use coordinates , only : nb, nbl, nbu, ne, nel, neu use coordinates , only : nb, nbm, nbp, ne, nem, nep
use coordinates , only : adx, ady, adz, advol, voli use coordinates , only : adx, ady, adz, advol, voli
use coordinates , only : periodic use coordinates , only : periodic
use coordinates , only : xmin, xmax use coordinates , only : xmin, xmax
@ -461,7 +463,7 @@ module integrals
! local parameters ! local parameters
! !
integer, parameter :: narr = 32 integer, parameter :: narr = 64
! local arrays ! local arrays
! !
@ -704,23 +706,21 @@ module integrals
! advection of magnetic energy through the lower X boundary ! advection of magnetic energy through the lower X boundary
! !
#if NDIMS == 3 #if NDIMS == 3
inarr(23) = inarr(23) + sum((pdata%q(iby,nbl:nb,nb:ne,nb:ne)**2 & inarr(23) = inarr(23) &
+ pdata%q(ibz,nbl:nb,nb:ne,nb:ne)**2) & + sum(sum(pdata%q(ibx:ibz,nbm:nb,nb:ne,nb:ne)**2,1) &
* pdata%q(ivx,nbl:nb,nb:ne,nb:ne) & * pdata%q(ivx ,nbm:nb,nb:ne,nb:ne)) * dyz
- (pdata%q(ivy,nbl:nb,nb:ne,nb:ne) & inarr(26) = inarr(26) &
* pdata%q(iby,nbl:nb,nb:ne,nb:ne) & - sum(sum(pdata%q(ivx:ivz,nbm:nb,nb:ne,nb:ne) &
+ pdata%q(ivz,nbl:nb,nb:ne,nb:ne) & * pdata%q(ibx:ibz,nbm:nb,nb:ne,nb:ne),1) &
* pdata%q(ibz,nbl:nb,nb:ne,nb:ne)) & * pdata%q(ibx ,nbm:nb,nb:ne,nb:ne)) * dyz
* pdata%q(ibx,nbl:nb,nb:ne,nb:ne)) * dyz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(23) = inarr(23) + sum((pdata%q(iby,nbl:nb,nb:ne, : )**2 & inarr(23) = inarr(23) &
+ pdata%q(ibz,nbl:nb,nb:ne, : )**2) & + sum(sum(pdata%q(ibx:ibz,nbm:nb,nb:ne, : )**2,1) &
* pdata%q(ivx,nbl:nb,nb:ne, : ) & * pdata%q(ivx ,nbm:nb,nb:ne, : )) * dyz
- (pdata%q(ivy,nbl:nb,nb:ne, : ) & inarr(26) = inarr(26) &
* pdata%q(iby,nbl:nb,nb:ne, : ) & - sum(sum(pdata%q(ivx:ivz,nbm:nb,nb:ne, : ) &
+ pdata%q(ivz,nbl:nb,nb:ne, : ) & * pdata%q(ibx:ibz,nbm:nb,nb:ne, : ),1) &
* pdata%q(ibz,nbl:nb,nb:ne, : )) & * pdata%q(ibx ,nbm:nb,nb:ne, : )) * dyz
* pdata%q(ibx,nbl:nb,nb:ne, : )) * dyz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
@ -728,40 +728,40 @@ module integrals
#if NDIMS == 3 #if NDIMS == 3
! the Y component of current density ! the Y component of current density
! !
tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,nb ,nb:ne,nbu:neu) & tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,nb ,nb:ne,nbp:nep) &
- pdata%q(ibx,nb ,nb:ne,nbl:nel)) / dh(3) & - pdata%q(ibx,nb ,nb:ne,nbm:nem)) / dh(3) &
- (pdata%q(ibz,nb ,nb:ne,nb :ne ) & - (pdata%q(ibz,nb ,nb:ne,nb :ne ) &
- pdata%q(ibz,nbl,nb:ne,nb :ne )) / dh(1) - pdata%q(ibz,nbm,nb:ne,nb :ne )) / dh(1)
! the Z component of current density ! the Z component of current density
! !
tmp(2,:,:) = (pdata%q(iby,nb ,nb :ne ,nb:ne) & tmp(2,:,:) = (pdata%q(iby,nb ,nb :ne ,nb:ne) &
- pdata%q(iby,nbl,nb :ne ,nb:ne)) / dh(1) & - pdata%q(iby,nbm,nb :ne ,nb:ne)) / dh(1) &
- 0.5d+00 * (pdata%q(ibx,nb ,nbu:neu,nb:ne) & - 0.5d+00 * (pdata%q(ibx,nb ,nbp:nep,nb:ne) &
- pdata%q(ibx,nb ,nbl:nel,nb:ne)) / dh(2) - pdata%q(ibx,nb ,nbm:nem,nb:ne)) / dh(2)
! diffusion of magnetic energy through the lower X boundary ! diffusion of magnetic energy through the lower X boundary
! !
inarr(26) = inarr(26) & inarr(29) = inarr(29) &
- sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne,nb:ne) & + sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne,nb:ne) &
- tmp(2,:,:) * pdata%q(iby,nb,nb:ne,nb:ne)) * dyz - tmp(2,:,:) * pdata%q(iby,nb,nb:ne,nb:ne)) * dyz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
! the Y component of current density ! the Y component of current density
! !
tmp(1,:,:) = (pdata%q(ibz,nbl,nb:ne, : ) & tmp(1,:,:) = - (pdata%q(ibz,nb ,nb:ne, : ) &
- pdata%q(ibz,nb ,nb:ne, : )) / dh(1) - pdata%q(ibz,nbm,nb:ne, : )) / dh(1)
! the Z component of current density ! the Z component of current density
! !
tmp(2,:,:) = (pdata%q(iby,nb ,nb :ne , : ) & tmp(2,:,:) = (pdata%q(iby,nb ,nb :ne , : ) &
- pdata%q(iby,nbl,nb :ne , : )) / dh(1) & - pdata%q(iby,nbm,nb :ne , : )) / dh(1) &
- 0.5d+00 * (pdata%q(ibx,nb ,nbu:neu, : ) & - 0.5d+00 * (pdata%q(ibx,nb ,nbp:nep, : ) &
- pdata%q(ibx,nb ,nbl:nel, : )) / dh(2) - pdata%q(ibx,nb ,nbm:nem, : )) / dh(2)
! diffusion of magnetic energy through the lower X boundary ! diffusion of magnetic energy through the lower X boundary
! !
inarr(26) = inarr(26) & inarr(29) = inarr(29) &
- sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne, : ) & + sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne, : ) &
- tmp(2,:,:) * pdata%q(iby,nb,nb:ne, : )) * dyz - tmp(2,:,:) * pdata%q(iby,nb,nb:ne, : )) * dyz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
@ -776,23 +776,21 @@ module integrals
! advection of magnetic energy through the upper X boundary ! advection of magnetic energy through the upper X boundary
! !
#if NDIMS == 3 #if NDIMS == 3
inarr(23) = inarr(23) + sum((pdata%q(iby,ne:neu,nb:ne,nb:ne)**2 & inarr(23) = inarr(23) &
+ pdata%q(ibz,ne:neu,nb:ne,nb:ne)**2) & - sum(sum(pdata%q(ibx:ibz,ne:nep,nb:ne,nb:ne)**2,1) &
* pdata%q(ivx,ne:neu,nb:ne,nb:ne) & * pdata%q(ivx ,ne:nep,nb:ne,nb:ne)) * dyz
- (pdata%q(ivy,ne:neu,nb:ne,nb:ne) & inarr(26) = inarr(26) &
* pdata%q(iby,ne:neu,nb:ne,nb:ne) & + sum(sum(pdata%q(ivx:ivz,ne:nep,nb:ne,nb:ne) &
+ pdata%q(ivz,ne:neu,nb:ne,nb:ne) & * pdata%q(ibx:ibz,ne:nep,nb:ne,nb:ne),1) &
* pdata%q(ibz,ne:neu,nb:ne,nb:ne)) & * pdata%q(ibx ,ne:nep,nb:ne,nb:ne)) * dyz
* pdata%q(ibx,ne:neu,nb:ne,nb:ne)) * dyz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(23) = inarr(23) + sum((pdata%q(iby,ne:neu,nb:ne, : )**2 & inarr(23) = inarr(23) &
+ pdata%q(ibz,ne:neu,nb:ne, : )**2) & - sum(sum(pdata%q(ibx:ibz,ne:nep,nb:ne, : )**2,1) &
* pdata%q(ivx,ne:neu,nb:ne, : ) & * pdata%q(ivx ,ne:nep,nb:ne, : )) * dyz
- (pdata%q(ivy,ne:neu,nb:ne, : ) & inarr(26) = inarr(26) &
* pdata%q(iby,ne:neu,nb:ne, : ) & + sum(sum(pdata%q(ivx:ivz,ne:nep,nb:ne, : ) &
+ pdata%q(ivz,ne:neu,nb:ne, : ) & * pdata%q(ibx:ibz,ne:nep,nb:ne, : ),1) &
* pdata%q(ibz,ne:neu,nb:ne, : )) & * pdata%q(ibx ,ne:nep,nb:ne, : )) * dyz
* pdata%q(ibx,ne:neu,nb:ne, : )) * dyz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
@ -800,40 +798,40 @@ module integrals
#if NDIMS == 3 #if NDIMS == 3
! the Y component of current density ! the Y component of current density
! !
tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,ne ,nb:ne,nbu:neu) & tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,ne ,nb:ne,nbp:nep) &
- pdata%q(ibx,ne ,nb:ne,nbl:nel)) / dh(3) & - pdata%q(ibx,ne ,nb:ne,nbm:nem)) / dh(3) &
- (pdata%q(ibz,neu,nb:ne,nb :ne ) & - (pdata%q(ibz,nep,nb:ne,nb :ne ) &
- pdata%q(ibz,ne ,nb:ne,nb :ne )) / dh(1) - pdata%q(ibz,ne ,nb:ne,nb :ne )) / dh(1)
! the Z component of current density ! the Z component of current density
! !
tmp(2,:,:) = (pdata%q(iby,neu,nb :ne ,nb:ne) & tmp(2,:,:) = (pdata%q(iby,nep,nb :ne ,nb:ne) &
- pdata%q(iby,ne ,nb :ne ,nb:ne)) / dh(1) & - pdata%q(iby,ne ,nb :ne ,nb:ne)) / dh(1) &
- 0.5d+00 * (pdata%q(ibx,ne ,nbu:neu,nb:ne) & - 0.5d+00 * (pdata%q(ibx,ne ,nbp:nep,nb:ne) &
- pdata%q(ibx,ne ,nbl:nel,nb:ne)) / dh(2) - pdata%q(ibx,ne ,nbm:nem,nb:ne)) / dh(2)
! diffusion of magnetic energy through the upper X boundary ! diffusion of magnetic energy through the upper X boundary
! !
inarr(26) = inarr(26) & inarr(29) = inarr(29) &
+ sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne,nb:ne) & - sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne,nb:ne) &
- tmp(2,:,:) * pdata%q(iby,ne,nb:ne,nb:ne)) * dyz - tmp(2,:,:) * pdata%q(iby,ne,nb:ne,nb:ne)) * dyz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
! the Y component of current density ! the Y component of current density
! !
tmp(1,:,:) = (pdata%q(ibz,neu,nb:ne, : ) & tmp(1,:,:) = (pdata%q(ibz,nep,nb:ne, : ) &
- pdata%q(ibz,ne ,nb:ne, : )) / dh(1) - pdata%q(ibz,ne ,nb:ne, : )) / dh(1)
! the Z component of current density ! the Z component of current density
! !
tmp(2,:,:) = (pdata%q(iby,neu,nb :ne , : ) & tmp(2,:,:) = (pdata%q(iby,nep,nb :ne , : ) &
- pdata%q(iby,ne ,nb :ne , : )) / dh(1) & - pdata%q(iby,ne ,nb :ne , : )) / dh(1) &
- 0.5d+00 * (pdata%q(ibx,ne ,nbu:neu, : ) & - 0.5d+00 * (pdata%q(ibx,ne ,nbp:nep, : ) &
- pdata%q(ibx,ne ,nbl:nel, : )) / dh(2) - pdata%q(ibx,ne ,nbm:nem, : )) / dh(2)
! diffusion of magnetic energy through the upper X boundary ! diffusion of magnetic energy through the upper X boundary
! !
inarr(26) = inarr(26) & inarr(29) = inarr(29) &
+ sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne, : ) & - sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne, : ) &
- tmp(2,:,:) * pdata%q(iby,ne,nb:ne, : )) * dyz - tmp(2,:,:) * pdata%q(iby,ne,nb:ne, : )) * dyz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
@ -860,25 +858,25 @@ module integrals
! advection of Bx along Y ! advection of Bx along Y
! !
#if NDIMS == 3 #if NDIMS == 3
inarr(13) = inarr(13) + sum(abs(pdata%q(ibx,nb:ne,nbl:nb,nb:ne)) & inarr(13) = inarr(13) + sum(abs(pdata%q(ibx,nb:ne,nbm:nb,nb:ne)) &
* pdata%q(ivy,nb:ne,nbl:nb,nb:ne)) * dxz * pdata%q(ivy,nb:ne,nbm:nb,nb:ne)) * dxz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(13) = inarr(13) + sum(abs(pdata%q(ibx,nb:ne,nbl:nb, : )) & inarr(13) = inarr(13) + sum(abs(pdata%q(ibx,nb:ne,nbm:nb, : )) &
* pdata%q(ivy,nb:ne,nbl:nb, : )) * dxz * pdata%q(ivy,nb:ne,nbm:nb, : )) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! shear of By along X ! shear of By along X
! !
#if NDIMS == 3 #if NDIMS == 3
inarr(15) = inarr(15) & inarr(15) = inarr(15) &
- sum(sign(pdata%q(iby,nb:ne,nbl:nb,nb:ne) & - sum(sign(pdata%q(iby,nb:ne,nbm:nb,nb:ne) &
* pdata%q(ivx,nb:ne,nbl:nb,nb:ne), & * pdata%q(ivx,nb:ne,nbm:nb,nb:ne), &
pdata%q(ibx,nb:ne,nbl:nb,nb:ne))) * dxz pdata%q(ibx,nb:ne,nbm:nb,nb:ne))) * dxz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(15) = inarr(15) & inarr(15) = inarr(15) &
- sum(sign(pdata%q(iby,nb:ne,nbl:nb, : ) & - sum(sign(pdata%q(iby,nb:ne,nbm:nb, : ) &
* pdata%q(ivx,nb:ne,nbl:nb, : ), & * pdata%q(ivx,nb:ne,nbm:nb, : ), &
pdata%q(ibx,nb:ne,nbl:nb, : ))) * dxz pdata%q(ibx,nb:ne,nbm:nb, : ))) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! mean magnetic energy at the lower Y boundary ! mean magnetic energy at the lower Y boundary
@ -892,23 +890,21 @@ module integrals
! advection of magnetic energy through the lower Y boundary ! advection of magnetic energy through the lower Y boundary
! !
#if NDIMS == 3 #if NDIMS == 3
inarr(24) = inarr(24) + sum((pdata%q(ibx,nb:ne,nbl:nb,nb:ne)**2 & inarr(24) = inarr(24) &
+ pdata%q(ibz,nb:ne,nbl:nb,nb:ne)**2) & + sum(sum(pdata%q(ibx:ibz,nb:ne,nbm:nb,nb:ne)**2,1) &
* pdata%q(ivy,nb:ne,nbl:nb,nb:ne) & * pdata%q(ivy ,nb:ne,nbm:nb,nb:ne)) * dxz
- (pdata%q(ivx,nb:ne,nbl:nb,nb:ne) & inarr(27) = inarr(27) &
* pdata%q(ibx,nb:ne,nbl:nb,nb:ne) & - sum(sum(pdata%q(ivx:ivz,nb:ne,nbm:nb,nb:ne) &
+ pdata%q(ivz,nb:ne,nbl:nb,nb:ne) & * pdata%q(ibx:ibz,nb:ne,nbm:nb,nb:ne),1) &
* pdata%q(ibz,nb:ne,nbl:nb,nb:ne)) & * pdata%q(iby ,nb:ne,nbm:nb,nb:ne)) * dxz
* pdata%q(iby,nb:ne,nbl:nb,nb:ne)) * dxz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(24) = inarr(24) + sum((pdata%q(ibx,nb:ne,nbl:nb, : )**2 & inarr(24) = inarr(24) &
+ pdata%q(ibz,nb:ne,nbl:nb, : )**2) & + sum(sum(pdata%q(ibx:ibz,nb:ne,nbm:nb, : )**2,1) &
* pdata%q(ivy,nb:ne,nbl:nb, : ) & * pdata%q(ivy ,nb:ne,nbm:nb, : )) * dxz
- (pdata%q(ivx,nb:ne,nbl:nb, : ) & inarr(27) = inarr(27) &
* pdata%q(ibx,nb:ne,nbl:nb, : ) & - sum(sum(pdata%q(ivx:ivz,nb:ne,nbm:nb, : ) &
+ pdata%q(ivz,nb:ne,nbl:nb, : ) & * pdata%q(ibx:ibz,nb:ne,nbm:nb, : ),1) &
* pdata%q(ibz,nb:ne,nbl:nb, : )) & * pdata%q(iby ,nb:ne,nbm:nb, : )) * dxz
* pdata%q(iby,nb:ne,nbl:nb, : )) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
@ -916,17 +912,17 @@ module integrals
#if NDIMS == 3 #if NDIMS == 3
! the Z component of current density ! the Z component of current density
! !
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb ,nb:ne) & tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbp:nep,nb ,nb:ne) &
- pdata%q(iby,nbl:nel,nb ,nb:ne)) / dh(1) & - pdata%q(iby,nbm:nem,nb ,nb:ne)) / dh(1) &
- (pdata%q(ibx,nb :ne ,nb ,nb:ne) & - (pdata%q(ibx,nb :ne ,nb ,nb:ne) &
- pdata%q(ibx,nb :ne ,nbl,nb:ne)) / dh(2) - pdata%q(ibx,nb :ne ,nbm,nb:ne)) / dh(2)
! the X component of current density ! the X component of current density
! !
tmp(:,2,:) = (pdata%q(ibz,nb:ne,nb ,nb :ne ) & tmp(:,2,:) = (pdata%q(ibz,nb:ne,nb ,nb :ne ) &
- pdata%q(ibz,nb:ne,nbl,nb :ne )) / dh(2) & - pdata%q(ibz,nb:ne,nbm,nb :ne )) / dh(2) &
- 0.5d+00 * (pdata%q(iby,nb:ne,nb ,nbu:neu) & - 0.5d+00 * (pdata%q(iby,nb:ne,nb ,nbp:nep) &
- pdata%q(iby,nb:ne,nb ,nbl:nel)) / dh(3) - pdata%q(iby,nb:ne,nb ,nbm:nem)) / dh(3)
! diffusion of Bx through the lower Y boundary ! diffusion of Bx through the lower Y boundary
! !
@ -935,21 +931,21 @@ module integrals
! diffusion of magnetic energy through the lower Y boundary ! diffusion of magnetic energy through the lower Y boundary
! !
inarr(27) = inarr(27) & inarr(30) = inarr(30) &
- sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb,nb:ne) & + sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb,nb:ne) &
- tmp(:,2,:) * pdata%q(ibz,nb:ne,nb,nb:ne)) * dxz - tmp(:,2,:) * pdata%q(ibz,nb:ne,nb,nb:ne)) * dxz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
! the Z component of current density ! the Z component of current density
! !
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb , : ) & tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbp:nep,nb , : ) &
- pdata%q(iby,nbl:nel,nb , : )) / dh(1) & - pdata%q(iby,nbm:nem,nb , : )) / dh(1) &
- (pdata%q(ibx,nb :ne ,nb , : ) & - (pdata%q(ibx,nb :ne ,nb , : ) &
- pdata%q(ibx,nb :ne ,nbl, : )) / dh(2) - pdata%q(ibx,nb :ne ,nbm, : )) / dh(2)
! the X component of current density ! the X component of current density
! !
tmp(:,2,:) = (pdata%q(ibz,nb:ne,nb , : ) & tmp(:,2,:) = (pdata%q(ibz,nb:ne,nb , : ) &
- pdata%q(ibz,nb:ne,nbl, : )) / dh(2) - pdata%q(ibz,nb:ne,nbm, : )) / dh(2)
! diffusion of Bx through the lower Y boundary ! diffusion of Bx through the lower Y boundary
! !
@ -958,8 +954,8 @@ module integrals
! diffusion of magnetic energy through the lower Y boundary ! diffusion of magnetic energy through the lower Y boundary
! !
inarr(27) = inarr(27) & inarr(30) = inarr(30) &
- sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb, : ) & + sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb, : ) &
- tmp(:,2,:) * pdata%q(ibz,nb:ne,nb, : )) * dxz - tmp(:,2,:) * pdata%q(ibz,nb:ne,nb, : )) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
@ -982,25 +978,25 @@ module integrals
! advection of Bx along Y ! advection of Bx along Y
! !
#if NDIMS == 3 #if NDIMS == 3
inarr(13) = inarr(13) - sum(abs(pdata%q(ibx,nb:ne,ne:neu,nb:ne)) & inarr(13) = inarr(13) - sum(abs(pdata%q(ibx,nb:ne,ne:nep,nb:ne)) &
* pdata%q(ivy,nb:ne,ne:neu,nb:ne)) * dxz * pdata%q(ivy,nb:ne,ne:nep,nb:ne)) * dxz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(13) = inarr(13) - sum(abs(pdata%q(ibx,nb:ne,ne:neu, : )) & inarr(13) = inarr(13) - sum(abs(pdata%q(ibx,nb:ne,ne:nep, : )) &
* pdata%q(ivy,nb:ne,ne:neu, : )) * dxz * pdata%q(ivy,nb:ne,ne:nep, : )) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! shear of By along X ! shear of By along X
! !
#if NDIMS == 3 #if NDIMS == 3
inarr(15) = inarr(15) & inarr(15) = inarr(15) &
+ sum(sign(pdata%q(iby,nb:ne,ne:neu,nb:ne) & + sum(sign(pdata%q(iby,nb:ne,ne:nep,nb:ne) &
* pdata%q(ivx,nb:ne,ne:neu,nb:ne), & * pdata%q(ivx,nb:ne,ne:nep,nb:ne), &
pdata%q(ibx,nb:ne,ne:neu,nb:ne))) * dxz pdata%q(ibx,nb:ne,ne:nep,nb:ne))) * dxz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(15) = inarr(15) & inarr(15) = inarr(15) &
+ sum(sign(pdata%q(iby,nb:ne,ne:neu, : ) & + sum(sign(pdata%q(iby,nb:ne,ne:nep, : ) &
* pdata%q(ivx,nb:ne,ne:neu, : ), & * pdata%q(ivx,nb:ne,ne:nep, : ), &
pdata%q(ibx,nb:ne,ne:neu, : ))) * dxz pdata%q(ibx,nb:ne,ne:nep, : ))) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! mean magnetic energy at the upper Y boundary ! mean magnetic energy at the upper Y boundary
@ -1014,23 +1010,21 @@ module integrals
! advection of magnetic energy through the upper Y boundary ! advection of magnetic energy through the upper Y boundary
! !
#if NDIMS == 3 #if NDIMS == 3
inarr(24) = inarr(24) - sum((pdata%q(ibx,nb:ne,ne:neu,nb:ne)**2 & inarr(24) = inarr(24) &
+ pdata%q(ibz,nb:ne,ne:neu,nb:ne)**2) & - sum(sum(pdata%q(ibx:ibz,nb:ne,ne:nep,nb:ne)**2,1) &
* pdata%q(ivy,nb:ne,ne:neu,nb:ne) & * pdata%q(ivy ,nb:ne,ne:nep,nb:ne)) * dxz
- (pdata%q(ivx,nb:ne,ne:neu,nb:ne) & inarr(27) = inarr(27) &
* pdata%q(ibx,nb:ne,ne:neu,nb:ne) & + sum(sum(pdata%q(ivx:ivz,nb:ne,ne:nep,nb:ne) &
+ pdata%q(ivz,nb:ne,ne:neu,nb:ne) & * pdata%q(ibx:ibz,nb:ne,ne:nep,nb:ne),1) &
* pdata%q(ibz,nb:ne,ne:neu,nb:ne)) & * pdata%q(iby ,nb:ne,ne:nep,nb:ne)) * dxz
* pdata%q(iby,nb:ne,ne:neu,nb:ne)) * dxz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(24) = inarr(24) - sum((pdata%q(ibx,nb:ne,ne:neu, : )**2 & inarr(24) = inarr(24) &
+ pdata%q(ibz,nb:ne,ne:neu, : )**2) & - sum(sum(pdata%q(ibx:ibz,nb:ne,ne:nep, : )**2,1) &
* pdata%q(ivy,nb:ne,ne:neu, : ) & * pdata%q(ivy ,nb:ne,ne:nep, : )) * dxz
- (pdata%q(ivx,nb:ne,ne:neu, : ) & inarr(27) = inarr(27) &
* pdata%q(ibx,nb:ne,ne:neu, : ) & + sum(sum(pdata%q(ivx:ivz,nb:ne,ne:nep, : ) &
+ pdata%q(ivz,nb:ne,ne:neu, : ) & * pdata%q(ibx:ibz,nb:ne,ne:nep, : ),1) &
* pdata%q(ibz,nb:ne,ne:neu, : )) & * pdata%q(iby ,nb:ne,ne:nep, : )) * dxz
* pdata%q(iby,nb:ne,ne:neu, : )) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
@ -1038,17 +1032,17 @@ module integrals
#if NDIMS == 3 #if NDIMS == 3
! the Z component of current density ! the Z component of current density
! !
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,ne ,nb:ne) & tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbp:nep,ne ,nb:ne) &
- pdata%q(iby,nbl:nel,ne ,nb:ne)) / dh(1) & - pdata%q(iby,nbm:nem,ne ,nb:ne)) / dh(1) &
- (pdata%q(ibx,nb :ne ,neu,nb:ne) & - (pdata%q(ibx,nb :ne ,nep,nb:ne) &
- pdata%q(ibx,nb :ne ,ne ,nb:ne)) / dh(2) - pdata%q(ibx,nb :ne ,ne ,nb:ne)) / dh(2)
! the X component of current density ! the X component of current density
tmp(:,2,:) = 0.5d+00 * (pdata%q(ibx,nb :ne ,ne,nbu:neu) & tmp(:,2,:) = 0.5d+00 * (pdata%q(ibx,nb :ne ,ne,nbp:nep) &
- pdata%q(ibx,nb :ne ,ne,nbl:nel)) / dh(3) & - pdata%q(ibx,nb :ne ,ne,nbm:nem)) / dh(3) &
- 0.5d+00 * (pdata%q(ibz,nbu:neu,ne,nb :ne ) & - 0.5d+00 * (pdata%q(ibz,nbp:nep,ne,nb :ne ) &
- pdata%q(ibz,nbl:nel,ne,nb :ne )) / dh(1) - pdata%q(ibz,nbm:nem,ne,nb :ne )) / dh(1)
! diffusion of Bx through the upper Y boundary ! diffusion of Bx through the upper Y boundary
! !
@ -1057,21 +1051,21 @@ module integrals
! diffusion of magnetic energy through the upper Y boundary ! diffusion of magnetic energy through the upper Y boundary
! !
inarr(27) = inarr(27) & inarr(30) = inarr(30) &
+ sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne,nb:ne) & - sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne,nb:ne) &
- tmp(:,2,:) * pdata%q(ibz,nb:ne,ne,nb:ne)) * dxz - tmp(:,2,:) * pdata%q(ibz,nb:ne,ne,nb:ne)) * dxz
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
! the Z component of current density ! the Z component of current density
! !
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,ne , : ) & tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbp:nep,ne , : ) &
- pdata%q(iby,nbl:nel,ne , : )) / dh(1) & - pdata%q(iby,nbm:nem,ne , : )) / dh(1) &
- (pdata%q(ibx,nb :ne ,neu, : ) & - (pdata%q(ibx,nb :ne ,nep, : ) &
- pdata%q(ibx,nb :ne ,ne , : )) / dh(2) - pdata%q(ibx,nb :ne ,ne , : )) / dh(2)
! the X component of current density ! the X component of current density
! !
tmp(:,2,:) = 0.5d+00 * (pdata%q(ibz,nbl:nel,ne, : ) & tmp(:,2,:) = 0.5d+00 * (pdata%q(ibz,nbm:nem,ne, : ) &
- pdata%q(ibz,nbu:neu,ne, : )) / dh(1) - pdata%q(ibz,nbp:nep,ne, : )) / dh(1)
! diffusion of Bx through the upper Y boundary ! diffusion of Bx through the upper Y boundary
! !
@ -1080,8 +1074,8 @@ module integrals
! diffusion of magnetic energy through the upper Y boundary ! diffusion of magnetic energy through the upper Y boundary
! !
inarr(27) = inarr(27) & inarr(30) = inarr(30) &
+ sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne, : ) & - sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne, : ) &
- tmp(:,2,:) * pdata%q(ibz,nb:ne,ne, : )) * dxz - tmp(:,2,:) * pdata%q(ibz,nb:ne,ne, : )) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
@ -1100,41 +1094,40 @@ module integrals
! advection of Bx along Z ! advection of Bx along Z
! !
inarr(14) = inarr(14) + sum(abs(pdata%q(ibx,nb:ne,nb:ne,nbl:nb)) & inarr(14) = inarr(14) + sum(abs(pdata%q(ibx,nb:ne,nb:ne,nbm:nb)) &
* pdata%q(ivz,nb:ne,nb:ne,nbl:nb)) * dxy * pdata%q(ivz,nb:ne,nb:ne,nbm:nb)) * dxy
! shear of Bz along X ! shear of Bz along X
! !
inarr(16) = inarr(16) - sum(sign(pdata%q(ibz,nb:ne,nb:ne,nbl:nb) & inarr(16) = inarr(16) - sum(sign(pdata%q(ibz,nb:ne,nb:ne,nbm:nb) &
* pdata%q(ivx,nb:ne,nb:ne,nbl:nb), & * pdata%q(ivx,nb:ne,nb:ne,nbm:nb), &
pdata%q(ibx,nb:ne,nb:ne,nbl:nb))) * dxy pdata%q(ibx,nb:ne,nb:ne,nbm:nb))) * dxy
! advection of magnetic energy through the lower Z boundary ! advection of magnetic energy through the lower Z boundary
! !
inarr(25) = inarr(25) + sum((pdata%q(ibx,nb:ne,nb:ne,nbl:nb)**2 & inarr(25) = inarr(25) &
+ pdata%q(iby,nb:ne,nb:ne,nbl:nb)**2) & + sum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nbm:nb)**2,1) &
* pdata%q(ivz,nb:ne,nb:ne,nbl:nb) & * pdata%q(ivz ,nb:ne,nb:ne,nbm:nb)) * dxy
- (pdata%q(ivx,nb:ne,nb:ne,nbl:nb) & inarr(28) = inarr(28) &
* pdata%q(ibx,nb:ne,nb:ne,nbl:nb) & - sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nbm:nb) &
+ pdata%q(ivy,nb:ne,nb:ne,nbl:nb) & * pdata%q(ibx:ibz,nb:ne,nb:ne,nbm:nb),1) &
* pdata%q(iby,nb:ne,nb:ne,nbl:nb)) & * pdata%q(ibz ,nb:ne,nb:ne,nbm:nb)) * dxy
* pdata%q(ibz,nb:ne,nb:ne,nbl:nb)) * dxy
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
! the X component of current density ! the X component of current density
! !
tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbu:neu,nb ) & tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbp:nep,nb ) &
- pdata%q(ibz,nb:ne,nbl:nel,nb )) / dh(2) & - pdata%q(ibz,nb:ne,nbm:nem,nb )) / dh(2) &
- (pdata%q(iby,nb:ne,nb :ne ,nb ) & - (pdata%q(iby,nb:ne,nb :ne ,nb ) &
- pdata%q(iby,nb:ne,nb :ne ,nbl)) / dh(3) - pdata%q(iby,nb:ne,nb :ne ,nbm)) / dh(3)
! the Y component of current density ! the Y component of current density
! !
tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,nb ) & tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,nb ) &
- pdata%q(ibx,nb :ne ,nb:ne,nbl)) / dh(3) & - pdata%q(ibx,nb :ne ,nb:ne,nbm)) / dh(3) &
- 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,nb ) & - 0.5d+00 * (pdata%q(ibz,nbp:nep,nb:ne,nb ) &
- pdata%q(ibz,nbl:nel,nb:ne,nb )) / dh(1) - pdata%q(ibz,nbm:nem,nb:ne,nb )) / dh(1)
! diffusion of Bx through the lower Z boundary ! diffusion of Bx through the lower Z boundary
! !
@ -1143,8 +1136,8 @@ module integrals
! diffusion of magnetic energy through the lower Z boundary ! diffusion of magnetic energy through the lower Z boundary
! !
inarr(28) = inarr(28) & inarr(31) = inarr(31) &
- sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,nb) & + sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,nb) &
- tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,nb)) * dxy - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,nb)) * dxy
end if ! resistivity end if ! resistivity
@ -1157,41 +1150,40 @@ module integrals
! advection of Bx along Z ! advection of Bx along Z
! !
inarr(14) = inarr(14) - sum(abs(pdata%q(ibx,nb:ne,nb:ne,ne:neu)) & inarr(14) = inarr(14) - sum(abs(pdata%q(ibx,nb:ne,nb:ne,ne:nep)) &
* pdata%q(ivz,nb:ne,nb:ne,ne:neu)) * dxy * pdata%q(ivz,nb:ne,nb:ne,ne:nep)) * dxy
! shear of Bz along X ! shear of Bz along X
! !
inarr(16) = inarr(16) + sum(sign(pdata%q(ibz,nb:ne,nb:ne,ne:neu) & inarr(16) = inarr(16) + sum(sign(pdata%q(ibz,nb:ne,nb:ne,ne:nep) &
* pdata%q(ivx,nb:ne,nb:ne,ne:neu), & * pdata%q(ivx,nb:ne,nb:ne,ne:nep), &
pdata%q(ibx,nb:ne,nb:ne,ne:neu))) * dxy pdata%q(ibx,nb:ne,nb:ne,ne:nep))) * dxy
! advection of magnetic energy through the upper Z boundary ! advection of magnetic energy through the upper Z boundary
! !
inarr(25) = inarr(25) + sum((pdata%q(ibx,nb:ne,nb:ne,ne:neu)**2 & inarr(25) = inarr(25) &
+ pdata%q(iby,nb:ne,nb:ne,ne:neu)**2) & - sum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,ne:nep)**2,1) &
* pdata%q(ivz,nb:ne,nb:ne,ne:neu) & * pdata%q(ivz ,nb:ne,nb:ne,ne:nep)) * dxy
- (pdata%q(ivx,nb:ne,nb:ne,ne:neu) & inarr(28) = inarr(28) &
* pdata%q(ibx,nb:ne,nb:ne,ne:neu) & + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,ne:nep) &
+ pdata%q(ivy,nb:ne,nb:ne,ne:neu) & * pdata%q(ibx:ibz,nb:ne,nb:ne,ne:nep),1) &
* pdata%q(iby,nb:ne,nb:ne,ne:neu)) & * pdata%q(ibz ,nb:ne,nb:ne,ne:nep)) * dxy
* pdata%q(ibz,nb:ne,nb:ne,ne:neu)) * dxy
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
! the X component of current density ! the X component of current density
! !
tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbu:neu,ne ) & tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbp:nep,ne ) &
- pdata%q(ibz,nb:ne,nbl:nel,ne )) / dh(2) & - pdata%q(ibz,nb:ne,nbm:nem,ne )) / dh(2) &
- (pdata%q(iby,nb:ne,nb :ne ,neu ) & - (pdata%q(iby,nb:ne,nb :ne ,nep) &
- pdata%q(iby,nb:ne,nb :ne ,ne )) / dh(3) - pdata%q(iby,nb:ne,nb :ne ,ne )) / dh(3)
! the Y component of current density ! the Y component of current density
! !
tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,neu) & tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,nep) &
- pdata%q(ibx,nb :ne ,nb:ne,ne )) / dh(3) & - pdata%q(ibx,nb :ne ,nb:ne,ne )) / dh(3) &
- 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,ne ) & - 0.5d+00 * (pdata%q(ibz,nbp:nep,nb:ne,ne ) &
- pdata%q(ibz,nbl:nel,nb:ne,ne )) / dh(1) - pdata%q(ibz,nbm:nem,nb:ne,ne )) / dh(1)
! diffusion of Bx through the upper Z boundary ! diffusion of Bx through the upper Z boundary
! !
@ -1200,8 +1192,8 @@ module integrals
! diffusion of magnetic energy through the upper Z boundary ! diffusion of magnetic energy through the upper Z boundary
! !
inarr(28) = inarr(28) & inarr(31) = inarr(31) &
+ sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,ne) & - sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,ne) &
- tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,ne)) * dxy - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,ne)) * dxy
end if ! resistivity end if ! resistivity
@ -1211,6 +1203,92 @@ module integrals
end if ! not periodic in Z end if ! not periodic in Z
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! get the conversion between the magnetic energy and kinetic and
! internal energies
!
#if NDIMS == 3
! the X component of current density
!
tmp(:,:,:) = 5.0d-01 * (pdata%q(ibz,nb:ne,nbp:nep,nb :ne ) &
- pdata%q(ibz,nb:ne,nbm:nem,nb :ne )) / dh(2) &
- 5.0d-01 * (pdata%q(iby,nb:ne,nb :ne ,nbp:nep) &
- pdata%q(iby,nb:ne,nb :ne ,nbm:nem)) / dh(3)
inarr(32) = inarr(32) &
+ sum((pdata%q(ivy,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibz,nb:ne,nb:ne,nb:ne) &
- pdata%q(ivz,nb:ne,nb:ne,nb:ne) &
* pdata%q(iby,nb:ne,nb:ne,nb:ne)) * tmp(:,:,:)) * dvol
#else /* NDIMS == 3 */
tmp(:,:,:) = 5.0d-01 * (pdata%q(ibz,nb:ne,nbp:nep, : ) &
- pdata%q(ibz,nb:ne,nbm:nem, : )) / dh(2)
inarr(32) = inarr(32) &
+ sum((pdata%q(ivy,nb:ne,nb:ne, : ) &
* pdata%q(ibz,nb:ne,nb:ne, : ) &
- pdata%q(ivz,nb:ne,nb:ne, : ) &
* pdata%q(iby,nb:ne,nb:ne, : )) * tmp(:,:,:)) * dvol
#endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) &
inarr(33) = inarr(33) - sum(tmp(:,:,:)**2) * dvol
! the Y component of current density
!
#if NDIMS == 3
tmp(:,:,:) = 5.0d-01 * (pdata%q(ibx,nb :ne ,nb:ne,nbp:nep) &
- pdata%q(ibx,nb :ne ,nb:ne,nbm:nem)) / dh(3) &
- 5.0d-01 * (pdata%q(ibz,nbp:nep,nb:ne,nb :ne ) &
- pdata%q(ibz,nbm:nem,nb:ne,nb :ne )) / dh(1)
inarr(32) = inarr(32) &
+ sum((pdata%q(ivz,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibx,nb:ne,nb:ne,nb:ne) &
- pdata%q(ivx,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibz,nb:ne,nb:ne,nb:ne)) * tmp(:,:,:)) * dvol
#else /* NDIMS == 3 */
tmp(:,:,:) = - 5.0d-01 * (pdata%q(ibz,nbp:nep,nb:ne, : ) &
- pdata%q(ibz,nbm:nem,nb:ne, : )) / dh(1)
inarr(32) = inarr(32) &
+ sum((pdata%q(ivz,nb:ne,nb:ne, : ) &
* pdata%q(ibx,nb:ne,nb:ne, : ) &
- pdata%q(ivx,nb:ne,nb:ne, : ) &
* pdata%q(ibz,nb:ne,nb:ne, : )) * tmp(:,:,:)) * dvol
#endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) &
inarr(33) = inarr(33) - sum(tmp(:,:,:)**2) * dvol
! the Z component of current density
!
#if NDIMS == 3
tmp(:,:,:) = (pdata%q(iby,nbp:nep,nb :ne ,nb:ne) &
- pdata%q(iby,nbm:nem,nb :ne ,nb:ne)) / dh(1) &
- (pdata%q(ibx,nb :ne ,nbp:nep,nb:ne) &
- pdata%q(ibx,nb :ne ,nbm:nem,nb:ne)) / dh(2)
inarr(32) = inarr(32) &
+ sum((pdata%q(ivx,nb:ne,nb:ne,nb:ne) &
* pdata%q(iby,nb:ne,nb:ne,nb:ne) &
- pdata%q(ivy,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibx,nb:ne,nb:ne,nb:ne)) * tmp(:,:,:)) * dvol
#else /* NDIMS == 3 */
tmp(:,:,:) = (pdata%q(iby,nbp:nep,nb :ne , : ) &
- pdata%q(iby,nbm:nem,nb :ne , : )) / dh(1) &
- (pdata%q(ibx,nb :ne ,nbp:nep, : ) &
- pdata%q(ibx,nb :ne ,nbm:nem, : )) / dh(2)
inarr(32) = inarr(32) &
+ sum((pdata%q(ivx,nb:ne,nb:ne, : ) &
* pdata%q(iby,nb:ne,nb:ne, : ) &
- pdata%q(ivy,nb:ne,nb:ne, : ) &
* pdata%q(ibx,nb:ne,nb:ne, : )) * tmp(:,:,:)) * dvol
#endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) &
inarr(33) = inarr(33) - sum(tmp(:,:,:)**2) * dvol
end if ! magnetized end if ! magnetized
! get the inflow speed ! get the inflow speed
@ -1262,8 +1340,10 @@ module integrals
inarr(13:16) = 5.0d-01 * inarr(13:16) inarr(13:16) = 5.0d-01 * inarr(13:16)
inarr(17:18) = resistivity * inarr(17:18) inarr(17:18) = resistivity * inarr(17:18)
inarr(22) = 2.5d-01 * inarr(22) / yarea inarr(22) = 2.5d-01 * inarr(22) / yarea
inarr(23:25) = 5.0d-01 * inarr(23:25) inarr(23:28) = 5.0d-01 * inarr(23:28)
inarr(26:28) = resistivity * inarr(26:28) inarr(29:31) = resistivity * inarr(29:31)
inarr(32) = 5.0d-01 * inarr(32)
inarr(33) = 2.5d-01 * resistivity * inarr(33)
! write down the integrals and statistics to appropriate files ! write down the integrals and statistics to appropriate files
! !
@ -1278,7 +1358,7 @@ module integrals
, avarr(6), mnarr(6), mxarr(6) & , avarr(6), mnarr(6), mxarr(6) &
, avarr(7), mnarr(7), mxarr(7) , avarr(7), mnarr(7), mxarr(7)
write(eunit,efmt) step, time, dtn, dte, maxval(errors(:)), errors(:) write(eunit,efmt) step, time, dtn, dte, maxval(errors(:)), errors(:)
write(runit,"(i9, 19es25.15e3)") step, time, inarr(11:28) write(runit,"(i9, 24es25.15e3)") step, time, inarr(11:33)
call flush_and_sync(funit) call flush_and_sync(funit)
call flush_and_sync(sunit) call flush_and_sync(sunit)