Index: /LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90	(revision 4652)
+++ /LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90	(revision 4653)
@@ -6,5 +6,5 @@
 
 subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
-                        wind_u,wind_v,temp,play,pinterf,cdrag_uv, &
+                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, &
                         tke,Km_out,Kh_out)
 
@@ -26,6 +26,6 @@
 
 
-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 : 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
 
@@ -43,4 +43,5 @@
 REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
 REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: qvap   ! specific humidity (kg/kg)
 REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
 REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
@@ -81,4 +82,6 @@
 REAL    :: netsource  ! net source term of tke
 REAL    :: ustar      ! friction velocity estimation
+REAL    :: invtau     
+REAL    :: rvap
 
 ! Initializations:
@@ -90,9 +93,9 @@
 END DO
 
-! Calculation of potential temperature: (if vapor -> todo virtual potential temperature)
+! Calculation of potential temperature: (if vapor -> virtual potential temperature)
 !=====================================
 
 preff=100000.
-! The result should not depend on the choice of preff
+! results should not depend on the choice of preff
 DO ilay=1,nlay
      DO igrid = 1, ngrid
@@ -101,4 +104,13 @@
 END DO
 
+! account for water vapor mass for buoyancy calculation
+IF (atke_ok_virtual) THEN
+  DO ilay=1,nlay
+     DO igrid = 1, ngrid
+        rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay)))
+        theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap)
+     END DO
+  END DO
+ENDIF
 
 
@@ -287,4 +299,29 @@
 
 
+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", &
@@ -312,5 +349,5 @@
 ! vertical diffusion of TKE 
 !==========================
-IF (ok_vdiff_tke) THEN
+IF (atke_ok_vdiff) THEN
    CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
 ENDIF
@@ -338,5 +375,5 @@
 
 ! routine that computes the vertical diffusion of TKE by the turbulence
-! using an implicit resolution (See note by Dufresne and Ghattas
+! using an implicit resolution (See note by Dufresne and Ghattas (2009))
 ! E Vignon, July 2023
 
Index: /LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4652)
+++ /LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4653)
@@ -13,6 +13,6 @@
   integer :: lunout,prt_level
   !$OMP THREADPRIVATE(lunout,prt_level)
-  real :: rg, rd, rpi, rcpd
-  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd)
+  real :: rg, rd, rpi, rcpd, rv
+  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv)
 
   real :: viscom, viscoh
@@ -22,15 +22,15 @@
   !$OMP THREADPRIVATE(lmin)
 
-  logical :: ok_vdiff_tke
-  !$OMP THREADPRIVATE(ok_vdiff_tke)
+  logical :: atke_ok_vdiff, atke_ok_virtual
+  !$OMP THREADPRIVATE(atke_ok_vdiff,atke_ok_virtual)
 
 CONTAINS
 
-SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in, rcpd_in)
+SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in, rcpd_in, rv_in)
 
    USE ioipsl_getin_p_mod, ONLY : getin_p
 
   integer, intent(in) :: lunout_in,prt_level_in
-  real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in
+  real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in
 
 
@@ -41,4 +41,5 @@
   rpi=rpi_in
   rcpd=rcpd_in
+  rv=rv_in
 
   viscom=1.46E-5
@@ -59,6 +60,12 @@
 
   ! activate vertical diffusion of TKE or not
-  ok_vdiff_tke=.false.
-  CALL getin_p('atke_ok_vdiff_tke',ok_vdiff_tke)
+  atke_ok_vdiff=.false.
+  CALL getin_p('atke_ok_vdiff',atke_ok_vdiff)
+
+
+  ! account for vapor for flottability
+  atke_ok_virtual=.true.
+  CALL getin_p('atke_ok_virtual',atke_ok_virtual)
+
 
   ! flag that controls the numerical treatment of diffusion coeffiient calculation
@@ -82,6 +89,13 @@
 
   ! constant for tke dissipation calculation
-  cepsilon=16.6/2./sqrt(2.) ! default value as in yamada4
+  cepsilon=5.87 ! default value as in yamada4
   CALL getin_p('atke_cepsilon',cepsilon)
+
+
+ ! coefficient for surface TKE 
+ ! following Lenderink & Holtslag 2004, ctkes=sqrt(cepsilon)
+ ! (provided by limit condition in neutral conditions)
+  ctkes=sqrt(cepsilon)
+
 
   ! slope of Pr=f(Ri) for stable conditions
