subroutine sources ! only in case of LES computation !-----------------------------------------------------------------| ! | !*** *sources* | ! calculates various terms from the subgrid TKE equation | ! | ! Hans Cuijpers I.M.A.U. 06/01/1995 | ! | ! purpose. | ! -------- | ! | ! Subroutine sources calculates all other terms in the | ! subgrid energy equation, except for the diffusion terms. | ! These terms are calculated in subroutine diff. | ! | !** interface. | ! ---------- | ! | ! *sources* is called from *program*. | ! | !-----------------------------------------------------------------| use modglobal, only : ib,ie,jb,je,kb,ke,dxi,delta,dy,dyi,dzfi,dzhi,grav,numol,prandtlmol,& dzh, delta use modfields, only : u0,v0,w0,e120,e12p,dthvdz,thl0,thvf use modsurfdata, only : dudz,dvdz,thvs ! use modmpi, only : myid implicit none real tdef2,prandtlmoli integer i,j,k,im,ip,jm,jp,km,kp prandtlmoli = 1./prandtlmol do k=kb+1,ke do j=jb,je do i=ib,ie kp=k+1 km=k-1 jp=j+1 jm=j-1 ip=i+1 im=i-1 tdef2 = 2. * ( & ((u0(ip,j,k)-u0(i,j,k)) * dxi )**2 & + ((v0(i,jp,k)-v0(i,j,k)) * dyi )**2 & + ((w0(i,j,kp)-w0(i,j,k)) * dzfi(k) )**2 ) tdef2 = tdef2 + 0.25 * ( & ((w0(i ,j,kp)-w0(im,j,kp)) * dxi + (u0(i ,j,kp)-u0(i ,j,k )) * dzhi(kp))**2 & + ((w0(i ,j,k )-w0(im,j,k )) * dxi + (u0(i ,j,k )-u0(i ,j,km)) * dzhi(k ))**2 & + ((w0(ip,j,k )-w0(i ,j,k )) * dxi + (u0(ip,j,k )-u0(ip,j,km)) * dzhi(k ))**2 & + ((w0(ip,j,kp)-w0(i,j,kp)) * dxi + (u0(ip,j,kp)-u0(ip,j,k )) * dzhi(kp))**2 ) tdef2 = tdef2 + 0.25 * ( & ((u0(i ,jp,k)-u0(i ,j ,k)) * dyi + (v0(i ,jp,k)-v0(im,jp,k)) * dxi)**2 & + ((u0(i ,j ,k)-u0(i ,jm,k)) * dyi + (v0(i ,j ,k)-v0(im,j ,k)) * dxi)**2 & + ((u0(ip,j ,k)-u0(ip,jm,k)) * dyi + (v0(ip,j ,k)-v0(i ,j ,k)) * dxi)**2 & + ((u0(ip,jp,k)-u0(ip,j ,k)) * dyi + (v0(ip,jp,k)-v0(i ,jp,k)) * dxi)**2 ) tdef2 = tdef2 + 0.25 * ( & ((v0(i,j ,kp)-v0(i,j ,k )) * dzhi(kp) + (w0(i,j ,kp)-w0(i,jm,kp)) * dyi)**2 & + ((v0(i,j ,k )-v0(i,j ,km)) * dzhi(k ) + (w0(i,j ,k )-w0(i,jm,k )) * dyi)**2 & + ((v0(i,jp,k )-v0(i,jp,km)) * dzhi(k ) + (w0(i,jp,k )-w0(i,j ,k )) * dyi)**2 & + ((v0(i,jp,kp)-v0(i,jp,k )) * dzhi(kp) + (w0(i,jp,kp)-w0(i,j ,kp)) * dyi)**2 ) ! sbshr(i,j,k) = ekm(i,j,k)*tdef2/ ( 2*e120(i,j,k)) ! sbbuo(i,j,k) = -ekh(i,j,k)*grav/thvs*dthvdz(i,j,k)/ ( 2*e120(i,j,k)) ! sbdiss(i,j,k) = - (ce1 + ce2*zlt(i,j,k)/delta(i,k)) * e120(i,j,k)**2 /(2.*zlt(i,j,k)) sbshr(i,j,k) = (ekm(i,j,k)-numol)*tdef2/ ( 2*e120(i,j,k)) ! subtract molecular viscosity ! sbbuo(i,j,k) = -(ekh(i,j,k)-numol*prandtlmoli)*grav/thvs*dthvdz(i,j,k)/ ( 2*e120(i,j,k)) ! subtract molecular diffusivity sbbuo(i,j,k) = -(ekh(i,j,k)-numol*prandtlmoli)*grav/thvs*dthvdz(i,j,k)/ ( 2*e120(i,j,k)) ! subtract molecular diffusivity and use thls instead of thvs (not defined) ! sbdiss(i,j,k) = - (ce1 + ce2*zlt(i,j,k)/delta(i,k)) * e120(i,j,k)**2 /(2.*damp*zlt(i,j,k)) ! add near-wall damping function ! added factor 2. for shear-driven flow sbdiss(i,j,k) = - 2. * (ce1 + ce2*zlt(i,j,k)/delta(i,k)) * e120(i,j,k)**2 /(2.*damp(i,j,k)*zlt(i,j,k)) ! add near-wall damping function !! added f end do end do end do ! ----------------------------------------------end i,j,k-loop ! special treatment for lowest level ! Don't do this - wall function at bottom ! do j=jb,je ! do i=ib,ie ! jp=j+1 ! jm=j-1 ! ip=i+1 ! im=i-1 ! ! tdef2 = 2. * ( & ! ((u0(ip,j,kb) - u0(i,j,kb))*dxi)**2 & ! + ((v0(i,jp,kb) - v0(i,j,kb))*dyi)**2 & ! + ((w0(i,j,kb+1) -w0(i,j,kb))*dzfi(kb))**2 & ! ) ! ! tdef2 = tdef2 + ( 0.25*(w0(i+1,j,kb+1)-w0(i-1,j,kb+1))*dxfi(i) + dudz(i,j))**2 ! ! tdef2 = tdef2 + 0.25 * ( & ! ((u0(i ,jp,kb) - u0(i ,j ,kb)) * dyi + (v0(i ,jp,kb) - v0(im,jp,kb)) * dxi)**2 & ! + ((u0(i ,j ,kb) - u0(i ,jm,kb)) * dyi + (v0(i ,j ,kb) - v0(im,j ,kb)) * dxi)**2 & ! + ((u0(ip,j ,kb) - u0(ip,jm,kb)) * dyi + (v0(ip,j ,kb) - v0(i ,j ,kb)) * dxi)**2 & ! + ((u0(ip,jp,kb) - u0(ip,j ,kb)) * dyi + (v0(ip,jp,kb) - v0(i ,jp,kb)) * dxi)**2 & ! ) ! ! tdef2 = tdef2 + ( 0.25 * (w0(i,jp,kb+1) - w0(i,jm,kb+1)) * dyi + dvdz(i,j))**2 ! ! ! ** Include shear and buoyancy production terms and dissipation ** ! ! sbshr(i,j,kb) = ekm(i,j,kb)*tdef2/ ( 2*e120(i,j,kb)) ! sbbuo(i,j,kb) = -ekh(i,j,kb)*grav/thvf(kb)*dthvdz(i,j,kb)/ ( 2*e120(i,j,kb)) ! sbdiss(i,j,kb) = - (ce1 + ce2*zlt(i,j,kb)/delta(i,kb)) * e120(i,j,kb)**2 /(2.*zlt(i,j,kb)) ! end do ! end do ! ------------------------------------------------ e12p(ib:ie,jb:je,kb:ke) = e12p(ib:ie,jb:je,kb:ke)+ & sbshr(ib:ie,jb:je,kb:ke)+sbbuo(ib:ie,jb:je,kb:ke)+sbdiss(ib:ie,jb:je,kb:ke) return end subroutine sources