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,delta,dxhi,dxfi,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
! Added by J. Tomas (thermodynamics routine is bypassed)
! thv does not exist (equals thl), therefore dthvdz is computed here
! IS+HJ: thermodynamics is back in now (which calculates dthvdz)
! do k=kb,ke
! do j=jb,je
! do i=ib,ie
! dthvdz(i,j,k) = (thl0(i,j,k+1)-thl0(i,j,k-1))/(dzh(k+1)+dzh(k))
! ! if (dthvdz(i,j,k) == 0.0) then
! ! write(6,*) 'dthvdz(i,j,k)=0,i,j,k',i,j,k
! ! end if
! end do
! end do
! end do
! End of addition by J. Tomas (thermodynamics routine is bypassed)
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)) * dxfi(i) )**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)) * dxhi(i) + &
(u0(i,j,kp)-u0(i,j,k)) * dzhi(kp) )**2 + &
!
((w0(i,j,k)-w0(im,j,k)) * dxhi(i) + &
(u0(i,j,k)-u0(i,j,km)) * dzhi(k) )**2 + &
!
((w0(ip,j,k)-w0(i,j,k)) * dxhi(ip) + &
(u0(ip,j,k)-u0(ip,j,km)) * dzhi(k) )**2 + &
!
((w0(ip,j,kp)-w0(i,j,kp)) * dxhi(ip) + &
(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)) * dxhi(i) )**2 + &
!
((u0(i,j,k)-u0(i,jm,k)) * dyi + &
(v0(i,j,k)-v0(im,j,k)) * dxhi(i) )**2 + &
!
((u0(ip,j,k)-u0(ip,jm,k)) * dyi + &
(v0(ip,j,k)-v0(i,j,k)) * dxhi(ip))**2 + &
!
((u0(ip,jp,k)-u0(ip,j,k)) * dyi + &
(v0(ip,jp,k)-v0(i,jp,k)) * dxhi(ip))**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
do j=jb,je
do i=ib,ie
jp=j+1
jm=j-1
tdef2= 2. * ( &
((u0(i+1,j,kb)-u0(i,j,kb))*dxfi(i))**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(i-1,jp,kb))*dxfi(i))**2 &
+((u0(i,j,kb)-u0(i,jm,kb))*dyi+(v0(i,j,kb)-v0(i-1,j,kb))*dxfi(i))**2 &
+((u0(i+1,j,kb)-u0(i+1,jm,kb))*dyi+(v0(i+1,j,kb)-v0(i,j,kb))*dxfi(i))**2 &
+((u0(i+1,jp,kb)-u0(i+1,j,kb))*dyi+(v0(i+1,jp,kb)-v0(i,jp,kb))*dxfi(i))**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)
!write(*,'(A,3(1pE9.2))') 'for check if called and e12p =', e12p(32,32,kb)
return
end subroutine sources