Index: /LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90	(revision 4653)
+++ /LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90	(revision 4654)
@@ -167,8 +167,10 @@
 !   iflag_pbl peut etre utilise comme longuer de melange
        IF (iflag_pbl.GE.31) THEN
-          CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, &
+          if (knon>0) then
+          CALL vdif_kcay(knon,klev,dtime,RG,RD,ypaprs,yt, &
                yzlev,yzlay,yu,yv,yteta, &
                ycdragm,yq2,q2diag,ykmm,ykmn,yustar, &
                iflag_pbl)
+          endif
        ELSE IF (iflag_pbl<20) THEN
           CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, &
Index: /LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90	(revision 4653)
+++ /LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90	(revision 4654)
@@ -1,5 +1,5 @@
-module lscp_ini_mod
-
-implicit none
+MODULE lscp_ini_mod
+
+IMPLICIT NONE
 
   ! PARAMETERS for lscp:
@@ -259,5 +259,7 @@
     CALL cloudth_ini(iflag_cloudth_vert,iflag_ratqs)
 
-end subroutine lscp_ini
-
-end module lscp_ini_mod
+RETURN
+
+END SUBROUTINE lscp_ini
+
+END MODULE lscp_ini_mod
Index: /LMDZ6/trunk/libf/phylmd/lscp_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4653)
+++ /LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4654)
@@ -7,5 +7,5 @@
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
-     paprs,pplay,t,q,ptconv,ratqs,                      &
+     paprs,pplay,temp,q,ptconv,ratqs,                      &
      d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri,  &
      pfraclr,pfracld,                                   &
@@ -91,4 +91,6 @@
 !
 USE print_control_mod, ONLY: prt_level, lunout
+
+! USE de modules contenant des fonctions.
 USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc 
 USE lscp_tools_mod, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat
@@ -96,4 +98,5 @@
 USE ice_sursat_mod, ONLY : ice_sursat
 
+! Use du module lscp_ini_mod contenant les constantes
 USE lscp_ini_mod, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
 USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
@@ -119,5 +122,5 @@
   REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t               ! temperature (K)
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q               ! total specific humidity (= vapor phase) [kg/kg] 
   INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
@@ -194,68 +197,75 @@
   !----------------
 
-  REAL qsl(klon), qsi(klon)
-  REAL zct, zcl,zexpo
-  REAL ctot(klon,klev)
-  REAL ctot_vol(klon,klev)
-  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
-  REAL zdqsdT_raw(klon)
-  REAL gammasat(klon),dgammasatdt(klon)                ! coefficient to make cold condensation at the correct RH and derivative wrt T
-  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
-  REAL cste
-  REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
-  REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
-  REAL erf
-  REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon)
-  REAL zrfl(klon), zrfln(klon), zqev, zqevt 
-  REAL zifl(klon), zifln(klon), ziflprev(klon),zqev0,zqevi, zqevti
-  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon)
-  REAL zoliql(klon), zoliqi(klon)
-  REAL zt(klon)
-  REAL zdz(klon),zrho(klon,klev),iwc(klon)
-  REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon)
-  REAL zmelt,zrain,zsnow,zprecip
-  REAL dzfice(klon)
-  REAL zsolid
-  REAL qtot(klon), qzero(klon)
-  REAL dqsl(klon),dqsi(klon)
-  REAL smallestreal
+  REAL,DIMENSION(klon) :: qsl, qsi
+  REAL :: zct, zcl,zexpo
+  REAL, DIMENSION(klon,klev) :: ctot
+  REAL, DIMENSION(klon,klev) :: ctot_vol
+  REAL, DIMENSION(klon) :: zqs, zdqs
+  REAL :: zdelta, zcor, zcvm5
+  REAL, DIMENSION(klon) :: zdqsdT_raw
+  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                ! coefficient to make cold condensation at the correct RH and derivative wrt T
+  REAL, DIMENSION(klon) :: Tbef,qlbef,DT
+  REAL :: num,denom
+  REAL :: cste
+  REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta
+  REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2
+  REAL :: erf
+  REAL, DIMENSION(klon,klev) :: icefrac_mpc
+  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
+  REAL, DIMENSION(klon) :: zrfl, zrfln
+  REAL :: zqev, zqevt 
+  REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
+  REAL :: zqev0,zqevi, zqevti
+  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
+  REAL, DIMENSION(klon) :: zoliql, zoliqi
+  REAL, DIMENSION(klon) :: zt
+  REAL, DIMENSION(klon,klev) :: zrho
+  REAL, DIMENSION(klon) :: zdz,iwc
+  REAL :: zchau,zfroi
+  REAL, DIMENSION(klon) :: zfice,zneb,znebprecip
+  REAL :: zmelt,zrain,zsnow,zprecip
+  REAL, DIMENSION(klon) :: dzfice
+  REAL :: zsolid
+  REAL, DIMENSION(klon) :: qtot, qzero
+  REAL, DIMENSION(klon) :: dqsl,dqsi
+  REAL :: smallestreal
   !  Variables for Bergeron process
-  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
-  REAL zqpreci(klon), zqprecl(klon)
+  REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
+  REAL, DIMENSION(klon) :: zqpreci, zqprecl
   ! Variables precipitation energy conservation
