subroutine closure
!-----------------------------------------------------------------|
! |
!*** *closure* calculates K-coefficients |
! |
! Hans Cuijpers I.M.A.U. 06/01/1995 |
! |
! purpose. |
! -------- |
! |
! All the K-closure factors are calculated. |
! |
! ekm(i,j,k) = k sub m : for velocity-closure |
! ekh(i,j,k) = k sub h : for temperture-closure |
! ekh(i,j,k) = k sub h = k sub c : for concentration-closure |
! |
! We will use the next model for these factors: |
! |
! k sub m = 0.12 * l * sqrt(E) |
! |
! k sub h = k sub c = ( 1 + (2*l)/D ) * k sub m |
! |
! where : l = mixing length ( in model = z2 ) |
! E = subgrid energy |
! D = grid-size distance |
! |
!** interface. |
! ---------- |
! |
! *closure* is called from *program*. |
! |
!-----------------------------------------------------------------|
use modglobal, only : ib,ie,jb,je,kb,ke,kh,ih,jh,jmax,delta,ekmin,grav, zf, fkar,jgb,jge,&
dxf,dxf2,dxhi,dxfi,dy2,dyi,dyiq,dzf,dzf2,dzfi,dzhi,rk3step,rslabs, &
numol,numoli,prandtlmoli,lles, rk3step,dxfiq,dzfiq,lbuoyancy,dzh
use modfields, only : dthvdz,e120,u0,v0,w0,thl0,mindist,wall,shear
use modsurfdata, only : dudz,dvdz,thvs,ustar
use modmpi, only : excjs, myid, nprocs, comm3d, mpierr, my_real,mpi_sum,slabsumi
use modboundary, only : closurebc
use modinletdata, only : utaui
implicit none
real, dimension(ib:ie) :: shearbot
real :: strain2,mlen,uhor,distplus,utaubot,a11,a12,a13, &
a21,a22,a23,a31,a32,a33,aa,b11,b12,b13,b21,b22, &
b23,b33,bb,const,const2
integer :: i,j,k,kp,km,jp,jm,im,ip,iw,jw,kw,c1,c2
! if (lles .and. rk3step == 1) then ! compute ekm and ekh only once in complete RK-cycle
if(lsmagorinsky) then
do k = kb,ke
kp=k+1
km=k-1
do i = ib,ie
ip=i+1
im=i-1
mlen = csz(i,k) * delta(i,k)
do j = jb,je
jp=j+1
jm=j-1
iw = wall(i,j,k,1) ! indices of closest wall
jw = wall(i,j,k,2)-myid*jmax ! indices of closest wall in local j-index
kw = wall(i,j,k,3)
c1 = wall(i,j,k,4) ! shear stress component
c2 = wall(i,j,k,5) ! shear stress component
if ((jw >= jb-1) .and. (jw <= je+1)) then ! check if jw is within the halo of this proc
!write(*,'(A,E9.2,A,E9.2,A,E9.2,A,E9.2)') 'component1:', c1, 'component2:', c2, 'shear c1:', shear(iw,jw,kw,c1), 'shear c2:', shear(iw,jw,kw,c2)
distplus = mindist(i,j,k)*sqrt(abs(shear(iw,jw,kw,c1))+abs(shear(iw,jw,kw,c2)))*numoli
damp(i,j,k) = sqrt(1. - exp((-distplus*0.04)**3.)) ! Wall-damping according to Piomelli
! write(*,'(A,2(1pE9.2))') 'damp, distplus', damp(i,j,k), distplus
else
damp(i,j,k) = 1.
end if
strain2 = ( &
((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 )
strain2 = strain2 + 0.125 * ( &
((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 )
strain2 = strain2 + 0.125 * ( &
((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 )
strain2 = strain2 + 0.125 * ( &
((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 )
ekm(i,j,k) = (mlen*damp(i,j,k)) ** 2. * sqrt(2. * strain2)
ekh(i,j,k) = ekm(i,j,k) * prandtli
end do
end do
end do
damp(:,:,:) = max(damp(:,:,:),dampmin)
ekm(:,:,:) = ekm(:,:,:) + numol ! add molecular viscosity
ekh(:,:,:) = ekh(:,:,:) + numol*prandtlmoli ! add molecular diffusivity
!write(*,'(A,3(1pE9.2))') 'strain2, ekm, ekh', strain2, ekm(10,10,10), ekh(10,10,10)
! ekm(:,:,:) = max(ekm(:,:,:),ekmin)
! ekh(:,:,:) = max(ekh(:,:,:),ekmin)
elseif(lvreman) then
if ((lbuoyancy) .and. (lbuoycorr)) then
const = prandtli*grav/(thvs*sqrt(2.*3.))
do k = kb,ke
kp = k+1
km = k-1
do j = jb,je
jp = j+1
jm = j-1
do i = ib,ie ! aij = du_i / dx_i
ip = i+1
im = i-1
a11 = (u0(ip,j,k) - u0(i,j,k)) * dxfi(i)
a12 = (((v0(ip,jp,k)+ v0(ip,j,k))*dxf(i) + &
(v0(i,jp,k) + v0(i,j,k)) *dxf(ip)) * dxhi(ip) &
- &
((v0(i,jp,k) + v0(i,j,k)) *dxf(im) + &
(v0(im,jp,k)+ v0(im,j,k))*dxf(i)) * dxhi(i)) * dxfiq(i)
a13 = (((w0(ip,j,kp)+ w0(ip,j,k))*dxf(i) + &
(w0(i,j,kp) + w0(i,j,k)) *dxf(ip)) * dxhi(ip) &
- &
((w0(i,j,kp) + w0(i,j,k)) *dxf(im) + &
(w0(im,j,kp)+ w0(im,j,k))*dxf(i)) * dxhi(i)) * dxfiq(i)
a21 = (u0(ip,jp,k) + u0(i,jp,k) - u0(ip,jm,k) - u0(i,jm,k) )*dyiq ! simplified after writing out interpolation.
a22 = (v0(i,jp,k) - v0(i,j,k)) * dyi
a23 = (w0(i,jp,kp) + w0(i,jp,k) - w0(i,jm,kp) - w0(i,jm,k) )*dyiq ! simplified after writing out interpolation.
a31 = (((u0(ip,j,kp) + u0(i,j,kp))*dzf(k) + &
(u0(ip,j,k) + u0(i,j,k)) *dzf(kp)) * dzhi(kp) &
- &
((u0(ip,j,k) + u0(i,j,k)) *dzf(km) + &
(u0(ip,j,km) + u0(i,j,km))*dzf(k)) * dzhi(k) )*dzfiq(k)
a32 = (((v0(i,jp,kp) + v0(i,j,kp))*dzf(k) + &
(v0(i,jp,k) + v0(i,j,k)) *dzf(kp)) * dzhi(kp) &
- &
((v0(i,jp,k) + v0(i,j,k)) *dzf(km) + &
(v0(i,jp,km) + v0(i,j,km))*dzf(k)) * dzhi(k) )*dzfiq(k)
a33 = (w0(i,j,kp) - w0(i,j,k)) * dzfi(k)
aa = a11*a11 + a21*a21 + a31*a31 + &
a12*a12 + a22*a22 + a32*a32 + &
a13*a13 + a23*a23 + a33*a33
b11 = dxf2(i)*a11*a11 + dy2*a21*a21 + dzf2(k)*a31*a31
b22 = dxf2(i)*a12*a12 + dy2*a22*a22 + dzf2(k)*a32*a32
b12 = dxf2(i)*a11*a12 + dy2*a21*a22 + dzf2(k)*a31*a32
b33 = dxf2(i)*a13*a13 + dy2*a23*a23 + dzf2(k)*a33*a33
b13 = dxf2(i)*a11*a13 + dy2*a21*a23 + dzf2(k)*a31*a33
b23 = dxf2(i)*a12*a13 + dy2*a22*a23 + dzf2(k)*a32*a33
bb = b11*b22 - b12*b12 + b11*b33 - b13*b13 + b22*b33 - b23*b23
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) then
const2=(bb/aa)
else
! write(*,*) "const",const
! write(*,*) "delta",delta
const2=(bb/aa)-(delta(i,k)**4)*dthvdz(i,j,k)*const
if (const2 <0.0) const2 = 0.0
end if
ekm(i,j,k)=c_vreman*sqrt(const2)
ekh(i,j,k)=ekm(i,j,k)*prandtli
end do
end do
end do
ekm(:,:,:) = ekm(:,:,:) + numol ! add molecular viscosity
ekh(:,:,:) = ekh(:,:,:) + numol*prandtlmoli ! add molecular diffusivity
else ! neutral case
do k = kb,ke
kp = k+1
km = k-1
do j = jb,je
jp = j+1
jm = j-1
do i = ib,ie ! aij = du_i / dx_i
ip = i+1
im = i-1
a11 = (u0(ip,j,k) - u0(i,j,k)) * dxfi(i)
a12 = (((v0(ip,jp,k)+ v0(ip,j,k))*dxf(i) + &
(v0(i,jp,k) + v0(i,j,k)) *dxf(ip)) * dxhi(ip) &
- &
((v0(i,jp,k) + v0(i,j,k)) *dxf(im) + &
(v0(im,jp,k)+ v0(im,j,k))*dxf(i)) * dxhi(i)) * dxfiq(i)
a13 = (((w0(ip,j,kp)+ w0(ip,j,k))*dxf(i) + &
(w0(i,j,kp) + w0(i,j,k)) *dxf(ip)) * dxhi(ip) &
- &
((w0(i,j,kp) + w0(i,j,k)) *dxf(im) + &
(w0(im,j,kp)+ w0(im,j,k))*dxf(i)) * dxhi(i)) * dxfiq(i)
a21 = (u0(ip,jp,k) + u0(i,jp,k) - u0(ip,jm,k) - u0(i,jm,k) )*dyiq !simplified after writing out interpolation.
a22 = (v0(i,jp,k) - v0(i,j,k)) * dyi
a23 = (w0(i,jp,kp) + w0(i,jp,k) - w0(i,jm,kp) - w0(i,jm,k) )*dyiq !simplified after writing out interpolation.
a31 = (((u0(ip,j,kp) + u0(i,j,kp))*dzf(k) + &
(u0(ip,j,k) + u0(i,j,k)) *dzf(kp)) * dzhi(kp) &
- &
((u0(ip,j,k) + u0(i,j,k)) *dzf(km) + &
(u0(ip,j,km) + u0(i,j,km))*dzf(k)) * dzhi(k) )*dzfiq(k)
a32 = (((v0(i,jp,kp) + v0(i,j,kp))*dzf(k) + &
(v0(i,jp,k) + v0(i,j,k)) *dzf(kp)) * dzhi(kp) &
- &
((v0(i,jp,k) + v0(i,j,k)) *dzf(km) + &
(v0(i,jp,km) + v0(i,j,km))*dzf(k)) * dzhi(k) )*dzfiq(k)
a33 = (w0(i,j,kp) - w0(i,j,k)) * dzfi(k)
aa = a11*a11 + a21*a21 + a31*a31 + &
a12*a12 + a22*a22 + a32*a32 + &
a13*a13 + a23*a23 + a33*a33
b11 = dxf2(i)*a11*a11 + dy2*a21*a21 + dzf2(k)*a31*a31
b22 = dxf2(i)*a12*a12 + dy2*a22*a22 + dzf2(k)*a32*a32
b12 = dxf2(i)*a11*a12 + dy2*a21*a22 + dzf2(k)*a31*a32
b33 = dxf2(i)*a13*a13 + dy2*a23*a23 + dzf2(k)*a33*a33
b13 = dxf2(i)*a11*a13 + dy2*a21*a23 + dzf2(k)*a31*a33
b23 = dxf2(i)*a12*a13 + dy2*a22*a23 + dzf2(k)*a32*a33
bb = b11*b22 - b12*b12 + b11*b33 - b13*b13 + b22*b33 - b23*b23
if (bb < 0.00000001) then
ekm(i,j,k) = 0.
ekh(i,j,k) = 0.
else
ekm(i,j,k) = c_vreman*sqrt(bb / aa)
ekh(i,j,k) = ekm(i,j,k)*prandtli
end if
end do
end do
end do
! ekm(:,:,:) = max(ekm(:,:,:),ekmin)
! ekh(:,:,:) = max(ekh(:,:,:),ekmin)
end if ! lbuoyancy
ekm(:,:,:) = ekm(:,:,:) + numol ! add molecular viscosity
ekh(:,:,:) = ekh(:,:,:) + numol*prandtlmoli ! add molecular diffusivity
!do TKE scheme
elseif (loneeqn ) then
do k=kb,ke
do j=jb,je
do i=ib,ie
iw = wall(i,j,k,1) ! indices of closest wall
jw = wall(i,j,k,2)-myid*jmax ! indices of closest wall in local j-index
kw = wall(i,j,k,3)
c1 = wall(i,j,k,4) ! shear stress component
c2 = wall(i,j,k,5) ! shear stress component
!ILS13 removed near-wall damping 25.06.2014
!if (jw >= jb-1 .and. jw <= je+1) then ! check if jw is within the halo of this proc
! distplus = mindist(i,j,k)*sqrt(abs(shear(iw,jw,kw,c1))+abs(shear(iw,jw,kw,c2)))*numoli
! damp(i,j,k) = sqrt(1. - exp((-distplus*0.04)**3.)) ! Wall-damping according to Piomelli
!else
damp(i,j,k) = 1.
!end if
if ((ldelta) .or. (dthvdz(i,j,k)<=0)) then
zlt(i,j,k) = delta(i,k)
ekm(i,j,k) = cm * zlt(i,j,k) *damp(i,j,k)* e120(i,j,k) !* 0.5! LES with near-wall damping !!! added factor 0.5 for shear-driven flow
ekh(i,j,k) = (ch1 + ch2) * ekm(i,j,k) ! maybe ekh should be calculated from (molecular) Prandtl number
ekm(i,j,k) = ekm(i,j,k) + numol ! add molecular viscosity
ekh(i,j,k) = ekh(i,j,k) + numol*prandtlmoli ! add molecular diffusivity
else
! zlt(i,j,k) = min(delta(i,k),cn*e120(i,j,k)/sqrt(grav/thvs*abs(dthvdz(i,j,k))))
zlt(i,j,k) = min(delta(i,k),cn*e120(i,j,k)/sqrt(grav/thvs*abs(dthvdz(i,j,k)))) !thls is used
ekm(i,j,k) = cm * zlt(i,j,k) *damp(i,j,k)* e120(i,j,k) !* 0.5 ! LES with near-wall damping !!! added factor 0.5 for shear-driven flow
ekh(i,j,k) = (ch1 + ch2 * zlt(i,j,k)/delta(i,k)) * ekm(i,j,k) ! needed in LES!
ekm(i,j,k) = ekm(i,j,k) + numol ! add molecular viscosity
ekh(i,j,k) = ekh(i,j,k) + numol*prandtlmoli ! add molecular diffusivity
endif
end do
end do
end do
! write(*,'(1A,F8.4)') 'ekh', ekh(2,3,3)
damp(:,:,:) = max(damp(:,:,:),dampmin)
! ekm(:,:,:) = max(ekm(:,:,:),ekmin)
! ekh(:,:,:) = max(ekh(:,:,:),ekmin)
else ! no subgrid model (DNS!)
ekm = numol
ekh = numol*prandtlmoli
end if
! write(*,'(A,3(1pE9.2))') 'strain2, ekm, ekh', strain2, ekm(10,10,10), ekh(10,10,10)
!do i=ib,ie
!write(*,'(A,1(1pE9.2))') 'damp, ekm', damp(i,32,kb), ekm(i,32,kb)
!end do
!*************************************************************
! Set boundary condition for K-closure factors. ! Also other BC's!!
!*************************************************************
call closurebc
return
end subroutine closure