subroutine press_stag_bldg(p_hat,dfdx,dfdy,divtz)   
! p_hat contains the physical space pressure on exit
!-------------------    
! Boundary Layer version with 4th order derivs in vertical.
!  04 December 1995
!	Mods.
!	12/6: added upper and lower boundary conditions.	
!    12/8: Corrected sign on imag. x and y deriv of RHS.
!	12/17 Forcing pressure equal zero at wall 
!			Removed forcing of <P>=0 for all z.
!		    Will need to change BC when hetero. surface stress.
!	12/17 Added ficticious node below wall (Note solver on Nz+1 x Nz+1
!          prev. version of this subroutine saved as p_bl4.old.12.17
!    12/18 Revised 2st deriv stencil at wall (avg of deriv at -1 and at 1)
!    12/21 Redid FDD to 2nd order accurate.
!    12/22 Now value of P(wall) diagnosed from prev P, to match gradient BC
!....1/13: major changes.
!		Broke out mean pressure for separate solution
!    1/20 back to tridag for natrix solution (same sol'n as LUDCMP...A Keeper!)
!....1/23 Staggered solution
!.........Using Nz+1 levels for computing P, but tossing out level below ground
!....4/1 Changed sign on Div T_iz at wall and lid (five places)
!-------------------          
use types,only:rprec
use param,only:ld,lh,nx,ny,nz,dt,dz
use sim_param,only:u,v,w,RHSx,RHSy,RHSz,RHSx_f,RHSy_f,RHSz_f
use fft
use immersedbc,only:fx,fy,fz
implicit none   
real(kind=rprec),dimension(ld,ny,nz),intent(in)::divtz
complex(kind=rprec),dimension(lh,ny,0:nz)::p_hat
real(kind=rprec),dimension(ld,ny,nz)::rH_x,rH_y,rH_z
complex(kind=rprec),dimension(lh,ny,nz)::H_x,H_y,H_z
equivalence (rH_x,H_x),(rH_y,H_y),(rH_z,H_z)
real(kind=rprec),dimension(ld,ny)::rtopw, rbottomw
complex(kind=rprec),dimension(lh,ny)::topw,bottomw
equivalence (rtopw,topw),(rbottomw,bottomw)
complex(kind=rprec),dimension(lh,ny,nz),intent(out)::dfdx,dfdy
complex(kind=rprec),dimension(nz+1)::RHS_col,p_col
real(kind=rprec),dimension(nz+1)::a1,b1,c1
real(kind=rprec)::const,ignore_me
! remove this stuff!
integer::jx,jy,jz,size,k

size=nz+1
p_hat(:,:,0)=0._rprec ! dont know for sure if this is right
!==========================================================================
! Get the right hand side ready 
! Loop over levels    
const=1._rprec/(nx*ny)
do jz=1,nz
! temp storage for sum of RHS terms.  normalized for fft
! sc: recall that the old timestep guys already contain the pressure
!   term
   rH_x(:,:,jz)=const*(RHSx(:,:,jz)-1._rprec/3._rprec*RHSx_f(:,:,jz)&
        +2._rprec/(3._rprec*dt)*u(:,:,jz)+2._rprec/3._rprec*fx(:,:,jz))
   rH_y(:,:,jz)=const*(RHSy(:,:,jz)-1._rprec/3._rprec*RHSy_f(:,:,jz)&
        +2._rprec/(3._rprec*dt)*v(:,:,jz)+2._rprec/3._rprec*fy(:,:,jz))
   rH_z(:,:,jz)=const*(RHSz(:,:,jz)-1._rprec/3._rprec*RHSz_f(:,:,jz)&
        +2._rprec/(3._rprec*dt)*w(:,:,jz)+2._rprec/3._rprec*fz(:,:,jz))
   call rfftwnd_f77_one_real_to_complex(forw,rH_x(:,:,jz),ignore_me)
   call rfftwnd_f77_one_real_to_complex(forw,rH_y(:,:,jz),ignore_me)
   call rfftwnd_f77_one_real_to_complex(forw,rH_z(:,:,jz),ignore_me)     
end do

! sc: could do out of place transform if wanted to...
rbottomw(:,:)=const*divtz(:,:,1)
rtopw(:,:)=const*divtz(:,:,nz)
call rfftwnd_f77_one_real_to_complex(forw,rbottomw(:,:),ignore_me)
call rfftwnd_f77_one_real_to_complex(forw,rtopw(:,:),ignore_me)

! set oddballs to 0
! probably can get rid of this if we're more careful below
H_x(lh,:,:)=0._rprec
H_y(lh,:,:)=0._rprec
H_z(lh,:,:)=0._rprec
H_x(:,ny/2+1,:)=0._rprec
H_y(:,ny/2+1,:)=0._rprec
H_z(:,ny/2+1,:)=0._rprec
topw(lh,:)=0._rprec
topw(:,ny/2+1)=0._rprec
bottomw(lh,:)=0._rprec
bottomw(:,ny/2+1)=0._rprec
!==========================================================================
! Loop over (Kx,Ky) to solve for Pressure amplitudes
do jy=1,ny
   if (jy.eq.ny/2+1) then
      p_hat(:,jy,:)=0._rprec
      goto 5111
   end if
! since this is all in Fourier space, don't we need to skip ny/2+1?
   do jx=1,lh-1
! the lh part should be zero, right?
! this is the k=0 part
! Zero wavenumber problem is non-unique for matrix solution
! So just integrate up from wall with arbitrary P=const starting point.
! why is the jz index the innermost loop?
      if (jx*jy.eq.1) then          
! what is in p_hat(:,:,0) here?
! does the (:,:,1) entry need to be real as in original?
! the 0. here is assuming that p(z=0)=0.
         p_hat(jx,jy,1)=0._rprec-dz*bottomw(jx,jy)
         do jz=2,nz
            p_hat(jx,jy,jz)=p_hat(jx,jy,jz-1)+H_z(jx,jy,jz)*dz
         end do
      else
! See end of wavenumber loops for endif
! Set up matrix solution for pressure Fourier coeff's
! Near Wall Nodes
         a1(1)=0._rprec
         b1(1)=-1._rprec
         c1(1)=1._rprec
         RHS_col(1)=-dz*bottomw(jx,jy) ! complex 
! Interior Nodes
         do jz=2,nz
            RHS_col(jz)=eye*(kx(jx,jy)*H_x(jx,jy,jz-1)+&
                 ky(jx,jy)*H_y(jx,jy,jz-1))&
                 +(H_z(jx,jy,jz)-H_z(jx,jy,jz-1))/dz
! this part OK
            a1(jz)=1._rprec/(dz**2)
            b1(jz)=-(kx(jx,jy)**2+ky(jx,jy)**2+2._rprec/(dz**2))
            c1(jz)=1._rprec/(dz**2)
         end do

! top nodes 
         RHS_col(nz+1)=-topw(jx,jy)*dz
! note: check to make sure the topw's aren't automatically 0
! this would save some time
         a1(nz+1)= -1._rprec
         b1(nz+1)=1._rprec
         c1(nz+1)=0._rprec
! Call Solver
! modified to handle complex
         call tridag (a1,b1,c1,RHS_col,p_col,size)
! Put Pressure amplitudes in Matrix
         do jz=0,nz
! if called above directly with a row of p_hat, we'd probably encounter
! some serious cache thrashing??? try it
            p_hat(jx,jy,jz)=p_col(jz+1)
         end do
! Ending IF from zero wavenumber solution
      end if
! End loop over kx,ky
   end do
   p_hat(lh,jy,:)=0._rprec
5111 continue
end do

!      do jz=0,nz
!        print*, p_hat(1,1,jz),p_hat(5,5,jz)
!      end do


!=========================================================================== 
!...Now need to get p_hat(wave,level) to physical p(jx,jy,jz)   
!.....Loop over height levels     
call rfftwnd_f77_one_complex_to_real(back,p_hat(:,:,0),ignore_me)
do jz=1,nz
do jy=1,ny
do jx=1,lh
! complex
   dfdx(jx,jy,jz)=eye*kx(jx,jy)*p_hat(jx,jy,jz)
   dfdy(jx,jy,jz)=eye*ky(jx,jy)*p_hat(jx,jy,jz)
! note the oddballs of p_hat are already 0, so we should be OK here
end do
end do
call rfftwnd_f77_one_complex_to_real(back,dfdx(:,:,jz),ignore_me)
call rfftwnd_f77_one_complex_to_real(back,dfdy(:,:,jz),ignore_me)
call rfftwnd_f77_one_complex_to_real(back,p_hat(:,:,jz),ignore_me)
end do
end subroutine press_stag_bldg
