Index: /trunk/LMDZ.VENUS/libf/phyvenus/clmain.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/clmain.F	(revision 2534)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/clmain.F	(revision 2535)
@@ -12,4 +12,5 @@
      .                  solsw, sollw, sollwdown, fder,
      .                  rlon, rlat, dx, dy, 
+     .                  q2,
      .                  debut, lafin, 
      .                  d_t,d_u,d_v,d_ts,
@@ -32,10 +33,15 @@
 cAA il faudra sortir ces memes champs en leur ajoutant une dimension, 
 cAA c'est a dire nbsrf (nbre de subsurface).
-      USE ioipsl
-      USE interface_surf
-      use dimphy
+!      USE ioipsl
+!      USE interface_surf
+      use dimphy, only: klon, klev
       use mod_grid_phy_lmdz, only: nbp_lev
       use cpdet_phy_mod, only: t2tpot
       use turb_mod, only :yustar
+      
+#ifdef CPP_XIOS      
+      use xios_output_mod, only: send_xios_field
+#endif      
+      
       IMPLICIT none
 c======================================================================
@@ -56,4 +62,6 @@
 c dy-----input-R- resolution des mailles en y (m)
 c
+c q2-----inoutput-R- $q^2$ TKE at inter-layers
+c
 c d_t------output-R- le changement pour "t"
 c d_u------output-R- le changement pour "u"
@@ -69,44 +77,50 @@
 c======================================================================
 c$$$ PB ajout pour soil
-#include "dimsoil.h"
-#include "iniprint.h"
-#include "clesphys.h"
-#include "compbl.h"
-c
-      REAL dtime
-      integer itap
-      REAL t(klon,klev)
-      REAL u(klon,klev), v(klon,klev)
-      REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
-! ADAPTATION GCM POUR CP(T)
-      real ppk(klon,klev)
-      REAL rlon(klon), rlat(klon), dx(klon), dy(klon)
-      REAL d_t(klon, klev)
-      REAL d_u(klon, klev), d_v(klon, klev)
-      REAL flux_t(klon,klev)
-      REAL dflux_t(klon)
-
-      REAL flux_u(klon,klev), flux_v(klon,klev)
-      REAL cdragh(klon), cdragm(klon)
-      real rmu0(klon)         ! cosinus de l'angle solaire zenithal
-      LOGICAL debut, lafin
-c
-      REAL ts(klon)
-      REAL d_ts(klon)
-      REAL albe(klon)
-C
-      REAL fder(klon)
-      REAL sollw(klon), solsw(klon), sollwdown(klon)
+      include "dimsoil.h"
+      include "iniprint.h"
+      include "clesphys.h"
+      include "compbl.h"
+c
+c Parametres d'entree
+      REAL, INTENT(IN) :: dtime
+      INTEGER, INTENT(IN) :: itap
+      REAL, INTENT(IN) :: t(klon,klev)
+      REAL, INTENT(IN) :: u(klon,klev), v(klon,klev)
+      REAL, INTENT(IN) :: paprs(klon,klev+1), pplay(klon,klev)
+! ADAPTATION GCM POUR CP(T)
+      REAL, INTENT(IN) :: ppk(klon,klev)
+      REAL, INTENT(IN) :: rlon(klon), rlat(klon), dx(klon), dy(klon)
+      REAL, INTENT(IN) :: rmu0(klon)         ! cosine of solar zenith angle
+      LOGICAL, INTENT(IN) :: debut, lafin
+      REAL, INTENT(IN) :: ts(klon)
+      REAL, INTENT(IN) :: sollw(klon), solsw(klon), sollwdown(klon)
+      
+c Paramètres IN/OUT
+      REAL, INTENT(INOUT) :: fder(klon)
+      REAL, INTENT(INOUT) :: flux_u(klon,klev), flux_v(klon,klev)
+      REAL, INTENT(INOUT) :: radsol(klon)
+      REAL, INTENT(INOUT) :: q2(klon,klev+1)
+      
+c Parametres de sorties
+      REAL, INTENT(OUT) :: d_t(klon, klev)
+      REAL, INTENT(OUT) :: d_u(klon, klev), d_v(klon, klev)
+      REAL, INTENT(OUT) :: flux_t(klon,klev)
+      REAL, INTENT(OUT) :: dflux_t(klon)
+      REAL, INTENT(OUT) :: cdragh(klon), cdragm(klon)
+      REAL, INTENT(OUT) :: d_ts(klon)
+      REAL, INTENT(OUT) :: albe(klon)
 cAA
