Index: LMDZ5/trunk/libf/phylmd/calltherm.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/calltherm.F90	(revision 1938)
+++ LMDZ5/trunk/libf/phylmd/calltherm.F90	(revision 1943)
@@ -174,10 +174,10 @@
          do isplit=1,nsplit_thermals
 
-          if (iflag_thermals.eq.1) then
-            CALL thermcell_2002(klon,klev,zdt   &
+          if (iflag_thermals>=1000) then
+            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
      &      ,pplay,paprs,pphi  &
      &      ,u_seri,v_seri,t_seri,q_seri  &
      &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
-     &      ,zfm_therm,zentr_therm  &
+     &      ,zfm_therm,zentr_therm,fraca,zw2  &
      &      ,r_aspect_thermals,30.,w2di_thermals  &
      &      ,tau_thermals)
@@ -273,8 +273,15 @@
        flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
 
-      if (iflag_thermals<=12) then
-         lmax=1
+! Calcul a posteriori du niveau max des thermiques pour les schémas qui
+! ne la sortent pas.
+      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
+         lmax(:)=1
          do k=1,klev-1
             zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
+         enddo
+         do k=1,klev-1
+            do i=1,klon
+               if (zfm_therm(i,k+1)>0.) lmax(i)=k
+            enddo
          enddo
       endif
Index: LMDZ5/trunk/libf/phylmd/hbtm.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/hbtm.F	(revision 1938)
+++ LMDZ5/trunk/libf/phylmd/hbtm.F	(revision 1943)
@@ -67,45 +67,27 @@
       INTEGER isommet
 cum      PARAMETER (isommet=klev) ! limite max sommet pbl
-      REAL vk
-      PARAMETER (vk=0.35)     ! Von Karman => passer a .41 ! cf U.Olgstrom
-      REAL ricr
-      PARAMETER (ricr=0.4)
-      REAL fak
-      PARAMETER (fak=8.5)     ! b calcul du Prandtl et de dTetas
-      REAL fakn
-      PARAMETER (fakn=7.2)    ! a
-      REAL onet
-      PARAMETER (onet=1.0/3.0)
-      REAL t_coup
-      PARAMETER(t_coup=273.15)
-      REAL zkmin
-      PARAMETER (zkmin=0.01)
-      REAL betam
-      PARAMETER (betam=15.0)  ! pour Phim / h dans la S.L stable
-      REAL betah
-      PARAMETER (betah=15.0)
-      REAL betas
-      PARAMETER (betas=5.0)   ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
-      REAL sffrac
-      PARAMETER (sffrac=0.1)  ! S.L. = z/h < .1
-      REAL binm
-      PARAMETER (binm=betam*sffrac)
-      REAL binh
-      PARAMETER (binh=betah*sffrac)
-      REAL ccon
-      PARAMETER (ccon=fak*sffrac*vk)
-c
+      REAL, PARAMETER :: vk=0.35  ! Von Karman => passer a .41 ! cf U.Olgstrom
+      REAL, PARAMETER :: ricr=0.4
+      REAL, PARAMETER :: fak=8.5     ! b calcul du Prandtl et de dTetas
+      REAL, PARAMETER :: fakn=7.2    ! a
+      REAL, PARAMETER :: onet=1.0/3.0
+      REAL, PARAMETER:: t_coup=273.15
+      REAL, PARAMETER :: zkmin=0.01
+      REAL, PARAMETER :: betam=15.0  ! pour Phim / h dans la S.L stable
+      REAL, PARAMETER :: betah=15.0
+      REAL, PARAMETER :: betas=5.0 ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
+      REAL, PARAMETER :: sffrac=0.1  ! S.L. = z/h < .1
+      REAL, PARAMETER :: usmin=1.e-12
+      REAL, PARAMETER :: binm=betam*sffrac
+      REAL, PARAMETER :: binh=betah*sffrac
+      REAL, PARAMETER :: ccon=fak*sffrac*vk
+      REAL, PARAMETER :: b1=70.,b2=20.
+      REAL, PARAMETER :: zref=2.  ! Niveau de ref a 2m peut eventuellement 
+c                              etre choisi a 10m
       REAL q_star,t_star
-      REAL b1,b2,b212,b2sr     ! Lambert correlations T' q' avec T* q*
-      PARAMETER (b1=70.,b2=20.)
+      REAL b212,b2sr     ! Lambert correlations T' q' avec T* q*
 c
       REAL z(klon,klev)
 cAM      REAL pcfm(klon,klev), pcfh(klon,klev)
-cAM
-      REAL zref
-      PARAMETER (zref=2.)    ! Niveau de ref a 2m peut eventuellement 
-c                              etre choisi a 10m
-cMA
-c
       INTEGER i, k, j
       REAL zxt
@@ -129,5 +111,5 @@
 cAM      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
 cAM      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
-      REAL obklen(klon)     ! Monin-Obukhov lengh
+      REAL unsobklen(klon)     ! Monin-Obukhov lengh
 cAM      REAL ztvd, ztvu, 
       REAL zdu2
@@ -345,5 +327,5 @@
          plcl(i) = 6000.
 c Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
-         obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))
+         unsobklen(i) = -RG*vk*heatv(i)/(t(i,1)*max(ustar(i),usmin)**3)
          trmb1(i)   = 0.
          trmb2(i)   = 0.
@@ -420,5 +402,5 @@
       DO i = 1, knon
       IF (check(i)) THEN
-        phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
+        phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet
 c ***************************************************
 c Wm ? et W* ? c'est la formule pour z/h < .1
@@ -621,6 +603,6 @@
           zxt=(Th_th(i)-zref*0.5*RG/RCPD/(1.+RVTMP2*qT_th(i)))
      .         *(1.+RETV*qT_th(i))
-          phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
-          phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
+          phiminv(i) = (1. - binm*pblh(i)*unsobklen(i))**onet
+          phihinv(i) = sqrt(1. - binh*pblh(i)*unsobklen(i))
           wm(i)      = ustar(i)*phiminv(i)
           fak2(i)    = wm(i)*pblh(i)*vk
@@ -657,5 +639,5 @@
 
             zh(i) = zmzp/pblh(i)
-            zl(i) = zmzp/obklen(i)
+            zl(i) = zmzp*unsobklen(i)
             zzh(i) = 0.
             IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
Index: LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 1938)
+++ LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 1943)
@@ -41,6 +41,4 @@
       REAL, SAVE, ALLOCATABLE :: d_u_ajs(:,:), d_v_ajs(:,:)
       !$OMP THREADPRIVATE(d_u_ajs, d_v_ajs)
-      REAL, SAVE, ALLOCATABLE :: d_u_ajsb(:,:), d_v_ajsb(:,:)
-      !$OMP THREADPRIVATE(d_u_ajsb, d_v_ajsb)
       REAL, SAVE, ALLOCATABLE :: d_t_eva(:,:),d_q_eva(:,:)
       !$OMP THREADPRIVATE(d_t_eva,d_q_eva)
@@ -303,5 +301,4 @@
       allocate(d_t_ajs(klon,klev),d_q_ajs(klon,klev))
       allocate(d_u_ajs(klon,klev),d_v_ajs(klon,klev))
-      allocate(d_u_ajsb(klon,klev),d_v_ajsb(klon,klev))
       allocate(d_t_eva(klon,klev),d_q_eva(klon,klev))
       allocate(d_t_lscst(klon,klev),d_q_lscst(klon,klev))
@@ -456,5 +453,4 @@
       deallocate(d_t_ajs,d_q_ajs)
       deallocate(d_u_ajs,d_v_ajs)
-      deallocate(d_u_ajsb,d_v_ajsb)
       deallocate(d_t_eva,d_q_eva)
       deallocate(d_t_lscst,d_q_lscst)
Index: LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 1938)
+++ LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 1943)
@@ -206,5 +206,5 @@
          rneb, rnebjn, zx_rh, frugs, agesno, d_t_dyn, d_q_dyn, &
          d_u_dyn, d_v_dyn, d_t_con, d_t_ajsb, d_t_ajs, &
-         d_u_ajsb, d_u_ajs, d_v_ajsb, d_v_ajs, &
+         d_u_ajs, d_v_ajs, &
          d_u_con, d_v_con, d_q_con, d_q_ajs, d_t_lsc, &
          d_t_eva, d_q_lsc, beta_prec, d_t_lscth, &
@@ -1046,11 +1046,9 @@
        CALL histwrite_phy(o_dtthe, zx_tmp_fi3d)
        IF (vars_defined) THEN
-          zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys - &
-               d_u_ajsb(1:klon,1:klev)/pdtphys
+          zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys
        ENDIF
        CALL histwrite_phy(o_duthe, zx_tmp_fi3d) 
        IF (vars_defined) THEN
-          zx_tmp_fi3d(1:klon,1:klev)=d_v_ajs(1:klon,1:klev)/pdtphys - &
-               d_v_ajsb(1:klon,1:klev)/pdtphys
+          zx_tmp_fi3d(1:klon,1:klev)=d_v_ajs(1:klon,1:klev)/pdtphys
        ENDIF
        CALL histwrite_phy(o_dvthe, zx_tmp_fi3d)
Index: LMDZ5/trunk/libf/phylmd/thermcell_dq.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/thermcell_dq.F90	(revision 1938)
+++ LMDZ5/trunk/libf/phylmd/thermcell_dq.F90	(revision 1943)
@@ -53,6 +53,6 @@
             cfl=max(cfl,fm(ig,k)/zzm)
             if (entr(ig,k).gt.zzm) then
-               print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
-               abort_message = ''
+               print*,'entr*dt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k)
+               abort_message = 'entr dt > m, 1st'
                CALL abort_gcm (modname,abort_message,1)
             endif
@@ -188,6 +188,6 @@
             cfl=max(cfl,fm(ig,k)/zzm)
             if (entr(ig,k).gt.zzm) then
-               print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
-               abort_message = ''
+               print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
+               abort_message = 'entr dt > m, 2nd'
                CALL abort_gcm (modname,abort_message,1)
             endif
Index: LMDZ5/trunk/libf/phylmd/thermcell_main.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/thermcell_main.F90	(revision 1938)
+++ LMDZ5/trunk/libf/phylmd/thermcell_main.F90	(revision 1943)
@@ -83,4 +83,5 @@
       REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
       real pphi(ngrid,nlay)
+      LOGICAL debut
 
 !   local:
@@ -178,5 +179,4 @@
 !FH/IM     save f0
       real zlevinter(klon)
