Index: LMDZ4/trunk/libf/phytherm/ajsec.F
===================================================================
--- LMDZ4/trunk/libf/phytherm/ajsec.F	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/ajsec.F	(revision 852)
@@ -2,5 +2,159 @@
 ! $Header$
 !
-      SUBROUTINE ajsec(paprs, pplay, t,q, d_t,d_q)
+      SUBROUTINE ajsec(paprs, pplay, t,q,limbas,d_t,d_q)
+      USE dimphy
+      IMPLICIT none
+c======================================================================
+c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
+c Objet: ajustement sec (adaptation du GCM du LMD)
+c======================================================================
+c Arguments:
+c t-------input-R- Temperature
+c
+c d_t-----output-R-Incrementation de la temperature
+c======================================================================
+cym#include "dimensions.h"
+cym#include "dimphy.h"
+#include "YOMCST.h"
+      REAL paprs(klon,klev+1), pplay(klon,klev)
+      REAL t(klon,klev), q(klon,klev)
+      REAL d_t(klon,klev), d_q(klon,klev)
+c
+      INTEGER limbas(klon), limhau ! les couches a ajuster
+c
+      LOGICAL mixq
+ccc      PARAMETER (mixq=.TRUE.)
+      PARAMETER (mixq=.FALSE.)
+c
+      REAL zh(klon,klev)
+      REAL zho(klon,klev)
+      REAL zq(klon,klev)
+      REAL zpk(klon,klev)
+      REAL zpkdp(klon,klev)
+      REAL hm, sm, qm
+      LOGICAL modif(klon), down
+      INTEGER i, k, k1, k2
+c
+c Initialisation:
+c
+cym
+      limhau=klev
+  
+      DO k = 1, klev
+      DO i = 1, klon
+         d_t(i,k) = 0.0
+         d_q(i,k) = 0.0
+      ENDDO
+      ENDDO
+c------------------------------------- detection des profils a modifier
+      DO k = 1, limhau
+      DO i = 1, klon
+         zpk(i,k) = pplay(i,k)**RKAPPA
+         zh(i,k) = RCPD * t(i,k)/ zpk(i,k)
+         zho(i,k) = zh(i,k)
+         zq(i,k) = q(i,k)
+      ENDDO
+      ENDDO
+c
+      DO k = 1, limhau
+      DO i = 1, klon
+         zpkdp(i,k) = zpk(i,k) * (paprs(i,k)-paprs(i,k+1))
+      ENDDO
+      ENDDO
+c
+      DO i = 1, klon
+         modif(i) = .FALSE.
+      ENDDO
+      DO k = 2, limhau
+      DO i = 1, klon
+      IF (.NOT.modif(i).and.k-1>limbas(i)) THEN
+         IF ( zh(i,k).LT.zh(i,k-1) ) modif(i) = .TRUE.
+      ENDIF
+      ENDDO
+      ENDDO
+c------------------------------------- correction des profils instables
+      DO 1080 i = 1, klon
+      IF (modif(i)) THEN
+          k2 = limbas(i)
+ 8000     CONTINUE
+            k2 = k2 + 1
+            IF (k2 .GT. limhau) goto 8001
+            IF (zh(i,k2) .LT. zh(i,k2-1)) THEN
+              k1 = k2 - 1
+              k = k1
+              sm = zpkdp(i,k2)
+              hm = zh(i,k2)
+              qm = zq(i,k2)
+ 8020         CONTINUE
+                sm = sm +zpkdp(i,k)
+                hm = hm +zpkdp(i,k) * (zh(i,k)-hm) / sm
+                qm = qm +zpkdp(i,k) * (zq(i,k)-qm) / sm
+                down = .FALSE.
+                IF (k1 .ne. limbas(i)) THEN
+                  IF (hm .LT. zh(i,k1-1)) down = .TRUE.
+                ENDIF
+                IF (down) THEN
+                  k1 = k1 - 1
+                  k = k1
+                ELSE
+                  IF ((k2 .EQ. limhau)) GOTO 8021
+                  IF ((zh(i,k2+1).GE.hm)) GOTO 8021
+                  k2 = k2 + 1
+                  k = k2
+                ENDIF
+              GOTO 8020
+ 8021         CONTINUE
+c------------ nouveau profil : constant (valeur moyenne)
+              DO k = k1, k2
+                zh(i,k) = hm
+                zq(i,k) = qm
+              ENDDO
+              k2 = k2 + 1
+            ENDIF
+          GOTO 8000
+ 8001     CONTINUE
+      ENDIF
+ 1080 CONTINUE
+c
+      DO k = 1, limhau
+      DO i = 1, klon
+         d_t(i,k) = (zh(i,k)-zho(i,k))*zpk(i,k)/RCPD
+         d_q(i,k) = zq(i,k) - q(i,k)
+      ENDDO
+      ENDDO
+c
+! FH : les d_q et d_t sont maintenant calcules de facon a valoir
+! effectivement 0. si on ne fait rien.
+!
+!     IF (limbas.GT.1) THEN
+!     DO k = 1, limbas-1
+!     DO i = 1, klon
+!        d_t(i,k) = 0.0
+!        d_q(i,k) = 0.0
+!     ENDDO
+!     ENDDO
+!     ENDIF
+c
+!     IF (limhau.LT.klev) THEN
+!     DO k = limhau+1, klev
+!     DO i = 1, klon
+!        d_t(i,k) = 0.0
+!        d_q(i,k) = 0.0
+!     ENDDO
+!     ENDDO
+!     ENDIF
+c
+      IF (.NOT.mixq) THEN
+      DO k = 1, klev
+      DO i = 1, klon
+         d_q(i,k) = 0.0
+      ENDDO
+      ENDDO
+      ENDIF
+c
+      RETURN
+      END
+
+      SUBROUTINE ajsec_convV2(paprs, pplay, t,q, d_t,d_q)
       USE dimphy
       IMPLICIT none
Index: LMDZ4/trunk/libf/phytherm/clouds_gno.F
===================================================================
--- LMDZ4/trunk/libf/phytherm/clouds_gno.F	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/clouds_gno.F	(revision 852)
@@ -219,5 +219,17 @@
 c -- test numerical convergence:
 
-c         print*,'avant test ',i,k,lconv(i),u(i),v(i)
+          if (beta(i).lt.1.e-10) then
+              print*,'avant test ',i,k,lconv(i),u(i),v(i),beta(i)
+              stop
+          endif
+          if (abs(fprime(i)).lt.1.e-11) then
+              print*,'avant test fprime<.e-11 '
+     s        ,i,k,lconv(i),u(i),v(i),beta(i),fprime(i)
+              print*,'klon,ND,R,RS,QSUB',
+     s        klon,ND,R(i,k),rs(i,k),qsub(i,k)
+              fprime(i)=sign(1.e-11,fprime(i))
+          endif
+
+
           if ( ABS(dist(i)/beta(i)) .LT. epsilon ) then 
 c           print*,'v-u **2',(v(i)-u(i))**2
Index: LMDZ4/trunk/libf/phytherm/coef_diff_turb_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phytherm/coef_diff_turb_mod.F90	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/coef_diff_turb_mod.F90	(revision 852)
@@ -18,5 +18,5 @@
   SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, &
        ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, &
-       ycoefm, ycoefh)
+       ycoefm, ycoefh ,yq2)
 ! 
 ! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the 
@@ -35,4 +35,5 @@
     REAL, DIMENSION(klon,klev), INTENT(IN)     :: yq, yt
     REAL, DIMENSION(klon), INTENT(IN)          :: yts, yrugos, yqsurf
+    REAL, DIMENSION(klon,klev+1)               :: yq2
   
 ! Output arguments
@@ -41,16 +42,9 @@
     REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefm
 
-! Local variables with attribute SAVE
-!****************************************************************************************
-    REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE  :: q2save
-    !$OMP THREADPRIVATE(q2save)
-    LOGICAL, SAVE                              :: first_call=.TRUE.
-    !$OMP THREADPRIVATE(first_call)    
-
 ! Other local variables
 !****************************************************************************************
     INTEGER                                    :: k, i, j
     REAL, DIMENSION(klon,klev)                 :: ycoefm0, ycoefh0, yzlay, yteta
-    REAL, DIMENSION(klon,klev+1)               :: yzlev, yq2, q2diag, ykmm, ykmn, ykmq
+    REAL, DIMENSION(klon,klev+1)               :: yzlev, q2diag, ykmm, ykmn, ykmq
     REAL, DIMENSION(klon)                      :: y_cd_h, y_cd_m
     REAL, DIMENSION(klon)                      :: yustar