-      REAL zcoefh(klon,klev)
-      REAL zu1(klon)
-      REAL zv1(klon)
+      REAL, INTENT(OUT) :: zcoefh(klon,klev)
+      REAL, INTENT(OUT) :: zu1(klon)
+      REAL, INTENT(OUT) :: zv1(klon)
 cAA
 c$$$ PB ajout pour soil
-      REAL ftsoil(klon,nsoilmx)
+      REAL, INTENT(INOUT) :: ftsoil(klon,nsoilmx)
+
+c======================================================================
+      EXTERNAL clqh, clvent, coefkz
+c======================================================================
+c Parametre locaux
       REAL ytsoil(klon,nsoilmx)
-c======================================================================
-      EXTERNAL clqh, clvent, coefkz
-c======================================================================
       REAL yts(klon)
       REAL yalb(klon)
@@ -135,10 +149,9 @@
       real y_cd_m(klon),y_cd_h(klon)
 c
-#include "YOMCST.h"
+      include "YOMCST.h"
       REAL u1lay(klon), v1lay(klon)
       REAL delp(klon,klev)
       INTEGER i, k
-      INTEGER ni(klon), knon, j
-      
+      INTEGER ni(klon), knon, j     
 c======================================================================
       REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
@@ -200,4 +213,6 @@
       y_flux_u = 0.0
       y_flux_v = 0.0
+      ycoefh=0.0
+      ycoefm=0.0
 C$$ PB
       y_dflux_t = 0.0
@@ -360,5 +375,5 @@
             call yamada4(knon,dtime,rg,rd,ypaprs,yt
      s      ,yzlev,yzlay,yu,yv,yteta
-     s      ,y_cd_m,ykmm,ykmn,ykmq,yustar,
+     s      ,y_cd_m,q2,ykmm,ykmn,ykmq,yustar,
      s      iflag_pbl)
          endif
@@ -383,5 +398,5 @@
       CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
      s            y_d_v,y_flux_v)
-
+ 
 c pour le couplage
       ytaux = y_flux_u(:,1)
@@ -447,15 +462,4 @@
       ENDDO
 
-c --------------------
-c TEST!!!!! PAS DE MELANGE PAR TURBULENCE !!!
-c       d_u = 0. 
-c       d_v = 0.
-c       flux_u = 0.
-c       flux_v = 0.
-c --------------------
-
-c     print*,"y_d_t apres clqh=",y_d_t(klon/2,:)
-
-      RETURN
       END
 
@@ -478,6 +482,6 @@
      s                flux_t, dflux_s)
 
-      USE interface_surf
-      use dimphy
+      use interface_surf, only: interfsurf_hq
+      use dimphy, only: klon, klev
       use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
       use cpdet_phy_mod, only: t2tpot,tpot2t,cpdet
@@ -488,36 +492,50 @@
 c Objet: diffusion verticale de "h"
 c======================================================================
-#include "YOMCST.h"
-#include "dimsoil.h"
-#include "iniprint.h"
+      include "YOMCST.h"
+      include "dimsoil.h"
+      include "iniprint.h"
 
 c Arguments:
-      INTEGER knon
-      REAL dtime              ! intervalle du temps (s)
-      REAL u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
-      REAL v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
-      REAL coef(klon,klev)    ! le coefficient d'echange (m**2/s)
+c Parametres d'entree
+      INTEGER, INTENT(IN) :: knon
+      REAL, INTENT(IN) :: dtime              ! intervalle du temps (s)
+      REAL, INTENT(IN) :: u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
+      REAL, INTENT(IN) :: v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
+      REAL, INTENT(IN) :: coef(klon,klev)    ! le coefficient d'echange (m**2/s)
 c                               multiplie par le cisaillement du 
 c                               vent (dV/dz); la premiere valeur
 c                               indique la valeur de Cdrag (sans unite)
