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

USER_PROBLEM: Implement reconnection rate measure.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-11-04 12:22:45 -03:00
parent 52d6933efb
commit 711fc23f03
1 changed files with 68 additions and 223 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

291
sources/user_problem.F90
View File

@ -153,10 +153,10 @@ module user_problem
call get_parameter("statistics_append", append)
select case(trim(append))
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
write(fname, "('energy_balance.dat')")
write(fname, "('reconnection.dat')")
inquire(file = fname, exist = flag)
case default
write(fname, "('energy_balance_',i2.2,'.dat')") rcount
write(fname, "('reconnection_',i2.2,'.dat')") rcount
flag = .false.
end select
@ -183,11 +183,11 @@ module user_problem
! write the integral file header
!
write(runit,'("#",a24,13a25)') 'time', 'Emag total', &
'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,'("#",a24,9a25)') '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'
write(runit,"('#')")
call print_section(verbose, "Parameters (* - set, otherwise calculated)")
@ -572,7 +572,7 @@ module user_problem
#ifdef MPI
use mpitools , only : reduce_sum
#endif /* MPI */
use operators , only : curl, gradient
use operators , only : curl, derivative_1st
use workspace , only : resize_workspace, work, work_in_use
implicit none
@ -681,108 +681,11 @@ module user_problem
!
rterms(1) = rterms(1) &
#if NDIMS == 3
+ sum(pdata%u(ibx:ibz,nb:ne,nb:ne,nb:ne)**2) * dvol
+ sum(abs(pdata%u(ibx,nb:ne,nb:ne,nb:ne))) * dvol
#else /* NDIMS == 3 */
+ sum(pdata%u(ibx:ibz,nb:ne,nb:ne, : )**2) * dvol
+ sum(abs(pdata%u(ibx,nb:ne,nb:ne, : ))) * dvol
#endif /* NDIMS == 3 */
if (.not. periodic(1)) then
if (pmeta%xmin < xmin) then
! interpolate velocity and magnetic field components
!
#if NDIMS == 3
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nbm:nb,nb:ne,nb:ne), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nbm:nb,nb:ne,nb:ne), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nbm:nb,nb:ne,nb:ne), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nbm:nb,nb:ne,nb:ne), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nbm:nb,nb:ne,nb:ne), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nbm:nb,nb:ne,nb:ne), 1)
#else /* NDIMS == 3 */
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nbm:nb,nb:ne, : ), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nbm:nb,nb:ne, : ), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nbm:nb,nb:ne, : ), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nbm:nb,nb:ne, : ), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nbm:nb,nb:ne, : ), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nbm:nb,nb:ne, : ), 1)
#endif /* NDIMS == 3 */
! advection of magnetic energy across the lower X-boundary
!
rterms(2) = rterms(2) + sum(sum(bb * bb, 1) * vv(1,:,:)) * dyz
! generation of magnetic energy at the lower X-boundary
!
rterms(5) = rterms(5) - sum(sum(vv * bb, 1) * bb(1,:,:)) * dyz
if (resistive) then
! interpolate the components of current density
!
#if NDIMS == 3
jj(2,:,:) = 5.0d-01 * sum(cur(2,nbm:nb,nb:ne,nb:ne), 1)
jj(3,:,:) = 5.0d-01 * sum(cur(3,nbm:nb,nb:ne,nb:ne), 1)
#else /* NDIMS == 3 */
jj(2,:,:) = 5.0d-01 * sum(cur(2,nbm:nb,nb:ne, : ), 1)
jj(3,:,:) = 5.0d-01 * sum(cur(3,nbm:nb,nb:ne, : ), 1)
#endif /* NDIMS == 3 */
! diffusion of magnetic energy across the lower X-boundary
!
rterms(8) = rterms(8) &
+ sum(jj(2,:,:) * bb(3,:,:) - jj(3,:,:) * bb(2,:,:)) * dyz
end if ! resistivity
end if
if (pmeta%xmax > xmax) then
#if NDIMS == 3
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,ne:nep,nb:ne,nb:ne), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,ne:nep,nb:ne,nb:ne), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,ne:nep,nb:ne,nb:ne), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,ne:nep,nb:ne,nb:ne), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,ne:nep,nb:ne,nb:ne), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,ne:nep,nb:ne,nb:ne), 1)
#else /* NDIMS == 3 */
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,ne:nep,nb:ne, : ), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,ne:nep,nb:ne, : ), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,ne:nep,nb:ne, : ), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,ne:nep,nb:ne, : ), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,ne:nep,nb:ne, : ), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,ne:nep,nb:ne, : ), 1)
#endif /* NDIMS == 3 */
! advection of magnetic energy across the upper X-boundary
!
rterms(2) = rterms(2) - sum(sum(bb * bb, 1) * vv(1,:,:)) * dyz
! generation of magnetic energy at the upper X-boundary
!
rterms(5) = rterms(5) + sum(sum(vv * bb, 1) * bb(1,:,:)) * dyz
if (resistive) then
! interpolate the components of current density
!
#if NDIMS == 3
jj(2,:,:) = 5.0d-01 * sum(cur(2,ne:nep,nb:ne,nb:ne), 1)
jj(3,:,:) = 5.0d-01 * sum(cur(3,ne:nep,nb:ne,nb:ne), 1)
#else /* NDIMS == 3 */
jj(2,:,:) = 5.0d-01 * sum(cur(2,ne:nep,nb:ne, : ), 1)
jj(3,:,:) = 5.0d-01 * sum(cur(3,ne:nep,nb:ne, : ), 1)
#endif /* NDIMS == 3 */
! diffusion of magnetic energy across the upper X-boundary
!
rterms(8) = rterms(8) &
- sum(jj(2,:,:) * bb(3,:,:) - jj(3,:,:) * bb(2,:,:)) * dyz
end if ! resistivity
end if
end if ! non-periodic along X
if (.not. periodic(2)) then
if (pmeta%ymin < ymin) then
@ -790,43 +693,43 @@ module user_problem
#if NDIMS == 3
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,nbm:nb,nb:ne), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,nbm:nb,nb:ne), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,nbm:nb,nb:ne), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,nbm:nb,nb:ne), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,nbm:nb,nb:ne), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nbm:nb,nb:ne), 1)
#else /* NDIMS == 3 */
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,nbm:nb, : ), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,nbm:nb, : ), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,nbm:nb, : ), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,nbm:nb, : ), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,nbm:nb, : ), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nbm:nb, : ), 1)
#endif /* NDIMS == 3 */
! advection of magnetic energy across the lower Y-boundary
! mean Bx at the Y-boundary
!
rterms(3) = rterms(3) + sum(sum(bb * bb, 1) * vv(2,:,:)) * dxz
rterms(2) = rterms(2) + sum(abs(bb(1,:,:))) * dxz
! generation of magnetic energy at the lower Y-boundary
! advection of Bx across the Y-boundary
!
rterms(6) = rterms(6) - sum(sum(vv * bb, 1) * bb(2,:,:)) * dxz
rterms(3) = rterms(3) &
+ sum(sign(bb(1,:,:) * vv(2,:,:), bb(1,:,:))) * dxz
! shear of By along the Y-boundary
!
rterms(5) = rterms(5) &
- sum(sign(bb(2,:,:) * vv(1,:,:), bb(1,:,:))) * dxz
if (resistive) then
! interpolate the components of current density
! interpolate Jz
!
#if NDIMS == 3
jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,nbm:nb,nb:ne), 1)
jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,nbm:nb,nb:ne), 1)
#else /* NDIMS == 3 */
jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,nbm:nb, : ), 1)
jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,nbm:nb, : ), 1)
#endif /* NDIMS == 3 */
! diffusion of magnetic energy through the lower Y boundary
! diffusion of Bx at the Y-boundary
!
rterms(9) = rterms(9) &
+ sum(jj(3,:,:) * bb(1,:,:) - jj(1,:,:) * bb(3,:,:)) * dxz
rterms(7) = rterms(7) &
+ sum(sign(jj(3,:,:), bb(1,:,:))) * dxz
end if ! resistivity
end if
@ -836,43 +739,43 @@ module user_problem
#if NDIMS == 3
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,ne:nep,nb:ne), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,ne:nep,nb:ne), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,ne:nep,nb:ne), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,ne:nep,nb:ne), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,ne:nep,nb:ne), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,ne:nep,nb:ne), 1)
#else /* NDIMS == 3 */
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,ne:nep, : ), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,ne:nep, : ), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,ne:nep, : ), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,ne:nep, : ), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,ne:nep, : ), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,ne:nep, : ), 1)
#endif /* NDIMS == 3 */
! advection of magnetic energy across the upper Y-boundary
! mean Bx at the boundary
!
rterms(3) = rterms(3) - sum(sum(bb * bb, 1) * vv(2,:,:)) * dxz
rterms(2) = rterms(2) + sum(abs(bb(1,:,:))) * dxz
! generation of magnetic energy at the upper Y-boundary
! advection of Bx across the Y-boundary
!
rterms(6) = rterms(6) + sum(sum(vv * bb, 1) * bb(2,:,:)) * dxz
rterms(3) = rterms(3) &
- sum(sign(bb(1,:,:) * vv(2,:,:), bb(1,:,:))) * dxz
! shear of By along the Y-boundary
!
rterms(5) = rterms(5) &
+ sum(sign(bb(2,:,:) * vv(1,:,:), bb(1,:,:))) * dxz
if (resistive) then
! interpolate the components of current density
! interpolate Jz
!
#if NDIMS == 3
jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,ne:nep,nb:ne), 1)
jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,ne:nep,nb:ne), 1)
#else /* NDIMS == 3 */
jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,ne:nep, : ), 1)
jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,ne:nep, : ), 1)
#endif /* NDIMS == 3 */
! diffusion of magnetic energy
! diffusion of Bx at the Y-boundary
!
rterms(9) = rterms(9) &
- sum(jj(3,:,:) * bb(1,:,:) - jj(1,:,:) * bb(3,:,:)) * dxz
rterms(7) = rterms(7) &
- sum(sign(jj(3,:,:), bb(1,:,:))) * dxz
end if ! resistivity
end if
@ -886,31 +789,30 @@ module user_problem
! interpolate velocity and magnetic field components
!
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,nb:ne,nbm:nb), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,nb:ne,nbm:nb), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,nb:ne,nbm:nb), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,nb:ne,nbm:nb), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,nb:ne,nbm:nb), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nb:ne,nbm:nb), 1)
! advection of magnetic energy across the lower Z-boundary
! advection of Bx across the Z-boundary
!
rterms(4) = rterms(4) + sum(sum(bb * bb, 1) * vv(3,:,:)) * dxy
rterms(4) = rterms(4) &
+ sum(sign(bb(1,:,:) * vv(3,:,:), bb(1,:,:))) * dxy
! generation of magnetic energy at the lower Z-boundary
! shear of Bz along the Z-boundary
!
rterms(7) = rterms(7) - sum(sum(vv * bb, 1) * bb(3,:,:)) * dxy
rterms(6) = rterms(6) &
- sum(sign(bb(3,:,:) * vv(1,:,:), bb(1,:,:))) * dxy
if (resistive) then
! interpolate the components of current density
! interpolate Jy
!
jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,nb:ne,nbm:nb), 1)
jj(2,:,:) = 5.0d-01 * sum(cur(2,nb:ne,nb:ne,nbm:nb), 1)
! diffusion of magnetic energy across the lower Z-boundary
! diffusion of Bx at the Z-boundary
!
rterms(10) = rterms(10) &
+ sum(jj(1,:,:) * bb(2,:,:) - jj(2,:,:) * bb(1,:,:)) * dxy
rterms(8) = rterms(8) &
+ sum(sign(jj(2,:,:), bb(1,:,:))) * dxy
end if ! resistivity
end if
@ -918,104 +820,47 @@ module user_problem
if (pmeta%zmax > zmax) then
vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,nb:ne,ne:nep), 1)
vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,nb:ne,ne:nep), 1)
vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,nb:ne,ne:nep), 1)
bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,nb:ne,ne:nep), 1)
bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,nb:ne,ne:nep), 1)
bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nb:ne,ne:nep), 1)
! advection of magnetic energy across the upper Z-boundary
! advection of Bx across the Z-boundary
!
rterms(4) = rterms(4) - sum(sum(bb * bb, 1) * vv(3,:,:)) * dxy
rterms(4) = rterms(4) &
- sum(sign(bb(1,:,:) * vv(3,:,:), bb(1,:,:))) * dxy
! generation of magnetic energy at the upper Z-boundary
! shear of Bz along the Z-boundary
!
rterms(7) = rterms(7) + sum(sum(vv * bb, 1) * bb(3,:,:)) * dxy
rterms(6) = rterms(6) &
+ sum(sign(bb(3,:,:) * vv(1,:,:), bb(1,:,:))) * dxy
if (resistive) then
! interpolate the components of current density
!
jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,nb:ne,ne:nep), 1)
jj(2,:,:) = 5.0d-01 * sum(cur(2,nb:ne,nb:ne,ne:nep), 1)
! diffusion of magnetic energy across the upper Z-boundary
! diffusion of Bx at the Z-boundary
!
rterms(10) = rterms(10) &
- sum(jj(1,:,:) * bb(2,:,:) - jj(2,:,:) * bb(1,:,:)) * dxy
rterms(8) = rterms(8) &
- sum(sign(jj(2,:,:), bb(1,:,:))) * dxy
end if ! resistivity
end if
end if ! non-periodic along Z
#endif /* NDIMS == 3 */
! conversion between magnetic and kinetic energies
! calculate X-derivative of ψ
!
rterms(11) = rterms(11) &
#if NDIMS == 3
+ 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 */
+ 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 */
rterms(11) = rterms(11) &
#if NDIMS == 3
+ 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 */
+ 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 */
rterms(11) = rterms(11) &
#if NDIMS == 3
+ 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 */
+ 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 */
call derivative_1st(1, dh(1), pdata%q(ibp,:,:,:), cur(1,:,:,:))
! conversion between magnetic energy and internal energies
!
if (resistive) then
rterms(12) = rterms(12) &
rterms(9) = rterms(9) &
#if NDIMS == 3
- sum(cur(1:3,nb:ne,nb:ne,nb:ne)**2) * dvol
- sum(sign(cur(1,nb:ne,nb:ne,nb:ne), &
pdata%q(ibx,nb:ne,nb:ne,nb:ne))) * dvol
#else /* NDIMS == 3 */
- 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,:,:,:))
rterms(13) = rterms(13) &
#if NDIMS == 3
- sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) &
* cur(1:3,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
- sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) &
* cur(1:3,nb:ne,nb:ne, : )) * dvol
- sum(sign(cur(1,nb:ne,nb:ne, : ), &
pdata%q(ibx,nb:ne,nb:ne, : ))) * dvol
#endif /* NDIMS == 3 */
pdata => pdata%next
@ -1027,11 +872,11 @@ module user_problem
call reduce_sum(rterms(:))
#endif /* MPI */
rterms(1) = 5.00d-01 * rterms(1)
rterms(8:10) = eta * rterms(8:10)
rterms(12) = eta * rterms(12)
rterms(2) = 5.00d-01 * rterms(2)
rterms(7) = eta * rterms(7)
rterms(8) = eta * rterms(8)
write(runit,"(14es25.15e3)") time, rterms(1:13)
write(runit,"(10es25.15e3)") time, rterms(1:9)
call flush_and_sync(runit)
100 continue