subroutine masscorr
!> correct the velocities to get prescribed flow rate
use modglobal, only : ib,ie,jb,je,ih,jh,kb,ke,dzf,dxf,dy,dt,rk3step,&
uflowrate,vflowrate,linoutflow,&
luoutflowr,lvoutflowr,luvolflowr,lvvolflowr
use modfields, only : um,up,vm,vp,uout,uouttot,udef,vout,vouttot,vdef,&
uoutarea,voutarea,fluidvol,IIu,IIv
use modmpi, only : myid,comm3d,mpierr,nprocs,MY_REAL,sumy_ibm
real, dimension(ib:ie, kb:ke) :: uvol
real, dimension(ib:ie, kb:ke) :: uvolold
real, dimension(ib:ie, kb:ke) :: vvol
real, dimension(ib:ie, kb:ke) :: vvolold
real, dimension(kb:ke) :: uoutold
real, dimension(kb:ke) :: voutold
real rk3coef,rk3coefi,&
uoutflow,voutflow,&
uflowrateold,vflowrateold
integer i,j,k
if ((.not.linoutflow) .and. (luoutflowr)) then
rk3coef = dt / (4. - dble(rk3step))
rk3coefi = 1 / rk3coef
udef = 0.
uout = 0.
uoutflow = 0.
uoutold = 0.
! integrate u fixed at outlet ie along y
call sumy_ibm(uout,up(ie,jb:je,kb:ke)*dy,ie,ie,jb,je,kb,ke,IIu(ie,jb:je,kb:ke)) ! u tendency at previous time step
call sumy_ibm(uoutold,um(ie,jb:je,kb:ke)*dy,ie,ie,jb,je,kb,ke,IIu(ie,jb:je,kb:ke)) ! u at previous time step
! integrate u in z
do k=kb,ke
uout(k) = rk3coef*uout(k)*dzf(k)
uoutold(k) = uoutold(k)*dzf(k)
end do
uoutflow = sum(uout(kb:ke))
uflowrateold = sum(uoutold(kb:ke))
! average over outflow area
uoutflow = uoutflow/uoutarea
uflowrateold = uflowrateold/uoutarea
! flow correction to match outflow rate
udef = uflowrate - (uoutflow + uflowrateold)
do k = kb,ke
do j = jb,je
do i = ib,ie
up(i,j,k) = up(i,j,k) + udef*rk3coefi
end do
end do
end do
! bss116 calculate uouttot which is used in modboundary.
! this really should be in the routine directly!
uouttot = sum(uout(kb:ke)) ! mass flow rate at outlet
elseif ((.not.linoutflow) .and. (luvolflowr)) then
rk3coef = dt / (4. - dble(rk3step))
rk3coefi = 1 / rk3coef
udef = 0.
uout = 0.
uoutflow = 0.
uoutold = 0.
uvol = 0.
uvolold = 0.
! integrate u in y
call sumy_ibm(uvol,up(ib:ie,jb:je,kb:ke)*dy,ib,ie,jb,je,kb,ke,IIu(ib:ie,jb:je,kb:ke)) ! u tendency at previous time step
call sumy_ibm(uvolold,um(ib:ie,jb:je,kb:ke)*dy,ib,ie,jb,je,kb,ke,IIu(ib:ie,jb:je,kb:ke)) ! u at previous time step
! integrate u in x
do k=kb,ke
uout(k) = sum(uvol(ib:ie,k)*dxf(ib:ie))
uoutold(k) = sum(uvolold(ib:ie,k)*dxf(ib:ie))
end do
! integrate u in z
do k=kb,ke
uout(k) = rk3coef*uout(k)*dzf(k)
uoutold(k) = uoutold(k)*dzf(k)
end do
uoutflow = sum(uout(kb:ke))
uflowrateold = sum(uoutold(kb:ke))
! average over fluid volume
uoutflow = uoutflow/fluidvol
uflowrateold = uflowrateold/fluidvol
! flow correction to match outflow rate
udef = uflowrate - (uoutflow + uflowrateold)
do k = kb,ke
do j = jb,je
do i = ib,ie
up(i,j,k) = up(i,j,k) + udef*rk3coefi
end do
end do
end do
end if
if ((.not.linoutflow) .and. (lvoutflowr)) then
rk3coef = dt / (4. - dble(rk3step))
rk3coefi = 1 / rk3coef
vdef = 0.
vout = 0.
voutflow = 0.
voutold = 0.
! integrate v fixed at outlet je along x
if (myid==nprocs-1) then
do k=kb,ke
vout(k) = sum(vp(ib:ie,je,k)*IIv(ib:ie,je,k)*dxf(ib:ie)) ! v tendency at previous time step
voutold(k) = sum(vm(ib:ie,je,k)*IIv(ib:ie,je,k)*dxf(ib:ie)) ! v at previous time step
end do
end if
call MPI_BCAST(vout,ke-kb+1,MY_REAL ,nprocs-1,comm3d,mpierr)
call MPI_BCAST(voutold,ke-kb+1,MY_REAL ,nprocs-1,comm3d,mpierr)
! integrate v in z
do k=kb,ke
vout(k) = rk3coef*vout(k)*dzf(k)
voutold(k) = voutold(k)*dzf(k)
end do
voutflow = sum(vout(kb:ke))
vflowrateold = sum(voutold(kb:ke))
! average over outflow area
voutflow = voutflow/voutarea
vflowrateold = vflowrateold/voutarea
! flow correction to match outflow rate
vdef = vflowrate - (voutflow + vflowrateold)
do k = kb,ke
do j = jb,je
do i = ib,ie
vp(i,j,k) = vp(i,j,k) + vdef*rk3coefi
end do
end do
end do
elseif ((.not.linoutflow) .and. (lvvolflowr)) then
rk3coef = dt / (4. - dble(rk3step))
rk3coefi = 1 / rk3coef
vdef = 0.
vout = 0.
voutflow = 0.
voutold = 0.
vvol = 0.
vvolold = 0.
! integrate v in y
call sumy_ibm(vvol,vp(ib:ie,jb:je,kb:ke)*dy,ib,ie,jb,je,kb,ke,IIv(ib:ie,jb:je,kb:ke)) ! v tendency at previous time step
call sumy_ibm(vvolold,vm(ib:ie,jb:je,kb:ke)*dy,ib,ie,jb,je,kb,ke,IIv(ib:ie,jb:je,kb:ke)) ! v at previous time step
! integrate v in x
do k=kb,ke
vout(k) = sum(vvol(ib:ie,k)*dxf(ib:ie))
voutold(k) = sum(vvolold(ib:ie,k)*dxf(ib:ie))
end do
! integrate v in z
do k=kb,ke
vout(k) = rk3coef*vout(k)*dzf(k)
voutold(k) = voutold(k)*dzf(k)
end do
voutflow = sum(vout(kb:ke))
vflowrateold = sum(voutold(kb:ke))
! average over fluid volume
voutflow = voutflow/fluidvol
vflowrateold = vflowrateold/fluidvol
! flow correction to match outflow rate
vdef = vflowrate - (voutflow + vflowrateold)
do k = kb,ke
do j = jb,je
do i = ib,ie
vp(i,j,k) = vp(i,j,k) + vdef*rk3coefi
end do
end do
end do
end if
end subroutine masscorr