Index: LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90	(revision 4449)
+++ LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90	(revision 4449)
@@ -0,0 +1,159 @@
+module atke_exchange_coeff_mod
+
+implicit none
+
+contains 
+
+subroutine atke_compute_km_kh(ngrid,nlay, &
+                        wind_u,wind_v,temp,play,pinterf, &
+                        tke,Km,Kh)
+
+!========================================================================
+! Routine that computes turbulent Km / Kh coefficients with a
+! 1.5 order closure scheme (TKE) with or without stationarity assumption
+!
+! This parameterization has been constructed in the framework of a
+! collective and collaborative workshop, 
+! the so-called 'Atelier TKE (ATKE)' with
+! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, 
+! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon
+!
+! Main assumptions of the model :
+! (1) dry atmosphere
+! (2) horizontal homogeneity (Dx=Dy=0.)
+!=======================================================================
+
+
+
+USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cn, cinf, rpi
+USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, rg, rd
+
+implicit none
+
+
+! Declarations:
+!=============
+
+INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
+INTEGER, INTENT(IN) :: nlay ! number of vertical index  
+
+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+1), INTENT(OUT)    :: Km   ! Exchange coefficient for momentum at interface between layers
+REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: Kh   ! Exchange coefficient for heat flux at interface between layers
+
+! Local variables
+REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange
+REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface
+REAL, DIMENSION(ngrid,nlay+1) :: dz_interf ! distance between two consecutive interfaces
+REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number
+REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number 
+REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum
+REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat
+
+INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
+REAL    :: Ri0        ! parameter for Sm stability function
+
+
+
+! Initializations:
+!================
+
+DO igrid=1,ngrid
+    z_interf(igrid,1) = 0.0
+    Ri(igrid,1) = 0.1! TO DO
+    Sm(igrid,1) = 0.0 ! TO DO
+END DO
+
+
+
+! Calculation of altitude of layers' bottom interfaces:
+!=====================================================
+
+DO ilay=2,nlay+1
+    DO igrid=1,ngrid
+        dz_interf(igrid,ilay) = rg * temp(igrid,ilay) / rg * log(play(igrid,ilay-1)/play(igrid,ilay))
+        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay)
+    ENDDO
+ENDDO
+
+
+! Computing the mixing length:
+!=============================
+
+DO ilay=2,nlay+1
+    DO igrid=1,ngrid
+        l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
+    ENDDO
+ENDDO
+
+! Computes the gradient Richardson's number and stability functions:
+!===================================================================
+
+DO ilay=2,nlay+1
+    DO igrid=1,ngrid
+        ! Warning: sure that gradient of temp and wind should be calculated with dz_interf and not with dz_lay?
+        Ri(igrid,ilay) = rg * (temp(igrid,ilay) - temp(igrid,ilay-1)) / dz_interf(igrid,ilay) / & 
+        (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 + & 
+        ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 ) 
+
+        ! calculation of Ri0 such that continuity in slope of Sm at Ri=0
+        Ri0=2.*ric*cinf/rpi/cn
+
+        IF (Ri(igrid,ilay) < 0.) THEN
+            Sm(igrid,ilay) = 2./rpi * cinf * atan(-Ri(igrid,ilay)/Ri0) + cn
+            Prandtl(igrid,ilay) = 1. + Ri(igrid,ilay) * pr_slope
+        ELSE
+            Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric))
+            Prandtl(igrid,ilay) = 2./rpi * atan(-Ri(igrid,ilay)) - 1. + pr_asym  ! exp(Ri(igrid,ilay))
+        END IF
+        
+        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
+
+    ENDDO
+ENDDO
+
+! Computing the TKE:
+!===================
+IF (iflag_atke == 0) THEN
+
+! stationary solution neglecting the vertical transport of TKE by turbulence
+    DO ilay=2,nlay+1
+        DO igrid=1,ngrid
+            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
+            (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 + & 
+            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 ) * &
+            (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
+        ENDDO
+    ENDDO
+
+ELSE ! TO DO
+    
+    print*, 'traitement non-stationnaire de la tke pas encore prevu'
+    stop
+
+END IF
+
+
+! Computing eddy diffusivity coefficients:
+!========================================
+DO ilay=2,nlay+1
+    DO igrid=1,ngrid
+        Km(igrid,ilay) = l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
+        Kh(igrid,ilay) = l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
+    END DO
+END DO
+
+
+
+end subroutine atke_compute_km_kh
+
+
+end module atke_exchange_coeff_mod
Index: LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4449)
+++ LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4449)
@@ -0,0 +1,72 @@
+MODULE atke_turbulence_ini_mod
+
+implicit none
+
+save
+
+  integer :: iflag_atke
+  !$OMP THREADPRIVATE(iflag_atke)
+  real :: kappa = 0.4 ! Von Karman constant
+  !$OMP THREADPRIVATE(kappa)
+  real :: l0, ric, ri0, cn, cinf, cepsilon, pr_slope, pr_asym
+  !$OMP THREADPRIVATE(l0, ric, cn, cinf, cepsilon, pr_slope, pr_asym)
+  integer :: lunout,prt_level
+  !$OMP THREADPRIVATE(lunout,prt_level)
+  real :: rg, rd, rpi
+  !$OMP THREADPRIVATE(rg, rd, rpi)
+
+
+CONTAINS
+
+SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_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
+
+
+  lunout=lunout_in
+  prt_level=prt_level_in
+  rd=rd_in
+  rg=rg_in
+  rpi=rpi_in
+
+  ! flag that controls options in atke_compute_km_kh
+  iflag_atke=0
+  CALL getin_p('iflag_atke',iflag_atke)
+
+  ! asymptotic mixing length [m]
+  l0=150.0 
+  CALL getin_p('atke_l0',l0)
+
+  ! critical Richardson number 
+  ric=0.25
+  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
+  CALL getin_p('atke_cinf',cinf)
+
+  ! constant for tke dissipation calculation
+  cepsilon=16.6 ! default value as in yamada4
+  CALL getin_p('atke_cepsilon',cepsilon)
+
+  ! slope of Pr=f(Ri) for stable conditions
+  pr_slope=5.0 ! default value from Zilitinkevich et al. 2005
+  CALL getin_p('atke_pr_slope',pr_slope)
+
+  ! asymptotic turbulent prandt number value for Ri=-Inf
+  pr_asym=0.4
+  CALL getin_p('atke_pr_asym',pr_asym)
+
+   
+ RETURN
+
+END SUBROUTINE atke_ini
+
+END MODULE  atke_turbulence_ini_mod
Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4448)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4449)
@@ -23,4 +23,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 ioipsl_getin_p_mod,  ONLY : getin_p
   USE cdrag_mod
