      SUBROUTINE psolver(ini,rhs,Ru,Rp,dphi,zp,zw)
      implicit none
c
c     RHS(IMAX,JMAX,KMAX)  on entry is the right hand side of the Poisson
c                          equation. On exit it contains the solution.
c
c     Ru(0:i1), Rp(0:i1)   arrays with the radial locations of velocity and
c                          pressure points
c
c     dphi                 azimuthal grid spacing
c
c     dz                   axial grid spacing
c
c
c

      include   'param.txt'
      REAL      RHS(IMAX,JMAX,KMAX),Ru(0:IMAX+1),Rp(0:IMAX+1),DPHI
      REAL      Y  (IMAX,KMAX)
      REAL      work(3*imax*kmax),pi,SW(10*imax*kmax),WW(10*imax*kmax)
c
c     Locals
c
      real      dz,dzi,bbb,z
      real      ar(imax),br(imax),cr(imax),dm(imax),brr(imax)
      real      az(kmax),bz(kmax),cz(kmax)
      real      yrt(jmax),wj(jmax+15)
      real      vfftj(imax*kmax,jmax)
      real      zp(0:kmax+1),zw(0:kmax+1)

c     WJ and WK are arrays with tables for Fourier transforms.
c

      integer   ipos,ier,ini
c
      call bvrffti(jmax,wj)

      pi = 4.*atan(1.)

c  J --> direction      (yrt)
      yrt(1    ) = 0.
      yrt(jmax ) = -4./(dphi*dphi)
      do j=3,jmax,2
      yrt(j-1)=(-4./(dphi*dphi))*(sin(float((j-1))*pi/(2.*jmax)))**2
      yrt(j  )= yrt(j-1)
      enddo

     
      do k=1,kmax
         do i=1,imax
         ipos = (k-1)*imax+i
           do j=1,jmax
             vfftj(ipos,j)=rhs(i,j,k)
           enddo
         enddo
      enddo

c     FFT in j direction

      call bvrfftf(imax*kmax,jmax,vfftj,rhs,imax*kmax,wj)
      do k=1,kmax
         do i=1,imax
         ipos = (k-1)*imax+i
            do j=1,jmax
               rhs(i,j,k) =vfftj(ipos,j)
            enddo
         enddo
      enddo
c
c
c
c  

      do j=1,jmax



	do k=1,kmax
	  do i=1,imax
	   Y(i,k) = rhs(i,j,k)
          enddo
        enddo

       do i=1,imax
       ar(i)=  Ru(I-1)/((Rp(I)-Rp(I-1))*Rp(I)*(Ru(I)-Ru(I-1)))
       br(i)=-(Ru(I)/(Rp(I+1)-Rp(I))+Ru(I-1)/(Rp(I)-Rp(I-1)))/
     $      (Rp(I)*(Ru(I)-Ru(I-1)))
       cr(i)=  Ru(I  ) /((Rp(I+1)-Rp(I))*Rp(I)*(Ru(I)-Ru(I-1)))
       dm(i)=1./((Rp(i)*dphi)**2)
       end do

       br(1)=   -(Ru(1)/(Rp(2)-Rp(1)))/(Rp(1)*(Ru(1)-Ru(0)))
       ar(1)=0.
       br(imax)=br(imax)-cr(imax)  ! dirichlet
c      br(imax)=br(imax)+cr(imax)   ! neumann
       cr(imax)=0.

       do k=1,kmax
           az(k)=1./( ( Zw(k)-Zw(k-1))*(Zp(k)-Zp(k-1)) )
           cz(k)=1./( ( Zw(k)-Zw(k-1))*(Zp(k+1)-Zp(k)) )
           bz(k)=-(az(k)+cz(k))
       enddo
       bz(1)=bz(1)+az(1)
       az(1)=0.
       bz(kmax)=bz(kmax)+cz(kmax)    ! orginal - neumann
c      bz(kmax)=bz(kmax)-cz(kmax)   ! dirichlet
       cz(kmax)=0.

      do i=1,imax
	brr(i)=br(i)+yrt(j)/(rp(i)**2)
      enddo 

      call blktri(0,1,kmax,az,bz,cz,1,
     $               imax,ar,brr,cr,imax,y,ier,work)
      call blktri(1,1,kmax,az,bz,cz,1,imax,ar,brr,cr,imax,y,
     $               ier,work)


c... these are crayfish routines which are on vn...

c      call S2GCIS(imax,ar,brr,cr,kmax,az,bz,cz,imax,ier,SW,WW)
c      call S2GCSS(imax,y,SW,WW) 

      do k=1,kmax
          do i=1,imax
          rhs(i,j,k)=Y(i,k)
          enddo
        enddo



      enddo  ! end j-loop


      do k=1,kmax
         do i=1,imax
         ipos = (k-1)*imax+i
           do j=1,jmax
             vfftj(ipos,j)=rhs(i,j,k)
           enddo
         enddo
      enddo
      call bvrfftb(imax*kmax,jmax,vfftj,rhs,imax*kmax,wj)
      do k=1,kmax
         do i=1,imax
         ipos = (k-1)*imax+i
            do j=1,jmax
               rhs(i,j,k) =vfftj(ipos,j)
            enddo
         enddo
      enddo
      end