@@ -94,5 +88,5 @@
             ycoefm0, ycoefh0)
 
-       DO k = 1, klev
+       DO k = 2, klev
           DO i = 1, knon
              ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
@@ -121,5 +115,5 @@
             ycoefm0,ycoefh0)
        
-       DO k = 1, klev
+       DO k = 2, klev
           DO i = 1, knon
              ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
@@ -137,12 +131,4 @@
 
     IF (iflag_pbl.GE.3) THEN
-
-       IF (first_call) THEN
-          PRINT*, "Using method MELLOR&YAMADA"
-          ! NB! q2save could/should be read and written from (re)startphy.nc
-          ALLOCATE(q2save(klon,klev+1,nbsrf))
-          q2save(:,:,:) = 1.e-8
-          first_call=.FALSE.
-       END IF
 
        yzlay(1:knon,1)= &
@@ -173,21 +159,13 @@
        END DO
 
-! Compresse q2save
-       DO k = 1, klev+1
-          DO j = 1, knon
-             i = ni(j)
-             yq2(j,k)=q2save(i,k,nsrf)
-          ENDDO
-       ENDDO
-         
-!   Bug introduit volontairement pour converger avec les resultats
-!   du papier sur les thermiques.
-       IF (1.EQ.1) THEN
-          y_cd_m(1:knon) = ycoefm(1:knon,1)
-          y_cd_h(1:knon) = ycoefh(1:knon,1)
-       ELSE
-          y_cd_h(1:knon) = ycoefm(1:knon,1)
-          y_cd_m(1:knon) = ycoefh(1:knon,1)
-       ENDIF
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un
+! bug sur les coefficients de surface :
+!          y_cd_h(1:knon) = ycoefm(1:knon,1)
+!          y_cd_m(1:knon) = ycoefh(1:knon,1)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+       y_cd_m(1:knon) = ycoefm(1:knon,1)
+       y_cd_h(1:knon) = ycoefh(1:knon,1)
        
        CALL ustarhb(knon,yu,yv,y_cd_m, yustar)
@@ -215,12 +193,4 @@
        ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
                 
-! update q2save with yq2
-       DO j=1,knon
-          DO k=1,klev+1
-             i=ni(j)
-             q2save(i,k,nsrf) = yq2(j,k)
-          END DO
-       END DO
-       
     ENDIF !(iflag_pbl.ge.3)
 
Index: LMDZ4/trunk/libf/phytherm/ini_histday.h
===================================================================
--- LMDZ4/trunk/libf/phytherm/ini_histday.h	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/ini_histday.h	(revision 852)
@@ -55,9 +55,9 @@
          CALL histdef(nid_day, "weakinv", "Weak inversion", "",
      s           iim,jjmp1,nhori, 1,1,1, -99, 32,
-     s           "ave(X)", zsto,zout)
+     s           "ave(X)", zstophy,zout)
 
          CALL histdef(nid_day, "dthmin", "dTheta mini", "K/m",
      s           iim,jjmp1,nhori, 1,1,1, -99, 32,
-     s           "ave(X)", zsto,zout)
+     s           "ave(X)", zstophy,zout)
 
 
@@ -547,4 +547,20 @@
      $         "ave(X)", zstophy,zout)
 C
+
+
+! FH Sorties specifiques pour Mellor et Yamada
+      if (iflag_pbl>1 .and. lev_histday.gt.10 ) then
+           call histdef(nid_day, "tke_"//clnsurf(nsrf), 
+     $         "Max Turb. Kinetic Energy "//clnsurf(nsrf), "-",  
+     $         iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     $         "ave(X)", zstophy,zout)
+
+           call histdef(nid_day, "tke_max_"//clnsurf(nsrf), 
+     $         "Max Turb. Kinetic Energy "//clnsurf(nsrf), "-",  
+     $         iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     $         "t_max(X)", zstophy,zout)
+      endif	
+
+C
          END DO 
 C           
Index: LMDZ4/trunk/libf/phytherm/ini_histmth.h
===================================================================
--- LMDZ4/trunk/libf/phytherm/ini_histmth.h	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/ini_histmth.h	(revision 852)
@@ -285,5 +285,5 @@
 c
          CALL histdef(nid_mth, "sens", "Sensible heat flux", "W/m2",
-     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                iim,jj_nb,nhori, 1,1,1, -99, 32, 
      .                "ave(X)", zstophy,zout)
 c
@@ -310,5 +310,5 @@
          CALL histdef(nid_mth, "fqfonte","Land ice melt",
      .                "kg/m2/s",iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "ave(X)", zsto,zout)
+     .                "ave(X)", zstophy,zout)
 
          DO nsrf = 1, nbsrf
@@ -700,4 +700,24 @@
       IF(lev_histmth.GE.4) THEN
 c
+!FH Sorties pour la couche limite
+      if (iflag_pbl>1) then
+         CALL histdef(nid_mth, "tke","TKE","m2/s2",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zstophy,zout)
+c
+         CALL histdef(nid_mth, "tke_max","TKE max","m2/s2",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "t_max(X)", zstophy,zout)
+      endif
+c
+         CALL histdef(nid_mth, "kz","Kz melange","m2/s",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zstophy,zout)
+c
+         CALL histdef(nid_mth, "kz_max","Kz melange max","m2/s",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "t_max(X)", zstophy,zout)
+
+
          CALL histdef(nid_mth, "clwcon", 
      .                "Convective Cloud Liquid water content"
@@ -730,4 +750,7 @@
      .                "ave(X)", zstophy,zout)
 c
+
+
+c
 c        CALL histdef(nid_mth, "ducon", "Convection du", "m/s2",
 c    .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
@@ -774,4 +797,12 @@
      .                "ave(X)", zstophy,zout)
 
+c
+         CALL histdef(nid_mth, "dtthe", "Dry adjust. dT", "K/s",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zstophy,zout)
+
+         CALL histdef(nid_mth,"dqthe","Dry adjust. dQ","(kg/kg)/s",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zstophy,zout)
 c
          CALL histdef(nid_mth, "dtajs", "Dry adjust. dT", "K/s",
@@ -914,4 +945,12 @@
 c
           CALL histdef(nid_mth, "dtcon","dtcon","K/s",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zstophy,zout)
+c
+         CALL histdef(nid_mth, "dtthe", "Dry adjust. dT", "K/s",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zstophy,zout)
+
+         CALL histdef(nid_mth,"dqthe","Dry adjust. dQ","(kg/kg)/s",
      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
      .                "ave(X)", zstophy,zout)
@@ -1622,4 +1661,26 @@
       IF(lev_histmth.GE.4) THEN
 c
+!FH Sorties pour la couche limite
+      if (iflag_pbl>1) then
+         CALL histdef(nid_mth, "tke","TKE","m2/s2",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zstophy,zout)
+c
+         CALL histdef(nid_mth, "tke_max","TKE max","m2/s2",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "t_max(X)", zstophy,zout)
+      endif
+c
+         CALL histdef(nid_mth, "kz","Kz melange","m2/s",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zstophy,zout)
+c
+         CALL histdef(nid_mth, "kz_max","Kz melange max","m2/s",
+     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
+     .                "t_max(X)", zstophy,zout)
+
+
+
+
          CALL histdef(nid_mth, "clwcon", 
      .                "Convective Cloud Liquid water content"
Index: LMDZ4/trunk/libf/phytherm/pbl_surface_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phytherm/pbl_surface_mod.F90	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/pbl_surface_mod.F90	(revision 852)
@@ -197,5 +197,5 @@
        wfbils,    wfbilo,    flux_t,   flux_u, flux_v,&
        dflux_t,   dflux_q,   zxsnow,                  &
-       zxfluxt,   zxfluxq,   q2m,      flux_q )
+       zxfluxt,   zxfluxq,   q2m,      flux_q, tke    )
 !****************************************************************************************
 ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
@@ -240,4 +240,5 @@
 ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
 !                    (orientation positive vers le bas)
+! tke---input/output-R- tke (kg/m**2/s)
 ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
 ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
@@ -348,4 +349,7 @@
     REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m
     REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q
+
+! Input/output
+    REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: tke
 
 
@@ -421,4 +425,5 @@
     REAL, DIMENSION(klon,klev)         :: delp
     REAL, DIMENSION(klon,klev+1)       :: ypaprs
+    REAL, DIMENSION(klon,klev+1)       :: ytke
     REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
     REAL, DIMENSION(klon,nbsrf)        :: pctsrf_pot