-      REAL t(klon,klev)       ! temperature (K)
-      REAL ts(klon)           ! temperature du sol (K)
-      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
-      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
-! ADAPTATION GCM POUR CP(T)
-      REAL ppk(klon,klev)     ! fonction d'Exner milieu de couche
-      REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
-      REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
-      REAL albedo(klon)       ! albedo de la surface
-      real rmu0(klon)         ! cosinus de l'angle solaire zenithal
-      real rlon(klon), rlat(klon), dx(klon), dy(klon)
-c
-      REAL d_t(klon,klev)     ! incrementation de "t"
-      REAL d_ts(klon)         ! incrementation de "ts"
-      REAL flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
+      REAL, INTENT(IN) :: t(klon,klev)       ! temperature (K)
+      REAL, INTENT(IN) :: ts(klon)           ! temperature du sol (K)
+      REAL, INTENT(IN) :: paprs(klon,klev+1) ! pression a inter-couche (Pa)
+      REAL, INTENT(IN) :: pplay(klon,klev)   ! pression au milieu de couche (Pa)
+! ADAPTATION GCM POUR CP(T)
+      REAL, INTENT(IN) :: ppk(klon,klev)     ! fonction d'Exner milieu de couche
+      REAL, INTENT(IN) :: delp(klon,klev)    ! epaisseur de couche en pression (Pa)
+      REAL, INTENT(INOUT) :: radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
+      REAL, INTENT(IN) :: rmu0(klon)         ! cosinus de l'angle solaire zenithal
+      REAL, INTENT(IN) :: rlon(klon), rlat(klon), dx(klon), dy(klon)
+      INTEGER, INTENT(IN) :: itime
+      LOGICAL, INTENT(IN) :: debut, lafin
+      REAL, INTENT(INOUT) :: fder(klon)
+      REAL, INTENT(IN) :: taux(klon), tauy(klon)
+      REAL, INTENT(IN) :: sollw(klon), sollwdown(klon)
+      REAL, INTENT(IN) :: swnet(klon)
+      
+      
+c$$$C PB ajout pour soil
+      LOGICAL, INTENT(IN) :: soil_model
+      REAL, INTENT(IN) :: tsoil(klon, nsoilmx)
+c
+c Parametre de sorties
+      REAL, INTENT(OUT) :: albedo(klon)       ! albedo de la surface
+      REAL, INTENT(OUT) :: d_t(klon,klev)     ! incrementation de "t"
+      REAL, INTENT(OUT) :: d_ts(klon)         ! incrementation de "ts"
+      REAL, INTENT(OUT) :: flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
 c                               sensible, flux de Cp*T, positif vers
 c                               le bas: j/(m**2 s) c.a.d.: W/m2
-      REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
-c======================================================================
+      REAL, INTENT(OUT) :: dflux_s(klon) ! derivee du flux sensible dF/dTs
+c======================================================================
+c Variables locales
       INTEGER i, k
       REAL zx_ch(klon,klev)
@@ -527,4 +545,5 @@
       REAL local_h(klon,klev) ! enthalpie potentielle
       REAL local_ts(klon)
+      REAL swdown(klon)
 c======================================================================
 c contre-gradient pour la chaleur sensible: Kelvin/metre
@@ -534,20 +553,13 @@
       REAL zdelz
 c======================================================================
-#include "compbl.h"
+      include "compbl.h"
 c======================================================================
 c Rajout pour l'interface
-      integer itime
-      logical debut, lafin
       real zlev1(klon)
-      real fder(klon), taux(klon), tauy(klon)
       real temp_air(klon)
       real epot_air(klon)
       real tq_cdrag(klon), petAcoef(klon)
       real petBcoef(klon)
-      real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
-      real p1lay(klon),pkh1(klon)
-c$$$C PB ajout pour soil
-      LOGICAL soil_model
-      REAL tsoil(klon, nsoilmx)
+      real p1lay(klon), pkh1(klon)
 
 ! Parametres de sortie
@@ -777,5 +789,4 @@
 c
 
-      RETURN
       END
       
@@ -791,5 +802,5 @@
      s                  d_ven,flux_v)
 
-      use dimphy
+      use dimphy, only: klon, klev
       IMPLICIT none
 c======================================================================