-  REAL zmqc(klon)
-  REAL zalpha_tr
-  REAL zfrac_lessi
-  REAL zprec_cond(klon)
-  REAL zmair(klon), zcpair, zcpeau
-  REAL zlh_solid(klon), zm_solid         ! for liquid -> solid conversion
-  REAL zrflclr(klon), zrflcld(klon)
-  REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)
-  REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon) 
-  REAL ziflclr(klon), ziflcld(klon)
-  REAL znebprecipclr(klon), znebprecipcld(klon)
-  REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)
-  REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
-  REAL velo(klon,klev), vr, ffallv
-  REAL qlmpc, qimpc, rnebmpc
-  REAL radocondi(klon,klev), radocondl(klon,klev)
-  REAL effective_zneb
-  REAL distcltop1D(klon), temp_cltop1D(klon)
-
-
-  INTEGER i, k, n, kk
-  INTEGER n_i(klon), iter
+  REAL, DIMENSION(klon) :: zmqc
+  REAL :: zalpha_tr
+  REAL :: zfrac_lessi
+  REAL, DIMENSION(klon) :: zprec_cond
+  REAL, DIMENSION(klon) :: zmair
+  REAL :: zcpair, zcpeau
+  REAL, DIMENSION(klon) :: zlh_solid
+  REAL :: zm_solid         ! for liquid -> solid conversion
+  REAL, DIMENSION(klon) :: zrflclr, zrflcld
+  REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
+  REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr 
+  REAL, DIMENSION(klon) :: ziflclr, ziflcld
+  REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
+  REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
+  REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
+  REAL, DIMENSION(klon,klev) :: velo
+  REAL :: vr, ffallv
+  REAL :: qlmpc, qimpc, rnebmpc
+  REAL, DIMENSION(klon,klev) :: radocondi, radocondl
+  REAL :: effective_zneb
+  REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D
+
+
+  INTEGER i, k, n, kk, iter
+  INTEGER, DIMENSION(klon) :: n_i
   INTEGER ncoreczq
-  INTEGER mpc_bl_points(klon,klev)
-
-  LOGICAL lognormale(klon)
-  LOGICAL keepgoing(klon)
-  LOGICAL,SAVE :: appel1er
-  !$OMP THREADPRIVATE(appel1er)
+  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
+
+  LOGICAL, DIMENSION(klon) :: lognormale
+  LOGICAL, DIMENSION(klon) :: keepgoing
 
   CHARACTER (len = 20) :: modname = 'lscp'
   CHARACTER (len = 80) :: abort_message
   
-  DATA appel1er /.TRUE./
 
 !===============================================================================
@@ -264,4 +274,5 @@
 
 ! Few initial checks
+
 
 IF (iflag_fisrtilp_qsat .LT. 0) THEN
@@ -332,4 +343,6 @@
 fcontrN(:,:) = 0.0
 fcontrP(:,:) = 0.0
+distcltop(:,:)=0.
+temp_cltop(:,:)=0.
 
 !c_iso: variable initialisation for iso 
@@ -346,5 +359,5 @@
     ! Initialisation temperature and specific humidity
     DO i = 1, klon
-        zt(i)=t(i,k)
+        zt(i)=temp(i,k)
         zq(i)=q(i,k)
         !c_iso init of iso
@@ -382,5 +395,5 @@
             zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
             ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
-            zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
+            zt(i) = ( (temp(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
              / (zcpair + zmqc(i)*zcpeau)
 
@@ -625,5 +638,5 @@
                          zq,zqta,fraca,                            &
                          qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
-                         ratqs,zqs,t,                              &
+                         ratqs,zqs,temp,                              &
                          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
 
@@ -637,5 +650,5 @@
                          zq,zqta,fraca,                                     &
                          qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
-                         ratqs,zqs,t, &
+                         ratqs,zqs,temp, &
                          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
 
@@ -648,5 +661,5 @@
                          zq,zqta,fraca,                                     &
                          qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
-                         ratqs,zqs,t, &
+                         ratqs,zqs,temp, &
                          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
 
@@ -658,5 +671,5 @@
 
                     CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, &
-                         t,ztv,zq,zqta,fraca, zpspsk,                      &
+                         temp,ztv,zq,zqta,fraca, zpspsk,                      &
                          paprs,pplay,ztla,zthl,ratqs,zqs,psfl,             &
                          qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol,    &
@@ -739,5 +752,5 @@
                   IF (iflag_t_glace.GE.4) THEN
                   ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
-                       CALL distance_to_cloud_top(klon,klev,k,t,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
+                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
                   ENDIF
                   CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:))
@@ -780,5 +793,5 @@
                         !------------------------------------
 
-                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), &
+                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), &
                              gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),                 &
                              rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                   & 
@@ -837,5 +850,5 @@
         ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
         IF (iflag_t_glace.GE.4) THEN
-            CALL distance_to_cloud_top(klon,klev,k,t,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
+            CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
             distcltop(:,k)=distcltop1D(:)
             temp_cltop(:,k)=temp_cltop1D(:)
@@ -1196,5 +1209,5 @@
         d_q(i,k) = zq(i) - q(i,k)
         ! c_iso: same for isotopes
-        d_t(i,k) = zt(i) - t(i,k)
+        d_t(i,k) = zt(i) - temp(i,k)
     ENDDO
 
@@ -1245,5 +1258,5 @@
         IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
 
-            IF (t(i,k) .GE. t_glace_min) THEN
+            IF (temp(i,k) .GE. t_glace_min) THEN
                 zalpha_tr = a_tr_sca(3)
             ELSE
@@ -1269,5 +1282,5 @@
             IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
 
-                IF (t(i,kk) .GE. t_glace_min) THEN
+                IF (temp(i,kk) .GE. t_glace_min) THEN
                     zalpha_tr = a_tr_sca(1)
                 ELSE
@@ -1318,4 +1331,7 @@
   ENDIF
 
+
+RETURN
+
 END SUBROUTINE lscp
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Index: /LMDZ6/trunk/libf/phylmd/physiqex_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiqex_mod.F90	(revision 4653)
+++ /LMDZ6/trunk/libf/phylmd/physiqex_mod.F90	(revision 4654)
@@ -9,18 +9,32 @@
      &            debut,lafin,pdtphys, &
      &            paprs,pplay,pphi,pphis,presnivs, &
-     &            u,v,rot,t,qx, &
+     &            u,v,rot,temp,qx, &
      &            flxmass_w, &
      &            d_u, d_v, d_t, d_qx, d_ps)
 
