subroutine diffu (putout)
use modglobal, only : ib,ie,ih,jb,je,jh,kb,ke,kh,kmax,dxhi,dxf,dxfi,lles,&
dzf,dzfi,dy,dyi,dy2i,dzhi,dzhiq,jmax,numol
use modfields, only : u0,v0,w0
use modsurfdata,only : ustar
use modmpi, only : myid
implicit none
real, intent(inout) :: putout(ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh)
real :: emmo,emom,emop,empo
real :: fu,dummy
real :: ucu, upcu
integer :: i,j,k,jm,jp,km,kp
if (lles) then
do k=kb,ke
kp=k+1
km=k-1
do j=jb,je
jp=j+1
jm=j-1
do i=ib,ie
emom = ( dzf(km) * ( ekm(i,j,k)*dxf(i-1) + ekm(i-1,j,k)*dxf(i) ) + & ! dx is non-equidistant
dzf(k) * ( ekm(i,j,km)*dxf(i-1) + ekm(i-1,j,km)*dxf(i) ) )*dxhi(i) * dzhiq(k)
emop = ( dzf(kp) * ( ekm(i,j,k)*dxf(i-1) + ekm(i-1,j,k)*dxf(i) ) + & ! dx is non-equidistant
dzf(k) * ( ekm(i,j,kp)*dxf(i-1) + ekm(i-1,j,kp)*dxf(i) ) )*dxhi(i) * dzhiq(kp)
empo = 0.25 * ( ( ekm(i,j,k)+ekm(i,jp,k))*dxf(i-1) + (ekm(i-1,j,k)+ekm(i-1,jp,k))*dxf(i) ) * dxhi(i) ! dx is non-equidistant
emmo = 0.25 * ( ( ekm(i,j,k)+ekm(i,jm,k))*dxf(i-1) +(ekm(i-1,jm,k)+ekm(i-1,j,k))*dxf(i) ) * dxhi(i) ! dx is non-equidistant
! Discretized diffusion term
putout(i,j,k) = putout(i,j,k) &
+ &
( ekm(i,j,k) * (u0(i+1,j,k)-u0(i,j,k))*dxfi(i) &
-ekm(i-1,j,k)* (u0(i,j,k)-u0(i-1,j,k))*dxfi(i-1) ) * 2. * dxhi(i) &
+ &
( empo * ( (u0(i,jp,k)-u0(i,j,k)) *dyi &
+(v0(i,jp,k)-v0(i-1,jp,k))*dxhi(i)) &
-emmo * ( (u0(i,j,k)-u0(i,jm,k)) *dyi &
+(v0(i,j,k)-v0(i-1,j,k)) *dxhi(i)) &
) * dyi &
+ &
( emop * ( (u0(i,j,kp)-u0(i,j,k)) *dzhi(kp) &
+(w0(i,j,kp)-w0(i-1,j,kp))*dxhi(i)) &
-emom * ( (u0(i,j,k)-u0(i,j,km)) *dzhi(k) &
+(w0(i,j,k)-w0(i-1,j,k)) *dxhi(i)) &
) *dzfi(k)
end do
end do
end do
else ! DNS
do k=kb,ke
kp=k+1
km=k-1
do j=jb,je
jp=j+1
jm=j-1
do i=ib,ie
! Discretized diffusion term
putout(i,j,k) = putout(i,j,k) &
+ &
( numol * (u0(i+1,j,k)-u0(i,j,k))*dxfi(i) &
-numol * (u0(i,j,k)-u0(i-1,j,k))*dxfi(i-1) ) * 2. * dxhi(i) &
+ &
( numol * ( (u0(i,jp,k)-u0(i,j,k)) *dyi &
+(v0(i,jp,k)-v0(i-1,jp,k))*dxhi(i)) &
-numol * ( (u0(i,j,k)-u0(i,jm,k)) *dyi &
+(v0(i,j,k)-v0(i-1,j,k)) *dxhi(i)) &
) * dyi &
+ &
( numol * ( (u0(i,j,kp)-u0(i,j,k)) *dzhi(kp) &
+(w0(i,j,kp)-w0(i-1,j,kp))*dxhi(i)) &
-numol * ( (u0(i,j,k)-u0(i,j,km)) *dzhi(k) &
+(w0(i,j,k)-w0(i-1,j,k)) *dxhi(i)) &
) *dzfi(k)
end do
end do
end do
end if ! lles
end subroutine diffu