@@ -814,16 +825,22 @@
 c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
 c======================================================================
-#include "iniprint.h"
-      INTEGER knon
-      REAL dtime
-      REAL u1lay(klon), v1lay(klon)
-      REAL coef(klon,klev)
-      REAL t(klon,klev), ven(klon,klev)
-      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
-      REAL d_ven(klon,klev)
-      REAL flux_v(klon,klev)
-c======================================================================
-#include "YOMCST.h"
-c======================================================================
+      include "iniprint.h"
+
+c Parametres d'entree
+      INTEGER, INTENT(IN) :: knon
+      REAL, INTENT(IN) :: dtime
+      REAL, INTENT(IN) :: u1lay(klon) , v1lay(klon)
+      REAL, INTENT(IN) :: coef(klon, klev)
+      REAL, INTENT(IN) :: t(klon, klev), ven(klon, klev)
+      REAL, INTENT(IN) :: paprs(klon, klev+1), pplay(klon, klev)
+      REAL, INTENT(IN) :: delp(klon, klev)
+      
+c Parametres de sorties
+      REAL, INTENT(OUT) :: d_ven(klon, klev)
+      REAL, INTENT(OUT) :: flux_v(klon, klev)
+c======================================================================
+      include "YOMCST.h"
+c======================================================================
+c Parametres locaux
       INTEGER i, k
       REAL zx_cv(klon,2:klev)
@@ -905,5 +922,4 @@
       ENDDO
 c
-      RETURN
       END
       
@@ -919,5 +935,5 @@
      .                  pcfm, pcfh)
 
-      use dimphy
+      use dimphy, only: klon, klev
       use cpdet_phy_mod, only: cpdet,t2tpot
 #ifdef CPP_XIOS      
@@ -944,47 +960,63 @@
 c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
 c======================================================================
-#include "YOMCST.h"
-#include "iniprint.h"
-#include "compbl.h"
-#include "clesphys.h"
+      include "YOMCST.h"
+      include "iniprint.h"
+      include "compbl.h"
+      include "clesphys.h"
 c
 c Arguments:
 c
-      INTEGER knon
-      REAL ts(klon)
-      REAL paprs(klon,klev+1), pplay(klon,klev)
-! ADAPTATION GCM POUR CP(T)
-      real ppk(klon,klev)
-      REAL u(klon,klev), v(klon,klev), t(klon,klev)
-c
-      REAL pcfm(klon,klev), pcfh(klon,klev)
-      INTEGER itop(klon)
+c Parametres d'entree
+      INTEGER, INTENT(IN) :: knon
+      REAL, INTENT(IN) :: ts(klon)
+      REAL, INTENT(IN) :: pplay(klon, klev)
+      REAL, INTENT(IN) ::paprs(klon, klev+1)
+! ADAPTATION GCM POUR CP(T)
+      REAL, INTENT(IN) :: ppk(klon, klev)
+      REAL, INTENT(IN) :: u(klon, klev), v(klon, klev), t(klon, klev)
+      
+c Parametre de sorties
+      REAL, INTENT(OUT) :: pcfm(klon, klev), pcfh(klon, klev)
 c
 c Quelques constantes et options:
 c
-      REAL cepdu2, ckap, cb, cc, cd, clam
+      REAL, PARAMETER :: cepdu2=((1.e-5)**2)
+c     REAL, PARAMETER :: cepdu2 =(0.1)**2
+      REAL, PARAMETER :: ckap=0.4
+      REAL, PARAMETER :: cb=5.0
+      REAL, PARAMETER :: cc=5.0
+      REAL, PARAMETER :: cd=5.0
+      REAL, PARAMETER :: clam=160.0
+      REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique
+      REAL, PARAMETER :: prandtl=0.4
+      INTEGER isommet ! le sommet de la couche limite
+      
+      LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante
+      LOGICAL, PARAMETER :: opt_ec=.FALSE. ! formule du Centre Europeen dans l'atmosphere
+      
+c      REAL cepdu2, ckap, cb, cc, cd, clam
 c TEST VENUS
 c     PARAMETER (cepdu2 =(0.1)**2)
