subroutine iolet
use modglobal, only:dxhi, dxhci, xh, zh, ib, ie, jb, je, ih, jh, kb, ke, kh, nsv, rk3step, dt, iinletgen, ltempeq, lmoist, ihc, idriver, dy, dzf, jtot, zh, lsdriver
use modfields, only:u0, um, v0, vm, w0, wm, e120, e12m, thl0, thlm, qt0, qtm, sv0, svm, uprof, vprof, e12prof, thlprof, &
qtprof, svprof, uouttot, wouttot
use modmpi, only:excjs, myid, slabsum
use modinletdata, only:u0inletbcold, v0inletbcold, w0inletbcold, uminletbc, vminletbc, wminletbc, totaluold, &
t0inletbcold, tminletbc, u0driver, v0driver, w0driver, e120driver, thl0driver, qt0driver, umdriver, vmdriver, wmdriver,&
e12mdriver, thlmdriver, qtmdriver, sv0driver, svmdriver
real rk3coef
real, dimension(kb:ke) :: uin
integer n, i, j, k, m
rk3coef = dt/(4.-dble(rk3step))
! Inlet boundary is located at ib (not ib-1)!
! Inlet
if ((iinletgen == 1) .or. (iinletgen == 2)) then
do j = jb, je
do k = kb, ke
u0(ib, j, k) = u0inletbcold(j, k)
um(ib, j, k) = uminletbc(j, k)
u0(ib - 1, j, k) = 2*u0(ib, j, k) - u0(ib + 1, j, k)
um(ib - 1, j, k) = 2*um(ib, j, k) - um(ib + 1, j, k)
v0(ib - 1, j, k) = v0inletbcold(j, k)
vm(ib - 1, j, k) = vminletbc(j, k)
! to be changed in the future: e12 should be taken from recycle plane!
e120(ib - 1, j, k) = e120(ib, j, k) ! extrapolate e12 from interior
e12m(ib - 1, j, k) = e12m(ib, j, k) ! extrapolate e12 from interior
do n = 1, nsv
do m = 1, ihc
sv0(ib - m, j, k, n) = 2*svprof(k, n) - sv0(ib + (m - 1), j, k, n)
svm(ib - m, j, k, n) = 2*svprof(k, n) - svm(ib + (m - 1), j, k, n)
enddo
enddo
end do
do k = kb, ke + 1
w0(ib - 1, j, k) = w0inletbcold(j, k)
wm(ib - 1, j, k) = wminletbc(j, k)
end do
end do
! Heat
if (ltempeq) then
do k = kb, ke
do j = jb, je
thl0(ib - 1, j, k) = t0inletbcold(j, k)
thlm(ib - 1, j, k) = tminletbc(j, k)
end do
end do
end if
if (lmoist) then
do k = kb, ke
do j = jb, je
qt0(ib - 1, j, k) = 2*qtprof(k) - qt0(ib, j, k) !watch!
qtm(ib - 1, j, k) = 2*qtprof(k) - qtm(ib, j, k)
end do
end do
end if
! Driver inlet
elseif (idriver == 2) then
do j=jb-1,je+1
do k=kb,ke !tg3315 removed +1 following above...
u0(ib,j,k)=u0driver(j,k) !max(0.,u0driver(j,k))
um(ib,j,k)=umdriver(j,k) !max(0.,umdriver(j,k))
u0(ib-1,j,k)= u0driver(j,k) !max(0.,2.*u0(ib,j,k)-u0(ib+1,j,k))
um(ib-1,j,k)= umdriver(j,k) !max(0.,2.*um(ib,j,k)-um(ib+1,j,k))
v0(ib,j,k) = v0driver(j,k) !max(0.,v0driver(j,k))
vm(ib,j,k) = vmdriver(j,k) !max(0.,vmdriver(j,k))
v0(ib-1,j,k) = v0driver(j,k) !max(0.,v0driver(j,k))
vm(ib-1,j,k) = vmdriver(j,k) !max(0.,vmdriver(j,k))
! to be changed in the future: e12 should be taken from recycle plane!
!e120(ib-1,j,k) = e120driver(j,k) ! extrapolate e12 from interior
!e12m(ib-1,j,k) = e12mdriver(j,k) ! extrapolate e12 from interior
if (lsdriver) then
do n=1,nsv
do m = 1,ihc
sv0(ib-m,j,k,n) = sv0driver(j,k,n)
svm(ib-m,j,k,n) = svmdriver(j,k,n)
!sv0(ib-m,j,k,n) = 2*svprof(k,n) - sv0(ib+(m-1),j,k,n)
!svm(ib-m,j,k,n) = 2*svprof(k,n) - svm(ib+(m-1),j,k,n)
enddo
sv0(ib,j,k,n) = sv0driver(j,k,n)
svm(ib,j,k,n) = svmdriver(j,k,n)
enddo
end if
end do
do k=kb,ke+1
w0(ib-1,j,k) = w0driver(j,k) !max(0.,w0driver(j,k))
wm(ib-1,j,k) = wmdriver(j,k) !max(0.,wmdriver(j,k))
w0(ib,j,k) = w0driver(j,k) !max(0.,w0driver(j,k))
wm(ib,j,k) = wmdriver(j,k) !max(0.,wmdriver(j,k))
end do
end do
! Heat
if (ltempeq ) then
do j=jb-1,je+1
do k=kb,ke+1
thl0(ib,j,k) = thl0driver(j,k)
thlm(ib,j,k) = thlmdriver(j,k)
thl0(ib-1,j,k) = thl0driver(j,k)
thlm(ib-1,j,k) = thlmdriver(j,k)
!thlm(ib-1,j,k) = 2*thlm(ib,j,k) - thlm(ib+1,j,k)
!thl0(ib-1,j,k) = 2*thl0(ib,j,k) - thl0(ib+1,j,k)
end do
end do
end if
if (lmoist ) then
do j=jb-1,je+1
do k=kb,ke+1
qt0(ib,j,k) = qt0driver(j,k)
! qt0(ib-1,j,k) = 2*qtprof(k) - qt0(ib,j,k)
qtm(ib,j,k) = qtmdriver(j,k)
! qtm(ib-1,j,k) = 2*qtprof(k) - qtm(ib,j,k)
qt0(ib-1,j,k) = qt0driver(j,k)
! qt0(ib-1,j,k) = 2*qtprof(k) - qt0(ib,j,k) !watch!
qtm(ib-1,j,k) = qtmdriver(j,k)
! qtm(ib-1,j,k) = 2*qtprof(k) - qtm(ib,j,k)
! qt0(ib-1,j,k) = 2*qt0(ib,j,k) - qt0(ib+1,j,k)
! qtm(ib-1,j,k) = 2*qtm(ib,j,k) - qtm(ib+1,j,k)
end do
end do
end if
else ! (if iinetgen==0)
do j = jb - 1, je + 1
do k = kb, ke + 1
! Momentum
u0(ib, j, k) = uprof(k)
um(ib, j, k) = uprof(k)
v0(ib - 1, j, k) = 2*vprof(k) - v0(ib, j, k) ! (v(ib)+v(ib-1))/2 = vprof
vm(ib - 1, j, k) = 2*vprof(k) - vm(ib, j, k) ! (v(ib)+v(ib-1))/2 = vprof
e120(ib - 1, j, k) = 2*e12prof(k) - e120(ib, j, k) ! (e12(ib)+e12(ib-1))/2=e12prof
e12m(ib - 1, j, k) = 2*e12prof(k) - e12m(ib, j, k) ! (e12(ib)+e12(ib-1))/2=e12prof
do n = 1, nsv
do m = 1, ihc
sv0(ib - m, j, k, n) = 2*svprof(k, n) - sv0(ib + (m - 1), j, k, n)
svm(ib - m, j, k, n) = 2*svprof(k, n) - svm(ib + (m - 1), j, k, n)
end do
end do
!end if
enddo
enddo
! Heat
if (ltempeq) then
do j = jb - 1, je + 1
do k = kb, ke + 1
thl0(ib - 1, j, k) = 2*thlprof(k) - thl0(ib, j, k)
thlm(ib - 1, j, k) = 2*thlprof(k) - thlm(ib, j, k)
end do
end do
end if
if (lmoist) then
do j = jb - 1, je + 1
do k = kb, ke + 1
qt0(ib - 1, j, k) = 2*qtprof(k) - qt0(ib, j, k)
qtm(ib - 1, j, k) = 2*qtprof(k) - qtm(ib, j, k)
end do
end do
end if
u0(ib - 1, :, :) = 2*u0(ib, :, :) - u0(ib + 1, :, :) ! (u(ib+1)+u(ib-1))/2 = u(ib)
um(ib - 1, :, :) = 2*um(ib, :, :) - um(ib + 1, :, :) ! (u(ib+1)+u(ib-1))/2 = u(ib)
w0(ib - 1, :, :) = -w0(ib, :, :) ! (w(ib)+w(ib-1))/2 = 0
wm(ib - 1, :, :) = -wm(ib, :, :)
end if ! iinletgen==1 .or. iinletgen==2
! tg3315 added to ensure that uouttot matches driven inflow regardless of forcing
! set up assuming we have a block at lowest cell kb
if (idriver==2) then
uin = 0.
uouttot = 0.
call slabsum(uin,kb,ke,um,ib-ih,ie+ih,jb-jh,je+jh,kb-kh,ke+kh,ie,ie,jb,je,kb,ke) ! determine horizontal (j) average outflow velocity old
do k=kb,ke
uin(k) = uin(k)*dzf(k)*dy ! flow rate through each slab at ib
end do
uouttot = sum(uin(kb+1:ke))/((zh(ke+1)-zh(kb+1))*jtot*dy) ! convective outflow velocity
end if
! Outlet
! Momentum
v0(ie + 1, :, :) = v0(ie, :, :) - (v0(ie + 1, :, :) - v0(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
w0(ie + 1, :, :) = w0(ie, :, :) - (w0(ie + 1, :, :) - w0(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
vm(ie + 1, :, :) = vm(ie, :, :) - (vm(ie + 1, :, :) - vm(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
wm(ie + 1, :, :) = wm(ie, :, :) - (wm(ie + 1, :, :) - wm(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
e120(ie + 1, :, :) = e120(ie, :, :) - (e120(ie + 1, :, :) - e120(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
e12m(ie + 1, :, :) = e12m(ie, :, :) - (e12m(ie + 1, :, :) - e12m(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
! Heat
if (ltempeq) then
thl0(ie + 1, :, :) = thl0(ie, :, :) - (thl0(ie + 1, :, :) - thl0(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
thlm(ie + 1, :, :) = thlm(ie, :, :) - (thlm(ie + 1, :, :) - thlm(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
end if
if (lmoist) then
qt0(ie + 1, :, :) = qt0(ie, :, :) - (qt0(ie + 1, :, :) - qt0(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
qtm(ie + 1, :, :) = qtm(ie, :, :) - (qtm(ie + 1, :, :) - qtm(ie, :, :))*dxhi(ie + 1)*rk3coef*uouttot
end if
! tg3315 !changed dxhi to dxhci!?
do n = 1, nsv
sv0(ie + 1, :, :, n) = sv0(ie, :, :, n) - (sv0(ie + 1, :, :, n) - sv0(ie, :, :, n))*dxhci(ie + 1)*rk3coef*uouttot
svm(ie + 1, :, :, n) = svm(ie, :, :, n) - (svm(ie + 1, :, :, n) - svm(ie, :, :, n))*dxhci(ie + 1)*rk3coef*uouttot
end do
return
end subroutine iolet