-      logical debut
        real seuil
       real csc(klon,klev)
@@ -256,42 +256,21 @@
 !
 
-      seuil=0.25
-
-      if (debut)  then
-!        call getin('dvdq',dvdq)
-!        call getin('dqimpl',dqimpl)
-
-         if (iflag_thermals==15.or.iflag_thermals==16) then
-            dvdq=0
-            dqimpl=-1
-         else
-            dvdq=1
-            dqimpl=1
-         endif
-
-         fm0=0.
-         entr0=0.
-         detr0=0.
-
-
-#undef wrgrads_thermcell
-#ifdef wrgrads_thermcell
-! Initialisation des sorties grads pour les thermiques.
-! Pour l'instant en 1D sur le point igout.
-! Utilise par thermcell_out3d.h
-         str10='therm'
-         call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
-     &   rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1.,   &
-     &   ptimestep,str10,'therm ')
-#endif
-
-
-
-      endif
-
-      fm=0. ; entr=0. ; detr=0.
-
-
-      icount=icount+1
+   seuil=0.25
+
+   if (debut) then
+      if (iflag_thermals==15.or.iflag_thermals==16) then
+         dvdq=0
+         dqimpl=-1
+      else
+         dvdq=1
+         dqimpl=1
+      endif
+
+      fm0=0.
+      entr0=0.
+      detr0=0.
+   endif
+   fm=0. ; entr=0. ; detr=0.
+   icount=icount+1
 
 !IM 090508 beg
Index: LMDZ5/trunk/libf/phylmd/thermcell_old.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/thermcell_old.F	(revision 1938)
+++ LMDZ5/trunk/libf/phylmd/thermcell_old.F	(revision 1943)
@@ -1,11 +1,11 @@
-      SUBROUTINE thermcell_2002(ngrid,nlay,ptimestep
+      SUBROUTINE thermcell_2002(ngrid,nlay,ptimestep,iflag_thermals
      s                  ,pplay,pplev,pphi
      s                  ,pu,pv,pt,po
      s                  ,pduadj,pdvadj,pdtadj,pdoadj
-     s                  ,fm0,entr0
-c    s                  ,pu_therm,pv_therm
+     s                  ,fm0,entr0,fraca,wa_moy
      s                  ,r_aspect,l_mix,w2di,tho)
 
       USE dimphy
+      USE write_field_phy
       IMPLICIT NONE
 
@@ -35,5 +35,4 @@
 
 #include "dimensions.h"
-cccc#include "dimphy.h"
 #include "YOMCST.h"
 
@@ -41,5 +40,5 @@
 c   ----------
 
-      INTEGER ngrid,nlay,w2di
+      INTEGER ngrid,nlay,w2di,iflag_thermals
       REAL tho
       real ptimestep,l_mix,r_aspect
@@ -50,15 +49,18 @@
       REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
       real pphi(ngrid,nlay)
-
-      integer idetr
-      save idetr
-      data idetr/3/
-c$OMP THREADPRIVATE(idetr)
+      real fraca(ngrid,nlay+1),zw2(ngrid,nlay+1)
+
+      integer,save :: idetr=3,lev_out=1
+c$OMP THREADPRIVATE(idetr,lev_out)
 
 c   local:
 c   ------
 
+      INTEGER, SAVE :: dvdq=0,flagdq=0,dqimpl=1
+      LOGICAL, SAVE :: debut=.true.
+!$OMP THREADPRIVATE(dvdq,flagdq,debut,dqimpl)
+
       INTEGER ig,k,l,lmax(klon,klev+1),lmaxa(klon),lmix(klon)
-      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
+      real zmax(klon),zw,zz,ztva(klon,klev),zzz
 
       real zlev(klon,klev+1),zlay(klon,klev)
@@ -79,5 +81,4 @@
       real zha(klon,klev)
       real wa_moy(klon,klev+1)
-      real fraca(klon,klev+1)
       real fracc(klon,klev+1)
       real zf,zf2
@@ -86,9 +87,4 @@
 
       real count_time
-      integer isplit,nsplit,ialt
-      parameter (nsplit=10)
-      data isplit/0/
-      save isplit
-c$OMP THREADPRIVATE(isplit)
 
       logical sorties
@@ -142,5 +138,5 @@
 c   ---------------------------------------------------
 
-      print*,'0 OK convect8'
+!     print*,'0 OK convect8'
 
       DO 1010 l=1,nlay
@@ -155,5 +151,5 @@
 1010  CONTINUE
 
-c     print*,'1 OK convect8'
+!     print*,'1 OK convect8'
 c                       --------------------
 c
@@ -177,4 +173,13 @@
 c-----------------------------------------------------------------------
 
+      if (debut)  then
+        flagdq=(iflag_thermals-1000)/100
+        dvdq=(iflag_thermals-(1000+flagdq*100))/10
+        if (flagdq==2) dqimpl=-1
+        if (flagdq==3) dqimpl=1
+        debut=.false.
+      endif
+      print*,'TH flag th ',iflag_thermals,flagdq,dvdq,dqimpl
+
       do l=2,nlay
          do ig=1,ngrid
@@ -192,5 +197,5 @@
       enddo
 
-c     print*,'2 OK convect8'
+!     print*,'2 OK convect8'
 c-----------------------------------------------------------------------
 c   Calcul des densites
@@ -217,5 +222,5 @@
       enddo
 
-c     print*,'3 OK convect8'
+!     print*,'3 OK convect8'
 c------------------------------------------------------------------
 c   Calcul de w2, quarre de w a partir de la cape
@@ -274,5 +279,5 @@
       enddo
 
-c     print*,'4 OK convect8'
+!     print*,'4 OK convect8'
 c Calcul de la couche correspondant a la hauteur du thermique
       do k=1,nlay-1
@@ -287,5 +292,5 @@
       enddo
 
-c     print*,'5 OK convect8'
+!     print*,'5 OK convect8'
 c   Calcule du w max du thermique
       do k=1,nlay
@@ -315,5 +320,5 @@
       enddo
 
-c     print*,'6 OK convect8'
+!     print*,'6 OK convect8'
 c   Longueur caracteristique correspondant a la hauteur des thermiques.
       do  ig=1,ngrid
@@ -330,4 +335,5 @@
 c      call dump2d(iim,jjm-1,zmax(2:ngrid-1),'ZMAX      ')
 
+!     print*,'OKl336'
 c   Calcul de l'entrainement.
 c   Le rapport d'aspect relie la largeur de l'ascendance a l'epaisseur
@@ -348,5 +354,6 @@
       enddo
 
-c     print*,'7 OK convect8'
+
+!     print*,'7 OK convect8'
       do k=1,klev+1
          do ig=1,ngrid
@@ -359,5 +366,5 @@
       enddo
 
-c     print*,'8 OK convect8'
+!     print*,'8 OK convect8'
       do ig=1,ngrid
          lmaxa(ig)=1
@@ -367,4 +374,5 @@
 
 
+!     print*,'OKl372'
       do l=1,nlay-2
          do ig=1,ngrid
@@ -422,5 +430,5 @@
       enddo
 
-c     print*,'9 OK convect8'
+!     print*,'9 OK convect8'
 c     print*,'WA1 ',wa_moy
 
@@ -433,4 +441,5 @@
 c   de la vitesse d'entrainement horizontal dans la couche alimentante.
 
+!     print*,'OKl439'
       do l=2,nlay
          do ig=1,ngrid
@@ -463,5 +472,5 @@
        enddo
 
-c     print*,'10 OK convect8'
+!     print*,'10 OK convect8'
 c     print*,'WA2 ',wa_moy
 c   calcul de la fraction de la maille concernée par l'ascendance en tenant
@@ -500,5 +509,5 @@
       enddo
 
-c     print*,'11 OK convect8'
+!     print*,'11 OK convect8'
 c     print*,'Ea3 ',wa_moy
 c------------------------------------------------------------------
@@ -535,5 +544,5 @@
       enddo
 
-c     print*,'12 OK convect8'
+!     print*,'12 OK convect8'
 c     print*,'WA4 ',wa_moy
 cc------------------------------------------------------------------
@@ -589,4 +598,5 @@
 
 4444   continue
+!     print*,'OK 444 '
 
       if (w2di.eq.1) then
@@ -598,26 +608,47 @@
       endif
 
-      if (1.eq.1) then
+      if (flagdq==0) then
          call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
      .    ,zh,zdhadj,zha)
          call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
      .    ,zo,pdoadj,zoa)
-      else
+         print*,'THERMALS OPT 1'
+      else if (flagdq==1) then
          call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
      .    ,zh,zdhadj,zha)
          call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
      .    ,zo,pdoadj,zoa)
+         print*,'THERMALS OPT 2'
+      else
+         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, 
+     .                    zh,zdhadj,zha,lev_out)
+         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,
+     .                   zo,pdoadj,zoa,lev_out)
+         print*,'THERMALS OPT 3',dqimpl
       endif
 
-      if (1.eq.0) then
-         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
-     .    ,fraca,zmax
-     .    ,zu,zv,pduadj,pdvadj,zua,zva)
-      else
+      print*,'TH VENT ',dvdq
+      if (dvdq==0) then
+!     print*,'TH VENT OK ',dvdq
          call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
      .    ,zu,pduadj,zua)
          call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
      .    ,zv,pdvadj,zva)
+      else if (dvdq==1) then
+         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
+     .    ,fraca,zmax
+     .    ,zu,zv,pduadj,pdvadj,zua,zva)
+      else if (dvdq==2) then
+         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse
+     &    ,fraca,zmax
+     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
+      else if (dvdq==3) then
+         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse 
+     &    ,zu,pduadj,zua,lev_out)
+         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse
+     &    ,zv,pdvadj,zva,lev_out)
       endif
+
+!     CALL writefield_phy('duadj',pduadj,klev)
 
       do l=1,nlay
@@ -632,5 +663,5 @@
 
 
-c     print*,'13 OK convect8'
+!     print*,'13 OK convect8'
 c     print*,'WA5 ',wa_moy
       do l=1,nlay
@@ -654,5 +685,5 @@
 c      enddo
 
-c     print*,'14 OK convect8'
+!     print*,'14 OK convect8'
 c------------------------------------------------------------------
 c   Calculs pour les sorties
@@ -679,123 +710,12 @@
          enddo
       enddo