-      PARAMETER (cepdu2 =(1.e-5)**2)
-
-      PARAMETER (CKAP=0.4)
-      PARAMETER (cb=5.0)
-      PARAMETER (cc=5.0)
-      PARAMETER (cd=5.0)
-      PARAMETER (clam=160.0)
-      REAL ric ! nombre de Richardson critique
-      PARAMETER(ric=0.4)
-      REAL prandtl
-      PARAMETER (prandtl=0.4)
-      INTEGER isommet ! le sommet de la couche limite
-
-      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
-      PARAMETER (tvirtu=.TRUE.)
-      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
-      PARAMETER (opt_ec=.FALSE.)
+c      PARAMETER (cepdu2 =(1.e-5)**2)
+c      PARAMETER (CKAP=0.4)
+c      PARAMETER (cb=5.0)
+c      PARAMETER (cc=5.0)
+c      PARAMETER (cd=5.0)
+c      PARAMETER (clam=160.0)
+c      REAL ric ! nombre de Richardson critique
+c      PARAMETER(ric=0.4)
+c      REAL prandtl
+c      PARAMETER (prandtl=0.4)
+c      INTEGER isommet ! le sommet de la couche limite
+
+c      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
+c      PARAMETER (tvirtu=.TRUE.)
+c      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
+c      PARAMETER (opt_ec=.FALSE.)
 
 c
 c Variables locales:
 c
+      INTEGER itop(klon)
       INTEGER i, k
       REAL zgeop(klon,klev)
@@ -1001,6 +1033,5 @@
       REAL z2geomf, zalh2, zalm2, zscfh, zscfm
 cIM
-      LOGICAL check
-      PARAMETER (check=.false.)
+      LOGICAL, PARAMETER :: check=.false.
 c
 c contre-gradient pour la chaleur sensible: Kelvin/metre
@@ -1008,16 +1039,12 @@
       REAL gamh(2:klev)
 c
-      LOGICAL appel1er
-      SAVE appel1er
+      LOGICAL, SAVE :: appel1er = .TRUE.
 c
 c Fonctions thermodynamiques et fonctions d'instabilite
       REAL fsta, fins, x
-      LOGICAL zxli ! utiliser un jeu de fonctions simples
-      PARAMETER (zxli=.FALSE.)
+      LOGICAL, PARAMETER :: zxli=.FALSE. ! utiliser un jeu de fonctions simples PARAMETER (zxli=.FALSE.)
 c
       fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
       fins(x) = SQRT(1.0-18.0*x)
-c
-      DATA appel1er /.TRUE./
 c
 c----------------------
@@ -1228,5 +1255,4 @@
 c#endif
 
-      RETURN
       END
 
@@ -1239,5 +1265,5 @@
      .                  pcfm, pcfh)
 
-      use dimphy
+      use dimphy, only: klon, klev
       use mod_grid_phy_lmdz, only: nbp_lev
       use cpdet_phy_mod, only: cpdet
@@ -1257,14 +1283,15 @@
 c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
 c======================================================================
-#include "YOMCST.h"
-#include "iniprint.h"
+      include "YOMCST.h"
+      include "iniprint.h"
 c
 c Arguments:
 c
-      INTEGER knon
-      REAL paprs(klon,klev+1), pplay(klon,klev)
-      REAL t(klon,klev)
-c
-      REAL pcfm(klon,klev), pcfh(klon,klev)
+      INTEGER,INTENT(IN) :: knon
+      REAL,INTENT(IN) :: paprs(klon,klev+1)
+      REAL,INTENT(IN) :: pplay(klon,klev)
+      REAL,INTENT(IN) :: t(klon,klev)
+c
+      REAL,INTENT(OUT) :: pcfm(klon,klev), pcfh(klon,klev)
 c
 c Variables locales:
@@ -1304,5 +1331,4 @@
       ENDDO
 c
-      RETURN
       END
 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 2534)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 2535)
@@ -1277,5 +1277,6 @@
      $            paprs,pplay,ppk,radsol,falbe,
      e            solsw, sollw, sollwdown, fder,
-     e            longitude_deg, latitude_deg, dx, dy,   
+     e            longitude_deg, latitude_deg, dx, dy,
+     &            q2, 
      e            debut, lafin,
      s            d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