@@ -451,4 +456,9 @@
 
 
+!jg+ temporary test
+    REAL, DIMENSION(klon,klev)         :: y_flux_u_old, y_flux_v_old
+    REAL, DIMENSION(klon,klev)         :: y_d_u_old, y_d_v_old
+!jg-
+    
 !****************************************************************************************
 ! End of declarations
@@ -514,5 +524,5 @@
     flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0      
     d_u = 0.0     ; d_v = 0.0        ; zcoefh = 0.0    ; yqsol = 0.0    
-    ytherm = 0.0
+    ytherm = 0.0  ; ytke=0.
      
     ytsoil = 999999. 
@@ -657,4 +667,5 @@
              ypplay(j,k) = pplay(i,k)
              ydelp(j,k) = delp(i,k)
+             ytke(j,k)=tke(i,k,nsrf)
              yu(j,k) = u(i,k)
              yv(j,k) = v(i,k)
@@ -687,5 +698,5 @@
        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf,  &
-            ycoefm, ycoefh)
+            ycoefm, ycoefh,ytke)
        
 !****************************************************************************************
@@ -839,4 +850,5 @@
 !****************************************************************************************
 
+       tke(:,:,nsrf)=0.
        DO k = 1, klev
           DO j = 1, knon
@@ -853,4 +865,7 @@
              flux_u(i,k,nsrf) = y_flux_u(j,k)
              flux_v(i,k,nsrf) = y_flux_v(j,k)
+
+             tke(i,k,nsrf)=ytke(j,k)
+
           ENDDO
        ENDDO
Index: LMDZ4/trunk/libf/phytherm/phyetat0.F
===================================================================
--- LMDZ4/trunk/libf/phytherm/phyetat0.F	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/phyetat0.F	(revision 852)
@@ -13,5 +13,5 @@
      .           rugsrel_p,tabcntr0,
      .           t_ancien_p,q_ancien_p,ancien_ok_p, rnebcon_p, ratqs_p,
-     .           clwcon_p)
+     .           clwcon_p,pbl_tke_p)
 
       USE dimphy
@@ -37,4 +37,6 @@
 #include "clesphys.h"
 #include "temps.h"
+#include "thermcell.h"
+#include "compbl.h"
 c======================================================================
       CHARACTER*(*) fichnom
@@ -45,4 +47,5 @@
       REAL solaire_etat0
       REAL tsol_p(klon,nbsrf)
+      REAL pbl_tke_p(klon,klev,nbsrf)
       REAL tsoil_p(klon,nsoilmx,nbsrf)
       REAL tslab_p(klon), seaice_p(klon)
@@ -83,4 +86,5 @@
       REAL rlat(klon_glo), rlon(klon_glo)
       REAL tsol(klon_glo,nbsrf)
+      REAL pbl_tke(klon_glo,klev,nbsrf)
       REAL tsoil(klon_glo,nsoilmx,nbsrf)
 cIM "slab" ocean
@@ -127,5 +131,5 @@
 c
       INTEGER nid, nvarid
-      INTEGER ierr, i, nsrf, isoil 
+      INTEGER ierr, i, nsrf, isoil ,k
       INTEGER length
       PARAMETER (length=100)
@@ -1492,4 +1496,46 @@
       xmax = MAXval(run_off_lic_0)
       PRINT*,'(ecart-type) run_off_lic_0:', xmin, xmax
+
+
+c Lecture de l'energie cinetique turbulente
+c
+
+      IF (iflag_pbl>1) then
+         PRINT*, 'phyetat0: Le champ <TKE> est absent'
+         PRINT*, '          Mais je vais essayer de lire TKE**'
+         DO nsrf = 1, nbsrf
+           IF (nsrf.GT.99) THEN
+             PRINT*, "Trop de sous-mailles"
+             CALL abort
+           ENDIF
+           WRITE(str2,'(i2.2)') nsrf
+           ierr = NF_INQ_VARID (nid, "TKE"//str2, nvarid)
+           IF (ierr.NE.NF_NOERR) THEN
+              PRINT*, "WARNING phyetat0: <TKE"//str2//"> est absent"
+              pbl_tke(:,:,nsrf)=1.e-8
+           ELSE
+#ifdef NC_DOUBLE
+              ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pbl_tke(1,1,nsrf))
+#else
+              ierr = NF_GET_VAR_REAL(nid, nvarid, pbl_tke(1,1,nsrf))
+#endif
+              IF (ierr.NE.NF_NOERR) THEN
+                PRINT*, "WARNING phyetat0: echec lect <TKE"//str2//">"
+                CALL abort
+              ENDIF
+           ENDIF
+
+           xmin = 1.0E+20
+           xmax = -1.0E+20
+           DO k = 1, klev
+           DO i = 1, klon_glo
+              xmin = MIN(pbl_tke(i,k,nsrf),xmin)
+              xmax = MAX(pbl_tke(i,k,nsrf),xmax)
+           ENDDO
+           ENDDO
+           PRINT*,'Temperature du sol TKE**:', nsrf, xmin, xmax
+         ENDDO
+      ENDIF
+
 c
 c Fermer le fichier:
@@ -1516,4 +1562,7 @@
       call Scatter( rlon,rlon_p)
       call Scatter( tsol,tsol_p)
+      IF (iflag_pbl>1) then
+         call Scatter( pbl_tke,pbl_tke_p)
+      endif
       call Scatter( tsoil,tsoil_p)
       call Scatter( tslab,tslab_p)
Index: LMDZ4/trunk/libf/phytherm/phyredem.F
===================================================================
--- LMDZ4/trunk/libf/phytherm/phyredem.F	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/phyredem.F	(revision 852)
@@ -9,5 +9,6 @@
      .           radsol_p,zmea_p,zstd_p,zsig_p,
      .           zgam_p,zthe_p,zpic_p,zval_p,rugsrel_p,
-     .           t_ancien_p, q_ancien_p, rnebcon_p, ratqs_p, clwcon_p)
+     .           t_ancien_p, q_ancien_p, rnebcon_p, ratqs_p, clwcon_p,
+     .           pbl_tke_p)
 
       USE dimphy
@@ -29,4 +30,6 @@
 #include "control.h"
 #include "temps.h"
+#include "thermcell.h"
+#include "compbl.h"
 c======================================================================
       CHARACTER*(*) fichnom
@@ -35,4 +38,5 @@
       REAL rlat_p(klon), rlon_p(klon)
       REAL tsol_p(klon,nbsrf)
+      REAL pbl_tke_p(klon,klev,nbsrf)
       REAL tsoil_p(klon,nsoilmx,nbsrf)
       CHARACTER*6 ocean
@@ -70,4 +74,5 @@
       REAL rlat(klon_glo), rlon(klon_glo)
       REAL tsol(klon_glo,nbsrf)
+      REAL pbl_tke(klon_glo,klev,nbsrf)
       REAL tsoil(klon_glo,nsoilmx,nbsrf)
       REAL tslab(klon_glo), seaice(klon_glo)
@@ -135,4 +140,5 @@
       call Gather( rlon_p,rlon)
       call Gather( tsol_p,tsol)
+      call Gather( pbl_tke_p,pbl_tke)
       call Gather( tsoil_p,tsoil)
       call Gather( tslab_p,tslab)
@@ -899,4 +905,35 @@
 c
 c
+!!!!!!!!!!!!!!!!!!!! DEB TKE PBL !!!!!!!!!!!!!!!!!!!!!!!!!
+c
+      IF (iflag_pbl>1) then
+      DO nsrf = 1, nbsrf
+        IF (nsrf.LE.99) THEN
+            WRITE(str2,'(i2.2)') nsrf
+            ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+            ierr = NF_DEF_VAR (nid,"TKE"//str2,NF_DOUBLE,1,idim3
+     $          ,nvarid)
+#else
+            ierr = NF_DEF_VAR (nid,"TKE"//str2,NF_FLOAT,1,idim3
+     $          ,nvarid)
+#endif
+            ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 15,
+     .                        "Energ. Cineti. Turb."//str2)
+            ierr = NF_ENDDEF(nid)
+        ELSE
+            PRINT*, "Trop de sous-mailles"
+            CALL abort
+        ENDIF
+#ifdef NC_DOUBLE
+        ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pbl_tke(:,:,nsrf))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,pbl_tke(:,:,nsrf))
+#endif
+      ENDDO
+      ENDIF
+
+!!!!!!!!!!!!!!!!!!!! FIN TKE PBL !!!!!!!!!!!!!!!!!!!!!!!!!
+c
       ierr = NF_CLOSE(nid)
 c