-
-c     print*,'15 OK convect8'
-
-      isplit=isplit+1
-
-
-c #define und
-	goto 123
-#ifdef und
-      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
-      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
-      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
-      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
-      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
-      CALL writeg1d(1,nlay,zla,'la      ','la      ')
-      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
-      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
-      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
-      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
-      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
-      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
-      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
-      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
-      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
-      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
-      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
-      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
-      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
-      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
-      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
-      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
-      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
-      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
-
-      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
-      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
-      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
-c   recalcul des flux en diagnostique...
-c     print*,'PAS DE TEMPS ',ptimestep
-       call dt2F(pplev,pplay,pt,pdtadj,wh)
-      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
-#endif
-123   continue
-! #define troisD
-#ifdef troisD
-c       if (sorties) then
-      print*,'Debut des wrgradsfi'
-
-c      print*,'16 OK convect8'
-         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
-         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
-         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
-         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
-         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
-         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
-c      print*,'WA6 ',wa_moy
-         call wrgradsfi(1,nlay,zla,'la        ','la        ')
-         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
-         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
-         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
-         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
-         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
-         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
-         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
-         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
-         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
-         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
-         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
-         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
-         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
-         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
-         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
-         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
-         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
-         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
-         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
-         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
-         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
-         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
-         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
-         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
-         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
-
-         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
-         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
-         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
-
-
-c      print*,'17 OK convect8'
-
-         do k=1,klev/10
-            write(str2,'(i2.2)') k
-            str10='wa'//str2
-            do l=1,nlay
-               do ig=1,ngrid
-                  zsortie(ig,l)=wa(ig,k,l)
-               enddo
-            enddo   
-            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
-            do l=1,nlay
-               do ig=1,ngrid
-                  zsortie(ig,l)=larg_part(ig,k,l)
-               enddo
-            enddo
-            str10='la'//str2
-            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
-         enddo
-
-
-c     print*,'18 OK convect8'
-c      endif
-      print*,'Fin des wrgradsfi'
-#endif
-
       endif
 
+!     print*,'15 OK convect8'
+
+
 c     if(wa_moy(1,4).gt.1.e-10) stop
 
-c     print*,'19 OK convect8'
+!     print*,'19 OK convect8'
       return
       end
@@ -925,9 +845,5 @@
       real ratqsdiff(klon,klev)
       real count_time
-      integer isplit,nsplit,ialt
-      parameter (nsplit=10)
-      data isplit/0/
-      save isplit
-c$OMP THREADPRIVATE(isplit)
+      integer ialt
 
       logical sorties
@@ -2197,5 +2113,5 @@
       enddo
 
-      print*,'12 OK convect8'
+!     print*,'12 OK convect8'
 c     print*,'WA4 ',wa_moy
 cc------------------------------------------------------------------
@@ -2347,5 +2263,5 @@
 c      enddo
 
-      print*,'14 OK convect8'
+!     print*,'14 OK convect8'
 c------------------------------------------------------------------
 c   Calculs pour les sorties
@@ -2433,183 +2349,11 @@
       enddo
       enddo     
-cdeja fait
-c      do l=1,nlay
-c         do ig=1,ngrid
-c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
-c            if (detr(ig,l).lt.0.) then
-c                entr(ig,l)=entr(ig,l)-detr(ig,l)
-c                detr(ig,l)=0.
-c     print*,'WARNING !!! detrainement negatif ',ig,l
-c            endif
-c         enddo
-c      enddo
-
-c     print*,'15 OK convect8'
-
-      isplit=isplit+1       
-
-
-c #define und
-	goto 123
-#ifdef und
-      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
-      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
-      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
-      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
-      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
-      CALL writeg1d(1,nlay,zla,'la      ','la      ')
-      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
-      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
-      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
-      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
-      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
-      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
-      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
-      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
-      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
-      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
-      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
-      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
-      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
-      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
-      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
-      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
-      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
-      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
-
-      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
-      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
-      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
-
-c   recalcul des flux en diagnostique...
-c     print*,'PAS DE TEMPS ',ptimestep
-       call dt2F(pplev,pplay,pt,pdtadj,wh)
-      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
-#endif
-123   continue
-! #define troisD
-#ifdef troisD
-c       if (sorties) then
-      print*,'Debut des wrgradsfi'
-
-c      print*,'16 OK convect8'
-c         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
-c         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
-         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
-         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
-c         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
-c         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
-c      print*,'WA6 ',wa_moy
-c         call wrgradsfi(1,nlay,zla,'la        ','la        ')
-c         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
-         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
-         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
-         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
-         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
-         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
-         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
-         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
-         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
-         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
-         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
-         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
-         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
-         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
-         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
-         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
-         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
-         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
-         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
-         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
-         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
-         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
-         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
-         call wrgradsfi(1,nlay,w_est,'w_est      ','w_est      ')
-con sort les moments
-         call wrgradsfi(1,nlay,thetath2,'zh2       ','zh2       ')
-         call wrgradsfi(1,nlay,wth2,'w2       ','w2       ')
-         call wrgradsfi(1,nlay,wth3,'w3       ','w3       ')
-         call wrgradsfi(1,nlay,q2,'q2       ','q2       ')
-         call wrgradsfi(1,nlay,dtheta,'dT       ','dT       ')
-c
-         call wrgradsfi(1,nlay,zw_sec,'zw_sec       ','zw_sec       ')
-         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
-         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
-         call wrgradsfi(1,nlay,nu,'nu       ','nu       ')
-
-         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
-         call wrgradsfi(1,nlay,zoa,'zoa        ','zoa        ')
-c         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
-c         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
-
-cAM:nouveaux diagnostiques
-         call wrgradsfi(1,nlay,zthl,'zthl        ','zthl        ')
-         call wrgradsfi(1,nlay,zta,'zta        ','zta        ')
-         call wrgradsfi(1,nlay,zl,'zl        ','zl        ')
-         call wrgradsfi(1,nlay,zdthladj,'zdthladj    ',
-     s        'zdthladj    ')
-         call wrgradsfi(1,nlay,ztla,'ztla      ','ztla      ')
-         call wrgradsfi(1,nlay,zqta,'zqta      ','zqta      ')
-         call wrgradsfi(1,nlay,zqla,'zqla      ','zqla      ')
-         call wrgradsfi(1,nlay,deltaz,'deltaz      ','deltaz      ')
-cCR:nouveaux diagnostiques
-      call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')
-      call wrgradsfi(1,nlay,detr_star  ,'detr_star   ','detr_star   ')     
-      call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
-      call wrgradsfi(1,nlay,zqsat    ,'zqsat   ','zqsat   ')
-      call wrgradsfi(1,nlay,zqsatth    ,'qsath   ','qsath   ')
-      call wrgradsfi(1,nlay,alim_star    ,'alim_star   ','alim_star   ')
-      call wrgradsfi(1,nlay,alim    ,'alim   ','alim   ')
-      call wrgradsfi(1,1,f,'f      ','f      ')
-      call wrgradsfi(1,1,alim_star_tot,'a_s_t      ','a_s_t      ')
-      call wrgradsfi(1,1,alim_star2,'a_2      ','a_2      ')
-      call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
-      call wrgradsfi(1,1,zmax_sec,'z_sec      ','z_sec      ')
-c      call wrgradsfi(1,1,zmax_sec2,'zz_se      ','zz_se      ')
-      call wrgradsfi(1,1,zmix,'zmix      ','zmix      ') 
-      call wrgradsfi(1,1,nivcon,'nivcon      ','nivcon      ')
-      call wrgradsfi(1,1,zcon,'zcon      ','zcon      ')
-      zsortie1d(:)=lmax(:)
-      call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
-      call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
-      call wrgradsfi(1,1,wmax_sec,'w_sec      ','w_sec      ')
-      zsortie1d(:)=lmix(:)
-      call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
-      zsortie1d(:)=lentr(:)
-      call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
-
-c      print*,'17 OK convect8'
-
-         do k=1,klev/10
-            write(str2,'(i2.2)') k
-            str10='wa'//str2
-            do l=1,nlay
-               do ig=1,ngrid
-                  zsortie(ig,l)=wa(ig,k,l)
-               enddo
-            enddo   
-            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
-            do l=1,nlay
-               do ig=1,ngrid
-                  zsortie(ig,l)=larg_part(ig,k,l)
-               enddo
-            enddo
-           str10='la'//str2
-            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
-         enddo
-
-
-c     print*,'18 OK convect8'
-c      endif
-      print*,'Fin des wrgradsfi'
-#endif
 
       endif
 
-c     if(wa_moy(1,4).gt.1.e-10) stop
-
-      print*,'19 OK convect8'
+!     print*,'19 OK convect8'
       return
       end