-      USE dimphy, only : klon,klev
-      USE infotrac_phy, only : nqtot
-      USE geometry_mod, only : latitude
-!      USE comcstphy, only : rg
-      USE ioipsl, only : ymds2ju
-      USE phys_state_var_mod, only : phys_state_var_init
-      USE phyetat0_mod, only: phyetat0
-      USE output_physiqex_mod, ONLY: output_physiqex
+
+USE dimphy, only : klon,klev
+USE infotrac_phy, only : nqtot
+USE geometry_mod, only : latitude
+USE ioipsl, only : ymds2ju
+USE phys_state_var_mod, only : phys_state_var_init
+USE phyetat0_mod, only: phyetat0
+USE output_physiqex_mod, ONLY: output_physiqex
+use vdif_ini, only : vdif_ini_
+USE lmdz_thermcell_ini, ONLY : thermcell_ini
+USE ioipsl_getin_p_mod, ONLY : getin_p
+USE wxios, ONLY: missing_val, using_xios
+USE lscp_mod, ONLY : lscp
+USE lscp_ini_mod, ONLY : lscp_ini
+USE add_phys_tend_mod, ONLY : add_phys_tend
+
 
       IMPLICIT none
+
+include "YOETHF.h"
+
+
+
+
 !
 ! Routine argument:
@@ -32,20 +46,110 @@
       logical,intent(in) :: lafin ! signals last call to physics
       real,intent(in) :: pdtphys ! physics time step (s)