@@ -111,7 +125,4 @@
   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
Index: /LMDZ6/trunk/libf/phylmd/call_atke_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/call_atke_mod.F90	(revision 4652)
+++ /LMDZ6/trunk/libf/phylmd/call_atke_mod.F90	(revision 4653)
@@ -9,5 +9,5 @@
 
 subroutine call_atke(dtime,ngrid,nlay,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
-                        wind_u,wind_v,temp,play,pinterf, &
+                        wind_u,wind_v,temp,qvap,play,pinterf, &
                         tke,Km_out,Kh_out)
 
@@ -38,4 +38,5 @@
 REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
 REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: qvap   ! specific humidity (kg/kg)
 REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
 REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
@@ -55,5 +56,5 @@
 
 call atke_compute_km_kh(ngrid,nlay,dtime,&
-                        wind_u,wind_v,temp,play,pinterf, cdrag_uv,&
+                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv,&
                         tke,Km_out,Kh_out)
 
@@ -74,5 +75,5 @@
 
    call atke_compute_km_kh(ngrid,nlay,dtime,&
-                        wind_u_predict,wind_v_predict,temp,play,pinterf,cdrag_uv, &
+                        wind_u_predict,wind_v_predict,temp,qvap,play,pinterf,cdrag_uv, &
                         tke,Km_out,Kh_out)
 
Index: /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4652)
+++ /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4653)
@@ -1704,5 +1704,5 @@
         IF (iflag_pbl>=50) THEN
 
-        CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt, &
+        CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt,yq, &
              ypplay,ypaprs,ytke,ycoefm, ycoefh)
 
@@ -1748,5 +1748,5 @@
         IF (iflag_pbl>=50) THEN
      
-        CALL call_atke(dtime,knon,klev,ycdragm_x,ycdragh_x,yus0,yvs0,yts_x,yu_x,yv_x,yt_x, &
+        CALL call_atke(dtime,knon,klev,ycdragm_x,ycdragh_x,yus0,yvs0,yts_x,yu_x,yv_x,yt_x,yq_x, &
              ypplay,ypaprs,ytke_x,ycoefm_x, ycoefh_x)
 
@@ -1787,5 +1787,5 @@
         IF (iflag_pbl>=50) THEN
         
-        CALL call_atke(dtime,knon,klev,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w,yu_w,yv_w,yt_w, &
+        CALL call_atke(dtime,knon,klev,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w,yu_w,yv_w,yt_w,yq_w, &
              ypplay,ypaprs,ytke_w,ycoefm_w, ycoefh_w)
 
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4652)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4653)
@@ -1805,5 +1805,5 @@
        CALL wake_ini(rg,rd,rv,prt_level)
        CALL yamada_ini(klon,lunout,prt_level)
-       CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD)
+       CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD, RV)
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
Index: /LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90	(revision 4652)
+++ /LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90	(revision 4653)
@@ -1981,5 +1981,5 @@
         IF (iflag_pbl>=50) THEN
 
-        CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt, &
+        CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt,yq, &
              ypplay,ypaprs,ytke,ycoefm, ycoefh)
 
@@ -2024,5 +2024,5 @@
         IF (iflag_pbl>=50) THEN
 
-        CALL call_atke(dtime,knon,klev,ycdragm_x,ycdragh_x,yus0,yvs0,yts_x,yu_x,yv_x,yt_x, &
+        CALL call_atke(dtime,knon,klev,ycdragm_x,ycdragh_x,yus0,yvs0,yts_x,yu_x,yv_x,yt_x,yq_x, &
              ypplay,ypaprs,ytke_x,ycoefm_x, ycoefh_x)
 
@@ -2064,5 +2064,5 @@
         IF (iflag_pbl>=50) THEN
 
-        CALL call_atke(dtime,knon,klev,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w,yu_w,yv_w,yt_w, &
+        CALL call_atke(dtime,knon,klev,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w,yu_w,yv_w,yt_w,yq_w, &
              ypplay,ypaprs,ytke_w,ycoefm_w, ycoefh_w)
 
Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4652)
+++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4653)
@@ -1917,5 +1917,5 @@
        CALL wake_ini(rg,rd,rv,prt_level)
        CALL yamada_ini(klon,lunout,prt_level)
-       CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD)
+       CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD, RV)
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
