Index: LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90	(revision 4479)
+++ LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90	(revision 4481)
@@ -26,5 +26,5 @@
 
 
-USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cn, cinf, rpi, rcpd
+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
@@ -66,5 +66,5 @@
 
 INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
-REAL    :: Ri0,Ri1    ! parameter for Sm stability function and Prandlt
+REAL    :: cn,Ri0,Ri1    ! parameter for Sm stability function and Prandlt
 REAL    :: preff      ! reference pressure for potential temperature calculations
 REAL    :: thetam     ! mean potential temperature at interface
@@ -77,10 +77,4 @@
     dz_interf(igrid,1) = 0.0
     z_interf(igrid,1) = 0.0
-!    Km(igrid,1) = 0.0
-!    Kh(igrid,1) = 0.0
-!    tke(igrid,1) = 0.0
-!    Km(igrid,nlay+1) = 0.0              
-!    Kh(igrid,nlay+1) = 0.0             
-!    tke(igrid,nlay+1) = 0.0
 END DO
 
@@ -129,4 +123,13 @@
 !===================================================================
 
+! calculation of cn = Sm value at Ri=0
+! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
+cn=(1./sqrt(cepsilon))**(2/3)
+! calculation of Ri0 such that continuity in slope of Sm at Ri=0
+Ri0=2./rpi*(cinf - cn)*ric/cn
+! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
+Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope
+
+
 DO ilay=2,nlay
     DO igrid=1,ngrid
@@ -136,15 +139,15 @@
         MAX(((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & 
         ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2,1E-10) 
-        ! calculation of Ri0 such that continuity in slope of Sm at Ri=0
-        Ri0=2./rpi*(cinf - cn)*ric/cn
-        ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
-        Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope
-
-        IF (Ri(igrid,ilay) < 0.) THEN
+        
+        IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
             Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn
             Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut 
-        ELSE
+        ELSE ! stable cases
             Sm(igrid,ilay) = max(0.,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
+               call abort_physic("atke_compute_km_kh", &
+               'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)
+            ENDIF
         END IF
         
Index: LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4479)
+++ LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4481)
@@ -9,6 +9,6 @@
   real :: kappa = 0.4 ! Von Karman constant
   !$OMP THREADPRIVATE(kappa)
-  real :: l0, ric, ri0, cn, cinf, cepsilon, pr_slope, pr_asym, pr_neut
-  !$OMP THREADPRIVATE(l0, ric, cn, cinf, cepsilon, pr_slope, pr_asym, pr_neut)
+  real :: l0, ric, ri0, cinf, cepsilon, pr_slope, pr_asym, pr_neut
+  !$OMP THREADPRIVATE(l0, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut)
   integer :: lunout,prt_level
   !$OMP THREADPRIVATE(lunout,prt_level)
@@ -53,8 +53,4 @@
   CALL getin_p('atke_ric',ric)
 
-  ! value of Sm stability function for neutral conditions (Ri=0) 
-  cn=1.0
-  CALL getin_p('atke_cn',cn)
-
   ! asymptotic value of Sm for Ri=-Inf
   cinf=1.5
@@ -74,6 +70,5 @@
 
   ! value of turbulent prandtl number in neutral conditions (Ri=0)
-  ! to guarantee tke>0 in stationary case, pr_neut has to be >= 1
-  pr_neut=1.0
+  pr_neut=0.8
   CALL getin_p('atke_pr_neut',pr_neut)
 
Index: LMDZ6/trunk/libf/phylmd/cdrag_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cdrag_mod.F90	(revision 4479)
+++ LMDZ6/trunk/libf/phylmd/cdrag_mod.F90	(revision 4481)
@@ -23,5 +23,5 @@
   USE print_control_mod, ONLY: lunout, prt_level
   USE ioipsl_getin_p_mod, ONLY : getin_p
-  USE atke_turbulence_ini_mod, ONLY : ric, cn, cinf, cepsilon, pr_slope, pr_asym, pr_neut
+  USE atke_turbulence_ini_mod, ONLY : ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut
 
   IMPLICIT NONE
@@ -128,5 +128,5 @@
   REAL zzzcd
   REAL, DIMENSION(klon) :: sm, prandtl              ! Stability function from atke turbulence scheme
-  REAL                  ri0, ri1                    ! to have dimensionless term in sm and prandtl
+  REAL                  ri0, ri1, cn                ! to have dimensionless term in sm and prandtl
 
   LOGICAL, SAVE :: firstcall = .TRUE.
@@ -168,4 +168,9 @@
 
   ALPHA=5.0
+
+! Consistent with atke scheme
+cn=(1./sqrt(cepsilon))**(2/3)
+ri0=2./rpi*(cinf - cn)*ric/cn
+ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope
 
 
@@ -246,4 +251,6 @@
 !*******************************************************
 
+
+
 !''''''''''''''
 ! Cas instables
@@ -291,6 +298,4 @@
          CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
          
-           ri0=2./rpi*(cinf - cn)*ric/cn
-           ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope
            sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
            prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut