subroutine wallfunmom(dir, rhs, bound_info)
use modglobal, only : ib, ie, ih, jb, je, jh, kb, ke, kh, xf, yf, zf, xh, yh, zh, &
eps1, fkar, dx, dy, dzf, iwallmom, xhat, yhat, zhat, vec0, nfcts, lwritefac, rk3step
use modfields, only : u0, v0, w0, thl0, tau_x, tau_y, tau_z
use initfac, only : facT, facz0, facz0h, facnorm, faca
use decomp_2d, only : zstart
use modmpi, only : comm3d, mpi_sum, mpierr, my_real
real, intent(in) :: dir(3)
real, intent(inout) :: rhs(ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh)
type(bound_info_type) :: bound_info
integer i, j, k, n, m, sec, pt, fac
real dist, stress, stress_dir, stress_aligned, area, vol, momvol, Tair, Tsurf, x, y, z, &
utan, udir, ctm, a, a_is, a_xn, a_yn, a_zn, stress_ix, stress_iy, stress_iz, xrec, yrec, zrec
real, dimension(3) :: uvec, norm, strm, span, stressvec
logical :: valid
real, dimension(1:nfcts) :: fac_tau_loc, fac_tau
!real, dimension(:), allocatable :: fac_tau, fac_pres
procedure(interp_velocity), pointer :: interp_velocity_ptr => null()
procedure(interp_temperature), pointer :: interp_temperature_ptr => null()
select case(alignment(dir))
case(1)
interp_velocity_ptr => interp_velocity_u
interp_temperature_ptr => interp_temperature_u
case(2)
interp_velocity_ptr => interp_velocity_v
interp_temperature_ptr => interp_temperature_v
case(3)
interp_velocity_ptr => interp_velocity_w
interp_temperature_ptr => interp_temperature_w
end select
fac_tau_loc = 0.
do sec = 1,bound_info%nfctsecsrank
!sec = bound_info%fctsecsrank(m) ! index of section
!n = bound_info%secbndptids(sec) ! index of boundary point
area = bound_info%secareas_loc(sec) ! area of section
fac = bound_info%secfacids_loc(sec) ! index of facet
norm = facnorm(fac,:) ! facet normal
if (bound_info%lskipsec_loc(sec)) cycle
!if (facz0(fac) < eps1) cycle
! i = bound_info%bndpts(n,1) - zstart(1) + 1
! j = bound_info%bndpts(n,2) - zstart(2) + 1
! k = bound_info%bndpts(n,3) - zstart(3) + 1
i = bound_info%secbndpts_loc(sec,1) - zstart(1) + 1 ! should be on this rank!
j = bound_info%secbndpts_loc(sec,2) - zstart(2) + 1 ! should be on this rank!
k = bound_info%secbndpts_loc(sec,3) - zstart(3) + 1 ! should be on this rank!
if ((i < ib) .or. (i > ie) .or. (j < jb) .or. (j > je)) then
write(*,*) "problem in wallfunmom", alignment(dir), bound_info%secbndpts_loc(sec,1), bound_info%secbndpts_loc(sec,2)
stop 1
end if
if (bound_info%lcomprec_loc(sec) .or. lnorec) then
uvec = interp_velocity_ptr(i, j, k)
if (iwallmom == 2) then
Tair = interp_temperature_ptr(i, j, k)
end if
dist = bound_info%bnddst_loc(sec)
else
xrec = bound_info%recpts_loc(sec,1)
yrec = bound_info%recpts_loc(sec,2)
zrec = bound_info%recpts_loc(sec,3)
uvec(1) = trilinear_interp_var(u0, bound_info%recids_u_loc(sec,:), xh, yf, zf, xrec, yrec, zrec)
uvec(2) = trilinear_interp_var(v0, bound_info%recids_v_loc(sec,:), xf, yh, zf, xrec, yrec, zrec)
uvec(3) = trilinear_interp_var(w0, bound_info%recids_w_loc(sec,:), xf, yf, zh, xrec, yrec, zrec)
if (iwallmom == 2) Tair = trilinear_interp_var(thl0, bound_info%recids_c_loc(sec,:), xf, yf, zf, xrec, yrec, zrec)
dist = bound_info%bnddst_loc(sec) + norm2((/xrec - xf(bound_info%secbndpts_loc(sec,1)), &
yrec - yf(bound_info%secbndpts_loc(sec,2)), &
zrec - zf(bound_info%secbndpts_loc(sec,3))/))
end if
if (log(dist/facz0(fac)) <= 1.) then
cycle ! ideally would set a value for dist that gives a resonable (large) flux
!dist = facz0(fac)+facz0h(fac)
end if
if (is_equal(uvec, vec0)) cycle
call local_coords(uvec, norm, span, strm, valid)
if (.not. valid) cycle
utan = dot_product(uvec, strm)
!utan = max(0.01, utan) ! uDALES 1
! calcualate momentum transfer coefficient
! make into interface somehow? because iwallmom doesn't change in the loop
if (iwallmom == 2) then ! stability included
ctm = mom_transfer_coef_stability(utan, dist, facz0(fac), facz0h(fac), Tair, facT(fac,1))
else if (iwallmom == 3) then ! neutral
ctm = mom_transfer_coef_neutral(dist, facz0(fac))
end if
stress = ctm * utan**2
if (bound_info%lcomprec_loc(sec)) then
a = dot_product(dir, strm)
stress_dir = a * stress
else
! Rotation from local (strm,span,norm) to global (xhat,yhat,zhat) basis
! \tau'_ij = a_ip a_jq \tau_pq
! \tau_pq in local coordinates is something like \tau \delta_13, because we only have \tau_{strm,norm})
a_is = dot_product(dir, strm)
a_xn = dot_product(xhat, norm)
a_yn = dot_product(yhat, norm)
a_zn = dot_product(zhat, norm)
stress_ix = a_is * a_xn * stress
stress_iy = a_is * a_yn * stress
stress_iz = a_is * a_zn * stress
stressvec(1) = stress_ix
stressvec(2) = stress_iy
stressvec(3) = stress_iz
stress_dir = norm2(stressvec)
end if
stress_dir = sign(stress_dir, dot_product(uvec, dir))
vol = dx*dy*dzf(k)
momvol = stress_dir * area / vol
rhs(i,j,k) = rhs(i,j,k) - momvol
fac_tau_loc(fac) = fac_tau_loc(fac) + stress_dir * area ! output stresses on facets
end do
if (lwritefac .and. rk3step==3) then
fac_tau_loc(1:nfcts) = fac_tau_loc(1:nfcts) / faca(1:nfcts)
call MPI_ALLREDUCE(fac_tau_loc(1:nfcts), fac_tau(1:nfcts), nfcts, MY_REAL, MPI_SUM, comm3d, mpierr)
select case(alignment(dir))
case(1)
fac_tau_x = fac_tau
case(2)
fac_tau_y = fac_tau
case(3)
fac_tau_z = fac_tau
end select
end if
! Do time-averaging like in modEB
end subroutine wallfunmom