-      real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
-      real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
-      real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer
-      real,intent(in) :: pphis(klon) ! surface geopotential
-      real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers
-      real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s)
-      real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s)
-      real,intent(in) :: rot(klon,klev) ! northward meridional wind (m/s)
-      real,intent(in) :: t(klon,klev) ! temperature (K)
-      real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air)
-      real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux
-      real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s)
-      real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s)
-      real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s)
-      real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
-      real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure
+      real,dimension(klon,klev+1),intent(in) :: paprs ! interlayer pressure (Pa)
+      real,dimension(klon,klev),intent(in) :: pplay ! mid-layer pressure (Pa)
+      real,dimension(klon,klev),intent(in) :: pphi ! geopotential at mid-layer
+      real,dimension(klon),intent(in) :: pphis ! surface geopotential
+      real,dimension(klev),intent(in) :: presnivs ! pseudo-pressure (Pa) of mid-layers
+      real,dimension(klon,klev),intent(in) :: u ! eastward zonal wind (m/s)
+      real,dimension(klon,klev),intent(in) :: v ! northward meridional wind (m/s)
+      real,dimension(klon,klev),intent(in) :: rot ! northward meridional wind (m/s)
+      real,dimension(klon,klev),intent(in) :: temp ! temperature (K)
+      real,dimension(klon,klev,nqtot),intent(in) :: qx  ! tracers (.../kg_air)
+      real,dimension(klon,klev),intent(in) :: flxmass_w ! vertical mass flux
+      real,dimension(klon,klev),intent(out) :: d_u ! physics tendency on u (m/s/s)
+      real,dimension(klon,klev),intent(out) :: d_v ! physics tendency on v (m/s/s)
+      real,dimension(klon,klev),intent(out) :: d_t ! physics tendency on t (K/s)
+      real,dimension(klon,klev,nqtot),intent(out) :: d_qx ! physics tendency on tracers
+      real,dimension(klon),intent(out) :: d_ps ! physics tendency on surface pressure
+
+      real, dimension(klon,klev) :: u_loc
+      real, dimension(klon,klev) :: v_loc
+      real, dimension(klon,klev) :: t_loc
+      real, dimension(klon,klev) :: h_loc
+      real, dimension(klon,klev) :: d_u_loc,d_v_loc,d_t_loc,d_h_loc
+
+      real, dimension(klon,klev) :: d_u_dyn,d_v_dyn,d_t_dyn
+      real, dimension(klon,klev,nqtot) :: d_q_dyn
+      real, allocatable, dimension(:,:), save :: u_prev,v_prev,t_prev
+      real, allocatable, dimension(:,:,:), save :: q_prev
+!$OMP THREADPRIVATE(u_prev,v_prev,t_prev,q_prev)
+
+
+
+      real, dimension(klon,klev) :: d_u_vdif,d_v_vdif,d_t_vdif,d_h_vdif
+      real, dimension(klon,klev) :: d_u_the,d_v_the,d_t_the
+      real, dimension(klon,klev,nqtot) :: q_loc,d_q_loc,d_q_vdif,d_q_the
+
+real, dimension(klon) :: capcal,z0m,z0h,dtsrf,emis,fluxsrf,cdh,cdv,tsrf_
+real, dimension(klon,klev) :: zzlay,masse,zpopsk
+real, dimension(klon,klev+1) :: zzlev,kz_v,kz_h,richardson
+
+real, save, allocatable, dimension(:) :: tsrf,f0,zmax0
+real, save, allocatable, dimension(:,:) :: q2
+!$OMP THREADPRIVATE(tsrf,q2,f0,zmax0)
+
+    real,save :: ratqsbas=0.002,ratqshaut=0.3,ratqsp0=50000.,ratqsdp=20000.
+    !$OMP THREADPRIVATE(ratqsbas,ratqshaut,ratqsp0,ratqsdp)
+
+
+real :: z1,z2,tau_thermals
+logical :: lwrite
+integer :: iflag_replay
+
+integer :: iflag_thermals=18
+
+!--------------------------------------------------------------
+! Declaration lscp
+!--------------------------------------------------------------
+  INTEGER                         :: iflag_cld_th    ! flag that determines the distribution of convective clouds    ! IN
+  INTEGER                         :: iflag_ice_thermo! flag to activate the ice thermodynamics                       ! IN
+  LOGICAL                         :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated           ! IN
+  LOGICAL, DIMENSION(klon,klev)   :: ptconv          ! grid points where deep convection scheme is active            ! IN
+  REAL, DIMENSION(klon,klev)      :: ztv             ! virtual potential temperature [K]                             ! IN
+  REAL, DIMENSION(klon,klev)      :: zqta            ! specific humidity within thermals [kg/kg]                     ! IN
+  REAL, DIMENSION(klon,klev+1)      :: frac_the,fm_the
+  REAL, DIMENSION(klon,klev)      :: zpspsk          ! exner potential (p/100000)**(R/cp)                            ! IN
+  REAL, DIMENSION(klon,klev)      :: ztla            ! liquid temperature within thermals [K]                        ! IN
+  REAL, DIMENSION(klon,klev)         :: zthl         ! liquid potential temperature [K]                              ! INOUT
+  REAL, DIMENSION(klon,klev)      :: ratqs            ! function of pressure that sets the large-scale               ! INOUT
+  REAL, DIMENSION(klon,klev)      :: beta             ! conversion rate of condensed water                           ! INOUT
+  REAL, DIMENSION(klon,klev)      :: rneb_seri        ! fraction nuageuse en memoire                                 ! INOUT
+  REAL, DIMENSION(klon,klev)      :: d_t_lscp         ! temperature increment [K]                                    ! OUT
+  REAL, DIMENSION(klon,klev)      :: d_q_lscp         ! specific humidity increment [kg/kg]                          ! OUT
+  REAL, DIMENSION(klon,klev)      :: d_ql_lscp        ! liquid water increment [kg/kg]                               ! OUT
+  REAL, DIMENSION(klon,klev)      :: d_qi_lscp        ! cloud ice mass increment [kg/kg]                             ! OUT
+  REAL, DIMENSION(klon,klev)      :: rneb             ! cloud fraction [-]                                           ! OUT
+  REAL, DIMENSION(klon,klev)      :: rneblsvol        ! cloud fraction per unit volume [-]                           ! OUT
+  REAL, DIMENSION(klon,klev)      :: pfraclr          ! precip fraction clear-sky part [-]                           ! OUT
+  REAL, DIMENSION(klon,klev)      :: pfracld          ! precip fraction cloudy part [-]                              ! OUT
+  REAL, DIMENSION(klon,klev)      :: radocond         ! condensed water used in the radiation scheme [kg/kg]         ! OUT
+  REAL, DIMENSION(klon,klev)      :: radicefrac       ! ice fraction of condensed water for radiation scheme         ! OUT
+  REAL, DIMENSION(klon,klev)      :: rhcl             ! clear-sky relative humidity [-]                              ! OUT
+  REAL, DIMENSION(klon)           :: rain             ! surface large-scale rainfall [kg/s/m2]                       ! OUT
+  REAL, DIMENSION(klon)           :: snow             ! surface large-scale snowfall [kg/s/m2]                       ! OUT
+  REAL, DIMENSION(klon,klev)      :: qsatl            ! saturation specific humidity wrt liquid [kg/kg]              ! OUT
+  REAL, DIMENSION(klon,klev)      :: qsats            ! saturation specific humidity wrt ice [kg/kg]                 ! OUT
+  REAL, DIMENSION(klon,klev+1)    :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]            ! OUT
+  REAL, DIMENSION(klon,klev+1)    :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]            ! OUT
+  REAL, DIMENSION(klon,klev)      :: distcltop        ! distance to cloud top [m]                                    ! OUT
+  REAL, DIMENSION(klon,klev)      :: temp_cltop       ! temperature of cloud top [K]                                 ! OUT
+  REAL, DIMENSION(klon,klev)      :: frac_impa        ! scavenging fraction due tu impaction [-]                     ! OUT
+  REAL, DIMENSION(klon,klev)      :: frac_nucl        ! scavenging fraction due tu nucleation [-]                    ! OUT
+  REAL, DIMENSION(klon,klev)      :: qclr             ! specific total water content in clear sky region [kg/kg]     ! OUT
+  REAL, DIMENSION(klon,klev)      :: qcld             ! specific total water content in cloudy region [kg/kg]        ! OUT
+  REAL, DIMENSION(klon,klev)      :: qss              ! specific total water content in supersat region [kg/kg]      ! OUT
+  REAL, DIMENSION(klon,klev)      :: qvc              ! specific vapor content in clouds [kg/kg]                     ! OUT
+  REAL, DIMENSION(klon,klev)      :: rnebclr          ! mesh fraction of clear sky [-]                               ! OUT
+  REAL, DIMENSION(klon,klev)      :: rnebss           ! mesh fraction of ISSR [-]                                    ! OUT
+  REAL, DIMENSION(klon,klev)      :: gamma_ss         ! coefficient governing the ice nucleation RHi threshold [-]   ! OUT   
+  REAL, DIMENSION(klon,klev)      :: Tcontr           ! threshold temperature for contrail formation [K]             ! OUT
+  REAL, DIMENSION(klon,klev)      :: qcontr           ! threshold humidity for contrail formation [kg/kg]            ! OUT
+  REAL, DIMENSION(klon,klev)      :: qcontr2          ! // (2nd expression more consistent with LMDZ expression of q)! OUT          
+  REAL, DIMENSION(klon,klev)      :: fcontrN          ! fraction of grid favourable to non-persistent contrails      ! OUT
+  REAL, DIMENSION(klon,klev)      :: fcontrP          ! fraction of grid favourable to persistent contrails          ! OUT
+!--------------------------------------------------------------
+
+  REAL, DIMENSION(klon,klev)      :: d_t_eva,d_q_eva,d_ql_eva,d_qi_eva
+include "YOMCST.h"
 
 !    include "clesphys.h"