@@ -1656,5 +1657,4 @@
 !!! jyg le 07/02/2012
        IF (iflag_split .eq.0) THEN
-!!!
 !!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
       IF (prt_level >=10) THEN
@@ -1669,6 +1669,13 @@
       print *,' args coef_diff_turb: ycdragh ', ycdragh
       print *,' args coef_diff_turb: ytke ', ytke
-
        ENDIF
+
+        IF (iflag_pbl>=50) THEN
+
+        CALL atke_compute_km_kh(knon,klev,yu,yv,yt, &
+             ypplay,ypaprs,ytke,ycoefm, ycoefh)
+
+        ELSE
+
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
@@ -1686,7 +1693,13 @@
        ENDDO
        ENDIF
+
+       ENDIF ! iflag_pbl >= 50
+
         IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh
-!
+
+
        ELSE  !(iflag_split .eq.0)
+
+      
       IF (prt_level >=10) THEN
       print *,' args coef_diff_turb: yu_x ',  yu_x  
@@ -1699,6 +1712,14 @@
       print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x
       print *,' args coef_diff_turb: ytke_x ', ytke_x
-
-       ENDIF
+      ENDIF
+
+
+        IF (iflag_pbl>=50) THEN
+     
+        CALL atke_compute_km_kh(knon,klev,yu_x,yv_x,yt_x, &
+             ypplay,ypaprs,ytke_x,ycoefm_x, ycoefh_x)
+
+        ELSE
+
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
@@ -1716,4 +1737,7 @@
        ENDDO
        ENDIF
+
+        ENDIF ! iflag_pbl >= 50
+
         IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x
 !
@@ -1728,5 +1752,13 @@
       print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w
       print *,' args coef_diff_turb: ytke_w ', ytke_w
-       ENDIF
+      ENDIF
+      
+        IF (iflag_pbl>=50) THEN
+        
+        CALL atke_compute_km_kh(knon,klev,yu_w,yv_w,yt_w, &
+             ypplay,ypaprs,ytke_w,ycoefm_w, ycoefh_w)
+
+        ELSE
+
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
@@ -1743,6 +1775,10 @@
        ENDDO
        ENDIF
+
+       ENDIF ! iflag_pbl >= 50
+
+
         IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w
-!
+
 !!!jyg le 10/04/2013
 !!   En attendant de traiter le transport des traceurs dans les poches froides, formule
@@ -1754,7 +1790,8 @@
         ENDDO
       ENDDO
-!!!
+
+
        ENDIF  ! (iflag_split .eq.0)
-!!!
+
        
 !****************************************************************************************
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4448)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4449)
@@ -82,5 +82,5 @@
     USE wake_ini_mod, ONLY : wake_ini
     USE yamada_ini_mod, ONLY : yamada_ini
-
+    USE atke_turbulence_ini_mod, ONLY : atke_ini
     USE thermcell_ini_mod, ONLY : thermcell_ini
     USE lscp_ini_mod, ONLY : lscp_ini
@@ -1748,4 +1748,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)
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