Index: LMDZ4/trunk/libf/phytherm/physiq.F
===================================================================
--- LMDZ4/trunk/libf/phytherm/physiq.F	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/physiq.F	(revision 852)
@@ -189,4 +189,5 @@
 
       integer lmax_th(klon)
+      integer limbas(klon)
       real ratqscth(klon,klev)
       real ratqsdiff(klon,klev)
@@ -945,4 +946,8 @@
       REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
       REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
+
+      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: pbl_tke ! turb kinetic energy
+c   !$OMP THREADPRIVATE(pbl_tke)
+
 c
       REAL zxfluxt(klon, klev)
@@ -1066,4 +1071,5 @@
 c$OMP THREADPRIVATE(d_u_con,d_v_con)
       REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
+      REAL d_t_ajsb(klon,klev), d_q_ajsb(klon,klev)
       REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
       REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
@@ -1554,4 +1560,16 @@
          itaprad = 0
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! Un petit travail à faire ici.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+         if (iflag_pbl>1) then
+             PRINT*, "Using method MELLOR&YAMADA" 
+         endif
+          ! NB! pbl_tke could/should be read and written from (re)startphy.nc
+          ALLOCATE(pbl_tke(klon,klev+1,nbsrf))  
+          pbl_tke(:,:,:) = 1.e-8 
+
+
          CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
      .       rlat,rlon,pctsrf, ftsol,
@@ -1561,5 +1579,12 @@
      .       radsol,clesphy0,
      .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
-     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon)
+     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
+     .       pbl_tke)
+
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 
        DO i=1,klon
@@ -1643,9 +1668,20 @@
          ENDIF
 
+           rugoro=0.
 c34EK
          IF (ok_orodr) THEN
-           DO i=1,klon
-             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
-           ENDDO
+
+           rugoro=0.
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
+! justement quand ok_orodr = false.
+! ce rugoro est utilise par la couche limite et fait double emploi
+! avec les paramétrisations spécifiques de Francois Lott.
+!           DO i=1,klon
+!             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
+!           ENDDO
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
            CALL SUGWD(klon,klev,paprs,pplay)
            DO i=1,klon
@@ -2061,5 +2097,5 @@
      d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
      -     dsens,     devap,     zxsnow,
-     -     zxfluxt,   zxfluxq,   q2m,     fluxq )
+     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
 c
 c
@@ -2233,4 +2269,8 @@
           DO i = 1, klon
             ema_pct(i)  = paprs(i,itop_con(i))
+            if (itop_con(i).gt.klev-3) then
+               print*,'La convection monte trop haut '
+               print*,'itop_con(,',i,',)=',itop_con(i)
+            endif
           ENDDO
           DO i = 1, klon
@@ -2321,4 +2361,6 @@
 
 
+      d_t_ajsb(:,:)=0.
+      d_q_ajsb(:,:)=0.
       d_t_ajs(:,:)=0.
       d_u_ajs(:,:)=0.
@@ -2336,13 +2378,8 @@
 c  ====
          IF(prt_level>9)WRITE(lunout,*)'pas de convection'
-      else if(iflag_thermals.eq.0) then
-
-c  Ajustement sec
-c  ==============
-         IF(prt_level>9)WRITE(lunout,*)'ajsec'
-         CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
-         t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
-         q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
+
+
       else
+
 c  Thermiques
 c  ==========
@@ -2351,4 +2388,7 @@
          print*,'JUSTE AVANT , iflag_thermals='
      s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
+
+
+         if (iflag_thermals.gt.1) then
          call calltherm(pdtphys
      s      ,pplay,paprs,pphi,weak_inversion
@@ -2357,4 +2397,5 @@
      s      ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth,
      s       ratqsdiff,zqsatth)
+         endif
 
 !        call calltherm(pdtphys
@@ -2363,4 +2404,42 @@
 !    s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
 !    s      ,fm_therm,entr_therm)
+
+c  Ajustement sec
+c  ==============
+
+! Dans le cas où on active les thermiques, on fait partir l'ajustement
+! a partir du sommet des thermiques.
+! Dans le cas contraire, on demarre au niveau 1.
+
+         if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
+
+         if(iflag_thermals.eq.0) then
+            IF(prt_level>9)WRITE(lunout,*)'ajsec'
+            limbas(:)=1
+         else
+            limbas(:)=lmax_th(:)
+         endif
+  
+! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
+! pour des test de convergence numerique.
+! Le nouveau ajsec est a priori mieux, meme pour le cas 
+! iflag_thermals = 0 (l'ancienne version peut faire des tendances
+! non nulles numeriquement pour des mailles non concernees.
+
+         if (iflag_thermals.eq.0) then
+            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
+     s      , d_t_ajsb, d_q_ajsb)
+         else
+            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
+     s      , d_t_ajsb, d_q_ajsb)
+         endif
+
+         t_seri(:,:) = t_seri(:,:) + d_t_ajsb(:,:)
+         q_seri(:,:) = q_seri(:,:) + d_q_ajsb(:,:)
+         d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
+         d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
+
+         endif
+
       endif
 c
@@ -2397,5 +2476,5 @@
          ptconvth(:,:)=.false.
          ratqsc(:,:)=0.
-         print*,'avant clouds_gno'
+         print*,'avant clouds_gno thermique'
          call clouds_gno
      s   (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
@@ -3321,5 +3400,6 @@
      .      radsol,
      .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
-     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon)
+     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon,
+     .      pbl_tke)
       ENDIF
       
Index: LMDZ4/trunk/libf/phytherm/thermcell_closure.F90
===================================================================
--- LMDZ4/trunk/libf/phytherm/thermcell_closure.F90	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/thermcell_closure.F90	(revision 852)
@@ -18,4 +18,5 @@
       REAL zmax(ngrid),zmax_sec(ngrid)
       REAL wmax(ngrid),wmax_sec(ngrid)
+      real zdenom
 
       REAL alim_star2(ngrid)
@@ -35,4 +36,16 @@
      &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
              enddo
+             zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig)
+             if (zdenom<1.e-14) then
+                print*,'ig=',ig
+                print*,'alim_star2',alim_star2(ig)
+                print*,'zmax',zmax(ig)
+                print*,'r_aspect',r_aspect
+                print*,'zdenom',zdenom
+                print*,'alim_star',alim_star(ig,:)
+                print*,'zmax_sec',zmax_sec(ig)
+                print*,'wmax_sec',wmax_sec(ig)
+                stop
+             endif
              if ((zmax_sec(ig).gt.1.e-10).and.(1.eq.1)) then 
              f(ig)=wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
@@ -41,5 +54,5 @@
      &                     zmax_sec(ig))*wmax_sec(ig))
              else
-             f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect*alim_star2(ig))
+             f(ig)=wmax(ig)/zdenom
             f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
      &                     zmax(ig))*wmax(ig))
Index: LMDZ4/trunk/libf/phytherm/thermcell_dry.F90
===================================================================
--- LMDZ4/trunk/libf/phytherm/thermcell_dry.F90	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/thermcell_dry.F90	(revision 852)
@@ -1,4 +1,4 @@
-       SUBROUTINE thermcell_dry(ngrid,nlay,zlev,pphi,ztv,entr_star,  &
-     &                            lentr,lmin,zmax,wmax,lev_out)
+       SUBROUTINE thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
+     &                            lalim,lmin,zmax,wmax,lev_out)
 
 !--------------------------------------------------------------------------
@@ -13,6 +13,6 @@
        REAL pphi(ngrid,nlay)
        REAl ztv(ngrid,nlay)
-       REAL entr_star(ngrid,nlay)
-       INTEGER lentr(ngrid)
+       REAL alim_star(ngrid,nlay)
+       INTEGER lalim(ngrid)
       integer lev_out                           ! niveau pour les print
 
@@ -33,5 +33,4 @@
           do l=1,nlay+1
              zw2(ig,l)=0.
-             f_star(ig,l)=0.
              wa_moy(ig,l)=0.
           enddo
