Index: LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4543)
+++ LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4545)
@@ -5,6 +5,6 @@
 save
 
-  integer :: iflag_atke
-  !$OMP THREADPRIVATE(iflag_atke)
+integer :: iflag_atke, iflag_num_atke
+!$OMP THREADPRIVATE(iflag_atke, iflag_num_atke)
   real :: kappa = 0.4 ! Von Karman constant
   !$OMP THREADPRIVATE(kappa)
@@ -45,4 +45,8 @@
   CALL getin_p('iflag_atke',iflag_atke)
 
+  ! flag that controls the numerical treatment of diffusion coeffiient calculation
+  iflag_num_atke=0
+  CALL getin_p('iflag_num_atke',iflag_num_atke)
+
   ! asymptotic mixing length [m]
   l0=150.0 
Index: LMDZ6/trunk/libf/phylmd/call_atke_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/call_atke_mod.F90	(revision 4545)
+++ LMDZ6/trunk/libf/phylmd/call_atke_mod.F90	(revision 4545)
@@ -0,0 +1,160 @@
+module call_atke_mod
+
+USE atke_exchange_coeff_mod, ONLY :  atke_compute_km_kh
+
+implicit none
+
+
+contains
+
+subroutine call_atke(dtime,ngrid,nlay,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
+                        wind_u,wind_v,temp,play,pinterf, &
+                        tke,Km_out,Kh_out)
+
+
+
+
+USE atke_turbulence_ini_mod, ONLY : iflag_num_atke, rg, rd
+
+implicit none
+
+
+! Declarations:
+!=============
+
+REAL, INTENT(IN)    :: dtime ! timestep (s)
+INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
+INTEGER, INTENT(IN) :: nlay ! number of vertical index  
+
+
+REAL, DIMENSION(ngrid), INTENT(IN)       :: cdrag_uv  ! drag coefficient for wind
+REAL, DIMENSION(ngrid), INTENT(IN)       :: cdrag_t   ! drag coefficient for temperature
+
+REAL, DIMENSION(ngrid), INTENT(IN)       :: u_surf    ! x wind velocity at the surface
+REAL, DIMENSION(ngrid), INTENT(IN)       :: v_surf    ! y wind velocity at the surface
+REAL, DIMENSION(ngrid), INTENT(IN)       :: temp_surf ! surface temperature
+
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_u   ! zonal velocity (m/s)
+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)       :: play   ! pressure (Pa)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
+
+
+REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers
+
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers
+
+
+REAL, DIMENSION(ngrid,nlay) :: wind_u_predict, wind_v_predict
+REAL, DIMENSION(ngrid) ::  wind1
+INTEGER i
+
+
+
+call atke_compute_km_kh(ngrid,nlay, &
+                        wind_u,wind_v,temp,play,pinterf, &
+                        tke,Km_out,Kh_out)
+
+if (iflag_num_atke .EQ. 1) then
+
+   !! pay attention that the treatment of the TKE
+   !! has to be adapted when solving the TKE with a prognostic equation
+
+   do i=1,ngrid
+      wind1(i)=sqrt(wind_u(i,1)**2+wind_v(i,1)**2)
+   enddo
+   call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_u,Km_out,u_surf,cdrag_uv,wind_u_predict)  
+   call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_v,Km_out,v_surf,cdrag_uv,wind_v_predict)
+
+
+   call atke_compute_km_kh(ngrid,nlay, &
+                        wind_u_predict,wind_v_predict,temp,play,pinterf, &
+                        tke,Km_out,Kh_out)
+
+end if
+
+end subroutine call_atke
+
+!----------------------------------------------------------------------------------------
+
+subroutine atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,x_in,K_in,x_surf,cdrag,x_predict)
+
+INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
+INTEGER, INTENT(IN) :: nlay ! number of vertical index  
+REAL, INTENT(IN) :: rg,rd,dtime ! gravity, R dry air and timestep
+REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: play      ! pressure middle of layers (Pa)
+REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: temp      ! temperature (K)
+REAL, DIMENSION(ngrid),        INTENT(IN)  :: wind1     ! wind speed first level (m/s)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)  :: pinterf ! pressure at interfaces(Pa)
+REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: x_in    ! variable at the beginning of timestep
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)  :: K_in    ! eddy diffusivity coef at the beginning of time step
+REAL, DIMENSION(ngrid),        INTENT(IN)  :: x_surf  ! surface variable
+REAL, DIMENSION(ngrid),        INTENT(IN)  :: cdrag  ! drag coefficient
+REAL, DIMENSION(ngrid,nlay),   INTENT(OUT) :: x_predict  ! variable at the end of time step after explicit prediction
+
+
+integer i,k
+real ml,F1,rho
+real, dimension(ngrid) :: play1,temp1
+real, dimension(ngrid,nlay+1) :: K_big
+
+! computation of K_big
+
+play1(:)=play(:,1)
+temp1(:)=temp(:,1)
+
+! "big K" calculation
+do  k=2,nlay-1
+   do i=1,ngrid
+      rho=pinterf(i,k)/rd/(0.5*(temp(i,k-1)+temp(i,k)))
+      K_big(i,k)=rg*K_in(i,k)/(play(i,k)-play(i,k+1))*(rho**2)
+   enddo
+enddo
+! speficic treatment for k=nlay
+do i=1,ngrid
+   rho=pinterf(i,nlay)/rd/temp(i,nlay)
+   K_big(i,nlay)=rg*K_in(i,nlay)/(2*(play(i,nlay)-pinterf(i,nlay+1)))*(rho**2)
+enddo
+
+
+
+! x_predict calculation for 2<=k<=nlay-1
+do  k=2,nlay-1
+   do i=1,ngrid
+      ml=(pinterf(i,k)-pinterf(i,k+1))/rg
+      x_predict(i,k)=x_in(i,k)-dtime/ml*(-K_big(i,k+1)*x_in(i,k+1) &
+                  + (K_big(i,k)+K_big(i,k+1))*x_in(i,k) &
+                  - K_big(i,k)*x_in(i,k-1))
+   enddo
+enddo
+
+! Specific treatment for k=1 
+do i=1,ngrid
+   ml=(pinterf(i,1)-pinterf(i,2))/rg
+   F1=-play1(i)/rd/temp1(i)*wind1(i)*cdrag(i)*(x_in(i,1)-x_surf(i)) ! attention convention sens du flux
+   x_predict(i,1)=x_in(i,1)-dtime/ml*(-K_big(i,2)*(x_in(i,2) - x_in(i,1))-F1)
+enddo
+
+! Specific treatment for k=nlay
+! flux at the top of the atmosphere=0
+do i=1,ngrid
+   ml=0.5*(pinterf(i,nlay)-pinterf(i,nlay+1))/rg
+   x_predict(i,nlay)=x_in(i,nlay)+dtime/ml*(K_big(i,nlay)*(x_in(i,nlay)-x_in(i,nlay-1)))
+enddo
+
+!K_big(:,1)=0.
+!do  k=1,nlay
+!   do i=1,ngrid
+!      print*, 'youhou', k, x_in(i,k), x_predict(i,k), K_big(i,k)
+!   end do
+!enddo
+
+
+
+end subroutine atke_explicit_prediction
+
+
+
+end module call_atke_mod
Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4543)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4545)
@@ -24,5 +24,5 @@
   USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
   USE coef_diff_turb_mod,  ONLY : coef_diff_turb
