subroutine bcpup(pup, pvp, pwp, rk3coef)
use modglobal, only : ib, ie, jb, je, ih, jh, kb, ke, kh, rk3step, dxi, dyi, dzhi, &
ibrank, ierank, jbrank, jerank, BCxm, BCym, BCtopm, &
BCtopm_freeslip, BCtopm_noslip, BCtopm_pressure, &
BCxm_periodic, BCxm_profile, BCxm_driver, &
BCym_periodic, BCym_profile
use modfields, only : pres0, up, vp, wp, um, vm, wm, w0, u0, v0, uouttot, vouttot, uinit, vinit, uprof, vprof, pres0, IIc, IIcs
use modmpi, only : excjs, excis, myid, avexy_ibm
use modinletdata, only : u0driver
use decomp_2d, only : exchange_halo_z
real, dimension(ib - ih:ie + ih, jb - jh:je + jh, kb:ke + kh), intent(inout) :: pup
real, dimension(ib - ih:ie + ih, jb - jh:je + jh, kb:ke + kh), intent(inout) :: pvp
real, dimension(ib - ih:ie + ih, jb - jh:je + jh, kb:ke + kh), intent(inout) :: pwp
real, dimension(kb:ke+kh) :: pres0ij
real, intent(in) :: rk3coef
real rk3coefi
integer i, j, k
rk3coefi = 1./rk3coef
! if (jbrank) write(*,*) "jb before exhange_halo ", pvp(ie/2,jb,ke)
! if (jerank) write(*,*) "je before exhange_halo ", pvp(ie/2,je+1,ke)
! Watch this communication as it is slightly different to normal -
! maybe safer to just resize to kb-kh:ke+kh
call exchange_halo_z(pup, opt_zlevel=(/ih,jh,0/))
call exchange_halo_z(pvp, opt_zlevel=(/ih,jh,0/))
call exchange_halo_z(pwp, opt_zlevel=(/ih,jh,0/))
! if (jbrank) write(*,*) "jb after exhange_halo ", pvp(ie/2,jb,ke)
! if (jerank) write(*,*) "je after exhange_halo ", pvp(ie/2,je+1,ke)
select case(BCtopm)
case(BCtopm_freeslip, BCtopm_noslip)
do j = jb, je
do i = ib, ie
pwp(i, j, kb) = 0.
pwp(i, j, ke + kh) = 0.
end do
end do
case(BCtopm_pressure)
call avexy_ibm(pres0ij(kb:ke+kh),pres0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+kh),IIcs(kb:ke+kh),.false.)
do j = jb, je
do i = ib, ie
pwp(i, j, kb) = 0.
!pwp(i, j, ke + 1) = wm(i, j, ke+1) * rk3coefi - (-pres0ij(ke) - pres0(i,j,ke)) * dzhi(ke+1) ! Doesn't work
pwp(i, j, ke + 1) = wm(i, j, ke+1) * rk3coefi + 2 * pres0ij(ke)*dzhi(ke+1)
wp(i, j, ke + 1) = pwp(i, j, ke+1) - wm(i,j,ke+1) * rk3coefi
end do
end do
end select !BCtopm
select case(BCxm)
case(BCxm_periodic)
if (ibrank .and. ierank) then ! not parallelised in x
do k = kb, ke
do j = jb, je
pup(ie+1, j, k) = pup(ib, j, k) ! cyclic
end do
end do
end if
case(BCxm_profile)
if (ibrank) then
do k=kb,ke
do j=jb-1,je+1
pup(ib, j, k) = uprof(k) * rk3coefi
up(ib, j, k) = 0.
end do
end do
end if
if (ierank) then
do k = kb+1, ke
do j = jb-1, je+1
! convective
pup(ie+1, j, k) = um(ie+1, j, k) * rk3coefi - (u0(ie+1, j, k) - u0(ie,j,k)) * dxi * uouttot !u0(ie,j,k) ! du/dt +u*du/dx=0 -> pup(i)=um(i)/rk3coef -um(i)*(um(i)-um(i-1))/dxf(i-1)
! Neumann
!pup(ie+1,j,k) = pup(ie,j,k)
up(ie+1, j, k) = pup(ie+1, j, k) - um(ie+1,j,k)*rk3coefi
end do
end do
! Neumann at bottom - performs better with no slip
pup(ie+1, :, kb) = pup(ie, :, kb)
up(ie+1, :, kb) = pup(ie+1,: , kb) - um(ie+1, :, kb) * rk3coefi
end if
case(BCxm_driver)
if (ibrank) then
do k = kb, ke
do j = jb - 1, je + 1
pup(ib, j, k) = u0driver(j, k) * rk3coefi
up(ib, j, k) = 0. ! u(ib) only evolves according to pressure correction
end do
end do
end if
if (ierank) then
do k = kb, ke
do j = jb-1, je+1
pup(ie+1, j, k) = um(ie+1, j, k) * rk3coefi - (u0(ie+1, j, k) - u0(ie, j, k)) * dxi * uouttot ! du/dt +u*du/dx=0 -> pup(i)=um(i)/rk3coef -um(i)*(um(i)-um(i-1))/dxf(i-1)
! !Neumann
!pup(ie+1,j,k) = pup(ie,j,k)
up(ie+1, j, k) = pup(ie+1, j, k) - um(ie+1, j, k) * rk3coefi
end do
end do
! Neumann at bottom - performs better with no slip
! pup(ie+1, :, kb) = pup(ie, :, kb)
! up(ie+1, :, kb) = pup(ie+1, :, kb) - um(ie+1, :, kb) * rk3coefi
end if
end select ! BCxm
select case(BCym)
case(BCym_periodic)
if (jbrank .and. jerank) then ! not parallelised in y
do k = kb, ke
do i = ib, ie
pvp(i, je+1, k) = pvp(i, jb, k) ! cyclic
end do
end do
end if
case(BCym_profile)
if (jbrank) then
do k = kb, ke
do i = ib-1, ie+1
pvp(i,jb,k) = vprof(k)*rk3coefi
vp(i,jb,k) = 0.
end do
end do
end if
if (jerank) then
do k = kb, ke
do i = ib-1, ie+1
! change to vouttot
pvp(i, je+1, k) = vm(i, je+1, k) * rk3coefi - (v0(i, je+1, k) - v0(i, je, k)) * dyi * vouttot
vp(i, je+1, k) = pvp(i, je+1, k) - vm(i,je+1,k)*rk3coefi
end do
end do
pvp(:, je+1, kb) = pvp(:, je, kb)
vp(:, je+1, kb) = pvp(:, je+1, kb) - vm(:, je+1, kb)*rk3coefi
end if
end select
end subroutine bcpup