@@ -47,18 +46,68 @@
        enddo
 !calcul de la vitesse a partir de la CAPE en melangeant thetav
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! A eliminer
+! Ce if complique etait fait pour reperer la premiere couche instable
+! Ici, c'est lmin.
+!
+!       do l=1,nlay-2
+!         do ig=1,ngrid
+!            if (ztv(ig,l).gt.ztv(ig,l+1)  &
+!     &         .and.alim_star(ig,l).gt.1.e-10  &
+!     &         .and.zw2(ig,l).lt.1e-10) then
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+! Calcul des F^*, integrale verticale de E^*
+       f_star(:,1)=0.
+       do l=1,nlay
+          f_star(:,l+1)=f_star(:,l)+alim_star(:,l)
+       enddo
+
+! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise
+       linter(:)=0.
+
+! couche la plus haute concernee par le thermique. 
+       lmax(:)=1
+
+! Le niveau linter est une variable continue qui se trouve dans la couche
+! lmax
+
        do l=1,nlay-2
          do ig=1,ngrid
-            if (ztv(ig,l).gt.ztv(ig,l+1)  &
-     &         .and.entr_star(ig,l).gt.1.e-10  &
-     &         .and.zw2(ig,l).lt.1e-10) then
-               f_star(ig,l+1)=entr_star(ig,l)
-!
+            if (l.eq.lmin(ig).and.lalim(ig).gt.1) then
+
+!------------------------------------------------------------------------
+!  Calcul de la vitesse en haut de la premiere couche instable.
+!  Premiere couche du panache thermique
+!------------------------------------------------------------------------
                zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
      &                     *(zlev(ig,l+1)-zlev(ig,l))  &
      &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
-            else if ((zw2(ig,l).ge.1e-10).and.  &
-     &               (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
-               f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
-               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)  &
+
+!------------------------------------------------------------------------
+! Tant que la vitesse en bas de la couche et la somme du flux de masse
+! et de l'entrainement (c'est a dire le flux de masse en haut) sont
+! positifs, on calcul
+! 1. le flux de masse en haut  f_star(ig,l+1)
+! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l)
+! 3. la vitesse au carré en haut zw2(ig,l+1)
+!------------------------------------------------------------------------
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  A eliminer : dans cette version, si zw2 est > 0 on a un therique.
+!  et donc, au dessus, f_star(ig,l+1) est forcement suffisamment 
+!  grand puisque on n'a pas de detrainement.
+!  f_star est une fonction croissante.
+!  c'est donc vraiment sur zw2 uniquement qu'il faut faire le test.
+!           else if ((zw2(ig,l).ge.1e-10).and.  &
+!    &               (f_star(ig,l)+alim_star(ig,l).gt.1.e-10)) then
+!              f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+            else if (zw2(ig,l).ge.1e-10) then
+
+               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+alim_star(ig,l)  &
      &                    *ztv(ig,l))/f_star(ig,l+1)
                zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+  &
@@ -67,11 +116,22 @@
             endif
 ! determination de zmax continu par interpolation lineaire
+!------------------------------------------------------------------------
+
+            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
+!               stop'On tombe sur le cas particulier de thermcell_dry'
+                print*,'On tombe sur le cas particulier de thermcell_dry'
+                zw2(ig,l+1)=0.
+                linter(ig)=l+1
+                lmax(ig)=l
+            endif
+
             if (zw2(ig,l+1).lt.0.) then
                linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
      &           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
                zw2(ig,l+1)=0.
-            else
+               lmax(ig)=l
+            endif
+
                wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
-            endif
 
             if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
@@ -84,15 +144,20 @@
        if (lev_out.ge.1) print*,'fin calcul zw2'
 !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! A eliminer :
+! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut
 ! Calcul de la couche correspondant a la hauteur du thermique
-      do ig=1,ngrid
-         lmax(ig)=lentr(ig)
-      enddo
-      do ig=1,ngrid
-         do l=nlay,lentr(ig)+1,-1
-            if (zw2(ig,l).le.1.e-10) then
-               lmax(ig)=l-1
-            endif
-         enddo
-      enddo
+!      do ig=1,ngrid
+!         lmax(ig)=lalim(ig)
+!      enddo
+!      do ig=1,ngrid
+!         do l=nlay,lalim(ig)+1,-1
+!            if (zw2(ig,l).le.1.e-10) then
+!               lmax(ig)=l-1
+!            endif
+!         enddo
+!      enddo
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 !    
 ! Determination de zw2 max
@@ -119,10 +184,25 @@
       do  ig=1,ngrid
 ! calcul de zlevinter
-          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
-     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
-     &    -zlev(ig,lmax(ig)))
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH A eliminer
+! Simplification
+!          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
+!     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
+!     &    -zlev(ig,lmax(ig)))
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+          zlevinter(ig)=zlev(ig,lmax(ig)) + &
+     &    (linter(ig)-lmax(ig))*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
            zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
       enddo
-!on stoppe après les calculs de zmax et wmax
+
+! Verification que lalim<=lmax
+      do ig=1,ngrid
+         if(lalim(ig)>lmax(ig)) then
+            print*,'WARNING thermcell_dry ig=',ig,'  lalim=',lalim(ig),'  lmax(ig)=',lmax(ig)
+            lmax(ig)=lalim(ig)
+         endif
+      enddo
       
       RETURN
Index: LMDZ4/trunk/libf/phytherm/thermcell_flux.F90
===================================================================
--- LMDZ4/trunk/libf/phytherm/thermcell_flux.F90	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/thermcell_flux.F90	(revision 852)
@@ -28,9 +28,10 @@
 
       integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
-      integer ncorecfm4,ncorecfm5,ncorecfm6
+      integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
       
 
       REAL entr(ngrid,klev),detr(ngrid,klev)
       REAL fm(ngrid,klev+1)
+      REAL zfm
 
       integer igout
@@ -40,8 +41,9 @@
       REAL f_old,ddd0,eee0,ddd,eee,zzz
 
-      REAL fracmax
-      save fracmax
-
-      fracmax=0.5
+      REAL fomass_max,alphamax
+      save fomass_max,alphamax
+
+      fomass_max=0.5
+      alphamax=0.7
 
       ncorecfm1=0
@@ -51,4 +53,6 @@
       ncorecfm5=0
       ncorecfm6=0
+      ncorecfm7=0
+      ncorecfm8=0
       ncorecalpha=0
 
@@ -56,4 +60,15 @@
       fm(:,:)=0.
       
+      if (lev_out.ge.10) then
+         write(lunout,*) 'Dans thermcell_flux 0'
+         write(lunout,*) 'flux base ',f(igout)
+         write(lunout,*) 'lmax ',lmax(igout)
+         write(lunout,*) 'lalim ',lalim(igout)
+         write(lunout,*) 'ig= ',igout
+         write(lunout,*) ' l E*    A*     D*  '
+         write(lunout,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) &
+     &    ,l=1,lmax(igout))
+      endif
+
 
 !-------------------------------------------------------------------------
@@ -66,5 +81,13 @@
             if (l.le.lmax(ig)) then
                if (entr_star(ig,l).gt.1.) then
-                    print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
+                    print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig)
+                    print*,'entr_star(ig,l)',entr_star(ig,l)
+                    print*,'alim_star(ig,l)',alim_star(ig,l)
+                    print*,'detr_star(ig,l)',detr_star(ig,l)
+!                   stop
+               endif
+            else
+               if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then
+                    print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig)
                     print*,'entr_star(ig,l)',entr_star(ig,l)
                     print*,'alim_star(ig,l)',alim_star(ig,l)
@@ -72,12 +95,4 @@
                     stop
                endif
-            else
-               if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0) then
-                    print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
-                    print*,'entr_star(ig,l)',entr_star(ig,l)
-                    print*,'alim_star(ig,l)',alim_star(ig,l)
-                    print*,'detr_star(ig,l)',detr_star(ig,l)
-                    stop
-               endif
             endif
          enddo
@@ -92,4 +107,15 @@
          detr(:,l)=f(:)*detr_star(:,l)
       enddo
+
+      if (lev_out.ge.10) then
+         write(lunout,*) 'Dans thermcell_flux 1'
+         write(lunout,*) 'flux base ',f(igout)
+         write(lunout,*) 'lmax ',lmax(igout)
+         write(lunout,*) 'lalim ',lalim(igout)
+         write(lunout,*) 'ig= ',igout
+         write(lunout,*) ' l   E    D     W2'
+         write(lunout,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) &
+     &    ,zw2(igout,l+1),l=1,lmax(igout))
+      endif
 
       fm(:,1)=0.
@@ -107,12 +133,22 @@
       enddo
 