-  USE atke_exchange_coeff_mod, ONLY :  atke_compute_km_kh
+  USE call_atke_mod,       ONLY :  call_atke
   USE ioipsl_getin_p_mod,  ONLY : getin_p
   USE cdrag_mod
@@ -790,4 +790,5 @@
     REAL, DIMENSION(klon)       :: uzon_w, vmer_w, speed_w, zri1_w, pref_w !speed_w, zri1_w, pref_w, added by Fuxing WANG, 04/03/2015
     REAL, DIMENSION(klon)       :: zgeo1_w, tair1_w, qair1_w, tairsol_w
+    REAL, DIMENSION(klon)       :: yus0, yvs0
 
 !!! jyg le 25/03/2013
@@ -993,4 +994,5 @@
  cdragh(:)=0. ; cdragm(:)=0.
  zu1(:)=0. ; zv1(:)=0.
+ yus0(:)=0. ; yvs0(:)=0.
 !albedo SB >>>
   alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0.
@@ -1699,5 +1701,5 @@
         IF (iflag_pbl>=50) THEN
 
-        CALL atke_compute_km_kh(knon,klev,yu,yv,yt, &
+        CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt, &
              ypplay,ypaprs,ytke,ycoefm, ycoefh)
 
@@ -1743,5 +1745,5 @@
         IF (iflag_pbl>=50) THEN
      
-        CALL atke_compute_km_kh(knon,klev,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, &
              ypplay,ypaprs,ytke_x,ycoefm_x, ycoefh_x)
 
@@ -1782,5 +1784,5 @@
         IF (iflag_pbl>=50) THEN
         
-        CALL atke_compute_km_kh(knon,klev,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, &
              ypplay,ypaprs,ytke_w,ycoefm_w, ycoefh_w)
 