Index: /trunk/LMDZ.VENUS/libf/phyvenus/yamada4.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/yamada4.F	(revision 2534)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/yamada4.F	(revision 2535)
@@ -3,10 +3,14 @@
 !
       SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
-c     s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
-     s   ,zlev,zlay,u,v,teta,cd,km,kn,kq,ustar
+     s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
      s   ,iflag_pbl)
 c.......................................................................
-      use dimphy
-      use turb_mod, only: q2,l0
+      use dimphy, only: klon, klev
+      use turb_mod, only: l0
+          
+#ifdef CPP_XIOS      
+      use xios_output_mod, only: send_xios_field
+#endif      
+    
       IMPLICIT NONE
 c.......................................................................
@@ -40,34 +44,36 @@
 
 c.......................................................................
-      REAL dt,g,rconst
-      real plev(klon,klev+1),temp(klon,klev)
-      real ustar(klon)
+      REAL, INTENT(IN) :: dt,g,rconst
+      REAL, INTENT(IN) :: plev(klon,klev+1),temp(klon,klev)
+      REAL, INTENT(IN) :: ustar(klon)
+      REAL, INTENT(INOUT) :: zlev(klon,klev+1)
+      REAL, INTENT(IN) :: zlay(klon,klev)
+      REAL, INTENT(IN) :: u(klon,klev)
+      REAL, INTENT(IN) :: v(klon,klev)
+      REAL, INTENT(IN) :: teta(klon,klev)
+      REAL, INTENT(IN) :: cd(klon)
+      INTEGER, INTENT(IN) :: iflag_pbl,ngrid
+      
+      REAL, INTENT(OUT) :: km(klon,klev+1)
+      REAL, INTENT(OUT) :: kn(klon,klev+1)
+      REAL, INTENT(OUT) :: kq(klon,klev+1)
+      
+      REAL, INTENT(INOUT) :: q2(klon,klev+1)
+c Variables locales:
       real kmin,qmin,pblhmin(klon),coriol(klon)
-      REAL zlev(klon,klev+1)
-      REAL zlay(klon,klev)
-      REAL u(klon,klev)
-      REAL v(klon,klev)
-      REAL teta(klon,klev)
-      REAL cd(klon)
       REAL qpre
       REAL unsdz(klon,klev)
       REAL unsdzdec(klon,klev+1)
 
-      REAL km(klon,klev+1)
       REAL kmpre(klon,klev+1),tmp2
       REAL mpre(klon,klev+1)
-      REAL kn(klon,klev+1)
-      REAL kq(klon,klev+1)
+      
       real ff(klon,klev+1),delta(klon,klev+1)
       real aa(klon,klev+1),aa0,aa1
-      integer iflag_pbl,ngrid
-
-
+      
       integer nlay,nlev
 
-      logical first
-      integer ipas
-      save first,ipas
-      data first,ipas/.true.,0/
+      logical,save :: first=.false. ! not neede any more
+      integer,save :: ipas=0
 
 
@@ -82,8 +88,4 @@
       real m2cstat,mcstat,kmcstat
       real l(klon,klev+1)
-      !real,save,allocatable :: l0(:)
-c  ATTENTION! mis ici car j'ai enlevé q2 des arguments...
-c   sinon, c'est au-dessus que ça se passe...
-      !REAL,save,allocatable :: q2(:,:)
 
       real sq(klon),sqz(klon),zz(klon,klev+1)
@@ -97,5 +99,5 @@
       real frif,falpha,fsm
       real fl,zzz,zl0,zq2,zn2
-
+      
 c     real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
 c    s  ,lyam(klon,klev),knyam(klon,klev)
@@ -109,4 +111,5 @@
      s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
 
+   
       if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
            stop'probleme de coherence dans appel a MY'
@@ -118,80 +121,72 @@
       nlev=klev+1
 
-      if (first) then
-!        IF (.not.ALLOCATED(l0)) allocate(l0(klon))
-!        IF (.not.ALLOCATED(q2)) allocate(q2(klon,klevp1))
-
-c (surtout pour k=1, à cause diagnostiques...)
-        q2 = 0.
-      endif
+c      if (first) then
+c        IF (.not.ALLOCATED(l0)) allocate(l0(klon))
+c        IF (.not.ALLOCATED(q2)) allocate(q2(klon,klev+1))
+c      endif
 c===================================