-      if (lev_out.ge.10) then
-         write(lunout,*) 'Dans thermcell_flux 1'
-         write(lunout,*) 'flux base ',f(igout)
-         write(lunout,*) 'lmax ',lmax(igout)
-         write(lunout,*) 'lalim ',lalim(igout)
-         write(lunout,'(i6,i4,3e15.5)') (ig,l,entr(igout,l),detr(igout,l) &
-     &    ,fm(igout,l+1),l=1,lmax(igout))
-      endif
+
+
+! Test provisoire : pour comprendre pourquoi on corrige plein de fois 
+! le cas fm6, on commence par regarder une premiere fois avant les
+! autres corrections.
+
+      do l=1,klev
+         do ig=1,ngrid
+            if (detr(ig,l).gt.fm(ig,l)) then
+               ncorecfm8=ncorecfm8+1
+!              igout=ig
+            endif
+         enddo
+      enddo
+
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'2  ')
 
 !-------------------------------------------------------------------------
@@ -130,4 +166,8 @@
          enddo
       enddo
+
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'3  ')
 
 !-------------------------------------------------------------------------
@@ -151,4 +191,8 @@
       enddo
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'4  ')
+
 
 !-------------------------------------------------------------------------
@@ -167,4 +211,8 @@
       enddo
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'5  ')
+
 !-------------------------------------------------------------------------
 !detr ne peut pas etre superieur a fm
@@ -172,4 +220,5 @@
 
       if(1.eq.1) then
+
       do l=1,klev
          do ig=1,ngrid
@@ -181,6 +230,23 @@
                ncorecfm6=ncorecfm6+1
                detr(ig,l)=fm(ig,l)
-               entr(ig,l)=fm(ig,l+1)
-            endif
+!              entr(ig,l)=fm(ig,l+1)
+
+! Dans le cas ou on est au dessus de la couche d'alimentation et que le
+! detrainement est plus fort que le flux de masse, on stope le thermique.
+               if (l.gt.lalim(ig)) then
+                  lmax(ig)=l
+                  fm(ig,l+1)=0.
+                  entr(ig,l)=0.
+               else
+                  ncorecfm7=ncorecfm7+1
+               endif
+            endif
+
+            if(l.gt.lmax(ig)) then
+               detr(ig,l)=0.
+               fm(ig,l+1)=0.
+               entr(ig,l)=0.
+            endif
+
             if (entr(ig,l).lt.0.) then
                print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
@@ -194,4 +260,8 @@
 
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'6  ')
+
 !-------------------------------------------------------------------------
 !fm ne peut pas etre negatif
@@ -207,5 +277,5 @@
             endif
             if (detr(ig,l).lt.0.) then
-               print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
+               print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig)
                print*,'detr(ig,l)',detr(ig,l)
                print*,'fm(ig,l)',fm(ig,l)
@@ -215,19 +285,41 @@
      enddo
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'7  ')
+
 !-----------------------------------------------------------------------
 !la fraction couverte ne peut pas etre superieure a 1            
 !-----------------------------------------------------------------------
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH Partie a revisiter.
+! Il semble qu'etaient codees ici deux optiques dans le cas
+! F/ (rho *w) > 1
+! soit limiter la hauteur du thermique en considerant que c'est 
+! la derniere chouche, soit limiter F a rho w.
+! Dans le second cas, il faut en fait limiter a un peu moins
+! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
+! dans thermcell_main et qu'il semble de toutes facons deraisonable
+! d'avoir des fractions de 1..
+! Ci dessous, et dans l'etat actuel, le premier des  deux if est
+! sans doute inutile.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+     if(1.eq.0) then
+
      do l=1,klev
         do ig=1,ngrid
            if (zw2(ig,l+1).gt.1.e-10) then