@@ -58,15 +162,22 @@
 
 
-real :: temp_newton(klon,klev)
-integer :: k
+real,dimension(klon,klev) :: temp_newton
+integer :: i,k,iq
+INTEGER, SAVE :: itap=0
+!$OMP THREADPRIVATE(itap)
+INTEGER, SAVE :: abortphy=0   ! Reprere si on doit arreter en fin de phys
+!$OMP THREADPRIVATE(abortphy)
+
+integer, save :: iflag_reevap=1,iflag_newton=0,iflag_vdif=1,iflag_lscp=1,iflag_cloudth_vert=3,iflag_ratqs=4
+!$OMP THREADPRIVATE(iflag_reevap,iflag_newton,iflag_vdif,iflag_lscp,iflag_cloudth_vert,iflag_ratqs)
+
 logical, save :: first=.true.
 !$OMP THREADPRIVATE(first)
-
-real,save :: rg=9.81
-!$OMP THREADPRIVATE(rg)
 
 ! For I/Os
 integer :: itau0
 real :: zjulian
+real,dimension(klon,klev) :: du0,dv0,dqbs0
+real,dimension(klon,klev) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
 
 
@@ -93,4 +204,43 @@
   ! Initialize IOIPSL output file
 #endif
+call suphel
+call vdif_ini_(klon,RCPD,RD,RG,RKAPPA)
+! Pourquoi ce tau_thermals en argument ??? AFAIRE
+tau_thermals=0.
+call getin_p('iflag_thermals',iflag_thermals)
+
+call getin_p('iflag_newton',iflag_newton)
+call getin_p('iflag_reevap',iflag_reevap)
+call getin_p('iflag_cloudth_vert',iflag_cloudth_vert)
+call getin_p('iflag_ratqs',iflag_ratqs)
+call getin_p('iflag_vdif',iflag_vdif)
+call getin_p('iflag_lscp',iflag_lscp)
+call getin_p('ratqsbas',ratqsbas)
+call getin_p('ratqshaut',ratqshaut)
+call getin_p('ratqsp0',ratqsp0)
+call getin_p('ratqsdp',ratqsdp)
+CALL thermcell_ini(iflag_thermals,0,tau_thermals,6, &
+   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
+CALL lscp_ini(pdtphys,.false.,iflag_ratqs, RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
+
+
+
+allocate(tsrf(klon),q2(klon,klev+1),f0(klon),zmax0(klon))
+allocate(u_prev(klon,klev),v_prev(klon,klev),t_prev(klon,klev),q_prev(klon,klev,nqtot))
+
+u_prev(:,:)=u(:,:)
+v_prev(:,:)=v(:,:)
+t_prev(:,:)=temp(:,:)
+q_prev(:,:,:)=qx(:,:,:)
+
+q2=1.e-10
+tsrf=temp(:,1)
+f0=0.
+zmax0=0.
+
+iflag_replay=0
+call getin_p('iflag_replay',iflag_replay)
+if ( iflag_replay >= 0 ) CALL iophys_ini(pdtphys)
+
 
 endif ! of if (debut)
@@ -100,4 +250,8 @@
 !------------------------------------------------------------
 
+d_u_dyn(:,:)=(u(:,:)-u_prev(:,:))/pdtphys
+d_v_dyn(:,:)=(v(:,:)-v_prev(:,:))/pdtphys
+d_t_dyn(:,:)=(temp(:,:)-t_prev(:,:))/pdtphys
+d_q_dyn(:,:,:)=(qx(:,:,:)-q_prev(:,:,:))/pdtphys
 
 ! set all tendencies to zero
@@ -108,18 +262,280 @@
 d_ps(1:klon)=0.
 
+u_loc(1:klon,1:klev)=u(1:klon,1:klev)
+v_loc(1:klon,1:klev)=v(1:klon,1:klev)
+t_loc(1:klon,1:klev)=temp(1:klon,1:klev)
+d_u_loc(1:klon,1:klev)=0.
+d_v_loc(1:klon,1:klev)=0.
+d_t_loc(1:klon,1:klev)=0.
+do iq=1,nqtot
+   do k=1,klev
+      do i=1,klon
+         q_loc(i,k,iq)=qx(i,k,iq)
+      enddo
+   enddo
+enddo
+
+du0(1:klon,1:klev)=0.
+dv0(1:klon,1:klev)=0.
+dqbs0(1:klon,1:klev)=0.
+
+
+
 !------------------------------------------------------------
 ! Calculs
 !------------------------------------------------------------
 
-! compute tendencies to return to the dynamics:
-! "friction" on the first layer
-d_u(1:klon,1)=-u(1:klon,1)/86400.
-d_v(1:klon,1)=-v(1:klon,1)/86400.
-! newtonian relaxation towards temp_newton()
+!------------------------------------------------------------
+! Rappel en temperature et frottement dans la premiere chouche
+!------------------------------------------------------------
+
+if ( iflag_newton == 1 ) then
+    ! compute tendencies to return to the dynamics:
+    ! "friction" on the first layer
+    d_u(1:klon,1)=-u(1:klon,1)/86400.
+    d_v(1:klon,1)=-v(1:klon,1)/86400.
+    ! newtonian relaxation towards temp_newton()
+    do k=1,klev
+       temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
+       d_t(1:klon,k)=(temp_newton(1:klon,k)-temp(1:klon,k))/5.e5
+    enddo
+else
+   temp_newton(:,:)=0.
+endif
+
+
+!------------------------------------------------------------
+! Reevaporation de la pluie
+!------------------------------------------------------------
+
+iflag_ice_thermo=1
+if ( iflag_reevap == 1 ) then
+      CALL reevap (klon,klev,iflag_ice_thermo,t_loc,q_loc(:,:,1),q_loc(:,:,2),q_loc(:,:,3), &
+&                  d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
+     do k=1,klev
+        do i=1,klon
+           t_loc(i,k)=t_loc(i,k)+d_t_eva(i,k)
+           q_loc(i,k,1)=q_loc(i,k,1)+d_q_eva(i,k)
+           q_loc(i,k,2)=q_loc(i,k,2)+d_ql_eva(i,k)
+           q_loc(i,k,3)=q_loc(i,k,3)+d_qi_eva(i,k)
+           q_loc(i,k,2)=0.
+           q_loc(i,k,3)=0.
+        enddo
+     enddo
+else
+     d_t_eva(:,:)=0.
+     d_q_eva(:,:)=0.
+     d_ql_eva(:,:)=0.
+     d_qi_eva(:,:)=0.
+endif
+
+
+
+!-----------------------------------------------------------------------
+! Variables intermédiaires (altitudes, temperature potentielle ...)
+!-----------------------------------------------------------------------
+
+DO k=1,klev
+   DO i=1,klon
+      zzlay(i,k)=pphi(i,k)/rg
+   ENDDO
+ENDDO
+DO i=1,klon
+   zzlev(i,1)=0.
+ENDDO
+DO k=2,klev
+   DO i=1,klon
+      z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
+      z2=(paprs(i,k)+pplay(i,k))/(paprs(i,k)-pplay(i,k))
+      zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2)
+   ENDDO
+ENDDO
+
+!   Transformation de la temperature en temperature potentielle
+DO k=1,klev
+   DO i=1,klon
+      zpopsk(i,k)=(pplay(i,k)/paprs(i,1))**rkappa
+      masse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
+   ENDDO
+ENDDO
+DO k=1,klev
+   DO i=1,klon
+      h_loc(i,k)=t_loc(i,k)/zpopsk(i,k)
+      d_h_loc(i,k)=d_t_loc(i,k)/zpopsk(i,k)
+      d_q_loc(i,k,1)=0.
+   ENDDO
+ENDDO
+
+!-----------------------------------------------------------------------
+! Diffusion verticale
+!-----------------------------------------------------------------------
+
+if ( iflag_vdif == 1 ) then
+      emis(:)=1.
+      !tsrf=300.
+      z0m=0.035
+      z0h=0.035
+      capcal=1e2
+      lwrite=.false.
+      print*,'lwrite ',lwrite
+            call vdif(klon,klev,                                               &
+     &                pdtphys,capcal,z0m,z0h,                                  &
+     &                pplay,paprs,zzlay,zzlev,                                 &
+     &                u_loc,v_loc,t_loc,h_loc,q_loc,tsrf,emis,                 &
+     &                d_u_loc,d_v_loc,d_h_loc,d_q_loc,fluxsrf,                 &
+     &                d_u_vdif,d_v_vdif,d_h_vdif,d_q_vdif,dtsrf,q2,kz_v,kz_h,  &
+     &                richardson,cdv,cdh,                                      &
+     &                lwrite)
+      do k=1,klev
+         do i=1,klon
+            d_t_vdif(i,k)=d_h_vdif(i,k)*zpopsk(i,k)
+            t_loc(i,k)=t_loc(i,k)+d_t_vdif(i,k)*pdtphys
+            u_loc(i,k)=u_loc(i,k)+d_u_vdif(i,k)*pdtphys
+            v_loc(i,k)=v_loc(i,k)+d_v_vdif(i,k)*pdtphys
+            q_loc(i,k,1)=q_loc(i,k,1)+d_q_vdif(i,k,1)*pdtphys
+         enddo
+      enddo
+      do i=1,klon
+         tsrf(i)=tsrf(i)+dtsrf(i)*pdtphys
+      enddo
+else
+      d_u_vdif(:,:)=0.
+      d_v_vdif(:,:)=0.
+      d_t_vdif(:,:)=0.
+      d_h_vdif(:,:)=0.
+      d_q_vdif(:,:,1)=0.
+      kz_v(:,:)=0.
+      kz_h(:,:)=0.
+      richardson(:,:)=0.
+endif
+
+!-----------------------------------------------------------------------
+! Thermiques
+!-----------------------------------------------------------------------
+
 do k=1,klev
-  temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
-  d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5
+   do i=1,klon
+      d_u_the(i,k)=0.
+      d_v_the(i,k)=0.
+      d_t_the(i,k)=0.
+      d_q_the(i,k,1)=0.
+   enddo
 enddo
 
+if ( iflag_thermals > 0 ) then
+
+
+          zqta(:,:)=q_loc(:,:,1)
+          call caltherm(pdtphys                            &
+         &      ,pplay,paprs,pphi                          &
+         &      ,u_loc,v_loc,t_loc,q_loc,debut             &
+         &      ,f0,zmax0,d_u_the,d_v_the,d_t_the,d_q_the  &
+         &     ,frac_the,fm_the,zqta,ztv,zpspsk,ztla,zthl  &
+         &   )
+    
+          do k=1,klev
+             do i=1,klon
+                t_loc(i,k)=t_loc(i,k)+d_t_the(i,k)
+                u_loc(i,k)=u_loc(i,k)+d_u_the(i,k)
+                v_loc(i,k)=v_loc(i,k)+d_v_the(i,k)
+                q_loc(i,k,1)=q_loc(i,k,1)+d_q_the(i,k,1)
+             enddo
+          enddo
+
+else
+          frac_the(:,:)=0.
+          fm_the(:,:)=0.
+          ztv(:,:)=t_loc(:,:)
+          zqta(:,:)=q_loc(:,:,1)
+          ztla(:,:)=0.
+          zthl(:,:)=0.
+          zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA
+ 
+endif
+
+!-----------------------------------------------------------------------
+! Condensation grande échelle
+!-----------------------------------------------------------------------
+
+iflag_cld_th=5
+ok_ice_sursat=.false.
+ptconv(:,:)=.false.
+distcltop=0.
+temp_cltop=0.
+beta(:,:)=1.
+rneb_seri(:,:)=0.
+do k=1,klev
+  ratqs(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
+  *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.)
+enddo
+
+
+if ( iflag_lscp == 1 ) then
+
+    call lscp(klon,klev,pdtphys,missing_val,            &
+     paprs,pplay,t_loc,q_loc,ptconv,ratqs,                      &
+     d_t_lscp, d_q_lscp, d_ql_lscp, d_qi_lscp, rneb, rneblsvol, rneb_seri,  &
+     pfraclr,pfracld,                                   &
+     radocond, radicefrac, rain, snow,                  &
+     frac_impa, frac_nucl, beta,                        &
+     prfl, psfl, rhcl, zqta, frac_the,                     &
+     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
+     iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
+     distcltop,temp_cltop,                               &
+     qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
+     Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
+     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+
+
+          do k=1,klev
+             do i=1,klon
+                t_loc(i,k)=t_loc(i,k)+d_t_lscp(i,k)
+                q_loc(i,k,1)=q_loc(i,k,1)+d_q_lscp(i,k)
+                q_loc(i,k,2)=q_loc(i,k,2)+d_ql_lscp(i,k)
+                q_loc(i,k,3)=q_loc(i,k,3)+d_qi_lscp(i,k)
+             enddo
+          enddo
+
+else
+     d_t_lscp(:,:)=0.
+     d_q_lscp(:,:)=0.
+     d_ql_lscp(:,:)=0.
+     d_qi_lscp(:,:)=0.
+     rneb(:,:)=0.
+     rneblsvol(:,:)=0.
+     pfraclr(:,:)=0.
+     pfracld(:,:)=0.
+     radocond(:,:)=0.
+     rain(:)=0.
+     snow(:)=0.
+     radicefrac(:,:)=0.
+     rhcl      (:,:)=0.
+     qsatl     (:,:)=0.
+     qsats     (:,:)=0.
+     prfl      (:,:)=0.
+     psfl      (:,:)=0.
+     distcltop (:,:)=0.
+     temp_cltop(:,:)=0.
+     frac_impa (:,:)=0.
+     frac_nucl (:,:)=0.
+     qclr      (:,:)=0.
+     qcld      (:,:)=0.
+     qss       (:,:)=0.
+     qvc       (:,:)=0.
+     rnebclr   (:,:)=0.
+     rnebss    (:,:)=0.
+     gamma_ss  (:,:)=0.
+     Tcontr    (:,:)=0.
+     qcontr    (:,:)=0.
+     qcontr2   (:,:)=0.
+     fcontrN   (:,:)=0.
+     fcontrP   (:,:)=0.
+endif
+
+
+d_u(:,:)=(u_loc(:,:)-u(:,:))/pdtphys
+d_v(:,:)=(v_loc(:,:)-v(:,:))/pdtphys
+d_t(:,:)=(t_loc(:,:)-temp(:,:))/pdtphys
+d_qx(:,:,:)=(q_loc(:,:,:)-qx(:,:,:))/pdtphys
 
 !------------------------------------------------------------
@@ -128,5 +544,15 @@
 
 
-call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,t,qx,0.*t,0.*t,0.*t,0.*t,0.*t,0.*t)
+tsrf_(:)=tsrf(:)
+if ( iflag_replay == -1 ) then
+   call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,temp,qx,0.*u,0.*u,0.*u,0.*u,q2,0.*u)
+else if (iflag_replay == 0 ) then
+     ! En mode replay, on sort aussi les variables de base
+! Les lignes qui suivent ont été générées automatiquement avec :
+! ( for i in `grep -i 'real.*::' physiqex_mod.F90 | sed -e '/^!/d' | grep '(klon,klev' | cut -d: -f3 | cut -d! -f1 | sed -e 's/,/ /g' -e '/rot/d'` ; do echo ' call iophys_ecrit("'$i'",klev,"","",'$i')' ; done ) > physiqex_out.h
+! ( for i in `grep -i 'real.*::' physiqex_mod.F90 | sed -e '/^!/d' | grep '(klon)' | cut -d: -f3 | cut -d! -f1 | sed -e 's/,/ /g' -e '/rot/d'` ; do echo ' call iophys_ecrit("'$i'",1,"","",'$i')' ; done ) >> physiqex_out.h
+include "physiqex_out.h"
+
+endif
 
 
@@ -136,4 +562,5 @@
 endif
 
+print*,'Fin physiqex'
 
 end subroutine physiqex
Index: /LMDZ6/trunk/libf/phylmd/vdif_kcay.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/vdif_kcay.F90	(revision 4653)
+++ /LMDZ6/trunk/libf/phylmd/vdif_kcay.F90	(revision 4654)
@@ -2,7 +2,7 @@
 ! $Header$
 
-SUBROUTINE vdif_kcay(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, &
+SUBROUTINE vdif_kcay(klon,klev, dt, g, rconst, plev, temp, zlev, zlay, u, v, &
     teta, cd, q2, q2diag, km, kn, ustar, l_mix)
-  USE dimphy
+
   IMPLICIT NONE
 
@@ -28,20 +28,26 @@
 
   ! .......................................................................
-  REAL dt, g, rconst
-  REAL plev(klon, klev+1), temp(klon, klev)
-  REAL ustar(klon), snstable
-  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 q2(klon, klev+1), q2s(klon, klev+1)
-  REAL q2diag(klon, klev+1)
-  REAL km(klon, klev+1)
-  REAL kn(klon, klev+1)
-  REAL sq(klon), sqz(klon), zz(klon, klev+1), zq, long0(klon)
-
-  INTEGER l_mix, iii
+  INTEGER, INTENT(IN) :: klon,klev
+  REAL,INTENT(IN) :: dt, g, rconst
+  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: plev
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: temp
+  REAL,DIMENSION(klon),INTENT(IN) :: ustar(klon)
+  REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: zlev
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: zlay
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: u
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: v
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: teta
+  REAL,DIMENSION(klon),INTENT(IN) :: cd(klon)
+  REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: q2
+  REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: q2diag
+  REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: km
+  REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: kn
+  INTEGER, INTENT(OUT) :: l_mix
+
+  REAL,DIMENSION(klon) :: sq, sqz, long0
+  REAL,DIMENSION(klon,klev+1) :: q2s,zz
+  REAL :: snstable,zq
+
+  INTEGER :: iii
   ! .......................................................................
 
@@ -56,8 +62,7 @@
 
   ! .......................................................................
-  INTEGER nlay, nlev, ngrid
-  REAL unsdz(klon, klev)
-  REAL unsdzdec(klon, klev+1)
-  REAL q(klon, klev+1)
+  INTEGER :: nlay, nlev, ngrid
+  REAL, DIMENSION(klon,klev) :: unsdz
+  REAL, DIMENSION(klon, klev+1) :: unsdzdec,q
 
   ! .......................................................................
@@ -70,8 +75,8 @@
 
   ! .......................................................................
-  REAL kmpre(klon, klev+1)
-  REAL qcstat
-  REAL q2cstat
-  REAL sss, sssq
+  REAL, DIMENSION(klon, klev+1) :: kmpre
+  REAL :: qcstat
+  REAL :: q2cstat
+  REAL :: sss, sssq
   ! .......................................................................
 
@@ -79,5 +84,5 @@
 
   ! .......................................................................
-  REAL long(klon, klev+1)
+  REAL, DIMENSION(klon, klev+1) :: long
   ! .......................................................................
 
@@ -97,13 +102,10 @@
 
   ! .......................................................................
-  REAL kmq3
-  REAL kmcstat
-  REAL knq3
-  REAL mcstat
-  REAL m2cstat
-  REAL m(klon, klev+1)
-  REAL mpre(klon, klev+1)
-  REAL m2(klon, klev+1)
-  REAL n2(klon, klev+1)
+  REAL :: kmq3
+  REAL :: kmcstat
+  REAL :: knq3
+  REAL :: mcstat
+  REAL :: m2cstat
+  REAL, DIMENSION(klon, klev+1) :: m,mpre,m2,n2
   ! .......................................................................
 
@@ -121,15 +123,10 @@
 
   ! .......................................................................
-  REAL gn
-  REAL gnmin
-  REAL gnmax
-  LOGICAL gninf
-  LOGICAL gnsup
-  REAL gm
-  ! REAL ri(klon,klev+1)
-  REAL sn(klon, klev+1)
-  REAL snq2(klon, klev+1)
-  REAL sm(klon, klev+1)
-  REAL smq2(klon, klev+1)
+  REAL :: gn,gm
+  REAL :: gnmin
+  REAL :: gnmax
+  LOGICAL :: gninf
+  LOGICAL :: gnsup
+  REAL, DIMENSION(klon, klev+1) :: sn, snq2, sm, smq2
   ! .......................................................................
 
@@ -142,9 +139,9 @@
 
   ! .......................................................................
-  REAL kappa
-  REAL long00
-  REAL a1, a2, b1, b2, c1
-  REAL cn1, cn2
-  REAL cm1, cm2, cm3, cm4
+  REAL :: kappa
+  REAL :: long00
+  REAL :: a1, a2, b1, b2, c1
+  REAL :: cn1, cn2
+  REAL :: cm1, cm2, cm3, cm4
   ! .......................................................................
 
@@ -155,8 +152,8 @@
 
   ! .......................................................................
-  REAL termq
-  REAL termq3
-  REAL termqm2
-  REAL termq3m2
+  REAL :: termq
+  REAL :: termq3
+  REAL :: termqm2
+  REAL :: termq3m2
   ! .......................................................................
 
@@ -165,15 +162,15 @@
 
   ! .......................................................................
-  REAL q2min
-  REAL q2max
+  REAL :: q2min
+  REAL :: q2max
   ! .......................................................................
   ! knmin : borne inferieure de kn
   ! kmmin : borne inferieure de km
   ! .......................................................................
-  REAL knmin
-  REAL kmmin
-  ! .......................................................................
-  INTEGER ilay, ilev, igrid
-  REAL tmp1, tmp2
+  REAL :: knmin
+  REAL :: kmmin
+  ! .......................................................................
+  INTEGER :: ilay, ilev, igrid
+  REAL :: tmp1, tmp2
   ! .......................................................................
   PARAMETER (kappa=0.4E+0)
@@ -202,5 +199,5 @@
   PARAMETER (cm4=-9.E+0*a1*a2)
 
-  LOGICAL first
+  LOGICAL :: first
   SAVE first
   DATA first/.TRUE./
@@ -211,13 +208,16 @@
 
   ! Initialisation de q2
+ngrid=klon
   nlay = klev
   nlev = klev + 1
 
-  CALL yamada(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, cd, &
-    q2diag, km, kn, ustar, l_mix)
-  IF (first .AND. 1==1) THEN
-    first = .FALSE.
-    q2 = q2diag
-  END IF
+! Initialisation avec un schema d'equilibre
+! CALL yamada(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, cd, &
+!   q2diag, km, kn, ustar, l_mix)
+! IF (first .AND. 1==1) THEN
+!   first = .FALSE.
+!   q2 = q2diag
+! END IF
+q2diag=0.
 
   DO ilev = 2, nlay
