! $Id$
!
! This module solves the Poisson equation
! (d^2/dx^2 + d^2/dy^2 + d^2/dz^2 - h) f = RHS(x,y,z)
! [which for h/=0 could also be called inhomogenous nonuniform Helmholtz
! equation] for the function f(x,y,z).
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lpoisson=.true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
!
!***************************************************************
module Poisson
!
use Cdata
use Cparam
use Fourier
use Messages
!
implicit none
!
real :: kmax=0.0
logical :: lrazor_thin=.false., lsemispectral=.false., lklimit_shear=.false.
logical :: lexpand_grid=.false.
logical :: lisoz=.false.
!
include 'poisson.h'
!
namelist /poisson_init_pars/ &
lsemispectral, kmax, lrazor_thin, lklimit_shear, lexpand_grid, lisoz
!
namelist /poisson_run_pars/ &
lsemispectral, kmax, lrazor_thin, lklimit_shear, lexpand_grid, lisoz
!
logical :: luse_fourier_transform = .false.
!
contains
!***********************************************************************
subroutine initialize_poisson
!
! Perform any post-parameter-read initialization i.e. calculate derived
! parameters.
!
! 18-oct-07/anders: adapted
!
if (lrazor_thin) then
if (nzgrid/=1) then
if (lroot) print*, 'inverse_laplacian_fft: razor-thin approximation only works with nzgrid==1'
call fatal_error('inverse_laplacian_fft','')
endif
endif
!
! Limit the wavenumber to the maximum circular region that is always available
! in k-space. The radial wavenumber kx changes with time due to shear as
!
! kx = kx0+qshear*Omega*t*ky
!
! Considering the available (kx,ky) space, it turns slowly from a square to a
! parallellogram (the hole for kx,ky<1 is ignored here):
!
! - - - - - - - -
! | | / /
! | | -> / /
! | | / /
! | | / /
! - - - - - - - -
!
! To make the gravity force isotropic at small scales one can limit k to
! the largest circular region that is present in both the square and the
! parallellogram. The circle has radius kmax=kNy/sqrt(2). See Gammie (2001).
!
if (lklimit_shear) kmax = kx_nyq/sqrt(2.0)
!
! Dimensionality
!
call decide_fourier_routine
!
endsubroutine initialize_poisson
!***********************************************************************
subroutine inverse_laplacian(phi)
!
! Dispatch solving the Poisson equation to inverse_laplacian_fft
! or inverse_laplacian_semispectral, based on the boundary conditions
!
! 17-jul-2007/wolf: coded wrapper
!
real, dimension(nx,ny,nz), intent(inout) :: phi
!
if (lcylindrical_coords) then
if (lroot) print*,'You are using cylindrical coordinates. '//&
'Use poisson_multigrid.f90 instead'
call fatal_error("inverse_laplacian","")
endif
!
if (lspherical_coords) then
if (lroot) then
print*,'There is no poisson solver for spherical '
print*,'coordinates yet. Please feel free to implement it. '
print*,'Many people will thank you for that.'
call fatal_error("inverse_laplacian","")
endif
endif
!
if (lsemispectral) then
call inverse_laplacian_semispectral(phi)
else if (lisoz) then
call inverse_laplacian_isoz(phi)
else if (lexpand_grid) then
if (lroot) then
print*,'The poisson solver for global disk was moved to an '
print*,'experimental module because it uses too much memory. '
print*,'If your memory requirements are modest ( < 500 proc) '
print*,'use POISSON=experimental/poisson_expand.f90 in your '
print*,'Makefile.local file'
call fatal_error("inverse_laplacian","")
endif
!call inverse_laplacian_expandgrid(phi)
else if (nprocx==1.and.(nxgrid/=1.and.nygrid/=1.and.nzgrid/=1.and.nxgrid/=nzgrid)) then
call inverse_laplacian_fft_z(phi)
else
call inverse_laplacian_fft(phi)
endif
!
endsubroutine inverse_laplacian
!***********************************************************************
subroutine inverse_laplacian_fft(phi)
!
! Solve the Poisson equation by Fourier transforming on a periodic grid.
! This method works both with and without shear.
!
! 15-may-2006/anders+jeff: coded
!
real, dimension (nx,ny,nz) :: phi
!
real, dimension (nx,ny,nz) :: b1
real :: k2
integer :: ikx, iky, ikz
!
! Identify version.
!
if (lroot .and. ip<10) call svn_id( &
"$Id$")
!
! The right-hand-side of the Poisson equation is purely real.
!
b1 = 0.0
!
! Forward transform (to k-space).
!
if (luse_fourier_transform) then
if (lshear) then
call fourier_transform_shear(phi,b1)
else
call fourier_transform(phi,b1)
endif
else
if (.not.(nxgrid/=1.and.nzgrid/=1.and.nxgrid/=nzgrid)) then
call fft_xyz_parallel(phi,b1)
else
call fatal_error("inverse_laplacian_fft",&
"not yet coded for the general case nxgrid/=nzgrid")
endif
endif
!
! Solve Poisson equation.
!
do ikz=1,nz; do iky=1,ny; do ikx=1,nx
if (.not. lrazor_thin) then
if (lshear) then
!
! Take into account that the Fourier transform has been done in shearing
! coordinates, and that the kx of each Fourier mode is different in the normal
! frame. The connection between kx0 (in the shearing frame) and the actual kx
! is
! kx = kx0 + qshear*Omega*t*ky
! Writing deltay/2 = qshear*Omega*Lx/2*t, one gets the expression
! kx = kx0 + deltay/Lx*ky
! The parallel Fourier transform returns the array in a tranposed state, so
! must be able to identify the x-direction in order to take shear into account.
! (see the subroutine transform_fftpack_shear in Mpicomm for details).
!
if (nzgrid/=1) then ! Order (kz,ky',kx)
k2 = (kx_fft(ikz+ipz*nz)+deltay/Lx*ky_fft(iky+ipy*ny))**2 + &
ky_fft(iky+ipy*ny)**2 + kz_fft(ikx+ipx*nx)**2
else ! Order (kx,ky',kz)
k2 = (kx_fft(ikx+ipx*nx)+deltay/Lx*ky_fft(iky+ipy*ny))**2 + &
ky_fft(iky+ipy*ny)**2 + kz_fft(ikz+ipz*nz)**2
endif
! The ordering of the array is not important here, because there is no shear!
else
k2 = kx_fft(ikx+ipx*nx)**2 + ky_fft(iky+ipy*ny)**2 + kz_fft(ikz+ipz*nz)**2
endif
!
! Solution of Poisson equation.
!
if (k2 > 0.0) then
phi(ikx,iky,ikz) = -phi(ikx,iky,ikz) / k2
b1(ikx,iky,ikz) = - b1(ikx,iky,ikz) / k2
else
phi(ikx,iky,ikz) = 0.0
b1(ikx,iky,ikz) = 0.0
endif
!
else
!
! Razor-thin approximation. Here we solve the equation
! del2Phi=4*pi*G*Sigma(x,y)*delta(z)
! The solution at scale k=(kx,ky) is
! Phi(x,y,z)=-(2*pi*G/|k|)*Sigma(x,y)*exp[i*(kx*x+ky*y)-|k|*|z|]
!
if (lshear) then
k2 = (kx_fft(ikx+ipx*nx)+deltay/Lx*ky_fft(iky+ipy*ny))**2+ky_fft(iky+ipy*ny)**2
else
k2 = kx_fft(ikx+ipx*nx)**2+ky_fft(iky+ipy*ny)**2
endif
!
if (k2 > 0.0) then
phi(ikx,iky,ikz) = -.5*phi(ikx,iky,ikz) / sqrt(k2)
b1(ikx,iky,ikz) = -.5* b1(ikx,iky,ikz) / sqrt(k2)
else
phi(ikx,iky,ikz) = 0.0
b1(ikx,iky,ikz) = 0.0
endif
endif
!
! Limit |k| < kmax
!
if (kmax>0.0) then
if (sqrt(k2)>=kmax) then
phi(ikx,iky,ikz) = 0.
b1(ikx,iky,ikz) = 0.
endif
endif
enddo; enddo; enddo
!
! Inverse transform (to real space).
!
if (luse_fourier_transform) then
if (lshear) then
call fourier_transform_shear(phi,b1,linv=.true.)
else
call fourier_transform(phi,b1,linv=.true.)
endif
else
call fft_xyz_parallel(phi,b1,linv=.true.)
endif
!
endsubroutine inverse_laplacian_fft
!***********************************************************************
subroutine inverse_laplacian_semispectral(phi)
!
! Solve the Poisson equation by Fourier transforming in the xy-plane and
! solving the discrete matrix equation in the z-direction.
!
! 19-dec-2006/anders: coded
!
use General, only: tridag
use Mpicomm, only: transp_xz, transp_zx
!
real, dimension (nx,ny,nz) :: phi
!
real, dimension (nx,ny,nz) :: b1
real, dimension (nzgrid,nx/nprocz) :: rhst, b1t
real, dimension (nzgrid) :: a_tri, b_tri, c_tri, r_tri, u_tri
real :: k2
integer :: ikx, iky
logical :: err
!
! Identify version.
!
if (lroot .and. ip<10) call svn_id( &
'$Id$')
!
! The right-hand-side of the Poisson equation is purely real.
!
b1 = 0.0
!
! Forward transform (to k-space).
!
if (lshear) then
call fourier_transform_shear_xy(phi,b1)
else
call fourier_transform_xy(phi,b1)
endif
!
! Solve for discrete z-direction with zero density above and below z-boundary.
!
do iky=1,ny
call transp_xz(phi(:,iky,:),rhst)
a_tri(:)=1.0/dz**2
c_tri(:)=1.0/dz**2
do ikx=1,nxgrid/nprocz
k2=kx_fft(ikx+nz*ipz)**2+ky_fft(iky)**2
b_tri=-2.0/dz**2-k2
!
if (k2==0.0) then
b_tri(1)=-2.0/dz**2-2*dz/xyz0(3)
c_tri(1)=1.0/dz**2+1.0
b_tri(nzgrid)=-2.0/dz**2-2*dz/xyz1(3)
a_tri(nzgrid)=1.0/dz**2+1.0
else
b_tri(1)=-2.0/dz**2-2*sqrt(k2)*dz
c_tri(1)=1.0/dz**2+1.0
b_tri(nzgrid)=-2.0/dz**2-2*sqrt(k2)*dz
a_tri(nzgrid)=1.0/dz**2+1.0
endif
!
r_tri=rhst(:,ikx)
call tridag(a_tri,b_tri,c_tri,r_tri,u_tri,err)
rhst(:,ikx)=u_tri
enddo
call transp_zx(rhst,phi(:,iky,:))
enddo
!
do iky=1,ny
call transp_xz(b1(:,iky,:),b1t)
a_tri(:)=1.0/dz**2
c_tri(:)=1.0/dz**2
do ikx=1,nxgrid/nprocz
k2=kx_fft(ikx+nz*ipz)**2+ky_fft(iky)**2
b_tri=-2.0/dz**2-k2
!
if (k2==0.0) then
b_tri(1)=-2.0/dz**2-2*dz/xyz0(3)
c_tri(1)=1.0/dz**2+1.0
b_tri(nzgrid)=-2.0/dz**2-2*dz/xyz1(3)
a_tri(nzgrid)=1.0/dz**2+1.0
else
b_tri(1)=-2.0/dz**2-2*sqrt(k2)*dz
c_tri(1)=1.0/dz**2+1.0
b_tri(nzgrid)=-2.0/dz**2-2*sqrt(k2)*dz
a_tri(nzgrid)=1.0/dz**2+1.0
endif
!
r_tri=b1t(:,ikx)
call tridag(a_tri,b_tri,c_tri,r_tri,u_tri,err)
b1t(:,ikx)=u_tri
enddo
call transp_zx(b1t,b1(:,iky,:))
enddo
!
! Inverse transform (to real space).
!
if (lshear) then
call fourier_transform_shear_xy(phi,b1,linv=.true.)
else
call fourier_transform_xy(phi,b1,linv=.true.)
endif
!
endsubroutine inverse_laplacian_semispectral
!***********************************************************************
subroutine inverse_laplacian_fft_z(phi)
!
! Solve the Poisson equation with nxgrid = nygrid /= nzgrid.
!
! 10-sep-2009/ccyang: coded
!
use Mpicomm, only: transp_xz, transp_zx
!
real, dimension(nx,ny,nz), intent(inout) :: phi
!
real, dimension(nx,ny,nz) :: b1
!
integer, parameter :: nxt = nx / nprocz
real, dimension(nzgrid,nxt) :: phirt, b1t
!
logical :: l1st = .true.
real, save, dimension(4*nzgrid+15) :: wsave
real, save, dimension(nzgrid) :: kz2
integer, save :: ikx0, iky0
!
complex, dimension(nzgrid) :: cz
!
integer :: ix, iy
real :: kx2, ky2, a0, a1
!
! Initialize the array wsave and other constants for future use.
!
if (l1st) then
call cffti(nzgrid, wsave)
kz2 = kz_fft**2
ikx0 = ipz * nxt
iky0 = ipy * ny
l1st = .false.
endif
!
! Forward transform in xy
!
b1 = 0.
if (nprocx==1) then
if (nxgrid==nygrid) then
if (lshear) then
call fourier_transform_shear_xy(phi, b1)
else
call fourier_transform_xy(phi, b1)
!#ccyang: Note that x and y are swapped in fourier_transform_xy.
endif
else
call fft_xy_parallel(phi,b1)
endif
else
call fatal_error("inverse_laplacian_fft_z",&
"Does not work for x-parallelization")
endif
!
! Convolution in z
!
if (lshear) a0 = deltay / Lx
do iy = 1, ny
!
call transp_xz(phi(:,iy,:), phirt)
call transp_xz(b1(:,iy,:), b1t)
!
if (lshear) then
ky2 = ky_fft(iky0+iy)**2
a1 = a0 * ky_fft(iky0+iy)
else
ky2 = kx_fft(iky0+iy)**2
endif
!
do ix = 1, nxt
!
if (lshear) then
kx2 = (kx_fft(ikx0+ix) + a1)**2
else
kx2 = ky_fft(ikx0+ix)**2
endif
!
cz = cmplx(phirt(:,ix), b1t(:,ix))
call cfftf (nzgrid, cz, wsave)
where (kx2 /= 0. .or. ky2 /= 0. .or. kz2 /= 0)
cz = -cz / (kx2 + ky2 + kz2) / nzgrid
elsewhere
cz = 0.
endwhere
call cfftb (nzgrid, cz, wsave)
!
phirt(:,ix) = real(cz)
b1t(:,ix) = aimag(cz)
!
enddo
!
call transp_zx(phirt, phi(:,iy,:))
call transp_zx(b1t, b1(:,iy,:))
!
enddo
!
! Inverse transform in xy
!
if (nprocx==1) then
if (nxgrid==nygrid) then
if (lshear) then
call fourier_transform_shear_xy(phi,b1,linv=.true.)
else
call fourier_transform_xy(phi,b1,linv=.true.)
!#ccyang: Note that x and y are swapped in fourier_transform_xy.
endif
else
call fft_xy_parallel(phi,b1,linv=.true.)
endif
else
call fatal_error("inverse_laplacian_fft_z",&
"Does not work for x-parallelization")
endif
!
return
!
endsubroutine inverse_laplacian_fft_z
!***********************************************************************
subroutine inverse_laplacian_isoz(phi)
!
! Solve the Poisson equation with isolated boundary condition in the z-
! direction and periodic in the x- and y- directions.
!
! 11-jan-2008/ccyang: coded
!
! Reference: Yang, Mac Low, & Menou 2011 (arXiv:1103.3268)
!
use Mpicomm, only: transp_xz, transp_zx
!
real, dimension(nx,ny,nz) :: phi
!
real, dimension(nx,ny,nz) :: b1
!
integer, parameter :: nzg2 = 2 * nzgrid
integer, parameter :: nxt = nx / nprocz
real, dimension(nzgrid,nxt) :: phirt, b1t
!
logical, save :: l1st = .true.
real, save, dimension(4*nzg2+15) :: wsave
real, save, dimension(nzg2) :: kz2
real, save :: dz2h
integer, save :: ikx0, iky0
!
complex, dimension(nzg2) :: cz
!
integer :: ix, iy, iz, iz1
real :: kx2, ky2, a0, a1
!
! Initialize the array wsave and other constants for future use.
!
if (l1st) then
call cffti(nzg2, wsave)
kz2 = (/(iz**2, iz = 0, nzgrid), &
((iz - nzgrid)**2, iz = 1, nzgrid - 1)/) * (pi / Lz)**2
dz2h = .5 * dz * dz
ikx0 = ipz * nxt
iky0 = ipy * ny
l1st = .false.
endif
!
! Forward transform in xy
!
b1 = 0.
if (lshear) then
call fourier_transform_shear_xy(phi, b1)
else
call fourier_transform_xy(phi, b1)
!#ccyang: Note that x and y are swapped in fourier_transform_xy.
endif
!
! Convolution in z
!
if (lshear) a0 = deltay / Lx
do iy = 1, ny
call transp_xz(phi(:,iy,:), phirt)
call transp_xz(b1(:,iy,:), b1t)
if (lshear) then
ky2 = ky_fft(iky0+iy)**2
a1 = a0 * ky_fft(iky0+iy)
else
ky2 = kx_fft(iky0+iy)**2
endif
do ix = 1, nxt
if (lshear) then
kx2 = (kx_fft(ikx0+ix) + a1)**2
else
kx2 = ky_fft(ikx0+ix)**2
endif
if (kx2 /= 0. .or. ky2 /= 0.) then
cz(1:nzgrid) = (/cmplx(phirt(:,ix), b1t(:,ix))/)
cz(nzgrid+1:nzg2) = 0.
call cfftf (nzg2, cz, wsave)
cz = -cz / (kx2 + ky2 + kz2) / nzg2
call cfftb (nzg2, cz, wsave)
else
cz = 0.
do iz = 1, nzgrid; do iz1 = 1, nzgrid
cz(iz) = cz(iz) + cmplx(phirt(iz1,ix), b1t(iz1,ix)) * &
abs(iz - iz1) * dz2h
enddo; enddo
endif
phirt(:,ix) = real(cz(1:nzgrid))
b1t(:,ix) = aimag(cz(1:nzgrid))
enddo
call transp_zx(phirt, phi(:,iy,:))
call transp_zx(b1t, b1(:,iy,:))
enddo
!
! Inverse transform in xy
!
if (lshear) then
call fourier_transform_shear_xy(phi,b1,linv=.true.)
else
call fourier_transform_xy(phi,b1,linv=.true.)
endif
!
return
!
endsubroutine inverse_laplacian_isoz
!***********************************************************************
subroutine inverse_laplacian_z_2nd_neumann(f)
!
! 19-mar-2018/MR: coded
! Second-order version in the vertical direction that uses tridag_neumann.
! On input: phi=div(U), on output: phi=potential of the irrotational flow.
!
use Fourier, only: fourier_transform_xy, kx_fft2, ky_fft2
use Mpicomm, only: transp_xz, transp_zx
!!!use Sub, only: tridag_neumann
!
real, dimension(:,:,:,:), intent(INOUT) :: f
real, dimension(nx,ny,0:nz+1) :: phi, b1
real, dimension(nx,ny), save :: k2dz2
real, dimension(nx,ny,1) :: uzhatlow_re, uzhatlow_im
real, dimension (nz-3) :: a, b
real, dimension (nz-2) :: c
!
integer, parameter :: nxt = nx / nprocz
logical, save :: l1st = .true.
real, save :: dz2
integer, save :: ikx0, iky0
!
complex, dimension(nx,ny,nz) :: rhs
!
integer :: ikx, iky, i1
real :: ky2
!
call fatal_error('inverse_laplacian_z_2nd_neumann',"don't call, not yet functioning")
!
! Initialize the array k2dz2 and other constants for future use.
!
if (l1st) then
dz2 = dz**2
ikx0 = ipz * nxt
iky0 = ipy * ny
do iky = 1, ny
ky2 = ky_fft2(iky0+iky)
do ikx = 1, nx
k2dz2(ikx,iky) = (kx_fft2(ikx0+ikx)+ky2)*dz2
enddo
enddo
l1st = .false.
endif
!
! Forward transform of div(U) and U on lower boundary in xy.
!
phi(:,:,1:nz)=f(l1:l2,m1:m2,n1:n2,iphiuu); b1 = 0.
call fourier_transform_xy(phi(:,:,1:nz), b1(:,:,1:nz)) ! Fourier transform of div(U) in cmplx(phi,b1)
print*, 'nach fourier(div)'
uzhatlow_re=f(l1:l2,m1:m2,n1,iuz:iuz); uzhatlow_im=0
call fourier_transform_xy(uzhatlow_re, uzhatlow_im) ! Fourier transform of U_z on lower boundary
if (lfirst_proc_z) then
rhs(:,:,1:1)=2.*dz*cmplx(uzhatlow_re,uzhatlow_im) ! auxiliary to compute phi on boundary
rhs(:,:,2)=3.*dz2*cmplx(phi(:,:,2),b1(:,:,2))+rhs(:,:,1)
i1=3
else
i1=1
endif
rhs(:,:,i1:nz)=dz2*cmplx(phi(:,:,i1:nz),b1(:,:,i1:nz))
!
!do iky = 1, ny
!call transp_xz(phi(:,iky,:), phit)
!call transp_xz(b1(:,iky,:), b1t)
!
! Solution of z dependent problem.
!
a=1.; b=-2.; c=1.; c(1)=2.
!!!call tridag_neumann(a,b,c,k2dz2,rhs,phi,b1)
!call transp_zx(phit, phi(:,iky,:))
!call transp_zx(b1t, b1(:,iky,:))
!
! Inverse transform in xy
!
call fourier_transform_xy(phi(:,:,1:nz),b1(:,:,1:nz),linv=.true.) ! potential in physical space in phi
f(l1:l2,m1:m2,n1:n2,iphiuu)=phi(:,:,1:nz)
!
endsubroutine inverse_laplacian_z_2nd_neumann
!***********************************************************************
subroutine decide_fourier_routine
!
! Decide, based on the dimensionality and on the geometry
! of the grid, which fourier transform routine is to be
! used. "fourier_transform" and "fourier_tranform_shear"
! are functional only without x-parallelization, and for
! cubic, 2D square, and 1D domains only. Everything else
! should use fft_parallel instead.
!
! 05-dec-2011/wlad: coded
!
logical :: lxpresent,lypresent,lzpresent
logical :: l3D,l2Dxy,l2Dyz,l2Dxz
logical :: l1Dx,l1Dy,l1Dz
!
! Check the dimensionality and store it in logicals
!
lxpresent=(nxgrid/=1)
lypresent=(nygrid/=1)
lzpresent=(nzgrid/=1)
!
l3D = lxpresent.and. lypresent.and. lzpresent
l2Dxy= lxpresent.and. lypresent.and..not.lzpresent
l2Dyz=.not.lxpresent.and. lypresent.and. lzpresent
l2Dxz= lxpresent.and..not.lypresent.and. lzpresent
l1Dx = lxpresent.and..not.lypresent.and..not.lzpresent
l1Dy =.not.lxpresent.and. lypresent.and..not.lzpresent
l1Dz =.not.lxpresent.and..not.lypresent.and. lzpresent
!
if (ldebug.and.lroot) then
if (l3D) print*,"This is a 3D simulation"
if (l2Dxy) print*,"This is a 2D xy simulation"
if (l2Dyz) print*,"This is a 2D yz simulation"
if (l2Dxz) print*,"This is a 2D xz simulation"
if (l1Dx) print*,"This is a 1D x simulation"
if (l1Dy) print*,"This is a 1D y simulation"
if (l1Dz) print*,"This is a 1D z simulation"
endif
!
! The subroutine "fourier_transform" should only be used
! for 1D, square 2D or cubic domains without x-parallelization.
! Everything else uses fft_parallel.
!
luse_fourier_transform=(nprocx==1.and.&
(l1dx .or.&
l1dy .or.&
l1dz .or.&
(l2dxy.and.nxgrid==nygrid) .or.&
(l2dxz.and.nxgrid==nzgrid) .or.&
(l2dyz.and.nygrid==nzgrid) .or.&
(l3d.and.nxgrid==nygrid.and.nygrid==nzgrid)))
!
endsubroutine decide_fourier_routine
!***********************************************************************
subroutine read_poisson_init_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=poisson_init_pars, IOSTAT=iostat)
!
endsubroutine read_poisson_init_pars
!***********************************************************************
subroutine write_poisson_init_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=poisson_init_pars)
!
endsubroutine write_poisson_init_pars
!***********************************************************************
subroutine read_poisson_run_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=poisson_run_pars, IOSTAT=iostat)
!
endsubroutine read_poisson_run_pars
!***********************************************************************
subroutine write_poisson_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=poisson_run_pars)
!
endsubroutine write_poisson_run_pars
!***********************************************************************
subroutine get_acceleration(acceleration)
!
use General, only: keep_compiler_quiet
!
real, dimension(nx,ny,nz,3), intent(out) :: acceleration !should I (CAN I?) make this allocatable?
!
call keep_compiler_quiet(acceleration)
!
endsubroutine get_acceleration
!***********************************************************************
endmodule Poisson