-           if ((((fm(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt.  &
-     &      1.)) then
+           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
+           if ( fm(ig,l+1) .gt. zfm) then
               f_old=fm(ig,l+1)
-              fm(ig,l+1)=rhobarz(ig,l+1)*zw2(ig,l+1)
-              zw2(ig,l+1)=0.
-              zqla(ig,l+1)=0.
+              fm(ig,l+1)=zfm
+!             zw2(ig,l+1)=0.
+!             zqla(ig,l+1)=0.
               detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
-              lmax(ig)=l+1
-              zmax(ig)=zlev(ig,lmax(ig))
+!             lmax(ig)=l+1
+!             zmax(ig)=zlev(ig,lmax(ig))
 !             print*,'alpha>1',l+1,lmax(ig)
               ncorecalpha=ncorecalpha+1
@@ -237,4 +329,9 @@
      enddo
 !
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'8  ')
+
+     endif
 
 !-----------------------------------------------------------------------
@@ -248,5 +345,5 @@
             eee0=entr(ig,l)
             ddd0=detr(ig,l)
-            eee=entr(ig,l)-masse(ig,l)*fracmax/ptimestep
+            eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
             ddd=detr(ig,l)-eee
             if (eee.gt.0.) then
@@ -272,6 +369,6 @@
                          print*,'detr',detr(ig,l)
                          print*,'masse',masse(ig,l)
-                         print*,'fracmax',fracmax
-                         print*,'masse(ig,l)*fracmax/ptimestep',masse(ig,l)*fracmax/ptimestep
+                         print*,'fomass_max',fomass_max
+                         print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
                          print*,'ptimestep',ptimestep
                          print*,'lmax(ig)',lmax(ig)
@@ -309,8 +406,57 @@
     &     ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
     &     ncorecfm6,'x fm6', &
+    &     ncorecfm7,'x fm7', &
+    &     ncorecfm8,'x fm8', &
     &     ncorecalpha,'x alpha'
       endif
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'fin')
+
 
       return
       end
+
+      subroutine printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,descr)
+
+     implicit none
+
+      integer ngrid,klev,lunout,igout,l,lm
+
+      integer lmax(klev),lalim(klev)
+      real ptimestep,masse(ngrid,klev),entr(ngrid,klev),detr(ngrid,klev)
+      real fm(ngrid,klev+1),f(ngrid)
+
+      character*3 descr
+
+      lm=lmax(igout)+5
+      if(lm.gt.klev) lm=klev
+
+      print*,'Impression jusque lm=',lm
+
+         write(lunout,*) 'Dans thermcell_flux '//descr
+         write(lunout,*) 'flux base ',f(igout)
+         write(lunout,*) 'lmax ',lmax(igout)
+         write(lunout,*) 'lalim ',lalim(igout)
+         write(lunout,*) 'ig= ',igout
+         write(lunout,'(a3,4a14)') 'l','M','E','D','F'
+         write(lunout,'(i4,4e14.4)') (l,masse(igout,l)/ptimestep, &
+     &     entr(igout,l),detr(igout,l) &
+     &    ,fm(igout,l+1),l=1,lm)
+
+
+      do l=lmax(igout)+1,klev
+          if (abs(entr(igout,l))+abs(detr(igout,l))+abs(fm(igout,l)).gt.0.) then
+          print*,'cas 1 : igout,l,lmax(igout)',igout,l,lmax(igout)
+          print*,'entr(igout,l)',entr(igout,l)
+          print*,'detr(igout,l)',detr(igout,l)
+          print*,'fm(igout,l)',fm(igout,l)
+          stop
+          endif
+      enddo
+
+      return
+      end
+
Index: LMDZ4/trunk/libf/phytherm/thermcell_main.F90
===================================================================
--- LMDZ4/trunk/libf/phytherm/thermcell_main.F90	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/thermcell_main.F90	(revision 852)
@@ -1,2 +1,5 @@
+!
+! $Header$
+!
       SUBROUTINE thermcell_main(ngrid,nlay,ptimestep  &
      &                  ,pplay,pplev,pphi,debut  &
@@ -55,7 +58,8 @@
 !   ------
 
-      integer,save :: igout=871
+!      integer,save :: igout=4259
+      integer,save :: igout=2928
       integer,save :: lunout=6
-      integer,save :: lev_out=0
+      integer,save :: lev_out=10
 
       INTEGER ig,k,l,ll
@@ -129,4 +133,5 @@
       real zlevinter(klon)
       logical debut
+       real seuil
 
 !
@@ -142,4 +147,7 @@
 !   ---------------
 !
+
+      seuil=0.25
+
       if (lev_out.ge.1) print*,'thermcell_main V4'
 
@@ -228,33 +236,59 @@
 
 !------------------------------------------------------------------
-!   Calcul de w2, carre de w a partir de la cape
-!
-!   Indicages:
-!   l'ascendance provenant du niveau k traverse l'interface l avec
-!   une vitesse wa(k,l).
-!
-!                       --------------------
-!
-!                       + + + + + + + + + + 
-!
-!  wa(k,l)   ----       --------------------    l
-!             /\
-!            /||\       + + + + + + + + + + 
-!             ||
-!             ||        --------------------
-!             ||
-!             ||        + + + + + + + + + + 
-!             ||
-!             ||        --------------------
-!             ||__
-!             |___      + + + + + + + + + +     k
-!
-!                       --------------------
-!
-!
-!
+!
+!             /|\
+!    --------  |  F_k+1 -------   
+!                              ----> D_k
+!             /|\              <---- E_k , A_k
+!    --------  |  F_k --------- 
+!                              ----> D_k-1
+!                              <---- E_k-1 , A_k-1
+!
+!
+!
+!
+!
+!    ---------------------------
+!
+!    ----- F_lmax+1=0 ----------         \
+!            lmax     (zmax)              |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------          |
+!                                         |  E
+!    ---------------------------          |  D
+!                                         |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------  \       |
+!            lalim                 |      |
+!    ---------------------------   |      |
+!                                  |      |
+!    ---------------------------   |      |
+!                                  | A    |
+!    ---------------------------   |      |
+!                                  |      |
+!    ---------------------------   |      |
+!    lmin  (=1 pour le moment)     |      |
+!    ----- F_lmin=0 ------------  /      /
+!
+!    ---------------------------
+!    //////////////////////////
+!
+!
+!=============================================================================
+!  Calculs initiaux ne faisant pas intervenir les changements de phase
+!=============================================================================
+
 !------------------------------------------------------------------
-!definition du profil d alimentation a partir de la flottabilite: 
-!alim_star, alim_star_tot, lalim et lmin
+!  1. alim_star est le profil vertical de l'alimentation à la base du
+!     panache thermique, calculé à partir de la flotabilité de l'air sec
+!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
 !------------------------------------------------------------------
 !
@@ -262,4 +296,8 @@
       CALL thermcell_init(ngrid,nlay,ztv,zlev,  &
      &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
+
+call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin  ')
+call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
+
 
       if (lev_out.ge.1) print*,'thermcell_main apres thermcell_init'
@@ -268,14 +306,31 @@
          write(lunout,*) 'lmin ',lmin(igout)
          write(lunout,*) 'lalim ',lalim(igout)
+         write(lunout,*) ' ig l alim_star thetav'
+         write(lunout,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
+     &   ,ztv(igout,l),l=1,lalim(igout)+4)
       endif
 
-!-----------------------------------------------------------------------------------
-!calcul des caracteristiques du thermique sec pour le calcul de detr et la fermeture
-!-----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
+!     panache sec conservatif (e=d=0) alimente selon alim_star 
+!     Il s'agit d'un calcul de type CAPE
+!     zmax_sec est utilisé pour déterminer la géométrie du thermique.
+!------------------------------------------------------------------------------
 !
       CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
      &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
 
+call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
+call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
+
       if (lev_out.ge.1) print*,'thermcell_main apres thermcell_dry'
+      if (lev_out.ge.10) then
+         write(lunout,*) 'Dans thermcell_main 1b'
+         write(lunout,*) 'lmin ',lmin(igout)
+         write(lunout,*) 'lalim ',lalim(igout)
+         write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
+         write(lunout,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
+     &    ,l=1,lalim(igout)+4)
+      endif
 
 
@@ -289,4 +344,6 @@
      &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
      &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
+      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
+      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
 
       if (lev_out.ge.1) print*,'thermcell_main apres thermcell_plume'
@@ -297,5 +354,5 @@
          write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
          write(lunout,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
-     &    ,f_star(igout,l+1),l=1,lmax(igout))
+     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
       endif
 
@@ -307,4 +364,10 @@
      &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
 
+
+      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
+      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
+      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
+      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
+
       if (lev_out.ge.1) print*,'thermcell_main apres thermcell_height'
 
@@ -322,5 +385,5 @@
 !-------------------------------------------------------------------------------
 
-      CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
+      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
      &       lalim,lmax,alim_star,  &
      &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
@@ -328,4 +391,6 @@
 
       if (lev_out.ge.1) print*,'thermcell_main apres thermcell_flux'
+      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
+      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
 
 !c------------------------------------------------------------------
@@ -391,4 +456,5 @@
       enddo
      
+      print*,'14a OK convect8'
 ! calcul du niveau de condensation
 ! initialisation
@@ -410,4 +476,5 @@
          enddo
       enddo
+      print*,'14b OK convect8'
       do k=nlay,1,-1
          do ig=1,ngrid
@@ -418,4 +485,5 @@
          enddo
       enddo
+      print*,'14c OK convect8'
 !calcul des moments
 !initialisation
@@ -429,4 +497,5 @@
          enddo
       enddo      
+      print*,'14d OK convect8'
       do l=1,nlay
          do ig=1,ngrid
@@ -440,8 +509,9 @@
             q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
 !test: on calcul q2/po=ratqsc
-            ratqscth(ig,l)=sqrt(q2(ig,l))/(po(ig,l)*1000.)
+            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
          enddo
       enddo
 !calcul du ratqscdiff
+      print*,'14e OK convect8'
       var=0.
       vardiff=0.
@@ -452,4 +522,5 @@
          enddo
       enddo
+      print*,'14f OK convect8'
       do ig=1,ngrid
           do l=1,lalim(ig)
@@ -462,4 +533,5 @@
           enddo
       enddo
+      print*,'14g OK convect8'
       do l=1,nlay
          do ig=1,ngrid
@@ -496,2 +568,34 @@
 
 !-----------------------------------------------------------------------------
+
+      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
+
+      integer klon,klev
+      real pplev(klon,klev+1),pplay(klon,klev)
+      real ztv(klon,klev)
+      real po(klon,klev)
+      real ztva(klon,klev)
+      real zqla(klon,klev)
+      real f_star(klon,klev)
+      real zw2(klon,klev)
+      integer long(klon)
+      real seuil
+      character*21 comment
+
+      print*,'TEST ',comment
+!  test sur la hauteur des thermiques ...
+         do i=1,klon
+            if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
+               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
+               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
+               do k=1,klev
+                  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
+               enddo
+!              stop
+            endif
+         enddo
+
+
+      return
+      end
+
Index: LMDZ4/trunk/libf/phytherm/thermcell_plume.F90
===================================================================
--- LMDZ4/trunk/libf/phytherm/thermcell_plume.F90	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/thermcell_plume.F90	(revision 852)
@@ -1,5 +1,5 @@
       SUBROUTINE thermcell_plume(ngrid,klev,ztv,zthl,po,zl,rhobarz,  &
      &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
-     &           lentr,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
+     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
      &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
 
@@ -29,5 +29,5 @@
       REAL l_mix
       REAL r_aspect
-      INTEGER lentr(ngrid)
+      INTEGER lalim(ngrid)
       integer lev_out                           ! niveau pour les print
 
@@ -133,10 +133,25 @@
 !tests sur la definition du detr
 !calcul de detr_star et entr_star
-                w_est(ig,3)=zw2(ig,2)*  &
-     &                   ((f_star(ig,2))**2)  &
-     &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
-     &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
-     &                   *(zlev(ig,3)-zlev(ig,2))
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH le test miraculeux de Catherine ? Le bout du tunel ?
+!               w_est(ig,3)=zw2(ig,2)*  &
+!    &                   ((f_star(ig,2))**2)  &
+!    &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
+!    &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
+!    &                   *(zlev(ig,3)-zlev(ig,2))
+!     if (l.gt.2) then
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
       if (l.gt.2) then
+          w_est(ig,3)=zw2(ig,2)* &
+     &      ((f_star(ig,2))**2) &
+     &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
+     &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
+     &      *(zlev(ig,3)-zlev(ig,2))
+
+
 !---------------------------------------------------------------------------
 !calcul de l entrainement et du detrainement lateral
@@ -217,23 +232,57 @@
 !calcul de entr_star
 !
-        if (detr_star(ig,l).gt.f_star(ig,l)) then
-           detr_star(ig,l)=f_star(ig,l)
-!a decommenter ou pas?
-!           entr_star(ig,l)=0.
-        endif
-
 ! Deplacement du calcul de entr_star pour eviter d'avoir aussi
 ! entr_star > fstar.
+! Redeplacer suite a la transformation du cas detr>f
 ! FH
-          entr_star(ig,l)=0.4*detr_star(ig,l)
-!
+        if(l.gt.lalim(ig)) then
+         entr_star(ig,l)=0.4*detr_star(ig,l)
+        else
+
+! FH :
+! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
+! en haut de la couche d'alimentation.
+! A remettre en questoin a la premiere occasion mais ca peut aider a 
+! ecrire un code robuste.
+! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
+! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
+! d'alimentation, ce qui n'est pas forcement heureux.
+         entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
+         detr_star(ig,l)=entr_star(ig,l)
+        endif
+
+!
+        if (detr_star(ig,l).gt.f_star(ig,l)) then
+
+!  Ce test est là entre autres parce qu'on passe par des valeurs
+!  delirantes de detr_star.
+!  ca vaut sans doute le coup de verifier pourquoi.
+
+           detr_star(ig,l)=f_star(ig,l)
+           if (l.gt.lalim(ig)+1) then
+               entr_star(ig,l)=0.
+               alim_star(ig,l)=0.
+! FH ajout pour forcer a stoper le thermique juste sous le sommet
+! de la couche (voir calcul de finter)
+               zw2(ig,l+1)=-1.e-10
+               linter(ig)=l+1
+            else
+               entr_star(ig,l)=detr_star(ig,l)
+            endif
+        endif
+
       else
          detr_star(ig,l)=0.
          entr_star(ig,l)=0.
       endif
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH inutile si on conserve comme on l'a fait plus haut entr=detr
+! dans la couche d'alimentation
 !pas d entrainement dans la couche alim
-      if ((l.lt.lentr(ig))) then
-           entr_star(ig,l)=0.
-      endif
+!      if ((l.le.lalim(ig))) then
+!           entr_star(ig,l)=0.
+!      endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !
 !prise en compte du detrainement et de l entrainement dans le calcul du flux
@@ -308,11 +357,21 @@
 !
 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max 
+
+            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
+!               stop'On tombe sur le cas particulier de thermcell_dry'
+                print*,'On tombe sur le cas particulier de thermcell_plume'
+                zw2(ig,l+1)=0.
+                linter(ig)=l+1
+            endif
+
+
         if (zw2(ig,l+1).lt.0.) then 
            linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
      &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
            zw2(ig,l+1)=0.
-        else
+        endif
+
            wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 
-        endif
+
         if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
 !   lmix est le niveau de la couche ou w (wa_moy) est maximum
Index: LMDZ4/trunk/libf/phytherm/write_histday.h
===================================================================
--- LMDZ4/trunk/libf/phytherm/write_histday.h	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/write_histday.h	(revision 852)
@@ -491,4 +491,14 @@
      $                     zx_tmp_fi2d) 
 C 
+! FH Sorties specifiques pour Mellor et Yamada
+      if (iflag_pbl>1 .and. lev_histday.gt.10 ) then
+
+        CALL histwrite_phy(nid_day,"tke_"//clnsurf(nsrf),itau_w,
+     $      pbl_tke(:,1:klev,nsrf))
+
+        CALL histwrite_phy(nid_day,"tke_max_"//clnsurf(nsrf),itau_w,
+     $      pbl_tke(:,1:klev,nsrf))
+	endif
+
       END DO  
 c=================================================================
Index: LMDZ4/trunk/libf/phytherm/write_histmth.h
===================================================================
--- LMDZ4/trunk/libf/phytherm/write_histmth.h	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/write_histmth.h	(revision 852)
@@ -681,4 +681,22 @@
       IF(lev_histmth.GE.4) THEN
 c
+c  FH Sorties pour la couche limite
+      CALL histwrite_phy(nid_mth,"kz",itau_w,ycoefh)
+      CALL histwrite_phy(nid_mth,"kz_max",itau_w,ycoefh)
+
+      if(iflag_pbl>1) then
+      zx_tmp_fi3d=0.
+      do nsrf=1,nbsrf
+         do k=1,klev
+          zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k)
+     ,    +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
+         enddo
+      enddo
+      CALL histwrite_phy(nid_mth,"tke",itau_w,zx_tmp_fi3d)
+      CALL histwrite_phy(nid_mth,"tke_max",itau_w,zx_tmp_fi3d)
+      endif
+
+
+
 cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, clwcon0, zx_tmp_3d)
       CALL histwrite_phy(nid_mth,"clwcon",itau_w,clwcon0)
