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,& dx,dxi,dxiq,dx2,dy2,dyi,dyiq,dzf,dzf2,dzfi,dzhi,rk3step,rslabs, & numol,numoli,prandtlmoli,lles, rk3step,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)) *dxi )**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))*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 ) strain2 = strain2 + 0.125 * ( & + ((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 ) 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 ! 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_j / dx_i ip = i+1 im = i-1 a11 = (u0(ip,j,k) - u0(i,j,k)) * dxi a12 = (v0(ip,jp,k) + v0(ip,j,k) - v0(im,jp,k) - v0(im,j,k) )*dxiq a13 = (w0(ip,j,kp) + w0(ip,j,k) - w0(im,j,kp) - w0(im,j,k) )*dxiq a21 = (u0(ip,jp,k) + u0(i,jp,k) - u0(ip,jm,k) - u0(i,jm,k) )*dyiq 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 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 = dx2*a11*a11 + dy2*a21*a21 + dzf2(k)*a31*a31 b22 = dx2*a12*a12 + dy2*a22*a22 + dzf2(k)*a32*a32 b12 = dx2*a11*a12 + dy2*a21*a22 + dzf2(k)*a31*a32 b33 = dx2*a13*a13 + dy2*a23*a23 + dzf2(k)*a33*a33 b13 = dx2*a11*a13 + dy2*a21*a23 + dzf2(k)*a31*a33 b23 = dx2*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 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_j / dx_i ip = i+1 im = i-1 a11 = (u0(ip,j,k) - u0(i,j,k)) * dxi a12 = (v0(ip,jp,k) + v0(ip,j,k) - v0(im,jp,k) - v0(im,j,k) )*dxiq a13 = (w0(ip,j,kp) + w0(ip,j,k) - w0(im,j,kp) - w0(im,j,k) )*dxiq a21 = (u0(ip,jp,k) + u0(i,jp,k) - u0(ip,jm,k) - u0(i,jm,k) )*dyiq 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 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 = dx2*a11*a11 + dy2*a21*a21 + dzf2(k)*a31*a31 b22 = dx2*a12*a12 + dy2*a22*a22 + dzf2(k)*a32*a32 b12 = dx2*a11*a12 + dy2*a21*a22 + dzf2(k)*a31*a32 b33 = dx2*a13*a13 + dy2*a23*a23 + dzf2(k)*a33*a33 b13 = dx2*a11*a13 + dy2*a21*a23 + dzf2(k)*a31*a33 b23 = dx2*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 ! 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 damp(:,:,:) = max(damp(:,:,:),dampmin) ! ekm(:,:,:) = max(ekm(:,:,:),ekmin) ! ekh(:,:,:) = max(ekh(:,:,:),ekmin) else ! no subgrid model (DNS!) ekm = numol ekh = numol*prandtlmoli end if !************************************************************* ! Set boundary condition for K-closure factors. ! Also other BC's!! !************************************************************* call closurebc return end subroutine closure