subroutine tderive
!-----------------------------------------------------------------|
! |
!*** *tderive* read input fields for initialisation |
! |
! Hans Cuijpers I.M.A.U. 06/01/1995 |
! |
! purpose. |
! -------- |
! |
! Refill array p with pressure values. The poisson-solver |
! produced a pressure array p in which the i-index varied |
! between 2 and i1, the j- and k-index between 1 and resp. |
! jmax and kmax. For our further calculations we'll change |
! the range for the j-index to vary between 2 and j1. |
! |
! Further we set cyclic boundary conditions for the pressure- |
! fluctuations in the x-y plane. |
! |
!** interface. |
! ---------- |
! |
! *tderive* is called from *program*. |
! |
!-----------------------------------------------------------------|
use modfields, only : up, vp, wp, pres0, IIc, IIcs
use modglobal, only : ib,ie,ih,jb,je,jh,kb,ke,kh,dxhi,dyi,dzhi,linoutflow,rslabs
use modmpi, only : myid,excj,slabsum,avexy_ibm
use modboundary,only : bcp
implicit none
integer i,j,k
real, dimension(kb-kh:ke+kh) :: pij
! logical, dimension(ib:ie, jb:je, kb:ke) :: pnan
real :: pijk
! integer :: ipnan
! Mathieu ATTTT: CHANGED!!! Loop removed!!!
! ** Boundary conditions **************
call bcp(p) ! boundary conditions for p.
!*****************************************************************
! ** Calculate time-derivative for the velocities with known ****
! ** pressure gradients. ***************************************
!*****************************************************************
! if (myid == 0) then
! write(*,*) "net ts", p(ib, jb, :)
! end if
!pnan = isnan(p(ib:ie, jb:je, kb:ke))
!ipnan = 0
!do k=kb,ke
!do j=jb,je
!do i=ib,ie
! if (pnan(i,j,k)) then
! ipnan = ipnan + 1
! end if
!end do
!end do
!end do
! write(*,*) "NaN", myid, ipnan
! write(*,*) "NaN", myid, dble(isnan(p)) !sum(real(isnan(p)))
do k=kb,ke
do j=jb,je
do i=ib,ie
vp(i,j,k) = vp(i,j,k)-(p(i,j,k)-p(i,j-1,k))*dyi
end do
end do
end do
if (linoutflow ) then
do k=kb,ke
do j=jb,je
do i=ib,ie+1
up(i,j,k) = up(i,j,k)-(p(i,j,k)-p(i-1,j,k))*dxhi(i) ! see equation 5.82 (u is computed from the mass conservation)
end do
end do
end do
else
do k=kb,ke
do j=jb,je
do i=ib,ie
up(i,j,k) = up(i,j,k)-(p(i,j,k)-p(i-1,j,k))*dxhi(i) ! see equation 5.82 (u is computed from the mass conservation)
end do
end do
end do
endif
do k=kb+1,ke
do j=jb,je
do i=ib,ie
wp(i,j,k) = wp(i,j,k)-(p(i,j,k)-p(i,j,k-1))*dzhi(k)
end do
end do
end do
! tg3315 02/02/2019
! account for pressure offset that results from ill-defined problem in pressure
! correction method when periodic horizontal BCs are applied. Arises within cyclic
! reduction scheme (called by BLKTRI) and due to the numerics in PRODP in
! cycred.f which define the BC in the periodic case. A linear offset existed in
! the pressure correction term (p) and this can be controlled by subtracting
! the volume averaged modified pressure from this value at all time steps.
! Periodic: p - <p>_ijk
! Makes no change on physical effect of modified pressure in code.
! tg3315 - update 24/06/19 -- there is a missing term in the application of the periodic BCs for pup, could this be part of the problem? Test with this to see if can avoid use of pijk below.
! refer to mvr for necessity of this
! useful refs:
! https://opensky.ucar.edu/islandora/object/technotes%3A98/datastream/PDF/download/citation.pdf
! https://epubs.siam.org/doi/pdf/10.1137/0711042
pij =0.; pijk=0.;
if (.not. linoutflow) then
call slabsum(pij(kb:ke),kb,ke,p(ib:ie,jb:je,kb:ke),ib,ie,jb,je,kb,ke,ib,ie,jb,je,kb,ke)
pij = pij/rslabs
pijk = sum(pij(kb:ke))/(ke-kb)
end if
do k=kb-1,ke+1
do j=jb-1,je+1
do i=ib-1,ie+1
pres0(i,j,k)=pres0(i,j,k)+p(i,j,k)-pijk ! update of the pressure: P_new = P_old + p
enddo
enddo
enddo
return
end subroutine tderive