Changeset 4663

Timestamp:

Sep 1, 2023, 4:01:08 PM (9 months ago)

Author:

evignon

Message:

ajout d'une longueur de melange sensible au cisaillement de vent dans atke

Location:

LMDZ6/trunk/libf/phylmd

Files:

• atke_exchange_coeff_mod.F90 (modified) (6 diffs)
• atke_turbulence_ini_mod.F90 (modified) (3 diffs)

LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90

@@ -28 +28 @@
 USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual
 USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd, rv, atke_ok_vdiff
-USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin, smmin
+USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin
 
 implicit none
@@ -68 +68 @@
 REAL, DIMENSION(ngrid,nlay+1) :: Sm          ! Stability function for momentum (at interface)
 REAL, DIMENSION(ngrid,nlay+1) :: Sh          ! Stability function for heat (at interface)
-LOGICAL, DIMENSION(ngrid,nlay+1) :: switch_num ! switch of numerical integration in very stable cases
 
 INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
@@ -174 +173 @@
 !==============================================================
 
-switch_num(:,:)=.false.
 
 IF (iflag_atke_lmix .EQ. 1 ) THEN
@@ -183 +181 @@
           IF (N2(igrid,ilay) .GT. 0.) THEN
              lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
-             IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
-                l_exchange(igrid,ilay)=max(lstrat,lmin)
-                switch_num(igrid,ilay)=.true.
-             ENDIF
+             lstrat=max(lstrat,lmin)
+             !Monin-Obukhov consistent interpolation, Van de Wiel et al. 2010
+             l_exchange(igrid,ilay)=(1./(2.*l_exchange(igrid,ilay))+sqrt(1./(4.*l_exchange(igrid,ilay) &
+                                   *l_exchange(igrid,ilay))+1./(2.*lstrat*lstrat)))**(-1.0)
           ENDIF 
       ENDDO
@@ -192 +190 @@
 
 ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
-! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms, eq 2
+! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms
 DO ilay=2,nlay
       DO igrid=1,ngrid
           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)),1.E-10)*(1.+sqrt(Ri(igrid,ilay))/2.))
-             IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
-                l_exchange(igrid,ilay)=max(lstrat,lmin)
-             ENDIF
+          IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
+             lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), &
+                    clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay)))
+             lstrat=max(lstrat,lmin)
+             !Monin-Obukhov consistent interpolation, Van de Wiel et al. 2010   
+             l_exchange(igrid,ilay)=(1./(2.*l_exchange(igrid,ilay))+sqrt(1./(4.*l_exchange(igrid,ilay) &
+                                   *l_exchange(igrid,ilay))+1./(2.*lstrat*lstrat)))**(-1.0)
+
           ENDIF
       ENDDO
@@ -299 +300 @@
 
 
-ELSE IF (iflag_atke == 5) THEN
-! semi implicit scheme from Arpege (V. Masson methodology with 
-! Taylor expansion of the dissipation term)
-! and implicit resolution when switch num (when we use the mixing length as a function of tke)
-    DO ilay=2,nlay
-        DO igrid=1,ngrid
-           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
-           if (switch_num(igrid,ilay) .and. N2(igrid,ilay)>0.) then
-             invtau=clmix*Sm(igrid,ilay)/sqrt(N2(igrid,ilay))*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))  & 
-                   -sqrt(N2(igrid,ilay))/(cepsilon*clmix)
-             !qq=qq/(1.-dtime*invtau)
-             !qq=qq*exp(dtime*invtau)
-             tke(igrid,ilay)=tke(igrid,ilay)/(1.-dtime*invtau)
-             tke(igrid,ilay)=max(0.,tke(igrid,ilay))
-           else
-             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)
-           endif
-        ENDDO
-    ENDDO
-
-
 ELSE
    call abort_physic("atke_compute_km_kh", &

LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90

@@ -9 +9 @@
   real :: kappa = 0.4 ! Von Karman constant
   !$OMP THREADPRIVATE(kappa)
-  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)
+  real :: l0,ric,ri0,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke
+  !$OMP THREADPRIVATE(l0,ric,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke)
   integer :: lunout,prt_level
   !$OMP THREADPRIVATE(lunout,prt_level)
@@ -94 +94 @@
 
  ! coefficient for surface TKE 
- ! following Lenderink & Holtslag 2004, ctkes=sqrt(cepsilon)
+ ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3)
  ! (provided by limit condition in neutral conditions)
-  ctkes=sqrt(cepsilon)
-
+  ctkes=(cepsilon)**(2./3.)
 
   ! slope of Pr=f(Ri) for stable conditions
@@ -118 +117 @@
   clmix=0.5
   CALL getin_p('atke_clmix',clmix)
-  
+ 
+  ! coefficient for mixing length depending on local wind shear
+  clmixshear=0.5
+  CALL getin_p('atke_clmixshear',clmixshear)
+
+
   ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri. 
   ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17