+
       SUBROUTINE thermcell_eau(ngrid,nlay,ptimestep
      s                  ,pplay,pplev,pphi
@@ -2710,9 +2454,5 @@
 
       real count_time
-      integer isplit,nsplit,ialt
-      parameter (nsplit=10)
-      data isplit/0/
-      save isplit
-c$OMP THREADPRIVATE(isplit)
+      integer ialt
 
       logical sorties
@@ -2861,5 +2601,5 @@
 c   ---------------------------------------------------
 
-      print*,'0 OK convect8'
+!     print*,'0 OK convect8'
 
       DO 1010 l=1,nlay
@@ -3568,15 +3308,756 @@
 c------------------------------------------------------------------
 
-      if(sorties) then
+      return
+      end
+
+      SUBROUTINE thermcell(ngrid,nlay,ptimestep
+     s                  ,pplay,pplev,pphi
+     s                  ,pu,pv,pt,po
+     s                  ,pduadj,pdvadj,pdtadj,pdoadj
+     s                  ,fm0,entr0
+c    s                  ,pu_therm,pv_therm
+     s                  ,r_aspect,l_mix,w2di,tho)
+
+      USE dimphy
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Calcul du transport verticale dans la couche limite en presence
+c   de "thermiques" explicitement representes
+c
+c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
+c
+c   le thermique est supposé homogène et dissipé par mélange avec
+c   son environnement. la longueur l_mix contrôle l'efficacité du
+c   mélange
+c
+c   Le calcul du transport des différentes espèces se fait en prenant
+c   en compte:
+c     1. un flux de masse montant
+c     2. un flux de masse descendant
+c     3. un entrainement
+c     4. un detrainement
+c
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+
+#include "dimensions.h"
+cccc#include "dimphy.h"
+#include "YOMCST.h"
+
+c   arguments:
+c   ----------
+
+      INTEGER ngrid,nlay,w2di
+      REAL tho
+      real ptimestep,l_mix,r_aspect
+      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
+      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
+      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
+      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
+      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
+      real pphi(ngrid,nlay)
+
+      integer idetr
+      save idetr
+      data idetr/3/
+c$OMP THREADPRIVATE(idetr)
+
+c   local:
+c   ------
+
+      INTEGER ig,k,l,lmaxa(klon),lmix(klon)
+      real zsortie1d(klon)
+c CR: on remplace lmax(klon,klev+1)
+      INTEGER lmax(klon),lmin(klon),lentr(klon)
+      real linter(klon)
+      real zmix(klon), fracazmix(klon) 
+c RC 
+      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
+
+      real zlev(klon,klev+1),zlay(klon,klev)
+      REAL zh(klon,klev),zdhadj(klon,klev)
+      REAL ztv(klon,klev)
+      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
+      REAL wh(klon,klev+1)
+      real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
+      real zla(klon,klev+1)
+      real zwa(klon,klev+1)
+      real zld(klon,klev+1)
+      real zwd(klon,klev+1)
+      real zsortie(klon,klev)
+      real zva(klon,klev)
+      real zua(klon,klev)
+      real zoa(klon,klev)
+
+      real zha(klon,klev)
+      real wa_moy(klon,klev+1)
+      real fraca(klon,klev+1)
+      real fracc(klon,klev+1)
+      real zf,zf2
+      real thetath2(klon,klev),wth2(klon,klev)
+!      common/comtherm/thetath2,wth2
+
+      real count_time
+      integer ialt
+
+      logical sorties
+      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
+      real zpspsk(klon,klev)
+
+c     real wmax(klon,klev),wmaxa(klon)
+      real wmax(klon),wmaxa(klon)
+      real wa(klon,klev,klev+1)
+      real wd(klon,klev+1)
+      real larg_part(klon,klev,klev+1)
+      real fracd(klon,klev+1)
+      real xxx(klon,klev+1)
+      real larg_cons(klon,klev+1)
+      real larg_detr(klon,klev+1)
+      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
+      real pu_therm(klon,klev),pv_therm(klon,klev)
+      real fm(klon,klev+1),entr(klon,klev)
+      real fmc(klon,klev+1)
+
+cCR:nouvelles variables
+      real f_star(klon,klev+1),entr_star(klon,klev)
+      real entr_star_tot(klon),entr_star2(klon)
+      real f(klon), f0(klon)
+      real zlevinter(klon)
+      logical first
+      data first /.false./
+      save first
+c$OMP THREADPRIVATE(first)
+cRC
+
+      character*2 str2
+      character*10 str10
+
+      character (len=20) :: modname='thermcell'
+      character (len=80) :: abort_message
+
+      LOGICAL vtest(klon),down
+
+      EXTERNAL SCOPY
+
+      integer ncorrec,ll
+      save ncorrec
+      data ncorrec/0/
+c$OMP THREADPRIVATE(ncorrec)
+      
+c
+c-----------------------------------------------------------------------
+c   initialisation:
+c   ---------------
+c
+       sorties=.true.
+      IF(ngrid.NE.klon) THEN
+         PRINT*
+         PRINT*,'STOP dans convadj'
+         PRINT*,'ngrid    =',ngrid
+         PRINT*,'klon  =',klon
+      ENDIF
+c
+c-----------------------------------------------------------------------
+c   incrementation eventuelle de tendances precedentes:
+c   ---------------------------------------------------
+
+!      print*,'0 OK convect8'
+
+      DO 1010 l=1,nlay
+         DO 1015 ig=1,ngrid
+            zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
+            zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
+            zu(ig,l)=pu(ig,l)
+            zv(ig,l)=pv(ig,l)
+            zo(ig,l)=po(ig,l)
+            ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
+1015     CONTINUE
+1010  CONTINUE
+
+!      print*,'1 OK convect8'
+c                       --------------------
+c
+c
+c                       + + + + + + + + + + +
+c
+c
+c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
+c  wh,wt,wo ...
+c
+c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
+c
+c
+c                       --------------------   zlev(1)
+c                       \\\\\\\\\\\\\\\\\\\\
+c
+c
+
+c-----------------------------------------------------------------------
+c   Calcul des altitudes des couches
+c-----------------------------------------------------------------------
+
+      do l=2,nlay
+         do ig=1,ngrid
+            zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
+         enddo
+      enddo
+      do ig=1,ngrid
+         zlev(ig,1)=0.
+         zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
+      enddo
       do l=1,nlay
          do ig=1,ngrid
-            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
-            zld(ig,l)=fracd(ig,l)*zmax(ig)
-            if(1.-fracd(ig,l).gt.1.e-10)
-     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
-         enddo
-      enddo
+            zlay(ig,l)=pphi(ig,l)/RG
+         enddo
+      enddo
+
+c      print*,'2 OK convect8'
+c-----------------------------------------------------------------------
+c   Calcul des densites
+c-----------------------------------------------------------------------
 
       do l=1,nlay
+         do ig=1,ngrid
+            rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
+         enddo
+      enddo
+
+      do l=2,nlay
+         do ig=1,ngrid
+            rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
+         enddo
+      enddo
+
+      do k=1,nlay
+         do l=1,nlay+1
+            do ig=1,ngrid
+               wa(ig,k,l)=0.
+            enddo
+         enddo
+      enddo
+
+c      print*,'3 OK convect8'
+c------------------------------------------------------------------
+c   Calcul de w2, quarre de w a partir de la cape
+c   a partir de w2, on calcule wa, vitesse de l'ascendance
+c
+c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
+c   w2 est stoke dans wa
+c
+c   ATTENTION: dans convect8, on n'utilise le calcule des wa
+c   independants par couches que pour calculer l'entrainement
+c   a la base et la hauteur max de l'ascendance.
+c
+c   Indicages:
+c   l'ascendance provenant du niveau k traverse l'interface l avec
+c   une vitesse wa(k,l).
+c
+c                       --------------------
+c
+c                       + + + + + + + + + + 
+c
+c  wa(k,l)   ----       --------------------    l
+c             /\
+c            /||\       + + + + + + + + + + 
+c             ||
+c             ||        --------------------
+c             ||
+c             ||        + + + + + + + + + + 
+c             ||
+c             ||        --------------------
+c             ||__
+c             |___      + + + + + + + + + +     k
+c
+c                       --------------------
+c
+c
+c
+c------------------------------------------------------------------
+
+cCR: ponderation entrainement des couches instables
+cdef des entr_star tels que entr=f*entr_star      
+      do l=1,klev
+         do ig=1,ngrid 
+            entr_star(ig,l)=0.
+         enddo
+      enddo
+c determination de la longueur de la couche d entrainement
+      do ig=1,ngrid
+         lentr(ig)=1
+      enddo
+
+con ne considere que les premieres couches instables
+      do k=nlay-2,1,-1
+         do ig=1,ngrid
+            if (ztv(ig,k).gt.ztv(ig,k+1).and.
+     s          ztv(ig,k+1).le.ztv(ig,k+2)) then
+               lentr(ig)=k
+            endif
+          enddo
+      enddo
+    
+c determination du lmin: couche d ou provient le thermique
+      do ig=1,ngrid
+         lmin(ig)=1
+      enddo
+      do ig=1,ngrid
+         do l=nlay,2,-1
+            if (ztv(ig,l-1).gt.ztv(ig,l)) then
+               lmin(ig)=l-1
+            endif
+         enddo
+      enddo
+c
+c definition de l'entrainement des couches
+      do l=1,klev-1
+         do ig=1,ngrid 
+            if (ztv(ig,l).gt.ztv(ig,l+1).and.
+     s          l.ge.lmin(ig).and.l.le.lentr(ig)) then 
+                 entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
+     s                           (zlev(ig,l+1)-zlev(ig,l))
+            endif
+         enddo
+      enddo
+c pas de thermique si couches 1->5 stables
+      do ig=1,ngrid
+         if (lmin(ig).gt.5) then
+            do l=1,klev
+               entr_star(ig,l)=0.
+            enddo
+         endif
+      enddo 
+c calcul de l entrainement total
+      do ig=1,ngrid
+         entr_star_tot(ig)=0.
+      enddo
+      do ig=1,ngrid
+         do k=1,klev
+            entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
+         enddo
+      enddo
+c
+      print*,'fin calcul entr_star'
+      do k=1,klev
+         do ig=1,ngrid 
+            ztva(ig,k)=ztv(ig,k)
+         enddo
+      enddo
+cRC
+c      print*,'7 OK convect8'
+      do k=1,klev+1
+         do ig=1,ngrid
+            zw2(ig,k)=0.
+            fmc(ig,k)=0.
+cCR
+            f_star(ig,k)=0.
+cRC
+            larg_cons(ig,k)=0.
+            larg_detr(ig,k)=0.
+            wa_moy(ig,k)=0.
+         enddo
+      enddo
+
+c      print*,'8 OK convect8'
+      do ig=1,ngrid
+         linter(ig)=1.
+         lmaxa(ig)=1
+         lmix(ig)=1
+         wmaxa(ig)=0.
+      enddo
+
+cCR: 
+      do l=1,nlay-2
+         do ig=1,ngrid
+            if (ztv(ig,l).gt.ztv(ig,l+1)
+     s         .and.entr_star(ig,l).gt.1.e-10
+     s         .and.zw2(ig,l).lt.1e-10) then
+               f_star(ig,l+1)=entr_star(ig,l)
+ctest:calcul de dteta
+               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
+     s                     *(zlev(ig,l+1)-zlev(ig,l))
+     s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
+               larg_detr(ig,l)=0.
+            else if ((zw2(ig,l).ge.1e-10).and.
+     s               (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)
+     s                    *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+
+     s                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
+     s                     *(zlev(ig,l+1)-zlev(ig,l))
+            endif
+c determination de zmax continu par interpolation lineaire
+            if (zw2(ig,l+1).lt.0.) then
+ctest
+               if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
+                  print*,'pb linter'
+               endif
+               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
+     s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
+               zw2(ig,l+1)=0.
+               lmaxa(ig)=l
+            else
+               if (zw2(ig,l+1).lt.0.) then
+                  print*,'pb1 zw2<0'
+               endif
+               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
+            endif
+            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
+c   lmix est le niveau de la couche ou w (wa_moy) est maximum
+               lmix(ig)=l+1
+               wmaxa(ig)=wa_moy(ig,l+1)
+            endif
+         enddo
+      enddo
+      print*,'fin calcul zw2'
+c
+c 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
+c pas de thermique si couches 1->5 stables
+      do ig=1,ngrid
+         if (lmin(ig).gt.5) then
+            lmax(ig)=1
+            lmin(ig)=1
+         endif
+      enddo 
+c    
+c Determination de zw2 max
+      do ig=1,ngrid
+         wmax(ig)=0.
+      enddo
+
+      do l=1,nlay
+         do ig=1,ngrid
+            if (l.le.lmax(ig)) then
+                if (zw2(ig,l).lt.0.)then
+                  print*,'pb2 zw2<0'
+                endif
+                zw2(ig,l)=sqrt(zw2(ig,l))
+                wmax(ig)=max(wmax(ig),zw2(ig,l))
+            else
+                 zw2(ig,l)=0.
+            endif
+          enddo
+      enddo
+
+c   Longueur caracteristique correspondant a la hauteur des thermiques.
+      do  ig=1,ngrid
+         zmax(ig)=0.
+         zlevinter(ig)=zlev(ig,1)
+      enddo
+      do  ig=1,ngrid
+c calcul de zlevinter
+          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
+     s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
+     s    -zlev(ig,lmax(ig)))
+       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
+      enddo
+
+      print*,'avant fermeture'
+c Fermeture,determination de f
+      do ig=1,ngrid
+         entr_star2(ig)=0.
+      enddo
+      do ig=1,ngrid
+         if (entr_star_tot(ig).LT.1.e-10) then
+            f(ig)=0.
+         else
+             do k=lmin(ig),lentr(ig)
+                entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
+     s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
+             enddo
+c Nouvelle fermeture
+             f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
+     s             *entr_star2(ig))*entr_star_tot(ig)
+ctest
+c             if (first) then
+c             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
+c     s             *wmax(ig))
+c             endif
+         endif
+c         f0(ig)=f(ig)
+c         first=.true.
+      enddo
+      print*,'apres fermeture'
+
+c Calcul de l'entrainement
+       do k=1,klev
+         do ig=1,ngrid 
+            entr(ig,k)=f(ig)*entr_star(ig,k)
+         enddo
+      enddo
+c Calcul des flux
+      do ig=1,ngrid
+         do l=1,lmax(ig)-1
+            fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
+         enddo
+      enddo
+
+cRC
+
+
+c      print*,'9 OK convect8'
+c     print*,'WA1 ',wa_moy
+
+c   determination de l'indice du debut de la mixed layer ou w decroit
+
+c   calcul de la largeur de chaque ascendance dans le cas conservatif.
+c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
+c   d'une couche est égale à la hauteur de la couche alimentante.
+c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
+c   de la vitesse d'entrainement horizontal dans la couche alimentante.
+
+      do l=2,nlay
+         do ig=1,ngrid
+            if (l.le.lmaxa(ig)) then
+               zw=max(wa_moy(ig,l),1.e-10)
+               larg_cons(ig,l)=zmax(ig)*r_aspect
+     s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
+            endif
+         enddo
+      enddo
+
+      do l=2,nlay
+         do ig=1,ngrid
+            if (l.le.lmaxa(ig)) then
+c              if (idetr.eq.0) then
+c  cette option est finalement en dur.
+                  if ((l_mix*zlev(ig,l)).lt.0.)then
+                   print*,'pb l_mix*zlev<0'
+                  endif
+                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
+c              else if (idetr.eq.1) then
+c                 larg_detr(ig,l)=larg_cons(ig,l)
+c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
+c              else if (idetr.eq.2) then
+c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
+c    s            *sqrt(wa_moy(ig,l))
+c              else if (idetr.eq.4) then
+c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
+c    s            *wa_moy(ig,l)
+c              endif
+            endif
+         enddo
+       enddo
+
+c      print*,'10 OK convect8'
+c     print*,'WA2 ',wa_moy
+c   calcul de la fraction de la maille concernée par l'ascendance en tenant
+c   compte de l'epluchage du thermique.
+c
+cCR def de  zmix continu (profil parabolique des vitesses)
+      do ig=1,ngrid
+           if (lmix(ig).gt.1.) then
+c test 
+              if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
+     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
+     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
+     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
+     s        then
+c             
+            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
+     s        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
+     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
+     s        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
+     s        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
+     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
+     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
+     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
+            else
+            zmix(ig)=zlev(ig,lmix(ig))
+            print*,'pb zmix'
+            endif
+         else 
+         zmix(ig)=0.
+         endif
+ctest
+         if ((zmax(ig)-zmix(ig)).lt.0.) then
+            zmix(ig)=0.99*zmax(ig)
+c            print*,'pb zmix>zmax'
+         endif
+      enddo
+c
+c calcul du nouveau lmix correspondant
+      do ig=1,ngrid
+         do l=1,klev
+            if (zmix(ig).ge.zlev(ig,l).and.
+     s          zmix(ig).lt.zlev(ig,l+1)) then
+              lmix(ig)=l
+             endif
+          enddo
+      enddo
+c
+      do l=2,nlay
+         do ig=1,ngrid
+            if(larg_cons(ig,l).gt.1.) then
+c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
+               fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
+     s            /(r_aspect*zmax(ig))
+c test
+               fraca(ig,l)=max(fraca(ig,l),0.)
+               fraca(ig,l)=min(fraca(ig,l),0.5)
+               fracd(ig,l)=1.-fraca(ig,l)
+               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
+            else
+c              wa_moy(ig,l)=0.
+               fraca(ig,l)=0.
+               fracc(ig,l)=0.
+               fracd(ig,l)=1.
+            endif
+         enddo
+      enddo                  
+cCR: calcul de fracazmix
+       do ig=1,ngrid
+          fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
+     s     (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
+     s    +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
+     s    -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
+       enddo
+c
+       do l=2,nlay
+          do ig=1,ngrid
+             if(larg_cons(ig,l).gt.1.) then
+               if (l.gt.lmix(ig)) then
+ctest
+                 if (zmax(ig)-zmix(ig).lt.1.e-10) then
+c                   print*,'pb xxx'
+                   xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
+                 else
+                 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
+                 endif
+           if (idetr.eq.0) then
+               fraca(ig,l)=fracazmix(ig)
+           else if (idetr.eq.1) then
+               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
+           else if (idetr.eq.2) then
+               fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
+           else
+               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
+           endif
+c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
+               fraca(ig,l)=max(fraca(ig,l),0.)
+               fraca(ig,l)=min(fraca(ig,l),0.5)
+               fracd(ig,l)=1.-fraca(ig,l)
+               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
+             endif
+            endif
+         enddo
+      enddo
+      
+      print*,'fin calcul fraca'
+c      print*,'11 OK convect8'
+c     print*,'Ea3 ',wa_moy
+c------------------------------------------------------------------
+c   Calcul de fracd, wd
+c   somme wa - wd = 0
+c------------------------------------------------------------------
+
+
+      do ig=1,ngrid
+         fm(ig,1)=0.
+         fm(ig,nlay+1)=0.
+      enddo
+
+      do l=2,nlay
+           do ig=1,ngrid
+              fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
+cCR:test
+              if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
+     s            .and.l.gt.lmix(ig)) then
+                 fm(ig,l)=fm(ig,l-1)
+c                 write(1,*)'ajustement fm, l',l
+              endif
+c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
+cRC
+           enddo
+         do ig=1,ngrid
+            if(fracd(ig,l).lt.0.1) then
+              abort_message = 'fracd trop petit'
+              CALL abort_gcm (modname,abort_message,1)
+            else
+c    vitesse descendante "diagnostique"
+               wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
+            endif
+         enddo
+      enddo
+
+      do l=1,nlay
+         do ig=1,ngrid
+c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
+            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
+         enddo
+      enddo
+
+!     print*,'12 OK convect8'
+c     print*,'WA4 ',wa_moy
+cc------------------------------------------------------------------
+c   calcul du transport vertical
+c------------------------------------------------------------------
+
+      go to 4444
+c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
+      do l=2,nlay-1
+         do ig=1,ngrid
+            if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
+     s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
+c     print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
+c    s         ,fm(ig,l+1)*ptimestep
+c    s         ,'   M=',masse(ig,l),masse(ig,l+1)
+             endif
+         enddo
+      enddo
+
+      do l=1,nlay
+         do ig=1,ngrid
+            if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
+c     print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
+c    s         ,entr(ig,l)*ptimestep
+c    s         ,'   M=',masse(ig,l)
+             endif
+         enddo
+      enddo
+
+      do l=1,nlay
+         do ig=1,ngrid
+            if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
+c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
+c    s         ,'   FM=',fm(ig,l)
+            endif
+            if(.not.masse(ig,l).ge.1.e-10
+     s         .or..not.masse(ig,l).le.1.e4) then
+c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
+c    s         ,'   M=',masse(ig,l)
+c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
+c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
+c     print*,'zlev(ig,l+1),zlev(ig,l)'
+c    s                ,zlev(ig,l+1),zlev(ig,l)
+c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
+c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
+            endif
+            if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
+c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
+c    s         ,'   E=',entr(ig,l)
+            endif
+         enddo
+      enddo
+
+4444   continue
+
+cCR:redefinition du entr
+       do l=1,nlay
          do ig=1,ngrid
             detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
