subroutine wallfunheat use modglobal, only : ib, ie, ih, jb, je, jh, kb, ke, kh, xf, yf, zf, xh, yh, zh, dx, dy, dzh, eps1, & xhat, yhat, zhat, vec0, fkar, ltempeq, lmoist, iwalltemp, iwallmoist, lEB, lwritefac, nfcts, rk3step, totheatflux, totqflux use modfields, only : u0, v0, w0, thl0, thlp, qt0, qtp, pres0 use initfac, only : facT, facz0, facz0h, facnorm, fachf, facef, facqsat, fachurel, facf, faclGR, faca use modmpi, only : comm3d, mpi_sum, mpierr, my_real use modsurfdata, only : z0, z0h use modibmdata, only : bctfxm, bctfxp, bctfym, bctfyp, bctfz use decomp_2d, only : zstart integer i, j, k, n, m, sec, fac real :: dist, flux, area, vol, tempvol, Tair, Tsurf, utan, cth, htc, cveg, hurel, qtair, qwall, resa, resc, ress, xrec, yrec, zrec real, dimension(3) :: uvec, norm, span, strm real, dimension(1:nfcts) :: fac_htc_loc, fac_cth_loc, fac_pres_loc logical :: valid fac_htc_loc = 0. fac_cth_loc = 0. fac_pres_loc = 0. do sec = 1,bound_info_c%nfctsecsrank ! sec = bound_info_c%fctsecsrank(m) ! index of section !n = bound_info_c%secbndptids(sec) ! index of boundary point fac = bound_info_c%secfacids_loc(sec) ! index of facet area = bound_info_c%secareas_loc(sec) ! area norm = facnorm(fac,:) ! i = bound_info_c%bndpts(n,1) - zstart(1) + 1 ! should be on this rank! ! j = bound_info_c%bndpts(n,2) - zstart(2) + 1 ! should be on this rank! ! k = bound_info_c%bndpts(n,3) - zstart(3) + 1 ! should be on this rank! i = bound_info_c%secbndpts_loc(sec,1) - zstart(1) + 1 ! should be on this rank! j = bound_info_c%secbndpts_loc(sec,2) - zstart(2) + 1 ! should be on this rank! k = bound_info_c%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 wallfunheat", i, j stop 1 end if fac_pres_loc(fac) = fac_pres_loc(fac) + pres0(i,j,k) * area ! output pressure on facets if (bound_info_c%lskipsec_loc(sec)) cycle !if (facz0(fac) < eps1) cycle if (bound_info_c%lcomprec_loc(sec) .or. lnorec) then ! section aligned with grid - use this cell's velocity uvec = interp_velocity_c(i, j, k) Tair = thl0(i,j,k) qtair = qt0(i,j,k) dist = bound_info_c%bnddst_loc(sec) else ! use velocity at reconstruction point xrec = bound_info_c%recpts_loc(sec,1) yrec = bound_info_c%recpts_loc(sec,2) zrec = bound_info_c%recpts_loc(sec,3) uvec(1) = trilinear_interp_var(u0, bound_info_c%recids_u_loc(sec,:), xh, yf, zf, xrec, yrec, zrec) uvec(2) = trilinear_interp_var(v0, bound_info_c%recids_v_loc(sec,:), xf, yh, zf, xrec, yrec, zrec) uvec(3) = trilinear_interp_var(w0, bound_info_c%recids_w_loc(sec,:), xf, yf, zh, xrec, yrec, zrec) Tair = trilinear_interp_var(thl0, bound_info_c%recids_c_loc(sec,:), xf, yf, zf, xrec, yrec, zrec) qtair = trilinear_interp_var( qt0, bound_info_c%recids_c_loc(sec,:), xf, yf, zf, xrec, yrec, zrec) ! dist = bound_info_c%bnddst(sec) + norm2((/xrec - xf(bound_info_c%bndpts(n,1)), & ! yrec - yf(bound_info_c%bndpts(n,2)), & ! zrec - zf(bound_info_c%bndpts(n,3))/)) dist = bound_info_c%bnddst_loc(sec) + norm2((/xrec - xf(bound_info_c%secbndpts_loc(sec,1)), & yrec - yf(bound_info_c%secbndpts_loc(sec,2)), & zrec - zf(bound_info_c%secbndpts_loc(sec,3))/)) end if if (log(dist/facz0(fac)) <= 1.) then cycle !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 ! Sensible heat if (ltempeq) then if (iwalltemp == 1) then ! probably remove this eventually, only relevant to grid-aligned facets !if (all(abs(norm - xhat) < eps1)) then if (is_equal(norm, xhat)) then flux = bctfxp !elseif (all(abs(norm + xhat) < eps1)) then elseif (is_equal(norm, -xhat)) then flux = bctfxm !elseif (all(abs(norm - yhat) < eps1)) then elseif (is_equal(norm, yhat)) then flux = bctfyp !elseif (all(abs(norm + yhat) < eps1)) then elseif (is_equal(norm, -yhat)) then flux = bctfxm !elseif (all(abs(norm - zhat) < eps1)) then elseif (is_equal(norm, zhat)) then flux = bctfz end if elseif (iwalltemp == 2) then call heat_transfer_coef_flux(utan, dist, facz0(fac), facz0h(fac), Tair, facT(fac, 1), cth, flux, htc) fac_cth_loc(fac) = fac_cth_loc(fac) + cth * area ! output heat transfer coefficients on facets fac_htc_loc(fac) = fac_htc_loc(fac) + htc * area ! output heat transfer coefficients on facets end if ! flux [Km/s] ! fluid volumetric sensible heat source/sink = flux * area / volume [K/s] ! facet sensible heat flux = volumetric heat capacity of air * flux * sectionarea / facetarea [W/m^2] thlp(i,j,k) = thlp(i,j,k) - flux * area / (dx*dy*dzh(k)) if (lEB) then totheatflux = totheatflux + flux*area ! [Km^3s^-1] This sums the flux over all facets fachf(fac) = fachf(fac) + flux * area ! [Km^2/s] (will be divided by facetarea(fac) in modEB) end if end if ! Latent heat if (lmoist .and. faclGR(fac)) then if (iwallmoist == 1) then ! probably remove this eventually, only relevant to grid-aligned facets if (is_equal(norm, xhat)) then flux = bcqfxp elseif (is_equal(norm, -xhat)) then flux = bcqfxm elseif (is_equal(norm, yhat)) then flux = bcqfyp elseif (is_equal(norm, -yhat)) then flux = bcqfym elseif (is_equal(norm, zhat)) then flux = bcqfz end if elseif (iwallmoist == 2) then if (abs(htc*abs(utan)) > 0.) then qwall = facqsat(fac) ! saturation humidity hurel = fachurel(fac) ! relative humidity resa = 1./(htc*abs(utan)) ! aerodynamic resistance resc = facf(fac,4) ! canopy resistance ress = facf(fac,5) ! soil resistance cveg = 0.8 ! vegetation fraction flux = moist_flux(cveg, resa, qtair, qwall, hurel, resc, ress) end if end if ! flux [kg/kg m/s] ! fluid volumetric latent heat source/sink = flux * area / volume [kg/kg / s] ! facet latent heat flux = volumetric heat capacity of air * flux * sectionarea / facetarea [W/m^2] totqflux = totqflux + flux*area ! [Km^3s^-1] This sums the flux over all facets qtp(i,j,k) = qtp(i,j,k) - flux * area / (dx*dy*dzh(k)) if (lEB) then facef(fac) = facef(fac) + flux * area ! [Km^2/s] (will be divided by facetarea(fac) in modEB) end if end if end do if (lwritefac .and. rk3step==3) then fac_cth_loc(1:nfcts) = fac_cth_loc(1:nfcts) / faca(1:nfcts) fac_htc_loc(1:nfcts) = fac_htc_loc(1:nfcts) / faca(1:nfcts) fac_pres_loc(1:nfcts) = fac_pres_loc(1:nfcts) / faca(1:nfcts) call MPI_ALLREDUCE(fac_cth_loc(1:nfcts), fac_cth(1:nfcts), nfcts, MY_REAL, MPI_SUM, comm3d, mpierr) call MPI_ALLREDUCE(fac_htc_loc(1:nfcts), fac_htc(1:nfcts), nfcts, MY_REAL, MPI_SUM, comm3d, mpierr) call MPI_ALLREDUCE(fac_pres_loc(1:nfcts), fac_pres(1:nfcts), nfcts, MY_REAL, MPI_SUM, comm3d, mpierr) end if end subroutine wallfunheat