Index: LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90 (revision 4643) +++ LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90 (revision 4644) @@ -6,5 +6,5 @@ subroutine atke_compute_km_kh(ngrid,nlay,dtime, & - wind_u,wind_v,temp,play,pinterf, & + wind_u,wind_v,temp,play,pinterf,cdrag_uv, & tke,Km_out,Kh_out) @@ -26,7 +26,7 @@ -USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd -USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd -USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin +USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, ok_vdiff_tke +USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd +USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin, smmin implicit none @@ -45,5 +45,5 @@ REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) - +REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag_uv ! surface drag coefficient for momentum REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers @@ -78,4 +78,7 @@ REAL :: lstrat ! mixing length depending on local stratification REAL :: taustrat ! caracteristic timescale for turbulence in very stable conditions +REAL :: netloss ! net loss term of tke +REAL :: netsource ! net source term of tke +REAL :: ustar ! friction velocity estimation ! Initializations: @@ -142,5 +145,5 @@ Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut ELSE ! stable cases - Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric)) + Sm(igrid,ilay) = max(smmin,cn*(1.-Ri(igrid,ilay)/Ric)) Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN @@ -182,5 +185,5 @@ l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) IF (N2(igrid,ilay) .GT. 0.) THEN - lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1E-10)*(1.+sqrt(Ri(igrid,ilay))/2.)) + lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1.E-10)*(1.+sqrt(Ri(igrid,ilay))/2.)) IF (lstrat .LT. l_exchange(igrid,ilay)) THEN l_exchange(igrid,ilay)=max(lstrat,lmin) @@ -204,7 +207,9 @@ -! Computing the TKE: -!=================== +! Computing the TKE k>=2: +!======================== IF (iflag_atke == 0) THEN + +! stationary solution (dtke/dt=0) DO ilay=2,nlay @@ -221,5 +226,5 @@ DO ilay=2,nlay DO igrid=1,ngrid - qq=sqrt(2.*tke(igrid,ilay)) + qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * & (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) & @@ -230,38 +235,52 @@ ENDDO + ELSE IF (iflag_atke == 2) THEN -! full implicit scheme resolved with a second order polynomial equation -! + exact resolution for very stable cases (iflag_atke_lmix=1) +! semi implicit scheme when l does not depend on tke +! positive-guaranteed if pr slope in stable condition >1 + + DO ilay=2,nlay + DO igrid=1,ngrid + qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) + qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.) & + *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & + /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) + tke(igrid,ilay)=0.5*(qq**2) + ENDDO + ENDDO + + +ELSE IF (iflag_atke == 3) THEN +! numerical resolution adapted from that in MAR (Deleersnijder 1992) +! positively defined by construction + DO ilay=2,nlay DO igrid=1,ngrid - qq=sqrt(2.*tke(igrid,ilay)) - IF (switch_num(igrid,ilay)) THEN - taustrat=clmix/2./sqrt(N2(igrid,ilay))*Sm(igrid,ilay)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) & - - sqrt(N2(igrid,ilay))/2./cepsilon/clmix - taustrat=-1./taustrat - qq=qq*exp(-dtime/taustrat) + qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) + IF (Ri(igrid,ilay) .LT. 0.) THEN + netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)) + netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) ELSE - delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * & - (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay) * shear2(igrid,ilay) & - *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) - qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay) + netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ & + l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay) + netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay) ENDIF + qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss) tke(igrid,ilay)=0.5*(qq**2) ENDDO ENDDO - - -ELSE IF (iflag_atke == 3) THEN - -! semi implicit scheme when l does not depend on tke -! positive-guaranteed if pr slope in stable condition >1 +ELSE IF (iflag_atke == 4) THEN +! semi implicit scheme from Arpege (V. Masson methodology with +! Taylor expansion of the dissipation term) DO ilay=2,nlay DO igrid=1,ngrid - qq=sqrt(2.*tke(igrid,ilay)) - qq=(qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & - /(1.+dtime*qq/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2))) - tke(igrid,ilay)=0.5*(qq**2) + qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) + qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) & + +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) & + /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) + qq=max(0.,qq) + tke(igrid,ilay)=0.5*(qq**2) ENDDO ENDDO @@ -273,4 +292,27 @@ END IF + +! We impose a 0 tke at nlay+1 +!============================== + +DO igrid=1,ngrid + tke(igrid,nlay+1)=0. +END DO + + +! Calculation of surface TKE (k=1) +!================================= +! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note) +DO igrid=1,ngrid + ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2)) + tke(igrid,1)=ctkes*(ustar**2) +END DO + + +! vertical diffusion of TKE +!========================== +IF (ok_vdiff_tke) THEN + CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke) +ENDIF @@ -292,4 +334,103 @@ end subroutine atke_compute_km_kh +!=============================================================================================== +subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke) + +! routine that computes the vertical diffusion of TKE by the turbulence +! using an implicit resolution (See note by Dufresne and Ghattas +! E Vignon, July 2023 + +USE atke_turbulence_ini_mod, ONLY : rd, cke, viscom + + +INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) +INTEGER, INTENT(IN) :: nlay ! number of vertical index + +REAL, INTENT(IN) :: dtime ! physics time step (s) +REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: z_lay ! altitude of mid-layers (m) +REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: z_interf ! altitude of bottom interfaces (m) +REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) +REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) +REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: l_exchange ! mixing length at interfaces between layers +REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: Sm ! stability function for eddy diffusivity for momentum at interface between layers + +REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers + + + +INTEGER :: igrid,ilay +REAL, DIMENSION(ngrid,nlay+1) :: Ke ! eddy diffusivity for TKE +REAL, DIMENSION(ngrid,nlay+1) :: dtke +REAL, DIMENSION(ngrid,nlay+1) :: ak, bk, ck, CCK, DDK +REAL :: gammak,Kem,KKb,KKt + + +! Few initialisations +CCK(:,:)=0. +DDK(:,:)=0. +dtke(:,:)=0. + + +! Eddy diffusivity for TKE + +DO ilay=2,nlay + DO igrid=1,ngrid + Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke + ENDDO +ENDDO +! at the top of the atmosphere set to 0 +Ke(:,nlay+1)=0. +! at the surface, set it equal to that at the first model level +Ke(:,1)=Ke(:,2) + + +! calculate intermediary variables + +DO ilay=2,nlay + DO igrid=1,ngrid + Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay)) + KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay)) + Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1)) + KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1)) + gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1)) + ak(igrid,ilay)=-gammak*dtime*KKb + ck(igrid,ilay)=-gammak*dtime*KKt + bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb) + ENDDO +ENDDO + +! calculate CCK and DDK coefficients +! downhill phase + +DO igrid=1,ngrid + CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay) + DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay) +ENDDO + + +DO ilay=nlay-1,2,-1 + DO igrid=1,ngrid + CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) & + / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) + DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) + ENDDO +ENDDO + +! calculate TKE +! uphill phase + +DO ilay=2,nlay+1 + DO igrid=1,ngrid + dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay) + ENDDO +ENDDO + +! update TKE +tke(:,:)=tke(:,:)+dtke(:,:) + + +end subroutine atke_vdiff_tke + + end module atke_exchange_coeff_mod Index: LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90 (revision 4643) +++ LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90 (revision 4644) @@ -9,6 +9,6 @@ real :: kappa = 0.4 ! Von Karman constant !$OMP THREADPRIVATE(kappa) - real :: l0, ric, ri0, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix - !$OMP THREADPRIVATE(l0, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix) + real :: l0, ric, ri0, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix, smmin, ctkes,cke + !$OMP THREADPRIVATE(l0,ric,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,smmin,ctkes,cke) integer :: lunout,prt_level !$OMP THREADPRIVATE(lunout,prt_level) @@ -19,7 +19,9 @@ !$OMP THREADPRIVATE(viscom,viscoh) - real :: lmin=0.001 ! minimum mixing length + real :: lmin=0.01 ! minimum mixing length !$OMP THREADPRIVATE(lmin) + logical :: ok_vdiff_tke + !$OMP THREADPRIVATE(ok_vdiff_tke) CONTAINS @@ -56,4 +58,7 @@ endif + ! activate vertical diffusion of TKE or not + ok_vdiff_tke=.false. + CALL getin_p('atke_ok_vdiff_tke',ok_vdiff_tke) ! flag that controls the numerical treatment of diffusion coeffiient calculation @@ -99,5 +104,20 @@ clmix=0.5 CALL getin_p('atke_clmix',clmix) - + + ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri. + ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17 + + smmin=0.17 + CALL getin_p('atke_smmin',smmin) + + ! coefficient for surface TKE (default value from Arpege, see E. Bazile note) + ctkes=3.75 + CALL getin_p('atke_ctkes',ctkes) + + ! ratio between the eddy diffusivity coeff for tke wrt that for momentum + ! default value from Lenderink et al. 2004 + cke=2. + CALL getin_p('atke_cke',cke) + RETURN Index: LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 (revision 4643) +++ LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 (revision 4644) @@ -55,9 +55,12 @@ call atke_compute_km_kh(ngrid,nlay,dtime,& - wind_u,wind_v,temp,play,pinterf, & + wind_u,wind_v,temp,play,pinterf, cdrag_uv,& tke,Km_out,Kh_out) + + if (iflag_num_atke .EQ. 1) then - + !! In this case, we make an explicit prediction of the wind shear to calculate the tke in a + !! forward backward way !! pay attention that the treatment of the TKE !! has to be adapted when solving the TKE with a prognostic equation @@ -71,8 +74,13 @@ call atke_compute_km_kh(ngrid,nlay,dtime,& - wind_u_predict,wind_v_predict,temp,play,pinterf, & + wind_u_predict,wind_v_predict,temp,play,pinterf,cdrag_uv, & tke,Km_out,Kh_out) end if + + + + + end subroutine call_atke