-      
+       
       ipas=ipas+1
-      if (0.eq.1.and.first) then
-      do ig=1,1000
-         ri=(ig-800.)/500.
-         if (ri.lt.ric) then
-            zrif=frif(ri)
-         else
-            zrif=rifc
-         endif
-         if(zrif.lt.0.16) then
-            zalpha=falpha(zrif)
-            zsm=fsm(zrif)
-         else
-            zalpha=1.12
-            zsm=0.085
-         endif
+!      if (0.eq.1.and.first) then
+!      do ig=1,1000
+!         ri=(ig-800.)/500.
+!         if (ri.lt.ric) then
+!            zrif=frif(ri)
+!         else
+!            zrif=rifc
+!         endif
+!         if(zrif.lt.0.16) then
+!            zalpha=falpha(zrif)
+!            zsm=fsm(zrif)
+!         else
+!            zalpha=1.12
+!            zsm=0.085
+!         endif
 c     print*,ri,rif,zalpha,zsm
-      enddo
-      endif
-
+c      enddo
+c      endif
+      
 c.......................................................................
 c  les increments verticaux
 c.......................................................................
 c
-c!!!!! allerte !!!!!c
-c!!!!! zlev n'est pas declare a nlev !!!!!c
-c!!!!! ---->
-                                                      DO ig=1,ngrid
-            zlev(ig,nlev)=zlay(ig,nlay)
-     &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )
-                                                      ENDDO
-c!!!!! <----
-c!!!!! allerte !!!!!c
-c
       DO k=1,nlay
-                                                      DO ig=1,ngrid
+       DO ig=1,ngrid
         unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
-                                                      ENDDO
+       ENDDO
       ENDDO
-                                                      DO ig=1,ngrid
+      DO ig=1,ngrid
       unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
-                                                      ENDDO
+      ENDDO
       DO k=2,nlay
-                                                      DO ig=1,ngrid
+       DO ig=1,ngrid
         unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
-                                                     ENDDO
+       ENDDO
       ENDDO
-                                                      DO ig=1,ngrid
+      DO ig=1,ngrid
       unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
-                                                     ENDDO
-c
-c.......................................................................
-
+      ENDDO
+c
+c.......................................................................
+ 
 c===================================
 c INITIALISATIONS (surtout pour k=1, à cause diagnostiques...)
-	dz = 0.
-	m2 = 0.
-	dtetadz = 0.
-	n2 = 0.
-	rif = 0.
-	alpha = 0.
-	sm = 0.
-	zz = 0.
-	l = 0.
+	dz(:,:) = 0.
+	m2(:,:) = 0.
+	dtetadz(:,:) = 0.
+	n2(:,:) = 0.
+	rif(:,:) = 0.
+	alpha(:,:) = 0.
+	sm(:,:) = 0.
+	zz(:,:) = 0.
+	l(:,:) = 0.
+        km(:,:) = 0.
+        kn(:,:) = 0.
+        kmpre(:,:) = 0.
+        mpre(:,:) = 0.
+
 c===================================
       do k=2,klev
-                                                          do ig=1,ngrid
+       do ig=1,ngrid
          dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
          m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
@@ -216,8 +211,6 @@
 c     print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)
 
-
-                                                          enddo
-      enddo
-
+       enddo
+      enddo
 
 c====================================================================
@@ -259,7 +252,6 @@
 
 c     print*,'Fin de l initialisation de q2 et l0'
-
-      endif ! first
-
+      endif ! of if (first.or.iflag_pbl.eq.6)
+    
 c====================================================================
 c  Calcul de la longueur de melange.
@@ -292,5 +284,4 @@
                                                           enddo
       enddo
-
 c====================================================================
 c   Yamada 2.0
@@ -402,7 +393,4 @@
 
       endif ! Fin du cas 8
-
-c     print*,'OK8'
-
 c====================================================================
 c   Calcul des coefficients de mélange
@@ -466,6 +454,4 @@
       endif
 
-c     print*,'YAMADA4 1'
-
 c   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
 
@@ -481,4 +467,4 @@
 c     print*,'OKFIN'
       first=.false.
-      return
+      
       end
