Index: LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90	(revision 4686)
+++ 	(revision )
@@ -1,451 +1,0 @@
-module atke_exchange_coeff_mod
-
-implicit none
-
-contains 
-
-subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
-                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, &
-                        tke,Km_out,Kh_out)
-
-!========================================================================
-! 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, T Dubos, 
-! 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, 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, clmixshear, iflag_atke_lmix, lmin, smmin
-
-implicit none
-
-
-! Declarations:
-!=============
-
-INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
-INTEGER, INTENT(IN) :: nlay  ! number of vertical index  
-
-REAL, INTENT(IN)    :: dtime ! physics time step (s)
-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)       :: 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)
-REAL, DIMENSION(ngrid), INTENT(IN)            :: cdrag_uv   ! surface drag coefficient for momentum
-
-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
-
-! Local variables
-REAL, DIMENSION(ngrid,nlay+1) :: Km          ! Exchange coefficient for momentum at interface between layers
-REAL, DIMENSION(ngrid,nlay+1) :: Kh          ! Exchange coefficient for heat flux at interface between layers
-REAL, DIMENSION(ngrid,nlay)   :: theta       ! Potential temperature
-REAL, DIMENSION(ngrid,nlay+1) :: l_exchange  ! Length of exchange (at interface)
-REAL, DIMENSION(ngrid,nlay+1) :: z_interf    ! Altitude at the interface
-REAL, DIMENSION(ngrid,nlay)   :: z_lay       ! Altitude of layers
-REAL, DIMENSION(ngrid,nlay)   :: dz_interf   ! distance between two consecutive interfaces
-REAL, DIMENSION(ngrid,nlay)   :: dz_lay      ! distance between two layer middles (NB: first and last are half layers)
-REAL, DIMENSION(ngrid,nlay+1) :: N2          ! square of Brunt Vaisala pulsation (at interface)
-REAL, DIMENSION(ngrid,nlay+1) :: shear2      ! square of wind shear (at interface)
-REAL, DIMENSION(ngrid,nlay+1) :: Ri          ! Richardson's number (at interface)
-REAL, DIMENSION(ngrid,nlay+1) :: Prandtl     ! Turbulent Prandtl's number (at interface)
-REAL, DIMENSION(ngrid,nlay+1) :: Sm          ! Stability function for momentum (at interface)
-REAL, DIMENSION(ngrid,nlay+1) :: Sh          ! Stability function for heat (at interface)
-
-INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
-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
-REAL    :: delta      ! discriminant of the second order polynomial
-REAL    :: qq         ! tke=qq**2/2
-REAL    :: shear      ! wind shear
-REAL    :: lstrat     ! mixing length depending on local stratification
-REAL    :: taustrat   ! caracteristic timescale for turbulence in very stable conditions
-REAL    :: netloss    ! net loss term of tke
-REAL    :: netsource  ! net source term of tke
-REAL    :: ustar      ! friction velocity estimation
-REAL    :: invtau     
-REAL    :: rvap
-
-! Initializations:
-!================
-
-DO igrid=1,ngrid
-    dz_interf(igrid,1) = 0.0
-    z_interf(igrid,1) = 0.0
-END DO
-
-! Calculation of potential temperature: (if vapor -> virtual potential temperature)
-!=====================================
-
-preff=100000.
-! results should not depend on the choice of preff
-DO ilay=1,nlay
-     DO igrid = 1, ngrid
-        theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
-     END DO
-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
-
-
-! Calculation of altitude of layers' middle and bottom interfaces:
-!=================================================================
-
-DO ilay=2,nlay+1
-    DO igrid=1,ngrid
-        dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
-        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
-    ENDDO
-ENDDO
-
-DO ilay=1,nlay
-   DO igrid=1,ngrid 
-      z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
-   ENDDO
-ENDDO
-
-
-! Computes the gradient Richardson's number and stability functions:
-!===================================================================
-
-! 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
-        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
-        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
-        N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay)
-        shear2(igrid,ilay)= (((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 )
-        Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10) 
-        
-        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 ! stable cases
-            Sm(igrid,ilay) = max(smmin,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
-        
-        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
-
-    ENDDO
-ENDDO
-
-
-! Computing the mixing length:
-!==============================================================
-
-
-IF (iflag_atke_lmix .EQ. 1 ) THEN
-
-   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))/sqrt(N2(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
-   ENDDO
-
-ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
-! 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. .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
-   ENDDO
-
-
-
-ELSE
-! default: neglect effect of local stratification and shear
-
-   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
-ENDIF
-
-
-! Computing the TKE k>=2:
-!========================
-IF (iflag_atke == 0) THEN
-
-! stationary solution (dtke/dt=0)
-
-   DO ilay=2,nlay
-        DO igrid=1,ngrid
-            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
-            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
-        ENDDO
-    ENDDO
-
-ELSE IF (iflag_atke == 1) THEN
-
-! full implicit scheme resolved with a second order polynomial equation
-
-    DO ilay=2,nlay
-        DO igrid=1,ngrid
-           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
-           delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. &
-                 +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + &
-                 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) &
-                 *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)))
-           qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2.
-           qq=max(0.,qq)
-           tke(igrid,ilay)=0.5*(qq**2) 
-        ENDDO
-    ENDDO
-
-
-ELSE IF (iflag_atke == 2) THEN
-
-! semi implicit scheme when l does not depend on tke
-! positive-guaranteed if pr slope in stable condition >1
-
-   DO ilay=2,nlay
-        DO igrid=1,ngrid
-           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
-           qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      & 
-               *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
-               /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) 
-           tke(igrid,ilay)=0.5*(qq**2) 
-        ENDDO
-    ENDDO
-
-
-ELSE IF (iflag_atke == 3) THEN
-! numerical resolution adapted from that in MAR (Deleersnijder 1992)
-! positively defined by construction
-
-    DO ilay=2,nlay
-        DO igrid=1,ngrid
-           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
-           IF (Ri(igrid,ilay) .LT. 0.) THEN
-              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
-              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
-           ELSE
-              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &
-                      l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)
-              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)
-           ENDIF
-           qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)
-           tke(igrid,ilay)=0.5*(qq**2)
-        ENDDO
-    ENDDO
-
-ELSE IF (iflag_atke == 4) THEN
-! semi implicit scheme from Arpege (V. Masson methodology with 
-! Taylor expansion of the dissipation term)
-    DO ilay=2,nlay
-        DO igrid=1,ngrid
-           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
-           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)
-        ENDDO
-    ENDDO
-
-
-ELSE
-   call abort_physic("atke_compute_km_kh", &
-        'numerical treatment of TKE not possible yet', 1)
-
-END IF
-
-! We impose a 0 tke at nlay+1
-!==============================
-
-DO igrid=1,ngrid
- tke(igrid,nlay+1)=0.
-END DO
-
-
-! Calculation of surface TKE (k=1)
-!=================================
-! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
-DO igrid=1,ngrid
- ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
- tke(igrid,1)=ctkes*(ustar**2)
-END DO
-
-
-! vertical diffusion of TKE 
-!==========================
-IF (atke_ok_vdiff) THEN
-   CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
-ENDIF
-
-
-! Computing eddy diffusivity coefficients:
-!========================================
-DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
-    DO igrid=1,ngrid
-        ! we add the molecular viscosity to Km,h 
-        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
-        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
-    END DO
-END DO
-
-! for output:
-!===========
-Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
-Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
-
-end subroutine atke_compute_km_kh
-
-!===============================================================================================
-subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
-
-! routine that computes the vertical diffusion of TKE by the turbulence
-! using an implicit resolution (See note by Dufresne and Ghattas (2009))
-! E Vignon, July 2023
-
-USE atke_turbulence_ini_mod, ONLY : rd, cke, viscom
-
-
-INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
-INTEGER, INTENT(IN) :: nlay  ! number of vertical index  
-
-REAL, INTENT(IN)    :: dtime ! physics time step (s)
-REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: z_lay   ! altitude of mid-layers (m)
-REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)       :: z_interf   ! altitude of bottom interfaces (m)
-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)     :: l_exchange     ! mixing length at interfaces between layers
-REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: Sm     ! stability function for eddy diffusivity for momentum at interface between layers
-
-REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke    ! turbulent kinetic energy at interface between layers
-
-
-
-INTEGER                                       :: igrid,ilay
-REAL, DIMENSION(ngrid,nlay+1)                 :: Ke     ! eddy diffusivity for TKE
-REAL, DIMENSION(ngrid,nlay+1)                 :: dtke
-REAL, DIMENSION(ngrid,nlay+1)                 :: ak, bk, ck, CCK, DDK
-REAL                                          :: gammak,Kem,KKb,KKt
-
-
-! Few initialisations
-CCK(:,:)=0.
-DDK(:,:)=0.
-dtke(:,:)=0.
-
-
-! Eddy diffusivity for TKE
-
-DO ilay=2,nlay
-    DO igrid=1,ngrid
-       Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke 
-    ENDDO
-ENDDO
-! at the top of the atmosphere set to 0
-Ke(:,nlay+1)=0.
-! at the surface, set it equal to that at the first model level
-Ke(:,1)=Ke(:,2)
-
-
-! calculate intermediary variables
-
-DO ilay=2,nlay 
-    DO igrid=1,ngrid
-    Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay))    
-    KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay))
-    Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1))
-    KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1))
-    gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1))
-    ak(igrid,ilay)=-gammak*dtime*KKb
-    ck(igrid,ilay)=-gammak*dtime*KKt
-    bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb)
-    ENDDO
-ENDDO
-
-! calculate CCK and DDK coefficients
-! downhill phase
-
-DO igrid=1,ngrid
-  CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
-  DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
-ENDDO
-
-
-DO ilay=nlay-1,2,-1
-    DO igrid=1,ngrid
-        CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
-                       / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
-        DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
-    ENDDO
-ENDDO
-
-! calculate TKE
-! uphill phase
-
-DO ilay=2,nlay+1
-    DO igrid=1,ngrid
-        dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay)
-    ENDDO
-ENDDO
-
-! update TKE
-tke(:,:)=tke(:,:)+dtke(:,:)
-
-
-end subroutine atke_vdiff_tke
-
-
-
-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 4686)
+++ 	(revision )
@@ -1,141 +1,0 @@
-MODULE atke_turbulence_ini_mod
-
-implicit none
-
-save
-
-integer :: iflag_atke, iflag_num_atke, iflag_atke_lmix
-  !$OMP THREADPRIVATE(iflag_atke, iflag_num_atke, iflag_atke_lmix)
-  real :: kappa = 0.4 ! Von Karman constant
-  !$OMP THREADPRIVATE(kappa)
-  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)
-  real :: rg, rd, rpi, rcpd, rv
-  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv)
-
-  real :: viscom, viscoh
-  !$OMP THREADPRIVATE(viscom,viscoh)
-
-  real :: lmin=0.01              ! minimum mixing length
-  !$OMP THREADPRIVATE(lmin)
-
-  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, 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, rv_in
-
-
-  lunout=lunout_in
-  prt_level=prt_level_in
-  rd=rd_in
-  rg=rg_in
-  rpi=rpi_in
-  rcpd=rcpd_in
-  rv=rv_in
-
-  viscom=1.46E-5
-  viscoh=2.06E-5
-
-  ! flag that controls options in atke_compute_km_kh
-  iflag_atke=0
-  CALL getin_p('iflag_atke',iflag_atke)
-
-  ! flag that controls the calculation of mixing length in atke
-  iflag_atke_lmix=0
-  CALL getin_p('iflag_atke_lmix',iflag_atke_lmix)
-
-  if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then
-        call abort_physic("atke_turbulence_ini", &
-        'stationary scheme must use mixing length formulation not depending on tke', 1)
-  endif
-
-  ! activate vertical diffusion of TKE or not
-  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
-  iflag_num_atke=0
-  CALL getin_p('iflag_num_atke',iflag_num_atke)
-
-  ! asymptotic mixing length in neutral conditions [m]
-  ! Sun et al 2011, JAMC 
-  ! between 10 and 40
-
-  l0=15.0 
-  CALL getin_p('atke_l0',l0)
-
-  ! critical Richardson number 
-  ric=0.25
-  CALL getin_p('atke_ric',ric)
-
-  ! asymptotic value of Sm for Ri=-Inf
-  cinf=1.5
-  CALL getin_p('atke_cinf',cinf)
-
-  ! constant for tke dissipation calculation
-  cepsilon=5.87 ! default value as in yamada4
-  CALL getin_p('atke_cepsilon',cepsilon)
-
-
- ! coefficient for surface TKE 
- ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3)
- ! (provided by limit condition in neutral conditions)
-  ctkes=(cepsilon)**(2./3.)
-
-  ! 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)
-  if (pr_slope .le. 1) then
-        call abort_physic("atke_turbulence_ini", &
-        'pr_slope has to be greater than 1 for consistency of the tke scheme', 1)
-  endif
-
-  ! asymptotic turbulent prandt number value for Ri=-Inf
-  pr_asym=0.4
-  CALL getin_p('atke_pr_asym',pr_asym)
-
-  ! value of turbulent prandtl number in neutral conditions (Ri=0)
-  pr_neut=0.8
-  CALL getin_p('atke_pr_neut',pr_neut)
-
-  ! coefficient for mixing length depending on local stratification
-  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
-  
-  smmin=0.17
-  CALL getin_p('atke_smmin',smmin)
-
-
-  ! ratio between the eddy diffusivity coeff for tke wrt that for momentum
-  ! default value from Lenderink et al. 2004
-  cke=2.
-  CALL getin_p('atke_cke',cke)
-
- RETURN
-
-END SUBROUTINE atke_ini
-
-END MODULE  atke_turbulence_ini_mod
Index: LMDZ6/trunk/libf/phylmd/call_atke_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/call_atke_mod.F90	(revision 4686)
+++ 	(revision )
@@ -1,169 +1,0 @@
-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,qvap,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)       :: 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)
-
-
-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,dtime,&
-                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv,&
-                        tke,Km_out,Kh_out)
-
-
-                               
-if (iflag_num_atke .EQ. 1) then
-   !! In this case, we make an explicit prediction of the wind shear to calculate the tke in a
-   !! forward backward way
-   !! 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,dtime,&
-                        wind_u_predict,wind_v_predict,temp,qvap,play,pinterf,cdrag_uv, &
-                        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/cdrag_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cdrag_mod.F90	(revision 4686)
+++ LMDZ6/trunk/libf/phylmd/cdrag_mod.F90	(revision 4687)
@@ -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, cinf, cepsilon, pr_slope, pr_asym, pr_neut
+  USE lmdz_atke_turbulence_ini, ONLY : ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut
 
   IMPLICIT NONE
