INTEGRALS: Remove reconnection rate terms from INTEGRALS.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
e4a7890959
commit
a3e0c0ccd2
@ -57,7 +57,6 @@ module integrals
|
|||||||
integer(kind=4), save :: funit = 11
|
integer(kind=4), save :: funit = 11
|
||||||
integer(kind=4), save :: sunit = 12
|
integer(kind=4), save :: sunit = 12
|
||||||
integer(kind=4), save :: eunit = 13
|
integer(kind=4), save :: eunit = 13
|
||||||
integer(kind=4), save :: runit = 14
|
|
||||||
integer(kind=4), save :: iintd = 1
|
integer(kind=4), save :: iintd = 1
|
||||||
|
|
||||||
! flag indicating if the files are actually written to disk
|
! flag indicating if the files are actually written to disk
|
||||||
@ -294,55 +293,6 @@ module integrals
|
|||||||
!
|
!
|
||||||
efmt = "(i9," // trim(adjustl(stmp)) // "(1x,1es18.8e3))"
|
efmt = "(i9," // trim(adjustl(stmp)) // "(1x,1es18.8e3))"
|
||||||
|
|
||||||
! depending on the append parameter, choose the right file initialization for
|
|
||||||
! the reconnection file
|
|
||||||
!
|
|
||||||
append = "off"
|
|
||||||
call get_parameter("reconnection_append", append)
|
|
||||||
select case(trim(append))
|
|
||||||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
|
||||||
write(fname, "('reconnection.dat')")
|
|
||||||
inquire(file = fname, exist = flag)
|
|
||||||
case default
|
|
||||||
write(fname, "('reconnection_',i2.2,'.dat')") irun
|
|
||||||
flag = .false.
|
|
||||||
end select
|
|
||||||
|
|
||||||
! check if the file exists; if not, create a new one, otherwise move to the end
|
|
||||||
!
|
|
||||||
if (flag .and. irun > 1) then
|
|
||||||
#ifdef __INTEL_COMPILER
|
|
||||||
open(newunit = runit, file = fname, form = 'formatted', status = 'old' &
|
|
||||||
, position = 'append', buffered = 'yes')
|
|
||||||
#else /* __INTEL_COMPILER */
|
|
||||||
open(newunit = runit, file = fname, form = 'formatted', status = 'old' &
|
|
||||||
, position = 'append')
|
|
||||||
#endif /* __INTEL_COMPILER */
|
|
||||||
write(runit,"('#')")
|
|
||||||
else
|
|
||||||
#ifdef __INTEL_COMPILER
|
|
||||||
open(newunit = runit, file = fname, form = 'formatted' &
|
|
||||||
, status = 'replace', buffered = 'yes')
|
|
||||||
#else /* __INTEL_COMPILER */
|
|
||||||
open(newunit = runit, file = fname, form = 'formatted' &
|
|
||||||
, status = 'replace')
|
|
||||||
#endif /* __INTEL_COMPILER */
|
|
||||||
end if
|
|
||||||
|
|
||||||
! write the integral file header
|
|
||||||
!
|
|
||||||
write(runit,'("#",a8,27a25)') &
|
|
||||||
'step', 'time', '|Bx| int', '|Bx| inf' &
|
|
||||||
, '|Bx| y-adv', '|Bx| z-adv', '|Bx| y-shr', '|Bx| z-shr' &
|
|
||||||
, '|Bx| y-dif', '|Bx| z-dif', '|Bx| Psi' &
|
|
||||||
, 'Vin lower' , 'Vin upper' &
|
|
||||||
, 'Emag', 'Emag inf' &
|
|
||||||
, 'Emag x-adv', 'Emag y-adv', 'Emag z-adv' &
|
|
||||||
, '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', 'Emag-Psi'
|
|
||||||
write(runit,"('#')")
|
|
||||||
|
|
||||||
end if ! store
|
end if ! store
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
@ -396,7 +346,6 @@ module integrals
|
|||||||
close(funit)
|
close(funit)
|
||||||
close(sunit)
|
close(sunit)
|
||||||
close(eunit)
|
close(eunit)
|
||||||
close(runit)
|
|
||||||
end if
|
end if
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
@ -425,16 +374,9 @@ module integrals
|
|||||||
! import external variables and subroutines
|
! import external variables and subroutines
|
||||||
!
|
!
|
||||||
use blocks , only : block_meta, block_data, list_data
|
use blocks , only : block_meta, block_data, list_data
|
||||||
use coordinates , only : ni => ncells, nn => bcells
|
use coordinates , only : ni => ncells, nb, ne
|
||||||
use coordinates , only : nb, nbm, nbp, ne, nem, nep
|
use coordinates , only : advol, voli
|
||||||
use coordinates , only : adx, ady, adz, advol, voli
|
use equations , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp
|
||||||
use coordinates , only : periodic
|
|
||||||
use coordinates , only : xmin, xmax
|
|
||||||
use coordinates , only : ymin, ymax, yarea
|
|
||||||
#if NDIMS == 3
|
|
||||||
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
|
use equations , only : ien, imx, imy, imz
|
||||||
use equations , only : magnetized, adiabatic_index, csnd
|
use equations , only : magnetized, adiabatic_index, csnd
|
||||||
use equations , only : errors
|
use equations , only : errors
|
||||||
@ -444,8 +386,6 @@ module integrals
|
|||||||
#ifdef MPI
|
#ifdef MPI
|
||||||
use mpitools , only : reduce_minimum, reduce_maximum, reduce_sum
|
use mpitools , only : reduce_minimum, reduce_maximum, reduce_sum
|
||||||
#endif /* MPI */
|
#endif /* MPI */
|
||||||
use operators , only : gradient, curl
|
|
||||||
use sources , only : resistivity
|
|
||||||
|
|
||||||
! local variables are not implicit by default
|
! local variables are not implicit by default
|
||||||
!
|
!
|
||||||
@ -453,11 +393,7 @@ module integrals
|
|||||||
|
|
||||||
! local variables
|
! local variables
|
||||||
!
|
!
|
||||||
real(kind=8) :: dvol, dvolh, dyz, dxz
|
real(kind=8) :: dvol, dvolh
|
||||||
#if NDIMS == 3
|
|
||||||
real(kind=8) :: dxy
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
real(kind=8), dimension(3) :: dh
|
|
||||||
|
|
||||||
! local pointers
|
! local pointers
|
||||||
!
|
!
|
||||||
@ -465,7 +401,7 @@ module integrals
|
|||||||
|
|
||||||
! local parameters
|
! local parameters
|
||||||
!
|
!
|
||||||
integer, parameter :: narr = 64
|
integer, parameter :: narr = 16
|
||||||
|
|
||||||
! local arrays
|
! local arrays
|
||||||
!
|
!
|
||||||
@ -475,11 +411,6 @@ module integrals
|
|||||||
#else /* NDIMS == 3 */
|
#else /* NDIMS == 3 */
|
||||||
real(kind=8), dimension(ni,ni, 1) :: vel, mag, sqd, tmp
|
real(kind=8), dimension(ni,ni, 1) :: vel, mag, sqd, tmp
|
||||||
#endif /* NDIMS == 3 */
|
#endif /* NDIMS == 3 */
|
||||||
#if NDIMS == 3
|
|
||||||
real(kind=8), dimension(3,nn,nn,nn) :: cur
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
real(kind=8), dimension(3,nn,nn, 1) :: cur
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! parameters
|
! parameters
|
||||||
!
|
!
|
||||||
@ -532,14 +463,6 @@ module integrals
|
|||||||
!
|
!
|
||||||
dvol = advol(pdata%meta%level)
|
dvol = advol(pdata%meta%level)
|
||||||
dvolh = 0.5d+00 * dvol
|
dvolh = 0.5d+00 * dvol
|
||||||
dyz = ady(pdata%meta%level) * adz(pdata%meta%level)
|
|
||||||
#if NDIMS == 3
|
|
||||||
dxy = adx(pdata%meta%level) * ady(pdata%meta%level)
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
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)
|
|
||||||
|
|
||||||
! sum up density and momenta components
|
! sum up density and momenta components
|
||||||
!
|
!
|
||||||
@ -687,546 +610,6 @@ module integrals
|
|||||||
avarr(7) = avarr(7) + sum(tmp(:,:,:)) * dvol
|
avarr(7) = avarr(7) + sum(tmp(:,:,:)) * dvol
|
||||||
mnarr(7) = min(mnarr(7), minval(tmp(:,:,:)))
|
mnarr(7) = min(mnarr(7), minval(tmp(:,:,:)))
|
||||||
mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:)))
|
mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:)))
|
||||||
|
|
||||||
!! RECONNECTION RATE MEASUREMENTS
|
|
||||||
!!
|
|
||||||
! calculate current density (J = ∇xB)
|
|
||||||
!
|
|
||||||
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), cur(1:3,:,:,:))
|
|
||||||
|
|
||||||
! the integral of |Bx|
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(11) = inarr(11) + sum(abs(pdata%q(ibx,nb:ne,nb:ne,nb:ne))) * dvol
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(11) = inarr(11) + sum(abs(pdata%q(ibx,nb:ne,nb:ne, : ))) * dvol
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! magnetic energy
|
|
||||||
!
|
|
||||||
inarr(22) = 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(24) = inarr(24) &
|
|
||||||
+ sum(sum(pdata%q(ibx:ibz,nbm:nb,nb:ne,nb:ne)**2,1) &
|
|
||||||
* pdata%q(ivx ,nbm:nb,nb:ne,nb:ne)) * dyz
|
|
||||||
inarr(27) = inarr(27) &
|
|
||||||
- sum(sum(pdata%q(ivx:ivz,nbm:nb,nb:ne,nb:ne) &
|
|
||||||
* pdata%q(ibx:ibz,nbm:nb,nb:ne,nb:ne),1) &
|
|
||||||
* pdata%q(ibx ,nbm:nb,nb:ne,nb:ne)) * dyz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(24) = inarr(24) &
|
|
||||||
+ sum(sum(pdata%q(ibx:ibz,nbm:nb,nb:ne, : )**2,1) &
|
|
||||||
* pdata%q(ivx ,nbm:nb,nb:ne, : )) * dyz
|
|
||||||
inarr(27) = inarr(27) &
|
|
||||||
- sum(sum(pdata%q(ivx:ivz,nbm:nb,nb:ne, : ) &
|
|
||||||
* pdata%q(ibx:ibz,nbm:nb,nb:ne, : ),1) &
|
|
||||||
* pdata%q(ibx ,nbm:nb,nb:ne, : )) * dyz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
if (resistivity > 0.0d+00) then
|
|
||||||
|
|
||||||
#if NDIMS == 3
|
|
||||||
! diffusion of magnetic energy through the lower X boundary
|
|
||||||
!
|
|
||||||
inarr(30) = inarr(30) &
|
|
||||||
+ sum(cur(2,nbm:nb,nb:ne,nb:ne) &
|
|
||||||
* pdata%q(ibz,nbm:nb,nb:ne,nb:ne) &
|
|
||||||
- cur(3,nbm:nb,nb:ne,nb:ne) &
|
|
||||||
* pdata%q(iby,nbm:nb,nb:ne,nb:ne)) * dyz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
! diffusion of magnetic energy through the lower X boundary
|
|
||||||
!
|
|
||||||
inarr(30) = inarr(30) &
|
|
||||||
+ sum(cur(2,nbm:nb,nb:ne, : ) &
|
|
||||||
* pdata%q(ibz,nbm:nb,nb:ne, : ) &
|
|
||||||
- cur(3,nbm:nb,nb:ne, : ) &
|
|
||||||
* pdata%q(iby,nbm: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(24) = inarr(24) &
|
|
||||||
- sum(sum(pdata%q(ibx:ibz,ne:nep,nb:ne,nb:ne)**2,1) &
|
|
||||||
* pdata%q(ivx ,ne:nep,nb:ne,nb:ne)) * dyz
|
|
||||||
inarr(27) = inarr(27) &
|
|
||||||
+ sum(sum(pdata%q(ivx:ivz,ne:nep,nb:ne,nb:ne) &
|
|
||||||
* pdata%q(ibx:ibz,ne:nep,nb:ne,nb:ne),1) &
|
|
||||||
* pdata%q(ibx ,ne:nep,nb:ne,nb:ne)) * dyz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(24) = inarr(24) &
|
|
||||||
- sum(sum(pdata%q(ibx:ibz,ne:nep,nb:ne, : )**2,1) &
|
|
||||||
* pdata%q(ivx ,ne:nep,nb:ne, : )) * dyz
|
|
||||||
inarr(27) = inarr(27) &
|
|
||||||
+ sum(sum(pdata%q(ivx:ivz,ne:nep,nb:ne, : ) &
|
|
||||||
* pdata%q(ibx:ibz,ne:nep,nb:ne, : ),1) &
|
|
||||||
* pdata%q(ibx ,ne:nep,nb:ne, : )) * dyz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
if (resistivity > 0.0d+00) then
|
|
||||||
|
|
||||||
#if NDIMS == 3
|
|
||||||
! diffusion of magnetic energy through the upper X boundary
|
|
||||||
!
|
|
||||||
|
|
||||||
inarr(30) = inarr(30) &
|
|
||||||
- sum(cur(2,ne:nep,nb:ne,nb:ne) &
|
|
||||||
* pdata%q(ibz,ne:nep,nb:ne,nb:ne) &
|
|
||||||
- cur(3,ne:nep,nb:ne,nb:ne) &
|
|
||||||
* pdata%q(iby,ne:nep,nb:ne,nb:ne)) * dyz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
! diffusion of magnetic energy through the upper X boundary
|
|
||||||
!
|
|
||||||
inarr(30) = inarr(30) &
|
|
||||||
- sum(cur(2,ne:nep,nb:ne, : ) &
|
|
||||||
* pdata%q(ibz,ne:nep,nb:ne, : ) &
|
|
||||||
- cur(3,ne:nep,nb:ne, : ) &
|
|
||||||
* pdata%q(iby,ne:nep,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
|
|
||||||
|
|
||||||
! mean Bx at boundary
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(12) = inarr(12) &
|
|
||||||
+ sum(abs(pdata%q(ibx,nb:ne,nbm:nb,nb:ne))) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(12) = inarr(12) &
|
|
||||||
+ sum(abs(pdata%q(ibx,nb:ne,nbm:nb, : ))) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! advection of Bx along Y
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(13) = inarr(13) &
|
|
||||||
+ sum(abs(pdata%q(ibx,nb:ne,nbm:nb,nb:ne)) &
|
|
||||||
* pdata%q(ivy,nb:ne,nbm:nb,nb:ne)) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(13) = inarr(13) &
|
|
||||||
+ sum(abs(pdata%q(ibx,nb:ne,nbm:nb, : )) &
|
|
||||||
* pdata%q(ivy,nb:ne,nbm:nb, : )) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! shear of By along X
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(15) = inarr(15) &
|
|
||||||
- sum(sign(pdata%q(iby,nb:ne,nbm:nb,nb:ne) &
|
|
||||||
* pdata%q(ivx,nb:ne,nbm:nb,nb:ne), &
|
|
||||||
pdata%q(ibx,nb:ne,nbm:nb,nb:ne))) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(15) = inarr(15) &
|
|
||||||
- sum(sign(pdata%q(iby,nb:ne,nbm:nb, : ) &
|
|
||||||
* pdata%q(ivx,nb:ne,nbm:nb, : ), &
|
|
||||||
pdata%q(ibx,nb:ne,nbm:nb, : ))) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! mean magnetic energy at the lower Y boundary
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(23) = inarr(23) &
|
|
||||||
+ sum(pdata%q(ibx:ibz,nb:ne,nbm:nb,nb:ne)**2) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(23) = inarr(23) &
|
|
||||||
+ sum(pdata%q(ibx:ibz,nb:ne,nbm:nb, : )**2) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! advection of magnetic energy through the lower Y boundary
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(25) = inarr(25) &
|
|
||||||
+ sum(sum(pdata%q(ibx:ibz,nb:ne,nbm:nb,nb:ne)**2,1) &
|
|
||||||
* pdata%q(ivy ,nb:ne,nbm:nb,nb:ne)) * dxz
|
|
||||||
inarr(28) = inarr(28) &
|
|
||||||
- sum(sum(pdata%q(ivx:ivz,nb:ne,nbm:nb,nb:ne) &
|
|
||||||
* pdata%q(ibx:ibz,nb:ne,nbm:nb,nb:ne),1) &
|
|
||||||
* pdata%q(iby ,nb:ne,nbm:nb,nb:ne)) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(25) = inarr(25) &
|
|
||||||
+ sum(sum(pdata%q(ibx:ibz,nb:ne,nbm:nb, : )**2,1) &
|
|
||||||
* pdata%q(ivy ,nb:ne,nbm:nb, : )) * dxz
|
|
||||||
inarr(28) = inarr(28) &
|
|
||||||
- sum(sum(pdata%q(ivx:ivz,nb:ne,nbm:nb, : ) &
|
|
||||||
* pdata%q(ibx:ibz,nb:ne,nbm:nb, : ),1) &
|
|
||||||
* pdata%q(iby ,nb:ne,nbm:nb, : )) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
if (resistivity > 0.0d+00) then
|
|
||||||
|
|
||||||
#if NDIMS == 3
|
|
||||||
! diffusion of Bx through the lower Y boundary
|
|
||||||
!
|
|
||||||
inarr(17) = inarr(17) &
|
|
||||||
+ sum(sign(cur(3,nb:ne,nbm:nb,nb:ne), &
|
|
||||||
pdata%q(ibx,nb:ne,nbm:nb,nb:ne))) * dxz
|
|
||||||
|
|
||||||
! diffusion of magnetic energy through the lower Y boundary
|
|
||||||
!
|
|
||||||
inarr(31) = inarr(31) &
|
|
||||||
+ sum(cur(3,nb:ne,nbm:nb,nb:ne) &
|
|
||||||
* pdata%q(ibx,nb:ne,nbm:nb,nb:ne) &
|
|
||||||
- cur(1,nb:ne,nbm:nb,nb:ne) &
|
|
||||||
* pdata%q(ibz,nb:ne,nbm:nb,nb:ne)) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
! diffusion of Bx through the lower Y boundary
|
|
||||||
!
|
|
||||||
inarr(17) = inarr(17) &
|
|
||||||
+ sum(sign(cur(3,nb:ne,nbm:nb, : ), &
|
|
||||||
pdata%q(ibx,nb:ne,nbm:nb, : ))) * dxz
|
|
||||||
|
|
||||||
! diffusion of magnetic energy through the lower Y boundary
|
|
||||||
!
|
|
||||||
inarr(31) = inarr(31) &
|
|
||||||
+ sum(cur(3,nb:ne,nbm:nb, : ) &
|
|
||||||
* pdata%q(ibx,nb:ne,nbm:nb, : ) &
|
|
||||||
- cur(1,nb:ne,nbm:nb, : ) &
|
|
||||||
* pdata%q(ibz,nb:ne,nbm:nb, : )) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
end if ! resistivity
|
|
||||||
|
|
||||||
end if ! ymin
|
|
||||||
|
|
||||||
! upper Y-boundary
|
|
||||||
!
|
|
||||||
if (pdata%meta%ymax >= ymax - dh(2)) then
|
|
||||||
|
|
||||||
! mean Bx at boundary
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(12) = inarr(12) &
|
|
||||||
+ sum(abs(pdata%q(ibx,nb:ne,ne:nep,nb:ne))) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(12) = inarr(12) &
|
|
||||||
+ sum(abs(pdata%q(ibx,nb:ne,ne:nep, : ))) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! advection of Bx along Y
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(13) = inarr(13) &
|
|
||||||
- sum(abs(pdata%q(ibx,nb:ne,ne:nep,nb:ne)) &
|
|
||||||
* pdata%q(ivy,nb:ne,ne:nep,nb:ne)) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(13) = inarr(13) &
|
|
||||||
- sum(abs(pdata%q(ibx,nb:ne,ne:nep, : )) &
|
|
||||||
* pdata%q(ivy,nb:ne,ne:nep, : )) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! shear of By along X
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(15) = inarr(15) &
|
|
||||||
+ sum(sign(pdata%q(iby,nb:ne,ne:nep,nb:ne) &
|
|
||||||
* pdata%q(ivx,nb:ne,ne:nep,nb:ne), &
|
|
||||||
pdata%q(ibx,nb:ne,ne:nep,nb:ne))) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(15) = inarr(15) &
|
|
||||||
+ sum(sign(pdata%q(iby,nb:ne,ne:nep, : ) &
|
|
||||||
* pdata%q(ivx,nb:ne,ne:nep, : ), &
|
|
||||||
pdata%q(ibx,nb:ne,ne:nep, : ))) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! mean magnetic energy at the upper Y boundary
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(23) = inarr(23) &
|
|
||||||
+ sum(pdata%q(ibx:ibz,nb:ne,ne:nep,nb:ne)**2) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(23) = inarr(23) &
|
|
||||||
+ sum(pdata%q(ibx:ibz,nb:ne,ne:nep, : )**2) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! advection of magnetic energy through the upper Y boundary
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(25) = inarr(25) &
|
|
||||||
- sum(sum(pdata%q(ibx:ibz,nb:ne,ne:nep,nb:ne)**2,1) &
|
|
||||||
* pdata%q(ivy ,nb:ne,ne:nep,nb:ne)) * dxz
|
|
||||||
inarr(28) = inarr(28) &
|
|
||||||
+ sum(sum(pdata%q(ivx:ivz,nb:ne,ne:nep,nb:ne) &
|
|
||||||
* pdata%q(ibx:ibz,nb:ne,ne:nep,nb:ne),1) &
|
|
||||||
* pdata%q(iby ,nb:ne,ne:nep,nb:ne)) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(25) = inarr(25) &
|
|
||||||
- sum(sum(pdata%q(ibx:ibz,nb:ne,ne:nep, : )**2,1) &
|
|
||||||
* pdata%q(ivy ,nb:ne,ne:nep, : )) * dxz
|
|
||||||
inarr(28) = inarr(28) &
|
|
||||||
+ sum(sum(pdata%q(ivx:ivz,nb:ne,ne:nep, : ) &
|
|
||||||
* pdata%q(ibx:ibz,nb:ne,ne:nep, : ),1) &
|
|
||||||
* pdata%q(iby ,nb:ne,ne:nep, : )) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
if (resistivity > 0.0d+00) then
|
|
||||||
|
|
||||||
#if NDIMS == 3
|
|
||||||
! diffusion of Bx through the upper Y boundary
|
|
||||||
!
|
|
||||||
inarr(17) = inarr(17) &
|
|
||||||
- sum(sign(cur(3,nb:ne,ne:nep,nb:ne), &
|
|
||||||
pdata%q(ibx,nb:ne,ne:nep,nb:ne))) * dxz
|
|
||||||
|
|
||||||
! diffusion of magnetic energy through the upper Y boundary
|
|
||||||
!
|
|
||||||
inarr(31) = inarr(31) &
|
|
||||||
- sum(cur(3,nb:ne,ne:nep,nb:ne) &
|
|
||||||
* pdata%q(ibx,nb:ne,ne:nep,nb:ne) &
|
|
||||||
- cur(1,nb:ne,ne:nep,nb:ne) &
|
|
||||||
* pdata%q(ibz,nb:ne,ne:nep,nb:ne)) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
! diffusion of Bx through the upper Y boundary
|
|
||||||
!
|
|
||||||
inarr(17) = inarr(17) &
|
|
||||||
- sum(sign(cur(3,nb:ne,ne:nep, : ), &
|
|
||||||
pdata%q(ibx,nb:ne,ne:nep, : ))) * dxz
|
|
||||||
|
|
||||||
! diffusion of magnetic energy through the upper Y boundary
|
|
||||||
!
|
|
||||||
inarr(31) = inarr(31) &
|
|
||||||
- sum(cur(3,nb:ne,ne:nep, : ) &
|
|
||||||
* pdata%q(ibx,nb:ne,ne:nep, : ) &
|
|
||||||
- cur(1,nb:ne,ne:nep, : ) &
|
|
||||||
* pdata%q(ibz,nb:ne,ne:nep, : )) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
end if ! resistivity
|
|
||||||
|
|
||||||
end if ! ymax
|
|
||||||
|
|
||||||
#if NDIMS == 3
|
|
||||||
! fluxes across the Z boundaries
|
|
||||||
!
|
|
||||||
if (.not. periodic(3)) then
|
|
||||||
|
|
||||||
! lower Z-boundary
|
|
||||||
!
|
|
||||||
if (pdata%meta%zmin <= zmin + dh(3)) then
|
|
||||||
|
|
||||||
! advection of Bx along Z
|
|
||||||
!
|
|
||||||
inarr(14) = inarr(14) &
|
|
||||||
+ sum(abs(pdata%q(ibx,nb:ne,nb:ne,nbm:nb)) &
|
|
||||||
* pdata%q(ivz,nb:ne,nb:ne,nbm:nb)) * dxy
|
|
||||||
|
|
||||||
! shear of Bz along X
|
|
||||||
!
|
|
||||||
inarr(16) = inarr(16) &
|
|
||||||
- sum(sign(pdata%q(ibz,nb:ne,nb:ne,nbm:nb) &
|
|
||||||
* pdata%q(ivx,nb:ne,nb:ne,nbm:nb), &
|
|
||||||
pdata%q(ibx,nb:ne,nb:ne,nbm:nb))) * dxy
|
|
||||||
|
|
||||||
! advection of magnetic energy through the lower Z boundary
|
|
||||||
!
|
|
||||||
inarr(26) = inarr(26) &
|
|
||||||
+ sum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nbm:nb)**2,1) &
|
|
||||||
* pdata%q(ivz ,nb:ne,nb:ne,nbm:nb)) * dxy
|
|
||||||
inarr(29) = inarr(29) &
|
|
||||||
- sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nbm:nb) &
|
|
||||||
* pdata%q(ibx:ibz,nb:ne,nb:ne,nbm:nb),1) &
|
|
||||||
* pdata%q(ibz ,nb:ne,nb:ne,nbm:nb)) * dxy
|
|
||||||
|
|
||||||
if (resistivity > 0.0d+00) then
|
|
||||||
|
|
||||||
! diffusion of Bx through the lower Z boundary
|
|
||||||
!
|
|
||||||
inarr(18) = inarr(18) &
|
|
||||||
+ sum(sign(cur(2,nb:ne,nb:ne,nbm:nb), &
|
|
||||||
pdata%q(ibx,nb:ne,nb:ne,nbm:nb))) * dxy
|
|
||||||
|
|
||||||
! diffusion of magnetic energy through the lower Z boundary
|
|
||||||
!
|
|
||||||
inarr(32) = inarr(32) &
|
|
||||||
+ sum(cur(1,nb:ne,nb:ne,nbm:nb) &
|
|
||||||
* pdata%q(iby,nb:ne,nb:ne,nbm:nb) &
|
|
||||||
- cur(2,nb:ne,nb:ne,nbm:nb) &
|
|
||||||
* pdata%q(ibx,nb:ne,nb:ne,nbm:nb)) * dxy
|
|
||||||
|
|
||||||
end if ! resistivity
|
|
||||||
|
|
||||||
end if ! zmin
|
|
||||||
|
|
||||||
! upper Z-boundary
|
|
||||||
!
|
|
||||||
if (pdata%meta%zmax >= zmax - dh(3)) then
|
|
||||||
|
|
||||||
! advection of Bx along Z
|
|
||||||
!
|
|
||||||
inarr(14) = inarr(14) &
|
|
||||||
- sum(abs(pdata%q(ibx,nb:ne,nb:ne,ne:nep)) &
|
|
||||||
* pdata%q(ivz,nb:ne,nb:ne,ne:nep)) * dxy
|
|
||||||
|
|
||||||
! shear of Bz along X
|
|
||||||
!
|
|
||||||
inarr(16) = inarr(16) &
|
|
||||||
+ sum(sign(pdata%q(ibz,nb:ne,nb:ne,ne:nep) &
|
|
||||||
* pdata%q(ivx,nb:ne,nb:ne,ne:nep), &
|
|
||||||
pdata%q(ibx,nb:ne,nb:ne,ne:nep))) * dxy
|
|
||||||
|
|
||||||
! advection of magnetic energy through the upper Z boundary
|
|
||||||
!
|
|
||||||
inarr(26) = inarr(26) &
|
|
||||||
- sum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,ne:nep)**2,1) &
|
|
||||||
* pdata%q(ivz ,nb:ne,nb:ne,ne:nep)) * dxy
|
|
||||||
inarr(29) = inarr(29) &
|
|
||||||
+ sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,ne:nep) &
|
|
||||||
* pdata%q(ibx:ibz,nb:ne,nb:ne,ne:nep),1) &
|
|
||||||
* pdata%q(ibz ,nb:ne,nb:ne,ne:nep)) * dxy
|
|
||||||
|
|
||||||
if (resistivity > 0.0d+00) then
|
|
||||||
|
|
||||||
! diffusion of Bx through the upper Z boundary
|
|
||||||
!
|
|
||||||
inarr(18) = inarr(18) &
|
|
||||||
- sum(sign(cur(2,nb:ne,nb:ne,ne:nep), &
|
|
||||||
pdata%q(ibx,nb:ne,nb:ne,ne:nep))) * dxy
|
|
||||||
|
|
||||||
! diffusion of magnetic energy through the upper Z boundary
|
|
||||||
!
|
|
||||||
inarr(32) = inarr(32) &
|
|
||||||
- sum(cur(1,nb:ne,nb:ne,ne:nep) &
|
|
||||||
* pdata%q(iby,nb:ne,nb:ne,ne:nep) &
|
|
||||||
- cur(2,nb:ne,nb:ne,ne:nep) &
|
|
||||||
* pdata%q(ibx,nb:ne,nb:ne,ne:nep)) * dxy
|
|
||||||
|
|
||||||
end if ! resistivity
|
|
||||||
|
|
||||||
end if ! zmax
|
|
||||||
|
|
||||||
end if ! not periodic in Z
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! get the conversion between the magnetic energy and kinetic and
|
|
||||||
! internal energies
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(33) = inarr(33) &
|
|
||||||
+ 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)) &
|
|
||||||
* cur(1,nb:ne,nb:ne,nb:ne)) * dvol
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(33) = inarr(33) &
|
|
||||||
+ 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, : )) &
|
|
||||||
* cur(1,nb:ne,nb:ne, : )) * dvol
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(33) = inarr(33) &
|
|
||||||
+ 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)) &
|
|
||||||
* cur(2,nb:ne,nb:ne,nb:ne)) * dvol
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(33) = inarr(33) &
|
|
||||||
+ 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, : )) &
|
|
||||||
* cur(2,nb:ne,nb:ne, : )) * dvol
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(33) = inarr(33) &
|
|
||||||
+ 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)) &
|
|
||||||
* cur(3,nb:ne,nb:ne,nb:ne)) * dvol
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(33) = inarr(33) &
|
|
||||||
+ 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, : )) &
|
|
||||||
* cur(3,nb:ne,nb:ne, : )) * dvol
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
if (resistivity > 0.0d+00) then
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(34) = inarr(34) - sum(cur(1:3,nb:ne,nb:ne,nb:ne)**2) * dvol
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(34) = inarr(34) - sum(cur(1:3,nb:ne,nb:ne, : )**2) * dvol
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
end if
|
|
||||||
|
|
||||||
! calculate gradient (∇ψ)
|
|
||||||
!
|
|
||||||
call gradient(dh(:), pdata%q(ibp,:,:,:), cur(1:3,:,:,:))
|
|
||||||
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(35) = inarr(35) &
|
|
||||||
- sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) &
|
|
||||||
* cur(1:3,nb:ne,nb:ne,nb:ne)) * dvol
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(35) = inarr(35) &
|
|
||||||
- sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) &
|
|
||||||
* cur(1:3,nb:ne,nb:ne, : )) * dvol
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
! contribution to |Bx| from ∇ψ
|
|
||||||
!
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(19) = inarr(19) &
|
|
||||||
- sum(sign(cur(1,nb:ne,nb:ne,nb:ne), &
|
|
||||||
pdata%q(ibx,nb:ne,nb:ne,nb:ne))) * dvol
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(19) = inarr(19) &
|
|
||||||
- sum(sign(cur(1,nb:ne,nb:ne, : ), &
|
|
||||||
pdata%q(ibx,nb:ne,nb:ne, : ))) * dvol
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
end if ! magnetized
|
|
||||||
|
|
||||||
! get the inflow speed
|
|
||||||
!
|
|
||||||
if (pdata%meta%ymin <= ymin + dh(2)) then
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(20) = inarr(20) + sum(pdata%q(ivy,nb:ne,nb,nb:ne)) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(20) = inarr(20) + sum(pdata%q(ivy,nb:ne,nb, : )) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
end if
|
|
||||||
if (pdata%meta%ymax >= ymax - dh(2)) then
|
|
||||||
#if NDIMS == 3
|
|
||||||
inarr(21) = inarr(21) - sum(pdata%q(ivy,nb:ne,ne,nb:ne)) * dxz
|
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
inarr(21) = inarr(21) - sum(pdata%q(ivy,nb:ne,ne, : )) * dxz
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
end if
|
end if
|
||||||
|
|
||||||
! associate the pointer with the next block on the list
|
! associate the pointer with the next block on the list
|
||||||
@ -1255,16 +638,6 @@ module integrals
|
|||||||
!
|
!
|
||||||
avarr(:) = avarr(:) * voli
|
avarr(:) = avarr(:) * voli
|
||||||
|
|
||||||
! apply factors to the reconnection rate terms
|
|
||||||
!
|
|
||||||
inarr(12) = 2.50d-01 * inarr(12) / yarea
|
|
||||||
inarr(13:16) = 5.00d-01 * inarr(13:16)
|
|
||||||
inarr(17:18) = 5.00d-01 * resistivity * inarr(17:18)
|
|
||||||
inarr(23) = 1.25d-01 * inarr(23) / yarea
|
|
||||||
inarr(24:29) = 5.00d-01 * inarr(24:29)
|
|
||||||
inarr(30:32) = 5.00d-01 * resistivity * inarr(30:32)
|
|
||||||
inarr(34) = resistivity * inarr(34)
|
|
||||||
|
|
||||||
! write down the integrals and statistics to appropriate files
|
! write down the integrals and statistics to appropriate files
|
||||||
!
|
!
|
||||||
if (stored) then
|
if (stored) then
|
||||||
@ -1278,12 +651,10 @@ 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, dth, dte, maxval(errors(:)), errors(:)
|
write(eunit,efmt) step, time, dth, dte, maxval(errors(:)), errors(:)
|
||||||
write(runit,"(i9, 26es25.15e3)") step, time, inarr(11:35)
|
|
||||||
|
|
||||||
call flush_and_sync(funit)
|
call flush_and_sync(funit)
|
||||||
call flush_and_sync(sunit)
|
call flush_and_sync(sunit)
|
||||||
call flush_and_sync(eunit)
|
call flush_and_sync(eunit)
|
||||||
call flush_and_sync(runit)
|
|
||||||
end if
|
end if
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
|
|||||||
@ -372,10 +372,10 @@ module user_problem
|
|||||||
call get_parameter("reconnection_append", append)
|
call get_parameter("reconnection_append", append)
|
||||||
select case(trim(append))
|
select case(trim(append))
|
||||||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||||||
write(fname, "('reconnection-new.dat')")
|
write(fname, "('reconnection.dat')")
|
||||||
inquire(file = fname, exist = flag)
|
inquire(file = fname, exist = flag)
|
||||||
case default
|
case default
|
||||||
write(fname, "('reconnection-new_',i2.2,'.dat')") rcount
|
write(fname, "('reconnection_',i2.2,'.dat')") rcount
|
||||||
flag = .false.
|
flag = .false.
|
||||||
end select
|
end select
|
||||||
|
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user