SUBROUTINE poisr(rhs,dx,dxh,dy,dz,dzh, &
ibc1,ibc2,kbc1,kbc2,ksen)
!
! CHANGES:
! includes jtot and ksen (calculated in poisson.f) in header
! help array rhst with dimensions (imax,jtot,ksen)
! ALL_ALL copy rhs to rhst and back for FFT's
!
! ibc?=1: neumann
! ibc?=2: periodic
! ibc?=3: dirichlet
!
! only FFT in j-direction, cyclic reduction in the others
! include'param.txt'
! include 'mpif.h'
! include 'mpi_cons.txt'
use modglobal, only : imax,isen,jmax,jsen,jtot,kmax,poisrcheck
use modfields, only : worksave
use modmpi, only : myid,comm3d,mpierr,my_real,nprocs, barrou,MPI_SUM
implicit none
! authors: m.j.b.m. pourquie, b.j. boersma
integer i, j, k, ksen
real rhs(0:imax+1,0:jmax+1,0:kmax+1),dx(0:IMAX+1),dxh(1:IMAX+1),dy
real, allocatable, dimension(:,:,:) :: rhs2
real, allocatable, dimension(:) :: work
integer ier, iperio, kperio
integer ibc1,ibc2,kbc1,kbc2
real dz(0:kmax+1),dzh(1:kmax+1),pi ! dz: kb-1:ke+1, dzh: kb:ke+1
real a(imax),b(imax),c(imax),bin(imax)
real az(kmax),bz(kmax),cz(kmax)
real yrt(jtot)
real, allocatable, dimension(:,:) :: vfftj
real, allocatable, dimension(:,:) :: y
real wj(jtot+15)
real angle, tst
integer ipos, jv
real suml, sum
!
! test write
! write(6,*)'POISR, kmax, ksen, nprocs,jmax,jtot ',kmax, ksen, nprocs ,jmax,jtot
allocate(rhs2(imax,jtot,ksen))
allocate(work(2*imax*jmax*kmax))
allocate(vfftj(imax*ksen,jtot))
allocate(y(imax,kmax))
pi=4.*atan(1.)
do i=1,imax
a(i) = 1./(dx(i)*dxh(i))
c(i) = 1./(dx(i)*dxh(i+1))
b(i) = - (a(i) + c(i))
enddo
if((ibc1).eq.1)then
! Neumann
b(1) = b(1) + a(1)
elseif(ibc1.eq.2)then
! periodic
b(1) = b(1)
elseif(ibc1.eq.3)then
! Dir
b(1) = b(1) - a(1)
endif
if((ibc2).eq.1)then
! Neumann
b(imax) = b(imax) + c(imax)
elseif((ibc2).eq.2)then
b(imax) = b(imax)
elseif((ibc2).eq.3)then
! b(imax) = b(imax) - c(kmax) ! Jasper T. : bug? Should be b(imax) - c(imax)?
b(imax) = b(imax) - c(imax)
endif
if(ibc1.ne.2)then
c(imax) = 0.
a(1) = 0.
endif
! do i=1,imax
! write(6,*)'i, abc(i)',i, a(i),b(i),c(i)
! enddo
!
! fill coefficients in k-direction
! do k=1,kmax ! Mathieu's version
! az(k) = 1./(dz(k)*dzh(k-1))
! cz(k) = 1./(dz(k)*dzh(k))
! bz(k) = - (az(k) + cz(k))
! enddo
do k=1,kmax
az(k) = 1./(dz(k)*dzh(k))
cz(k) = 1./(dz(k)*dzh(k+1))
bz(k) = - (az(k) + cz(k))
enddo
!
! BC:
!
!periodic bz(1) = bz(1) - az(1)
!periodic az(1) = 0.
!periodic bz(kmax) = bz(kmax) - az(kmax)
!periodic az(kmax) = 0.
if((kbc1).eq.1)then
! Neumann
bz(1) = bz(1) + az(1)
elseif(kbc1.eq.2)then
! periodic
bz(1) = bz(1)
endif
if((kbc2).eq.1)then
! Neumann
bz(kmax) = bz(kmax) + cz(kmax)
elseif((kbc2).eq.2)then
bz(kmax) = bz(kmax)
elseif((kbc2).eq.3)then
!
! p = 0.
bz(kmax) = bz(kmax) - cz(kmax)
endif
if(kbc1.ne.2)then
cz(kmax) = 0.
az(1) = 0.
endif
! do k=1,kmax
! write(6,*)'k, abc(k)',k, az(k),bz(k),cz(k)
! enddo
!
! initialise for FFT
! PAR
! call vrffti(jmax,wj)
call vrffti(jtot,wj)
!
yrt(1)=0.
! PAR
! yrt(jmax)=-4./(dy*dy)
yrt(jtot)=-4./(dy*dy)
do j=3,jtot,2
!
! 2.*(cos2(alpha) -1) = 2.*(-2.*sin(alpha)**2
yrt(j-1)=(-4./(dy*dy))*(sin(float((j-1))*pi/(2.*jtot)))**2
yrt(j)= yrt(j-1)
enddo
! help = rhs
!
! sum check (comment out if not deeded)
!
!! suml = 0.
!! do k=1,kmax
!! do j=1,jmax
!! do i=1,imax
!! suml = suml + rhs(i,j,k)*dx(i)*dy*dz(k)
!! enddo
!! enddo
!! enddo
!! call MPI_ALLREDUCE(suml, sum, 1, MY_REAL, MPI_SUM, comm3d,mpierr)
!! if(myid.eq.0) write(6,*)'solver sum = ', sum
!
! end sum check
!
call barrou()
call ALL_ALL_j2(rhs,rhs2,0,ksen)
call barrou()
do k=1,ksen
do i=1,imax
ipos=(k-1)*imax+i
do j=1,jtot
vfftj(ipos,j)=rhs2(i,j,k)
! suml = suml + rhs2(i,j,k)*dy*dz(k)*dx(i)
enddo
enddo
enddo
! call MPI_ALLREDUCE(suml, sum, 1, MY_REAL, MPI_SUM, comm3d,mpierr)
! if(myid.eq.0) write(6,*)'solver sum = ', sum
!
call vrfftf(imax*ksen,jtot,vfftj,rhs2,imax*ksen,wj)
!
! goto 1234
do k=1,ksen
do i=1,imax
ipos=(k-1)*imax+i
do j=1,jtot
rhs2(i,j,k) = vfftj(ipos,j)
enddo
enddo
enddo
call barrou()
call ALL_ALL_j2(rhs,rhs2,1,ksen)
call barrou()
! call sumchk3(rhs,imax,jmax,kmax,8,1)
! PAR CHECK!!!
do j=1,jmax ! begin loop over angles
! do j=1,jtot ! begin loop over angles
!
! add proper part to diagonal
jv = j + myid*jmax
do i=1,imax
! PAR CHECK!!
bin(i)=b(i) + yrt(jv)!!!!!!!!!/(Rp(i)**2)
enddo
!ATTT do k=1,ksen
do k=1,kmax
do i=1,imax
ipos=(k-1)*imax+i
!PAR CHECK
y(i,k) = rhs(i,j,k)!vfftj(ipos,j)
enddo
enddo
iperio=1
kperio=1
if(ibc1.eq.2)iperio=0 ! i periodic
if(kbc1.eq.2)kperio=0 ! k periodic
if(poisrcheck.eq.0) then
poisrcheck=1
CALL BLKTRI(0,kperio,kmax,az,bz,cz,iperio,imax,a,bin,c,imax,y &
!
! ^ 0 for periodic BC
,ier,work)
! write(6,*)'ier = ', ier, iperio, kperio
! if(ier.ne.0)stop 'IER'
worksave = work
write(6,*) 'First time step in POISR, poisrcheck=', poisrcheck
end if
CALL BLKTRI(1,kperio,kmax,az,bz,cz,iperio,imax,a,bin,c,imax,y &
,ier,worksave)
!
! write(6,*)'myid, ier = ', myid, ier, j
! if(ier.ne.0)stop 'IER'
do k=1,kmax
do i=1,imax
ipos=(k-1)*imax+i
!PAR CHECK
! vfftj(ipos,j) = y(i,k)
rhs(i,j,k) = y(i,k)
enddo
enddo
enddo ! end loop over angles
1234 continue
! call sumchk3(rhs,imax,jmax,kmax,9,1)
!
call barrou()
call ALL_ALL_j2(rhs,rhs2,0,ksen)
call barrou()
!
do k=1,ksen
do i=1,imax
ipos=(k-1)*imax+i
do j=1,jtot
vfftj(ipos,j)=rhs2(i,j,k)
enddo
enddo
enddo
call vrfftb(imax*ksen,jtot,vfftj,rhs2,imax*ksen,wj)
!!ATTTT dok=1,kmax
do k=1,ksen
do i=1,imax
ipos=(k-1)*imax+i
do j=1,jtot
rhs2(i,j,k)=vfftj(ipos,j)
! if(abs(rhs(i,j,k)-help(i,j,k)).gt.1.e-13)then
! write(6,*)'ERRR', i, j, k, rhs(i,j,k), help(i,j,k)
! endif
enddo
enddo
enddo
call barrou()
call ALL_ALL_j2(rhs,rhs2,1,ksen)
call barrou()
! call sumchk3(rhs,imax,jmax,kmax,2,1)
deallocate(rhs2)
deallocate(work)
deallocate(vfftj)
deallocate(y)
return
end subroutine poisr