@@ -757,11 +775,20 @@
 c
 cIM: 101003 : K/30min ==> K/s
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys
+cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
+      CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d)
+c
+      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys
+cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
+      CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d)
+c
+cIM: 101003 : K/30min ==> K/s
       zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys
 cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
-      CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d)
+      CALL histwrite_phy(nid_mth,"dtthe",itau_w,zx_tmp_fi3d)
 c
       zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys
 cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
-      CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d)
+      CALL histwrite_phy(nid_mth,"dqthe",itau_w,zx_tmp_fi3d)
 c
 cIM: 101003 : K/day ==> K/s
@@ -933,13 +960,26 @@
       CALL histwrite_phy(nid_mth,"dtcon",itau_w,zx_tmp_fi3d)
 c
-c temperature tendency due to dry convective processes
-       DO l=1, klev
-       DO i=1, klon
-        zx_tmp_fi3d(i,l)=d_t_ajs(i,l)/pdtphys
-       ENDDO !i
-       ENDDO !l
-c
-cym      CALL gr_fi_ecrit(klev, klon,iim,jjmp1, zx_tmp_fi3d,zx_tmp_3d)
+
+
+
+cIM: 101003 : K/30min ==> K/s
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys
+cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d)
+c
+      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys
+cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
+      CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d)
+c
+cIM: 101003 : K/30min ==> K/s
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys
+cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
+      CALL histwrite_phy(nid_mth,"dtthe",itau_w,zx_tmp_fi3d)
+c
+      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys
+cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
+      CALL histwrite_phy(nid_mth,"dqthe",itau_w,zx_tmp_fi3d)
+c
+
 c
 c  temperature tendency due to large scale precipitation
@@ -1672,4 +1712,22 @@
       IF(lev_histmth.GE.4) THEN
 c
+c  FH Sorties pour la couche limite
+      CALL histwrite_phy(nid_mth,"kz",itau_w,ycoefh)
+      CALL histwrite_phy(nid_mth,"kz_max",itau_w,ycoefh)
+
+      if(iflag_pbl>1) then
+      zx_tmp_fi3d=0.
+      do nsrf=1,nbsrf
+         do k=1,klev
+          zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k)
+     ,    +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
+         enddo
+      enddo
+      CALL histwrite_phy(nid_mth,"tke",itau_w,zx_tmp_fi3d)
+      CALL histwrite_phy(nid_mth,"tke_max",itau_w,zx_tmp_fi3d)
+      endif
+
+
+
 cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, clwcon0, zx_tmp_3d)
       CALL histwrite_phy(nid_mth,"clwcon",itau_w,clwcon0)
@@ -1739,8 +1797,4 @@
       CALL histwrite_phy(nid_mth,"ratqs",itau_w,ratqs)
 c
-      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys
-cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
-      CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d)
-c
 cIM: 101003 : K/day ==> K/s
       zx_tmp_fi3d(1:klon,1:klev)=heat(1:klon,1:klev)/RDAY 
Index: LMDZ4/trunk/libf/phytherm/yamada4.F
===================================================================
--- LMDZ4/trunk/libf/phytherm/yamada4.F	(revision 849)
+++ LMDZ4/trunk/libf/phytherm/yamada4.F	(revision 852)
@@ -458,20 +458,28 @@
 c   Diagnostique pour stokage
 
+      if(1.eq.0)then
       rino=rif
-      smyam(:,1:klev)=sm(:,1:klev)
-      styam=sm(:,1:klev)*alpha(:,1:klev)
-      lyam(1:klon,1:klev)=l(:,1:klev)
-      knyam(1:klon,1:klev)=kn(:,1:klev)
+      smyam(1:ngrid,1)=0.
+      styam(1:ngrid,1)=0.
+      lyam(1:ngrid,1)=0.
+      knyam(1:ngrid,1)=0.
+      w2yam(1:ngrid,1)=0.
+      t2yam(1:ngrid,1)=0.
+
+      smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
+      styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
+      lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
+      knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
 
 c   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
 
-        if(1.eq.0)then
-      w2yam=q2(:,1:klev)*0.24
-     s    +lyam(:,1:klev)*5.17*kn(:,1:klev)*n2(:,1:klev)
-     s   /sqrt(q2(:,1:klev))
-
-      t2yam=9.1*kn(:,1:klev)*dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))
-     s  *lyam(:,1:klev)
-	endif
+      w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24
+     s    +lyam(1:ngrid,2:klev)*5.17*kn(1:ngrid,2:klev)
+     s    *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
+
+      t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev)
+     s    *dtetadz(1:ngrid,2:klev)**2
+     s    /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
+      endif
 
 c     print*,'OKFIN'
