subroutine trees use modglobal, only : ib,ie,jb,je,kb,ke,ih,jh,dzf,xh,dxf,numol,prandtlmol,rlv,cp,tree,& ntrees,ltrees,itot,jtot,cd,ud,lmoist,nsv,dxf,dy,dzf,dzfi,zf,dy,zh,& ltempeq,pref0,r_s,lad,rhoa,ntree_max,lsize,Qstar,dQdt,BCxs,tr_A,& rslabs,kmax,rv,rd use modfields, only : um,vm,wm,thlm,qtm,svp,up,vp,wp,thlp,qtp,svm,Rn,qc,qa,thlm,clai,& ladzf,ladzh,IIcs,tr_omega,& tr_u,tr_v,tr_w,tr_qt,tr_qtR,tr_qtA,tr_thl,tr_sv,thlpcar use modmpi, only : myidx,myidy,mpi_sum,mpierr,comm3d,mpierr,my_real,nprocx,nprocy use modsurfdata,only : wtsurf,wttop,wqtop use modibmdata, only : bctfz use modsubgriddata, only : ekh implicit none integer :: i,j,k,n,m,il,iu,jl,ju,kl,ku real :: e_sat,e_vap,qh,qe,shade,r_a,s,D,numoli,gam,Rq,qhmin,qhmax,omega real :: Vq_dum,Vq_dum2,Vq_dum3,VT_dum,VT_dum2,VT_dum3,wTbot_dum,wTbot_dum2,wTbot,wTbot_dum3 numoli = 1/numol gam = (cp*pref0*rv)/(rlv*rd) qhmin=0.;qhmax=0. if (ltrees .eqv. .false.) return ! dummy variables used for steady-state BC calculation VT_dum2 = 0. Vq_dum2 = 0. wTbot_dum2 = 0. ! loop over all tree canopies do n = 1,ntrees ! drag in z-direction il = tree(n,1) - myidx*itot/nprocx iu = tree(n,2) - myidx*itot/nprocx kl = tree(n,5) ku = tree(n,6) + 1 jl = tree(n,3) - myidy*jtot/nprocy ju = tree(n,4) - myidy*jtot/nprocy if (iu < ib .or. il > ie .or. ju < jb .or. jl > je) then cycle else if (iu > ie) iu=ie if (il < ib) il=ib if (ju > je) ju=je if (jl < jb) jl=jb do k = kl,ku do j = jl,ju do i = il,iu tr_w(i,j,k) = - cd * ladzh(ntree_max-(ku-k)+1) * wm(i,j,k) * & sqrt( wm(i,j,k)**2 & + (0.25*(um(i,j,k) + um(i+1,j,k) + um(i,j,k-1) + um(i+1,j,k-1)))**2 & + (0.25*(vm(i,j,k) + vm(i,j+1,k) + vm(i,j,k-1) + vm(i,j+1,k-1)))**2 ) end do end do end do end if ! drag in y-direction il = tree(n,1) - myidx*itot/nprocx iu = tree(n,2) - myidx*itot/nprocx kl = tree(n,5) ku = tree(n,6) jl = tree(n,3) - myidy*jtot/nprocy ju = tree(n,4) + 1 - myidy*jtot/nprocy if (iu < ib .or. il > ie .or. ju < jb .or. jl > je) then cycle else if (iu > ie) iu=ie if (il < ib) il=ib if (ju > je) ju=je if (jl < jb) jl=jb do k = kl,ku do j = jl,ju do i = il,iu tr_v(i,j,k) = - cd * ladzf(ntree_max-(ku-k)) * vm(i,j,k) * & sqrt( vm(i,j,k)**2 & + (0.25*(um(i,j,k) + um(i+1,j,k) + um(i,j-1,k) + um(i+1,j-1,k)))**2 & + (0.25*(wm(i,j,k) + wm(i,j,k+1) + wm(i,j-1,k) + wm(i,j-1,k+1)))**2 ) end do end do end do end if ! drag in x-direction il = tree(n,1) - myidx*itot/nprocx iu = tree(n,2) + 1 - myidx*itot/nprocx kl = tree(n,5) ku = tree(n,6) jl = tree(n,3) - myidy*jtot/nprocy ju = tree(n,4) - myidy*jtot/nprocy if (iu < ib .or. il > ie .or. ju < jb .or. jl > je) then cycle else if (iu > ie) iu=ie if (il < ib) il=ib if (ju > je) ju=je if (jl < jb) jl=jb do k = kl,ku do j = jl,ju do i = il,iu tr_u(i,j,k) = - cd * ladzf(ntree_max-(ku-k)) * um(i,j,k) * & sqrt( um(i,j,k)**2 & + (0.25*(vm(i,j,k) + vm(i,j+1,k) + vm(i-1,j,k) + vm(i-1,j+1,k)))**2 & + (0.25*(wm(i,j,k) + wm(i,j,k+1) + wm(i-1,j,k) + wm(i-1,j,k+1)))**2 ) end do end do end do end if ! scalar volumetric sources/sinks ! Canopy EB if (lmoist .and. ltempeq) then il = tree(n,1) - myidx*itot/nprocx iu = tree(n,2) - myidx*itot/nprocx kl = tree(n,5) ku = tree(n,6) jl = tree(n,3) - myidy*jtot/nprocy ju = tree(n,4) - myidy*jtot/nprocy if (iu < ib .or. il > ie .or. ju < jb .or. jl > je) then cycle else if (iu > ie) iu=ie if (il < ib) il=ib if (ju > je) ju=je if (jl < jb) jl=jb do k = kl,ku do j = jl,ju do i = il,iu ! psychometrics ! saturation vapour pressure pressure e_sat = 610.8*exp((17.27*(thlm(i,j,k)-273.15))/(thlm(i,j,k)-35.85)) ! water vapour partial pressure e_vap = (qtm(i,j,k) * pref0) / (0.378 * qtm(i,j,k) + 0.622) ! qtcheck.m and https://earthscience.stackexchange.com/questions/2360/how-do-i-convert-specific-humidity-to-relative-humidity ! vapour pressure deficit D = max(e_sat - e_vap,0.) ! output warning - cannot physically occur if (e_sat<e_vap) then write(*,*) 'D, e_sat, e_vap, thlm, qtm', D, e_sat, e_vap, thlm(i,j,k), qtm(i,j,k) end if ! slope of the curve relating saturation vapour pressure to temperature s = (4098*e_sat)/((thlm(i,j,k)-35.85)**2) ! aerodynamic resistance r_a = 130*sqrt(lsize/(sqrt((0.5*(um(i,j,k)+um(i+1,j,k)))**2+(0.5*(vm(i,j,k)+vm(i,j+1,k)))**2+(0.5*(wm(i,j,k)+wm(i,j,k+1)))**2))) ! decoupling factor omega = 1/(1 + 2*(gam/(s+2*gam)) * ((r_s)/r_a) ) ! latent heat qe = omega*(s/(s+2*gam))*(qa(ntree_max-(ku-k))/(dzf(k)*ladzf(ntree_max-(ku-k)))) + (1-omega)*(1/(gam*(r_s)))*rhoa*cp*D ! sensible heat qh = qa(ntree_max-(ku-k))/(dzf(k)*ladzf(ntree_max-(ku-k))) - qe ! volumetric sinks/source of specific humidity and temp tr_qt(i,j,k) = ladzf(ntree_max-(ku-k))*qe/(rhoa*rlv) tr_qtR(i,j,k) = ladzf(ntree_max-(ku-k))*( omega*(s/(s+2*gam))*(qa(ntree_max-(ku-k))/(dzf(k)*ladzf(ntree_max-(ku-k)))))/(rhoa*rlv) tr_qtA(i,j,k) = ladzf(ntree_max-(ku-k))*((1-omega)*(1/(gam*(r_s)))*rhoa*cp*D)/(rhoa*rlv) tr_thl(i,j,k) = ladzf(ntree_max-(ku-k))*qh/(rhoa*cp) tr_omega(i,j,k) = omega ! add to rhs qtp(i,j,k) = qtp(i,j,k) + tr_qt(i,j,k) !ladz(ntree_max-(ku-k))*qe/(rhoa*rlv) thlp(i,j,k) = thlp(i,j,k) + tr_thl(i,j,k) !ladz(ntree_max-(ku-k))*qh/(rhoa*cp) ! fill dummy variables for steady-state BCs VT_dum2 = VT_dum2 - tr_thl(i,j,k)*dzf(k)*dxf(i)*dy Vq_dum2 = Vq_dum2 - tr_qt(i,j,k)*dzf(k)*dxf(i)*dy end do end do end do ! fraction of total net radiation reaching the bottom of the tree Rq = Rn(ntree_max+1-(ku+1-kl)) / Qstar ! sensible heat flux from shaded surfaces shade = ( (1 - 0.7 ) * Rn(ntree_max+1-(ku+1-kl) ) - 0.33*Rq*dQdt + Rq*38 )/ (rhoa*cp) ! dummy var for steady-state BC wTbot_dum2 = wTbot_dum2 - shade*(xh(iu+1)-xh(il))*dy*(ju-jl+1) ! overwrite standard heat flux defined in modibm in shaded positions thlp(il:iu,jl:ju,kb+1) = thlp(il:iu,jl:ju,kb+1) + (bctfz + shade) * dzfi(kb+1) end if end if !lmoist and ltempeq ! scalar deposition if (nsv>0) then ! specific to infinite canyon study (tg3315 thesis - Chapter 5) if (BCxs==2) then ! no deposition effect from first or last tree to avoid interaction with BCs if (n==1) cycle if (n==ntrees) cycle end if il = tree(n,1) - myidx*itot/nprocx iu = tree(n,2) - myidx*itot/nprocx kl = tree(n,5) ku = tree(n,6) jl = tree(n,3) - myidy*jtot/nprocy ju = tree(n,4) - myidy*jtot/nprocy if (iu < ib .or. il > ie .or. ju < jb .or. jl > je) then cycle else if (iu > ie) iu=ie if (il < ib) il=ib if (ju > je) ju=je if (jl < jb) jl=jb do m = 1,nsv ! define this for scalar variables that are deposited do k = kl,ku do j = jl,ju do i = il,iu tr_sv(i,j,k,m) = - svm(i,j,k,m) * ladzf(ntree_max-(ku-k)) * ud svp(i,j,k,m) = svp(i,j,k,m) + tr_sv(i,j,k,m) end do end do end do end do end if end if !nsv end do ! ntrees ! define these outside loop to avoid double counting cell faces wp(ib:ie,jb:je,kb:ke) = wp(ib:ie,jb:je,kb:ke) + tr_w(ib:ie,jb:je,kb:ke) vp(ib:ie,jb:je,kb:ke) = vp(ib:ie,jb:je,kb:ke) + tr_v(ib:ie,jb:je,kb:ke) up(ib:ie,jb:je,kb:ke) = up(ib:ie,jb:je,kb:ke) + tr_u(ib:ie,jb:je,kb:ke) ! steady-state BC calc. if (lmoist .and. ltempeq) then ! dummy vars - not all used? wTbot_dum3 = 0. wTbot_dum = 0. wTbot = 0. VT_dum3 = 0. Vq_dum3 = 0. VT_dum = 0. Vq_dum = 0. ! sum dummy variables assigned in main loop call MPI_ALLREDUCE(wTbot_dum2, wTbot_dum3,1,MY_REAL,MPI_SUM,comm3d,mpierr) call MPI_ALLREDUCE(VT_dum2, VT_dum3 ,1,MY_REAL,MPI_SUM,comm3d,mpierr) call MPI_ALLREDUCE(Vq_dum2, Vq_dum3 ,1,MY_REAL,MPI_SUM,comm3d,mpierr) ! not necessary wTbot_dum = wTbot_dum3 Vq_dum = Vq_dum3 VT_dum = VT_dum3 ! total flux of temp from ground and roofs accounting for shading wTbot = wTbot_dum + ((xh(ie+1)-xh(ib))*jtot*dy - tr_A)*bctfz ! temp flux at top equal to half of the total flux of heat (ground, roofs and canopies) - consistent with tg3315 thesis Chapter 4 - method iii wttop = 0.5*(wTbot + VT_dum)/((xh(ie+1)-xh(ib))*jtot*dy) ! volumetric sink term throughout domain to account for the remaining heating - "" thlpcar = wttop/((sum(real(IIcs(kb+1:ke))/rslabs)/real(kmax-1.))*(zh(ke+1)-zh(kb+1))) ! top flux of qt equal to total from canopies - therefore uniform flux profile wqtop = (Vq_dum ) / ((xh(ie+1)-xh(ib))*jtot*dy) end if end subroutine trees