@@ -3588,12 +4069,101 @@
          enddo
       enddo
+cRC
+      if (w2di.eq.1) then
+         fm0=fm0+ptimestep*(fm-fm0)/tho
+         entr0=entr0+ptimestep*(entr-entr0)/tho
+      else
+         fm0=fm
+         entr0=entr
+      endif
+
+      if (1.eq.1) then
+         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
+     .    ,zh,zdhadj,zha)
+         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
+     .    ,zo,pdoadj,zoa)
+      else
+         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
+     .    ,zh,zdhadj,zha)
+         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
+     .    ,zo,pdoadj,zoa)
+      endif
+
+      if (1.eq.0) then
+         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
+     .    ,fraca,zmax
+     .    ,zu,zv,pduadj,pdvadj,zua,zva)
+      else
+         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
+     .    ,zu,pduadj,zua)
+         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
+     .    ,zv,pdvadj,zva)
+      endif
+
+      do l=1,nlay
+         do ig=1,ngrid
+            zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
+            zf2=zf/(1.-zf)
+            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
+            wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
+         enddo
+      enddo
+
+
+
+c     print*,'13 OK convect8'
+c     print*,'WA5 ',wa_moy
+      do l=1,nlay
+         do ig=1,ngrid
+            pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
+         enddo
+      enddo
+
+
+c     do l=1,nlay
+c        do ig=1,ngrid
+c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
+c     print*,'WARN!!! ig=',ig,'  l=',l
+c    s         ,'   pdtadj=',pdtadj(ig,l)
+c           endif
+c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
+c     print*,'WARN!!! ig=',ig,'  l=',l
+c    s         ,'   pdoadj=',pdoadj(ig,l)
+c           endif
+c        enddo
+c      enddo
+
+!     print*,'14 OK convect8'
+c------------------------------------------------------------------
+c   Calculs pour les sorties
+c------------------------------------------------------------------
+
+      if(sorties) then
+      do l=1,nlay
+         do ig=1,ngrid
+            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
+            zld(ig,l)=fracd(ig,l)*zmax(ig)
+            if(1.-fracd(ig,l).gt.1.e-10)
+     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
+         enddo
+      enddo
+
+cdeja fait
+c      do l=1,nlay
+c         do ig=1,ngrid
+c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
+c            if (detr(ig,l).lt.0.) then
+c                entr(ig,l)=entr(ig,l)-detr(ig,l)
+c                detr(ig,l)=0.
+c     print*,'WARNING !!! detrainement negatif ',ig,l
+c            endif
+c         enddo
+c      enddo
 
 c     print*,'15 OK convect8'
 
