Changeset 4663 for LMDZ6/trunk
- Timestamp:
- Sep 1, 2023, 4:01:08 PM (14 months ago)
- Location:
- LMDZ6/trunk/libf/phylmd
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90
r4653 r4663 28 28 USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual 29 29 USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd, rv, atke_ok_vdiff 30 USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin, smmin30 USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin 31 31 32 32 implicit none … … 68 68 REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum (at interface) 69 69 REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat (at interface) 70 LOGICAL, DIMENSION(ngrid,nlay+1) :: switch_num ! switch of numerical integration in very stable cases71 70 72 71 INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) … … 174 173 !============================================================== 175 174 176 switch_num(:,:)=.false.177 175 178 176 IF (iflag_atke_lmix .EQ. 1 ) THEN … … 183 181 IF (N2(igrid,ilay) .GT. 0.) THEN 184 182 lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)) 185 IF (lstrat .LT. l_exchange(igrid,ilay)) THEN186 l_exchange(igrid,ilay)=max(lstrat,lmin)187 switch_num(igrid,ilay)=.true.188 ENDIF183 lstrat=max(lstrat,lmin) 184 !Monin-Obukhov consistent interpolation, Van de Wiel et al. 2010 185 l_exchange(igrid,ilay)=(1./(2.*l_exchange(igrid,ilay))+sqrt(1./(4.*l_exchange(igrid,ilay) & 186 *l_exchange(igrid,ilay))+1./(2.*lstrat*lstrat)))**(-1.0) 189 187 ENDIF 190 188 ENDDO … … 192 190 193 191 ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN 194 ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms , eq 2192 ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms 195 193 DO ilay=2,nlay 196 194 DO igrid=1,ngrid 197 195 l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) 198 IF (N2(igrid,ilay) .GT. 0.) THEN 199 lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1.E-10)*(1.+sqrt(Ri(igrid,ilay))/2.)) 200 IF (lstrat .LT. l_exchange(igrid,ilay)) THEN 201 l_exchange(igrid,ilay)=max(lstrat,lmin) 202 ENDIF 196 IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN 197 lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), & 198 clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))) 199 lstrat=max(lstrat,lmin) 200 !Monin-Obukhov consistent interpolation, Van de Wiel et al. 2010 201 l_exchange(igrid,ilay)=(1./(2.*l_exchange(igrid,ilay))+sqrt(1./(4.*l_exchange(igrid,ilay) & 202 *l_exchange(igrid,ilay))+1./(2.*lstrat*lstrat)))**(-1.0) 203 203 204 ENDIF 204 205 ENDDO … … 299 300 300 301 301 ELSE IF (iflag_atke == 5) THEN302 ! semi implicit scheme from Arpege (V. Masson methodology with303 ! Taylor expansion of the dissipation term)304 ! and implicit resolution when switch num (when we use the mixing length as a function of tke)305 DO ilay=2,nlay306 DO igrid=1,ngrid307 qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)308 if (switch_num(igrid,ilay) .and. N2(igrid,ilay)>0.) then309 invtau=clmix*Sm(igrid,ilay)/sqrt(N2(igrid,ilay))*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &310 -sqrt(N2(igrid,ilay))/(cepsilon*clmix)311 !qq=qq/(1.-dtime*invtau)312 !qq=qq*exp(dtime*invtau)313 tke(igrid,ilay)=tke(igrid,ilay)/(1.-dtime*invtau)314 tke(igrid,ilay)=max(0.,tke(igrid,ilay))315 else316 qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &317 +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &318 /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))319 qq=max(0.,qq)320 tke(igrid,ilay)=0.5*(qq**2)321 endif322 ENDDO323 ENDDO324 325 326 302 ELSE 327 303 call abort_physic("atke_compute_km_kh", & -
LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90
r4653 r4663 9 9 real :: kappa = 0.4 ! Von Karman constant 10 10 !$OMP THREADPRIVATE(kappa) 11 real :: l0, ric, ri0, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix, smmin,ctkes,cke12 !$OMP THREADPRIVATE(l0,ric,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix, smmin,ctkes,cke)11 real :: l0,ric,ri0,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke 12 !$OMP THREADPRIVATE(l0,ric,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke) 13 13 integer :: lunout,prt_level 14 14 !$OMP THREADPRIVATE(lunout,prt_level) … … 94 94 95 95 ! coefficient for surface TKE 96 ! following Lenderink & Holtslag 2004, ctkes= sqrt(cepsilon)96 ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3) 97 97 ! (provided by limit condition in neutral conditions) 98 ctkes=sqrt(cepsilon) 99 98 ctkes=(cepsilon)**(2./3.) 100 99 101 100 ! slope of Pr=f(Ri) for stable conditions … … 118 117 clmix=0.5 119 118 CALL getin_p('atke_clmix',clmix) 120 119 120 ! coefficient for mixing length depending on local wind shear 121 clmixshear=0.5 122 CALL getin_p('atke_clmixshear',clmixshear) 123 124 121 125 ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri. 122 126 ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17
Note: See TracChangeset
for help on using the changeset viewer.