Index: LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90	(revision 4687)
+++ LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90	(revision 4687)
@@ -0,0 +1,448 @@
+module lmdz_atke_exchange_coeff
+
+implicit none
+
+contains 
+
+subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
+                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, &
+                        tke,Km_out,Kh_out)
+
+!========================================================================
+! 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, T Dubos, 
+! 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 lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual
+USE lmdz_atke_turbulence_ini, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd, rv, atke_ok_vdiff
+USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin
+
+implicit none
+
+
+! Declarations:
+!=============
+
+INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
+INTEGER, INTENT(IN) :: nlay  ! number of vertical index  
+
+REAL, INTENT(IN)    :: dtime ! physics time step (s)
+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)       :: 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)
+REAL, DIMENSION(ngrid), INTENT(IN)            :: cdrag_uv   ! surface drag coefficient for momentum
+
+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
+
+! Local variables
+REAL, DIMENSION(ngrid,nlay+1) :: Km          ! Exchange coefficient for momentum at interface between layers
+REAL, DIMENSION(ngrid,nlay+1) :: Kh          ! Exchange coefficient for heat flux at interface between layers
+REAL, DIMENSION(ngrid,nlay)   :: theta       ! Potential temperature
+REAL, DIMENSION(ngrid,nlay+1) :: l_exchange  ! Length of exchange (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: z_interf    ! Altitude at the interface
+REAL, DIMENSION(ngrid,nlay)   :: z_lay       ! Altitude of layers
+REAL, DIMENSION(ngrid,nlay)   :: dz_interf   ! distance between two consecutive interfaces
+REAL, DIMENSION(ngrid,nlay)   :: dz_lay      ! distance between two layer middles (NB: first and last are half layers)
+REAL, DIMENSION(ngrid,nlay+1) :: N2          ! square of Brunt Vaisala pulsation (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: shear2      ! square of wind shear (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Ri          ! Richardson's number (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Prandtl     ! Turbulent Prandtl's number (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Sm          ! Stability function for momentum (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Sh          ! Stability function for heat (at interface)
+
+INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
+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
+REAL    :: delta      ! discriminant of the second order polynomial
+REAL    :: qq         ! tke=qq**2/2
+REAL    :: shear      ! wind shear
+REAL    :: lstrat     ! mixing length depending on local stratification
+REAL    :: taustrat   ! caracteristic timescale for turbulence in very stable conditions
+REAL    :: netloss    ! net loss term of tke
+REAL    :: netsource  ! net source term of tke
+REAL    :: ustar      ! friction velocity estimation
+REAL    :: invtau     
+REAL    :: rvap
+
+! Initializations:
+!================
+
+DO igrid=1,ngrid
+    dz_interf(igrid,1) = 0.0
+    z_interf(igrid,1) = 0.0
+END DO
+
+! Calculation of potential temperature: (if vapor -> virtual potential temperature)
+!=====================================
+
+preff=100000.
+! results should not depend on the choice of preff
+DO ilay=1,nlay
+     DO igrid = 1, ngrid
+        theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
+     END DO
+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
+
+
+! Calculation of altitude of layers' middle and bottom interfaces:
+!=================================================================
+
+DO ilay=2,nlay+1
+    DO igrid=1,ngrid
+        dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
+        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
+    ENDDO
+ENDDO
+
+DO ilay=1,nlay
+   DO igrid=1,ngrid 
+      z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
+   ENDDO
+ENDDO
+
+
+! Computes the gradient Richardson's number and stability functions:
+!===================================================================
+
+! 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
+        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
+        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
+        N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay)
+        shear2(igrid,ilay)= (((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 )
+        Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10) 
+        
+        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 ! stable cases
+            Sm(igrid,ilay) = max(smmin,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
+        
+        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
+
+    ENDDO
+ENDDO
+
+
+! Computing the mixing length:
+!==============================================================
+
+
+IF (iflag_atke_lmix .EQ. 1 ) THEN
+
+   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))/sqrt(N2(igrid,ilay))
+             lstrat=max(lstrat,lmin)
+             !Inverse quadratic interpolation, Van de Wiel et al. 2010
+             l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
+          ENDIF 
+      ENDDO
+   ENDDO
+
+ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
+! 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. .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)
+             !Inverse quadratic interpolation, Van de Wiel et al. 2010   
+             l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
+          ENDIF
+      ENDDO
+   ENDDO
+
+
+
+ELSE
+! default: neglect effect of local stratification and shear
+
+   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
+ENDIF
+
+
+! Computing the TKE k>=2:
+!========================
+IF (iflag_atke == 0) THEN
+
+! stationary solution (dtke/dt=0)
+
+   DO ilay=2,nlay
+        DO igrid=1,ngrid
+            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
+            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
+        ENDDO
+    ENDDO
+
+ELSE IF (iflag_atke == 1) THEN
+
+! full implicit scheme resolved with a second order polynomial equation
+
+    DO ilay=2,nlay
+        DO igrid=1,ngrid
+           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+           delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. &
+                 +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + &
+                 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) &
+                 *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)))
+           qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2.
+           qq=max(0.,qq)
+           tke(igrid,ilay)=0.5*(qq**2) 
+        ENDDO
+    ENDDO
+
+
+ELSE IF (iflag_atke == 2) THEN
+
+! semi implicit scheme when l does not depend on tke
+! positive-guaranteed if pr slope in stable condition >1
+
+   DO ilay=2,nlay
+        DO igrid=1,ngrid
+           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+           qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      & 
+               *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
+               /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) 
+           tke(igrid,ilay)=0.5*(qq**2) 
+        ENDDO
+    ENDDO
+
+
+ELSE IF (iflag_atke == 3) THEN
+! numerical resolution adapted from that in MAR (Deleersnijder 1992)
+! positively defined by construction
+
+    DO ilay=2,nlay
+        DO igrid=1,ngrid
+           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+           IF (Ri(igrid,ilay) .LT. 0.) THEN
+              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
+              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
+           ELSE
+              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &
+                      l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)
+              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)
+           ENDIF
+           qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)
+           tke(igrid,ilay)=0.5*(qq**2)
+        ENDDO
+    ENDDO
+
+ELSE IF (iflag_atke == 4) THEN
+! semi implicit scheme from Arpege (V. Masson methodology with 
+! Taylor expansion of the dissipation term)
+    DO ilay=2,nlay
+        DO igrid=1,ngrid
+           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+           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)
+        ENDDO
+    ENDDO
+
+
+ELSE
+   call abort_physic("atke_compute_km_kh", &
+        'numerical treatment of TKE not possible yet', 1)
+
+END IF
+
+! We impose a 0 tke at nlay+1
+!==============================
+
+DO igrid=1,ngrid
+ tke(igrid,nlay+1)=0.
+END DO
+
+
+! Calculation of surface TKE (k=1)
+!=================================
+! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
+DO igrid=1,ngrid
+ ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
+ tke(igrid,1)=ctkes*(ustar**2)
+END DO
+
+
+! vertical diffusion of TKE 
+!==========================
+IF (atke_ok_vdiff) THEN
+   CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
+ENDIF
+
+
+! Computing eddy diffusivity coefficients:
+!========================================
+DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
+    DO igrid=1,ngrid
+        ! we add the molecular viscosity to Km,h 
+        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
+        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
+    END DO
+END DO
+
+! for output:
+!===========
+Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
+Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
+
+end subroutine atke_compute_km_kh
+
+!===============================================================================================
+subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
+
+! routine that computes the vertical diffusion of TKE by the turbulence
+! using an implicit resolution (See note by Dufresne and Ghattas (2009))
+! E Vignon, July 2023
+
+USE lmdz_atke_turbulence_ini, ONLY : rd, cke, viscom
+
+
+INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
+INTEGER, INTENT(IN) :: nlay  ! number of vertical index  
+
+REAL, INTENT(IN)    :: dtime ! physics time step (s)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: z_lay   ! altitude of mid-layers (m)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)       :: z_interf   ! altitude of bottom interfaces (m)
+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)     :: l_exchange     ! mixing length at interfaces between layers
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: Sm     ! stability function for eddy diffusivity for momentum at interface between layers
+
+REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke    ! turbulent kinetic energy at interface between layers
+
+
+
+INTEGER                                       :: igrid,ilay
+REAL, DIMENSION(ngrid,nlay+1)                 :: Ke     ! eddy diffusivity for TKE
+REAL, DIMENSION(ngrid,nlay+1)                 :: dtke
+REAL, DIMENSION(ngrid,nlay+1)                 :: ak, bk, ck, CCK, DDK
+REAL                                          :: gammak,Kem,KKb,KKt
+
+
+! Few initialisations
+CCK(:,:)=0.
+DDK(:,:)=0.
+dtke(:,:)=0.
+
+
+! Eddy diffusivity for TKE
+
+DO ilay=2,nlay
+    DO igrid=1,ngrid
+       Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke 
+    ENDDO
+ENDDO
+! at the top of the atmosphere set to 0
+Ke(:,nlay+1)=0.
+! at the surface, set it equal to that at the first model level
+Ke(:,1)=Ke(:,2)
+
+
+! calculate intermediary variables
+
+DO ilay=2,nlay 
+    DO igrid=1,ngrid
+    Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay))    
+    KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay))
+    Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1))
+    KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1))
+    gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1))
+    ak(igrid,ilay)=-gammak*dtime*KKb
+    ck(igrid,ilay)=-gammak*dtime*KKt
+    bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb)
+    ENDDO
+ENDDO
+
+! calculate CCK and DDK coefficients
+! downhill phase
+
+DO igrid=1,ngrid
+  CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
+  DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
+ENDDO
+
+
+DO ilay=nlay-1,2,-1
+    DO igrid=1,ngrid
+        CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
+                       / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
+        DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
+    ENDDO
+ENDDO
+
+! calculate TKE
+! uphill phase
+
+DO ilay=2,nlay+1
+    DO igrid=1,ngrid
+        dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay)
+    ENDDO
+ENDDO
+
+! update TKE
+tke(:,:)=tke(:,:)+dtke(:,:)
+
+
+end subroutine atke_vdiff_tke
+
+
+
+end module lmdz_atke_exchange_coeff
Index: LMDZ6/trunk/libf/phylmd/lmdz_atke_turbulence_ini.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_atke_turbulence_ini.F90	(revision 4687)
+++ LMDZ6/trunk/libf/phylmd/lmdz_atke_turbulence_ini.F90	(revision 4687)
@@ -0,0 +1,141 @@
+MODULE lmdz_atke_turbulence_ini
+
+implicit none
+
+save
+
+integer :: iflag_atke, iflag_num_atke, iflag_atke_lmix
+  !$OMP THREADPRIVATE(iflag_atke, iflag_num_atke, iflag_atke_lmix)
+  real :: kappa = 0.4 ! Von Karman constant
+  !$OMP THREADPRIVATE(kappa)
+  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)
+  real :: rg, rd, rpi, rcpd, rv
+  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv)
+
+  real :: viscom, viscoh
+  !$OMP THREADPRIVATE(viscom,viscoh)
+
+  real :: lmin=0.01              ! minimum mixing length
+  !$OMP THREADPRIVATE(lmin)
+
+  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, 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, rv_in
+
+
+  lunout=lunout_in
+  prt_level=prt_level_in
+  rd=rd_in
+  rg=rg_in
+  rpi=rpi_in
+  rcpd=rcpd_in
+  rv=rv_in
+
+  viscom=1.46E-5
+  viscoh=2.06E-5
+
+  ! flag that controls options in atke_compute_km_kh
+  iflag_atke=0
+  CALL getin_p('iflag_atke',iflag_atke)
+
+  ! flag that controls the calculation of mixing length in atke
+  iflag_atke_lmix=0
+  CALL getin_p('iflag_atke_lmix',iflag_atke_lmix)
+
+  if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then
+        call abort_physic("atke_turbulence_ini", &
+        'stationary scheme must use mixing length formulation not depending on tke', 1)
+  endif
+
+  ! activate vertical diffusion of TKE or not
+  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
+  iflag_num_atke=0
+  CALL getin_p('iflag_num_atke',iflag_num_atke)
+
+  ! asymptotic mixing length in neutral conditions [m]
+  ! Sun et al 2011, JAMC 
+  ! between 10 and 40
+
+  l0=15.0 
+  CALL getin_p('atke_l0',l0)
+
+  ! critical Richardson number 
+  ric=0.25
+  CALL getin_p('atke_ric',ric)
+
+  ! asymptotic value of Sm for Ri=-Inf
+  cinf=1.5
+  CALL getin_p('atke_cinf',cinf)
+
+  ! constant for tke dissipation calculation
+  cepsilon=5.87 ! default value as in yamada4
+  CALL getin_p('atke_cepsilon',cepsilon)
+
+
+ ! coefficient for surface TKE 
+ ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3)
+ ! (provided by limit condition in neutral conditions)
+  ctkes=(cepsilon)**(2./3.)
+
+  ! 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)
+  if (pr_slope .le. 1) then
+        call abort_physic("atke_turbulence_ini", &
+        'pr_slope has to be greater than 1 for consistency of the tke scheme', 1)
+  endif
+
+  ! asymptotic turbulent prandt number value for Ri=-Inf
+  pr_asym=0.4
+  CALL getin_p('atke_pr_asym',pr_asym)
+
+  ! value of turbulent prandtl number in neutral conditions (Ri=0)
+  pr_neut=0.8
+  CALL getin_p('atke_pr_neut',pr_neut)
+
+  ! coefficient for mixing length depending on local stratification
+  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
+  
+  smmin=0.17
+  CALL getin_p('atke_smmin',smmin)
+
+
+  ! ratio between the eddy diffusivity coeff for tke wrt that for momentum
+  ! default value from Lenderink et al. 2004
+  cke=2.
+  CALL getin_p('atke_cke',cke)
+
+ RETURN
+
+END SUBROUTINE atke_ini
+
+END MODULE  lmdz_atke_turbulence_ini
Index: LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90	(revision 4687)
+++ LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90	(revision 4687)
@@ -0,0 +1,169 @@
+module lmdz_call_atke
+
+USE lmdz_atke_exchange_coeff, 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,qvap,play,pinterf, &
+                        tke,Km_out,Kh_out)
+
+
+
+
+USE lmdz_atke_turbulence_ini, 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)       :: 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)
+
+
+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,dtime,&
+                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv,&
+                        tke,Km_out,Kh_out)
+
+
+                               
+if (iflag_num_atke .EQ. 1) then
+   !! In this case, we make an explicit prediction of the wind shear to calculate the tke in a
+   !! forward backward way
+   !! 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,dtime,&
+                        wind_u_predict,wind_v_predict,temp,qvap,play,pinterf,cdrag_uv, &
+                        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 lmdz_call_atke
Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4686)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4687)
@@ -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 call_atke_mod,       ONLY :  call_atke
+  USE lmdz_call_atke,      ONLY : call_atke
   USE ioipsl_getin_p_mod,  ONLY : getin_p
   USE cdrag_mod
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4686)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4687)
@@ -83,5 +83,5 @@
     USE lmdz_wake_ini, ONLY : wake_ini
     USE yamada_ini_mod, ONLY : yamada_ini
-    USE atke_turbulence_ini_mod, ONLY : atke_ini
+    USE lmdz_atke_turbulence_ini, ONLY : atke_ini
     USE lmdz_thermcell_ini, ONLY : thermcell_ini
     USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