-      isplit=isplit+1
-
 
 c #define und
-	goto 123
+      goto 123
 #ifdef und
       CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
@@ -3632,5 +4202,5 @@
 #endif
 123   continue
-! #define troisD
+#define troisD
 #ifdef troisD
 c       if (sorties) then
@@ -3676,13 +4246,4 @@
          call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
 
-cAM:nouveaux diagnostiques
-         call wrgradsfi(1,nlay,zthl,'zthl        ','zthl        ')
-         call wrgradsfi(1,nlay,zta,'zta        ','zta        ')
-         call wrgradsfi(1,nlay,zl,'zl        ','zl        ')
-         call wrgradsfi(1,nlay,zdthladj,'zdthladj    ',
-     s        'zdthladj    ')
-         call wrgradsfi(1,nlay,ztla,'ztla      ','ztla      ')
-         call wrgradsfi(1,nlay,zqta,'zqta      ','zqta      ')
-         call wrgradsfi(1,nlay,zqla,'zqla      ','zqla      ')
 cCR:nouveaux diagnostiques
       call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
@@ -3728,993 +4289,5 @@
 c     if(wa_moy(1,4).gt.1.e-10) stop
 
-c     print*,'19 OK convect8'
-      return
-      end
-
-      SUBROUTINE thermcell(ngrid,nlay,ptimestep
-     s                  ,pplay,pplev,pphi
-     s                  ,pu,pv,pt,po
-     s                  ,pduadj,pdvadj,pdtadj,pdoadj
-     s                  ,fm0,entr0
-c    s                  ,pu_therm,pv_therm
-     s                  ,r_aspect,l_mix,w2di,tho)
-
-      USE dimphy
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Calcul du transport verticale dans la couche limite en presence
-c   de "thermiques" explicitement representes
-c
-c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
-c
-c   le thermique est supposé homogène et dissipé par mélange avec
-c   son environnement. la longueur l_mix contrôle l'efficacité du
-c   mélange
-c
-c   Le calcul du transport des différentes espèces se fait en prenant
-c   en compte:
-c     1. un flux de masse montant
-c     2. un flux de masse descendant
-c     3. un entrainement
-c     4. un detrainement
-c
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-#include "dimensions.h"
-cccc#include "dimphy.h"
-#include "YOMCST.h"
-
-c   arguments:
-c   ----------
-
-      INTEGER ngrid,nlay,w2di
-      REAL tho
-      real ptimestep,l_mix,r_aspect
-      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
-      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
-      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
-      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
-      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
-      real pphi(ngrid,nlay)
-
-      integer idetr
-      save idetr
-      data idetr/3/
-c$OMP THREADPRIVATE(idetr)
-
-c   local:
-c   ------
-
-      INTEGER ig,k,l,lmaxa(klon),lmix(klon)
-      real zsortie1d(klon)
-c CR: on remplace lmax(klon,klev+1)
-      INTEGER lmax(klon),lmin(klon),lentr(klon)
-      real linter(klon)
-      real zmix(klon), fracazmix(klon) 
-c RC 
-      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
-
-      real zlev(klon,klev+1),zlay(klon,klev)
-      REAL zh(klon,klev),zdhadj(klon,klev)
-      REAL ztv(klon,klev)
-      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
-      REAL wh(klon,klev+1)
-      real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
-      real zla(klon,klev+1)
-      real zwa(klon,klev+1)
-      real zld(klon,klev+1)
-      real zwd(klon,klev+1)
-      real zsortie(klon,klev)
-      real zva(klon,klev)
-      real zua(klon,klev)
-      real zoa(klon,klev)
-
-      real zha(klon,klev)
-      real wa_moy(klon,klev+1)
-      real fraca(klon,klev+1)
-      real fracc(klon,klev+1)
-      real zf,zf2
-      real thetath2(klon,klev),wth2(klon,klev)
-!      common/comtherm/thetath2,wth2
-
-      real count_time
-      integer isplit,nsplit,ialt
-      parameter (nsplit=10)
-      data isplit/0/
-      save isplit
-c$OMP THREADPRIVATE(isplit)
-
-      logical sorties
-      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
-      real zpspsk(klon,klev)
-
-c     real wmax(klon,klev),wmaxa(klon)
-      real wmax(klon),wmaxa(klon)
-      real wa(klon,klev,klev+1)
-      real wd(klon,klev+1)
-      real larg_part(klon,klev,klev+1)
-      real fracd(klon,klev+1)
-      real xxx(klon,klev+1)
-      real larg_cons(klon,klev+1)
-      real larg_detr(klon,klev+1)
-      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
-      real pu_therm(klon,klev),pv_therm(klon,klev)
-      real fm(klon,klev+1),entr(klon,klev)
-      real fmc(klon,klev+1)
-
-cCR:nouvelles variables
-      real f_star(klon,klev+1),entr_star(klon,klev)
-      real entr_star_tot(klon),entr_star2(klon)
-      real f(klon), f0(klon)
-      real zlevinter(klon)
-      logical first
-      data first /.false./
-      save first
-c$OMP THREADPRIVATE(first)
-cRC
-
-      character*2 str2
-      character*10 str10
-
-      character (len=20) :: modname='thermcell'
-      character (len=80) :: abort_message
-
-      LOGICAL vtest(klon),down
-
-      EXTERNAL SCOPY
-
-      integer ncorrec,ll
-      save ncorrec
-      data ncorrec/0/
-c$OMP THREADPRIVATE(ncorrec)
-      
-c
-c-----------------------------------------------------------------------
-c   initialisation:
-c   ---------------
-c
-       sorties=.true.
-      IF(ngrid.NE.klon) THEN
-         PRINT*
-         PRINT*,'STOP dans convadj'
-         PRINT*,'ngrid    =',ngrid
-         PRINT*,'klon  =',klon
-      ENDIF
-c
-c-----------------------------------------------------------------------
-c   incrementation eventuelle de tendances precedentes:
-c   ---------------------------------------------------
-
-       print*,'0 OK convect8'
-
-      DO 1010 l=1,nlay
-         DO 1015 ig=1,ngrid
-            zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
-            zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
-            zu(ig,l)=pu(ig,l)
-            zv(ig,l)=pv(ig,l)
-            zo(ig,l)=po(ig,l)
-            ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
-1015     CONTINUE
-1010  CONTINUE
-
-       print*,'1 OK convect8'
-c                       --------------------
-c
-c
-c                       + + + + + + + + + + +
-c
-c
-c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
-c  wh,wt,wo ...
-c
-c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
-c
-c
-c                       --------------------   zlev(1)
-c                       \\\\\\\\\\\\\\\\\\\\
-c
-c
-
-c-----------------------------------------------------------------------
-c   Calcul des altitudes des couches
-c-----------------------------------------------------------------------
-
-      do l=2,nlay
-         do ig=1,ngrid
-            zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
-         enddo
-      enddo
-      do ig=1,ngrid
-         zlev(ig,1)=0.
-         zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
-      enddo
-      do l=1,nlay
-         do ig=1,ngrid
-            zlay(ig,l)=pphi(ig,l)/RG
-         enddo
-      enddo
-
-c      print*,'2 OK convect8'
-c-----------------------------------------------------------------------
-c   Calcul des densites
-c-----------------------------------------------------------------------
-
-      do l=1,nlay
-         do ig=1,ngrid
-            rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
-         enddo
-      enddo
-
-      do l=2,nlay
-         do ig=1,ngrid
-            rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
-         enddo
-      enddo
-
-      do k=1,nlay
-         do l=1,nlay+1
-            do ig=1,ngrid
-               wa(ig,k,l)=0.
-            enddo
-         enddo
-      enddo
-
-c      print*,'3 OK convect8'
-c------------------------------------------------------------------
-c   Calcul de w2, quarre de w a partir de la cape
-c   a partir de w2, on calcule wa, vitesse de l'ascendance
-c
-c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
-c   w2 est stoke dans wa
-c
-c   ATTENTION: dans convect8, on n'utilise le calcule des wa
-c   independants par couches que pour calculer l'entrainement
-c   a la base et la hauteur max de l'ascendance.
-c
-c   Indicages:
-c   l'ascendance provenant du niveau k traverse l'interface l avec
-c   une vitesse wa(k,l).
-c
-c                       --------------------
-c
-c                       + + + + + + + + + + 
-c
-c  wa(k,l)   ----       --------------------    l
-c             /\
-c            /||\       + + + + + + + + + + 
-c             ||
-c             ||        --------------------
-c             ||
-c             ||        + + + + + + + + + + 
-c             ||
-c             ||        --------------------
-c             ||__
-c             |___      + + + + + + + + + +     k
-c
-c                       --------------------
-c
-c
-c
-c------------------------------------------------------------------
-
-cCR: ponderation entrainement des couches instables
-cdef des entr_star tels que entr=f*entr_star      
-      do l=1,klev
-         do ig=1,ngrid 
-            entr_star(ig,l)=0.
-         enddo
-      enddo
-c determination de la longueur de la couche d entrainement
-      do ig=1,ngrid
-         lentr(ig)=1
-      enddo
-
-con ne considere que les premieres couches instables
-      do k=nlay-2,1,-1
-         do ig=1,ngrid
-            if (ztv(ig,k).gt.ztv(ig,k+1).and.
-     s          ztv(ig,k+1).le.ztv(ig,k+2)) then
-               lentr(ig)=k
-            endif
-          enddo
-      enddo
-    
-c determination du lmin: couche d ou provient le thermique
-      do ig=1,ngrid
-         lmin(ig)=1
-      enddo
-      do ig=1,ngrid
-         do l=nlay,2,-1
-            if (ztv(ig,l-1).gt.ztv(ig,l)) then
-               lmin(ig)=l-1
-            endif
-         enddo
-      enddo
-c
-c definition de l'entrainement des couches
-      do l=1,klev-1
-         do ig=1,ngrid 
-            if (ztv(ig,l).gt.ztv(ig,l+1).and.
-     s          l.ge.lmin(ig).and.l.le.lentr(ig)) then 
-                 entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
-     s                           (zlev(ig,l+1)-zlev(ig,l))
-            endif
-         enddo
-      enddo
-c pas de thermique si couches 1->5 stables
-      do ig=1,ngrid
-         if (lmin(ig).gt.5) then
-            do l=1,klev
-               entr_star(ig,l)=0.
-            enddo
-         endif
-      enddo 
-c calcul de l entrainement total
-      do ig=1,ngrid
-         entr_star_tot(ig)=0.
-      enddo
-      do ig=1,ngrid
-         do k=1,klev
-            entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
-         enddo
-      enddo
-c
-      print*,'fin calcul entr_star'
-      do k=1,klev
-         do ig=1,ngrid 
-            ztva(ig,k)=ztv(ig,k)
-         enddo
-      enddo
-cRC
-c      print*,'7 OK convect8'
-      do k=1,klev+1
-         do ig=1,ngrid
-            zw2(ig,k)=0.
-            fmc(ig,k)=0.
-cCR
-            f_star(ig,k)=0.
-cRC
-            larg_cons(ig,k)=0.
-            larg_detr(ig,k)=0.
-            wa_moy(ig,k)=0.
-         enddo
-      enddo
-
-c      print*,'8 OK convect8'
-      do ig=1,ngrid
-         linter(ig)=1.
-         lmaxa(ig)=1
-         lmix(ig)=1
-         wmaxa(ig)=0.
-      enddo
-
-cCR: 
-      do l=1,nlay-2
-         do ig=1,ngrid
-            if (ztv(ig,l).gt.ztv(ig,l+1)
-     s         .and.entr_star(ig,l).gt.1.e-10
-     s         .and.zw2(ig,l).lt.1e-10) then
-               f_star(ig,l+1)=entr_star(ig,l)
-ctest:calcul de dteta
-               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
-     s                     *(zlev(ig,l+1)-zlev(ig,l))
-     s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
-               larg_detr(ig,l)=0.
-            else if ((zw2(ig,l).ge.1e-10).and.
-     s               (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)
-     s                    *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+
-     s                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
-     s                     *(zlev(ig,l+1)-zlev(ig,l))
-            endif
-c determination de zmax continu par interpolation lineaire
-            if (zw2(ig,l+1).lt.0.) then
-ctest
-               if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
-                  print*,'pb linter'
-               endif
-               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
-     s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
-               zw2(ig,l+1)=0.
-               lmaxa(ig)=l
-            else
-               if (zw2(ig,l+1).lt.0.) then
-                  print*,'pb1 zw2<0'
-               endif
-               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
-            endif
-            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
-c   lmix est le niveau de la couche ou w (wa_moy) est maximum
-               lmix(ig)=l+1
-               wmaxa(ig)=wa_moy(ig,l+1)
-            endif
-         enddo
-      enddo
-      print*,'fin calcul zw2'
-c
-c 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
-c pas de thermique si couches 1->5 stables
-      do ig=1,ngrid
-         if (lmin(ig).gt.5) then
-            lmax(ig)=1
-            lmin(ig)=1
-         endif
-      enddo 
-c    
-c Determination de zw2 max
-      do ig=1,ngrid
-         wmax(ig)=0.
-      enddo
-
-      do l=1,nlay
-         do ig=1,ngrid
-            if (l.le.lmax(ig)) then
-                if (zw2(ig,l).lt.0.)then
-                  print*,'pb2 zw2<0'
-                endif
-                zw2(ig,l)=sqrt(zw2(ig,l))
-                wmax(ig)=max(wmax(ig),zw2(ig,l))
-            else
-                 zw2(ig,l)=0.
-            endif
-          enddo
-      enddo
-
-c   Longueur caracteristique correspondant a la hauteur des thermiques.
-      do  ig=1,ngrid
-         zmax(ig)=0.
-         zlevinter(ig)=zlev(ig,1)
-      enddo
-      do  ig=1,ngrid
-c calcul de zlevinter
-          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
-     s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
-     s    -zlev(ig,lmax(ig)))
-       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
-      enddo
-
-      print*,'avant fermeture'
-c Fermeture,determination de f
-      do ig=1,ngrid
-         entr_star2(ig)=0.
-      enddo
-      do ig=1,ngrid
-         if (entr_star_tot(ig).LT.1.e-10) then
-            f(ig)=0.
-         else
-             do k=lmin(ig),lentr(ig)
-                entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
-     s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
-             enddo
-c Nouvelle fermeture
-             f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
-     s             *entr_star2(ig))*entr_star_tot(ig)
-ctest
-c             if (first) then
-c             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
-c     s             *wmax(ig))
-c             endif
-         endif
-c         f0(ig)=f(ig)
-c         first=.true.
-      enddo
-      print*,'apres fermeture'
-
-c Calcul de l'entrainement
-       do k=1,klev
-         do ig=1,ngrid 
-            entr(ig,k)=f(ig)*entr_star(ig,k)
-         enddo
-      enddo
-c Calcul des flux
-      do ig=1,ngrid
-         do l=1,lmax(ig)-1
-            fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
-         enddo
-      enddo
-
-cRC
-
-
-c      print*,'9 OK convect8'
-c     print*,'WA1 ',wa_moy
-
-c   determination de l'indice du debut de la mixed layer ou w decroit
-
-c   calcul de la largeur de chaque ascendance dans le cas conservatif.
-c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
-c   d'une couche est égale à la hauteur de la couche alimentante.
-c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
-c   de la vitesse d'entrainement horizontal dans la couche alimentante.
-
-      do l=2,nlay
-         do ig=1,ngrid
-            if (l.le.lmaxa(ig)) then
-               zw=max(wa_moy(ig,l),1.e-10)
-               larg_cons(ig,l)=zmax(ig)*r_aspect
-     s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
-            endif
-         enddo
-      enddo
-
-      do l=2,nlay
-         do ig=1,ngrid
-            if (l.le.lmaxa(ig)) then
-c              if (idetr.eq.0) then
-c  cette option est finalement en dur.
-                  if ((l_mix*zlev(ig,l)).lt.0.)then
-                   print*,'pb l_mix*zlev<0'
-                  endif
-                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
-c              else if (idetr.eq.1) then
-c                 larg_detr(ig,l)=larg_cons(ig,l)
-c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
-c              else if (idetr.eq.2) then
-c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
-c    s            *sqrt(wa_moy(ig,l))
-c              else if (idetr.eq.4) then
-c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
-c    s            *wa_moy(ig,l)
-c              endif
-            endif
-         enddo
-       enddo
-
-c      print*,'10 OK convect8'
-c     print*,'WA2 ',wa_moy
-c   calcul de la fraction de la maille concernée par l'ascendance en tenant
-c   compte de l'epluchage du thermique.
-c
-cCR def de  zmix continu (profil parabolique des vitesses)
-      do ig=1,ngrid
-           if (lmix(ig).gt.1.) then
-c test 
-              if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
-     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
-     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
-     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
-     s        then
-c             
-            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
-     s        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
-     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
-     s        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
-     s        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
-     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
-     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
-     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
-            else
-            zmix(ig)=zlev(ig,lmix(ig))
-            print*,'pb zmix'
-            endif
-         else 
-         zmix(ig)=0.
-         endif
-ctest
-         if ((zmax(ig)-zmix(ig)).lt.0.) then
-            zmix(ig)=0.99*zmax(ig)
-c            print*,'pb zmix>zmax'
-         endif
-      enddo
-c
-c calcul du nouveau lmix correspondant
-      do ig=1,ngrid
-         do l=1,klev
-            if (zmix(ig).ge.zlev(ig,l).and.
-     s          zmix(ig).lt.zlev(ig,l+1)) then
-              lmix(ig)=l
-             endif
-          enddo
-      enddo
-c
-      do l=2,nlay
-         do ig=1,ngrid
-            if(larg_cons(ig,l).gt.1.) then
-c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
-               fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
-     s            /(r_aspect*zmax(ig))
-c test
-               fraca(ig,l)=max(fraca(ig,l),0.)
-               fraca(ig,l)=min(fraca(ig,l),0.5)
-               fracd(ig,l)=1.-fraca(ig,l)
-               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
-            else
-c              wa_moy(ig,l)=0.
-               fraca(ig,l)=0.
-               fracc(ig,l)=0.
-               fracd(ig,l)=1.
-            endif
-         enddo
-      enddo                  
-cCR: calcul de fracazmix
-       do ig=1,ngrid
-          fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
-     s     (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
-     s    +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
-     s    -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
-       enddo
-c
-       do l=2,nlay
-          do ig=1,ngrid
-             if(larg_cons(ig,l).gt.1.) then
-               if (l.gt.lmix(ig)) then
-ctest
-                 if (zmax(ig)-zmix(ig).lt.1.e-10) then
-c                   print*,'pb xxx'
-                   xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
-                 else
-                 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
-                 endif
-           if (idetr.eq.0) then
-               fraca(ig,l)=fracazmix(ig)
-           else if (idetr.eq.1) then
-               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
-           else if (idetr.eq.2) then
-               fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
-           else
-               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
-           endif
-c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
-               fraca(ig,l)=max(fraca(ig,l),0.)
-               fraca(ig,l)=min(fraca(ig,l),0.5)
-               fracd(ig,l)=1.-fraca(ig,l)
-               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
-             endif
-            endif
-         enddo
-      enddo
-      
-      print*,'fin calcul fraca'
-c      print*,'11 OK convect8'
-c     print*,'Ea3 ',wa_moy
-c------------------------------------------------------------------
-c   Calcul de fracd, wd
-c   somme wa - wd = 0
-c------------------------------------------------------------------
-
-
-      do ig=1,ngrid
-         fm(ig,1)=0.
-         fm(ig,nlay+1)=0.
-      enddo
-
-      do l=2,nlay
-           do ig=1,ngrid
-              fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
-cCR:test
-              if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
-     s            .and.l.gt.lmix(ig)) then
-                 fm(ig,l)=fm(ig,l-1)
-c                 write(1,*)'ajustement fm, l',l
-              endif
-c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
-cRC
-           enddo
-         do ig=1,ngrid
-            if(fracd(ig,l).lt.0.1) then
-              abort_message = 'fracd trop petit'
-              CALL abort_gcm (modname,abort_message,1)
-            else
-c    vitesse descendante "diagnostique"
-               wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
-            endif
-         enddo
-      enddo
-
-      do l=1,nlay
-         do ig=1,ngrid
-c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
-            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
-         enddo
-      enddo
-
-      print*,'12 OK convect8'
-c     print*,'WA4 ',wa_moy
-cc------------------------------------------------------------------
-c   calcul du transport vertical
-c------------------------------------------------------------------
-
-      go to 4444
-c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
-      do l=2,nlay-1
-         do ig=1,ngrid
-            if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
-     s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
-c     print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
-c    s         ,fm(ig,l+1)*ptimestep
-c    s         ,'   M=',masse(ig,l),masse(ig,l+1)
-             endif
-         enddo
-      enddo
-
-      do l=1,nlay
-         do ig=1,ngrid
-            if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
-c     print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
-c    s         ,entr(ig,l)*ptimestep
-c    s         ,'   M=',masse(ig,l)
-             endif
-         enddo
-      enddo
-
-      do l=1,nlay
-         do ig=1,ngrid
-            if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
-c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
-c    s         ,'   FM=',fm(ig,l)
-            endif
-            if(.not.masse(ig,l).ge.1.e-10
-     s         .or..not.masse(ig,l).le.1.e4) then
-c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
-c    s         ,'   M=',masse(ig,l)
-c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
-c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
-c     print*,'zlev(ig,l+1),zlev(ig,l)'
-c    s                ,zlev(ig,l+1),zlev(ig,l)
-c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
-c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
-            endif
-            if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
-c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
-c    s         ,'   E=',entr(ig,l)
-            endif
-         enddo
-      enddo
-
-4444   continue
-
-cCR:redefinition du entr
-       do l=1,nlay
-         do ig=1,ngrid
-            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
-            if (detr(ig,l).lt.0.) then
-                entr(ig,l)=entr(ig,l)-detr(ig,l)
-                detr(ig,l)=0.
-c     print*,'WARNING !!! detrainement negatif ',ig,l
-            endif
-         enddo
-      enddo
-cRC
-      if (w2di.eq.1) then
-         fm0=fm0+ptimestep*(fm-fm0)/tho
-         entr0=entr0+ptimestep*(entr-entr0)/tho
-      else
-         fm0=fm
-         entr0=entr
-      endif
-
-      if (1.eq.1) then
-         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
-     .    ,zh,zdhadj,zha)
-         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
-     .    ,zo,pdoadj,zoa)
-      else
-         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
-     .    ,zh,zdhadj,zha)
-         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
-     .    ,zo,pdoadj,zoa)
-      endif
-
-      if (1.eq.0) then
-         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
-     .    ,fraca,zmax
-     .    ,zu,zv,pduadj,pdvadj,zua,zva)
-      else
-         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
-     .    ,zu,pduadj,zua)
-         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
-     .    ,zv,pdvadj,zva)
-      endif
-
-      do l=1,nlay
-         do ig=1,ngrid
-            zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
-            zf2=zf/(1.-zf)
-            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
-            wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
-         enddo
-      enddo
-
-
-
-c     print*,'13 OK convect8'
-c     print*,'WA5 ',wa_moy
-      do l=1,nlay
-         do ig=1,ngrid
-            pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
-         enddo
-      enddo
-
-
-c     do l=1,nlay
-c        do ig=1,ngrid
-c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
-c     print*,'WARN!!! ig=',ig,'  l=',l
-c    s         ,'   pdtadj=',pdtadj(ig,l)
-c           endif
-c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
-c     print*,'WARN!!! ig=',ig,'  l=',l
-c    s         ,'   pdoadj=',pdoadj(ig,l)
-c           endif
-c        enddo
-c      enddo
-
-      print*,'14 OK convect8'
-c------------------------------------------------------------------
-c   Calculs pour les sorties
-c------------------------------------------------------------------
-
-      if(sorties) then
-      do l=1,nlay
-         do ig=1,ngrid
-            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
-            zld(ig,l)=fracd(ig,l)*zmax(ig)
-            if(1.-fracd(ig,l).gt.1.e-10)
-     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
-         enddo
-      enddo
-
-cdeja fait
-c      do l=1,nlay
-c         do ig=1,ngrid
-c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
-c            if (detr(ig,l).lt.0.) then
-c                entr(ig,l)=entr(ig,l)-detr(ig,l)
-c                detr(ig,l)=0.
-c     print*,'WARNING !!! detrainement negatif ',ig,l
-c            endif
-c         enddo
-c      enddo
-
-c     print*,'15 OK convect8'
-
-      isplit=isplit+1
-
-
-c #define und
-	goto 123
-#ifdef und
-      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
-      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
-      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
-      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
-      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
-      CALL writeg1d(1,nlay,zla,'la      ','la      ')
-      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
-      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
-      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
-      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
-      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
-      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
-      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
-      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
-      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
-      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
-      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
-      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
-      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
-      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
-      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
-      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
-      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
-      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
-
-      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
-      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
-      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
-
-c   recalcul des flux en diagnostique...
-c     print*,'PAS DE TEMPS ',ptimestep
-       call dt2F(pplev,pplay,pt,pdtadj,wh)
-      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
-#endif
-123   continue
-#define troisD
-#ifdef troisD
-c       if (sorties) then
-      print*,'Debut des wrgradsfi'
-
-c      print*,'16 OK convect8'
-         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
-         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
-         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
-         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
-         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
-         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
-c      print*,'WA6 ',wa_moy
-         call wrgradsfi(1,nlay,zla,'la        ','la        ')
-         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
-         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
-         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
-         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
-         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
-         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
-         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
-         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
-         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
-         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
-         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
-         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
-         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
-         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
-         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
-         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
-         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
-         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
-         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
-         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
-         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
-         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
-         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
-         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
-         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
-
-         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
-         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
-         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
-
-cCR:nouveaux diagnostiques
-      call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
-      call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
-      call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
-      call wrgradsfi(1,1,zmix,'zmix      ','zmix      ') 
-      zsortie1d(:)=lmax(:)
-      call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
-      call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
-      zsortie1d(:)=lmix(:)
-      call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
-      zsortie1d(:)=lentr(:)
-      call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
-
-c      print*,'17 OK convect8'
-
-         do k=1,klev/10
-            write(str2,'(i2.2)') k
-            str10='wa'//str2
-            do l=1,nlay
-               do ig=1,ngrid
-                  zsortie(ig,l)=wa(ig,k,l)
-               enddo
-            enddo   
-            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
-            do l=1,nlay
-               do ig=1,ngrid
-                  zsortie(ig,l)=larg_part(ig,k,l)
-               enddo
-            enddo
-            str10='la'//str2
-            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
-         enddo
-
-
-c     print*,'18 OK convect8'
-c      endif
-      print*,'Fin des wrgradsfi'
-#endif
-
-      endif
-
-c     if(wa_moy(1,4).gt.1.e-10) stop
-
-      print*,'19 OK convect8'
+!     print*,'19 OK convect8'
       return
       end
