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,du0dz,dv0dz,Rig 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 ! Vreman (2004) model for eddy viscosity; 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 ekm(i,j,k) = c_vreman*sqrt( max( (bb/aa), 0.0 ) ) ! Eddy viscosity end do end do end do ! ekm(:,:,:) = max(ekm(:,:,:),ekmin) ! Optional buoyancy correction on Vreman model for (non neutral) stable stratification ! Huusko et al., 2025. Large eddy simulation of canonical atmospheric boundary layer flows with the spectral element method in Nek5000. JAMES, 17(10), p.e2025MS005233 if ((lbuoyancy) .and. (lbuoycorr)) then do k = kb,ke kp = k+1 km = k-1 do j = jb,je jp = j+1 do i = ib,ie ip = i+1 !! No need to calculate dthvdz here as it anyway gets calculated in thermodynamics module in calthv subroutine. Avoid redundant computation here. !dthvdz(i,j,k) = ( thl0(i,j,kp) - thl0(i,j,km) )/(dzh(kp)+dzh(k)) du0dz = 0.5*( (u0(i,j,kp)+u0(ip,j,kp)) - (u0(i,j,km)+u0(ip,j,km)) )/(dzh(kp)+dzh(k)) dv0dz = 0.5*( (v0(i,j,kp)+v0(i,jp,kp)) - (v0(i,j,km)+v0(i,jp,km)) )/(dzh(kp)+dzh(k)) Rig = ( (grav/thl0(i,j,k)) * dthvdz(i,j,k) ) / (du0dz**2 + dv0dz**2 + 1.e-10) ! Eq. B4 ekm(i,j,k) = ekm(i,j,k) * sqrt( 1.0 - min(max(Rig,0.0),Rigc) / Rigc ) ! Eq. B3 ! Bouyancy correction applied only for stable stratification Rig > 0 ! For Rig >= Rigc, ekm = 0 end do end do end do end if ekh(:,:,:) = ekm(:,:,:)*prandtli ! Eddy diffusivity from eddy viscosity, Kh = Km/Pr ! ekh(:,:,:) = max(ekh(:,:,:),ekmin) 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