pure subroutine update_reservoir_zonal_east(hTr, h_layer, mass_flux_x, &
tres_e, T_data, &
nx_total, ny_total, nz, &
i_e, j0, j1, &
L_out, L_in, dt, it, n_tr)
!! Update reservoir for the east open edge.
!! Wall-face index for mass flux = i_e+1 (east face of the last physical cell i_e).
!! Interior cell = i_e.
!! Outward normal at east = +x: u_n = +mass_flux_x(i_e+1,j,k)/max(h,hmin).
integer, intent(in) :: nx_total, ny_total, nz, n_tr, it
integer, intent(in) :: i_e, j0, j1
real(wp), intent(in) :: hTr(nx_total, ny_total, nz)
real(wp), intent(in) :: h_layer(nx_total, ny_total, nz)
real(wp), intent(in) :: mass_flux_x(nx_total + 1, ny_total, nz)
real(wp), intent(inout) :: tres_e(ny_total, nz, n_tr)
real(wp), intent(in) :: T_data
real(wp), intent(in) :: L_out, L_in, dt
integer :: j, k
real(wp) :: h_int, T_int, u_n, c_out, c_in, tres_old, denom
do concurrent(j=j0:j1, k=1:nz) local(h_int, T_int, u_n, c_out, c_in, tres_old, denom)
h_int = max(h_layer(i_e, j, k), RES_H_MIN)
T_int = hTr(i_e, j, k)/h_int
! Outward normal at east = +x.
u_n = mass_flux_x(i_e + 1, j, k)/h_int
tres_old = tres_e(j, k, it)
if (L_out == 0.0_wp .and. u_n > 0.0_wp) then
tres_e(j, k, it) = T_int
else if (L_in == 0.0_wp .and. u_n < 0.0_wp) then
tres_e(j, k, it) = T_data
else
c_out = merge(max(0.0_wp, u_n)*dt/L_out, 0.0_wp, L_out > 0.0_wp)
c_in = merge(max(0.0_wp, -u_n)*dt/L_in, 0.0_wp, L_in > 0.0_wp)
denom = 1.0_wp + c_out + c_in
tres_e(j, k, it) = (tres_old + c_out*T_int + c_in*T_data)/denom
end if
end do
end subroutine update_reservoir_zonal_east