@@ -5239,9 +4812,5 @@
 
       real count_time
-      integer isplit,nsplit,ialt
-      parameter (nsplit=10)
-      data isplit/0/
-      save isplit
-c$OMP THREADPRIVATE(isplit)
+      integer ialt
 
       logical sorties
@@ -6007,117 +5576,5 @@
       enddo
 
-cdeja fait
-c      do l=1,nlay
-c         do ig=1,ngrid
-c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
-c            if (detr(ig,l).lt.0.) then
-c                entr(ig,l)=entr(ig,l)-detr(ig,l)
-c                detr(ig,l)=0.
-c     print*,'WARNING !!! detrainement negatif ',ig,l
-c            endif
-c         enddo
-c      enddo
-
-c     print*,'15 OK convect8'
-
-      isplit=isplit+1
-
-
-c #define und
-	goto 123
-#ifdef und
-      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
-      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
-      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
-      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
-      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
-      CALL writeg1d(1,nlay,zla,'la      ','la      ')
-      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
-      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
-      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
-      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
-      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
-      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
-      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
-      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
-      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
-      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
-      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
-      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
-      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
-      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
-      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
-      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
-      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
-      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
-
-      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
-      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
-      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
-
-c   recalcul des flux en diagnostique...
-c     print*,'PAS DE TEMPS ',ptimestep
-       call dt2F(pplev,pplay,pt,pdtadj,wh)
-      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
-#endif
-123   continue
-! #define troisD
-#ifdef troisD
-c       if (sorties) then
-      print*,'Debut des wrgradsfi'
-
-c      print*,'16 OK convect8'
-         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
-         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
-         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
-         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
-         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
-         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
-c      print*,'WA6 ',wa_moy
-         call wrgradsfi(1,nlay,zla,'la        ','la        ')
-         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
-         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
-         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
-         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
-         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
-         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
-         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
-         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
-         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
-         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
-         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
-         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
-         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
-         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
-         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
-         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
-         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
-         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
-         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
-         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
-         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
-         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
-         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
-         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
-         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
-
-         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
-         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
-         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
-
-cCR:nouveaux diagnostiques
-      call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
-      call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
-      call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
-      call wrgradsfi(1,1,zmix,'zmix      ','zmix      ') 
-      zsortie1d(:)=lmax(:)
-      call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
-      call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
-      zsortie1d(:)=lmix(:)
-      call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
-      zsortie1d(:)=lentr(:)
-      call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
-
-c      print*,'17 OK convect8'
+
 
          do k=1,klev/10
@@ -6139,15 +5596,8 @@
          enddo
 
-
-c     print*,'18 OK convect8'
-c      endif
-      print*,'Fin des wrgradsfi'
-#endif
-
       endif
 
-c     if(wa_moy(1,4).gt.1.e-10) stop
-
-c      print*,'19 OK convect8'
+
+
       return
       end
