Index: /LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- /LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 4089)
@@ -130,5 +130,5 @@
 
   INCLUDE "compbl.h"
-  INCLUDE "thermcell.h"
+  INCLUDE "alpale.h"
   
   deg2rad= pi/180.0
Index: /LMDZ6/trunk/libf/phylmd/alpale.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/alpale.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/alpale.F90	(revision 4089)
@@ -48,5 +48,5 @@
   REAL, DIMENSION(klon), INTENT(OUT)                         :: Ale_wake, Alp_wake
 
-  include "thermcell.h"
+  include "alpale.h"
   include "YOMCST.h"
   include "YOETHF.h"
Index: /LMDZ6/trunk/libf/phylmd/alpale_th.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/alpale_th.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/alpale_th.F90	(revision 4089)
@@ -47,5 +47,5 @@
   REAL, DIMENSION(klon), INTENT(OUT)                         :: birth_rate
 
-  include "thermcell.h"
+  include "alpale.h"
 
 ! Local variables
Index: /LMDZ6/trunk/libf/phylmd/calltherm.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/calltherm.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/calltherm.F90	(revision 4089)
@@ -7,5 +7,5 @@
      &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs  &
      &      ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,  &
-     &       ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, &
+     &       ratqsdiff,zqsatth,ale_bl,alp_bl,lalim_conv,wght_th, &
      &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl &
 !!! nrlmd le 10/04/2012
@@ -17,12 +17,28 @@
      &      ,alp_bl_conv,alp_bl_stat &
 !!! fin nrlmd le 10/04/2012
-     &      ,zqla,ztva )
+     &      ,zqla,ztva &
+#ifdef ISO         
+     &      ,xt_seri,d_xt_ajs &
+#ifdef DIAGISO         
+     &      ,q_the,xt_the &
+#endif
+#endif         
+     &   )
 
       USE dimphy
       USE indice_sol_mod
       USE print_control_mod, ONLY: prt_level,lunout
+#ifdef ISO
+      use infotrac_phy, ONLY: ntraciso
+#ifdef ISOVERIF
+      USE isotopes_mod, ONLY: iso_eau,iso_HDO
+      USE isotopes_verif_mod, ONLY: iso_verif_aberrant_enc_vect2D, &
+        iso_verif_egalite_vect2D
+#endif   
+#endif
 
       implicit none
-      include "thermcell.h"
+      include "clesphys.h"
+      include "thermcell_old.h"
 
 
@@ -79,8 +95,8 @@
       real zqsatth(klon,klev)  
 !nouvelles variables pour la convection
-      real Ale_bl(klon)
-      real Alp_bl(klon)
-      real Ale(klon)
-      real Alp(klon)
+      real ale_bl(klon)
+      real alp_bl(klon)
+      real ale(klon)
+      real alp(klon)
 !RC
       !on garde le zmax du pas de temps precedent
@@ -102,4 +118,13 @@
 !********************************************************
 
+      real, dimension(klon) :: pcon
+      real, dimension(klon,klev) :: rhobarz,wth3
+      integer,dimension(klon) :: lalim
+      real, dimension(klon,klev+1) :: fm
+      real, dimension(klon,klev) :: alim_star
+      real, dimension(klon) :: zmax
+
+
+
 
 ! variables locales
@@ -115,6 +140,21 @@
       character (len=80) :: abort_message
 
-      integer i,k
+      integer i,k,isplit
       logical, save :: first=.true.
+      logical :: new_thermcell
+
+#ifdef ISO
+      REAL xt_seri(ntraciso,klon,klev),xtmemoire(ntraciso,klon,klev)
+      REAL d_xt_ajs(ntraciso,klon,klev)
+      real d_xt_the(ntraciso,klon,klev)
+#ifdef DIAGISO
+      real q_the(klon,klev)
+      real xt_the(ntraciso,klon,klev)
+#endif
+      real qprec(klon,klev)
+      integer ixt
+#endif
+
+
 !$OMP THREADPRIVATE(first)
 !********************************************************
@@ -144,6 +184,6 @@
          detr_therm(:,:)=0.
 
-         Ale_bl(:)=0.
-         Alp_bl(:)=0.
+         ale_bl(:)=0.
+         alp_bl(:)=0.
          if (prt_level.ge.10) then
           print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
@@ -159,6 +199,15 @@
                 logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
                if (logexpr2(i,k)) then
+#ifdef ISO
+                qprec(i,k)=q_seri(i,k)
+#endif
                 q_seri(i,k)=1.e-15
                 nbptspb=nbptspb+1
+#ifdef ISO
+                do ixt=1,ntraciso
+                  xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k))
+                  ! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt)))
+                enddo
+#endif
                endif
 !               if (logexpr0) &
@@ -169,5 +218,21 @@
          if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
 
+
+         new_thermcell=iflag_thermals>=15.and.iflag_thermals<=18
+#ifdef ISO
+      if (.not.new_thermcell) then
+           CALL abort_gcm('calltherm 234','isos pas prevus ici',1)
+      endif
+#ifdef ISOVERIF
+      if (iso_eau.gt.0) then
+       call iso_verif_egalite_vect2D( &
+     &           xt_seri,q_seri, &
+     &           'calltherm 174',ntraciso,klon,klev) 
+      endif !if (iso_eau.gt.0) then
+#endif    
+#endif
          zdt=dtime/REAL(nsplit_thermals)
+
+
          do isplit=1,nsplit_thermals
 
@@ -207,12 +272,4 @@
               abort_message = 'cas non prevu dans calltherm'
               CALL abort_physic (modname,abort_message,1)
-
-!           CALL thermcell_pluie(klon,klev,zdt  &
-!   &      ,pplay,paprs,pphi,zlev  &
-!    &      ,u_seri,v_seri,t_seri,q_seri  &
-!    &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
-!    &      ,zfm_therm,zentr_therm,zqla  &
-!    &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
-!    &      ,tau_thermals,3)
           else if (iflag_thermals.eq.12) then
             CALL calcul_sec(klon,klev,zdt  &
@@ -223,16 +280,7 @@
      &      ,tau_thermals)
           else if (iflag_thermals==13.or.iflag_thermals==14) then
-            CALL thermcellV0_main(itap,klon,klev,zdt  &
-     &      ,pplay,paprs,pphi,debut  &
-     &      ,u_seri,v_seri,t_seri,q_seri  &
-     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
-     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
-     &      ,ratqscth,ratqsdiff,zqsatth  &
-     &      ,r_aspect_thermals,l_mix_thermals  &
-     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
-     &      ,zmax0,f0,zw2,fraca)
-          else if (iflag_thermals>=15.and.iflag_thermals<=18) then
-
-!            print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
+              abort_message = 'thermcellV0_main enleve svn>2084'
+              CALL abort_physic (modname,abort_message,1)
+          else if (new_thermcell) then
             CALL thermcell_main(itap,klon,klev,zdt  &
      &      ,pplay,paprs,pphi,debut  &
@@ -241,18 +289,26 @@
      &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
      &      ,ratqscth,ratqsdiff,zqsatth  &
-!    &      ,r_aspect_thermals,l_mix_thermals &
-!    &      ,tau_thermals,iflag_thermals_ed,iflag_coupl &
-     &      ,Ale,Alp,lalim_conv,wght_th &
      &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
-     &      ,ztla,zthl &
-!!! nrlmd le 10/04/2012
-     &      ,pbl_tke,pctsrf,omega,airephy &
-     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
-     &      ,n2,s2,ale_bl_stat &
-     &      ,therm_tke_max,env_tke_max &
-     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
-     &      ,alp_bl_conv,alp_bl_stat &
-!!! fin nrlmd le 10/04/2012
-     &      ,ztva ) 
+     &      ,ztla,zthl,ztva &
+     &      ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
+#ifdef ISO         
+     &      ,xt_seri,d_xt_the &
+#endif         
+     &   ) 
+
+            CALL thermcell_alp(klon,klev,zdt  &                      ! in
+     &        ,pplay,paprs  &                                        ! in
+     &        ,zfm_therm,zentr_therm,lmax  &                         ! in
+     &        ,pbl_tke,pctsrf,omega,airephy &                        ! in
+     &        ,zw2,fraca &                                           ! in
+     &        ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
+     &        ,ale,alp,lalim_conv,wght_th &                          ! out
+     &        ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out
+     &        ,n2,s2,ale_bl_stat &                                   ! out
+     &        ,therm_tke_max,env_tke_max &                           ! out
+     &        ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
+     &        ,alp_bl_conv,alp_bl_stat &                             ! out
+     &        )
+
            if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
          else
@@ -303,4 +359,9 @@
             detr_therm(:,k)=detr_therm(:,k)  &
      &       +zdetr_therm(:,k)*fact(:)
+#ifdef ISO
+            do ixt=1,ntraciso
+              d_xt_the(ixt,:,k)=d_xt_the(ixt,:,k)*dtime*fact(:)
+            enddo
+#endif
       ENDDO
        fm_therm(:,klev+1)=0.
@@ -313,4 +374,7 @@
             d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
             d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
+#ifdef ISO
+            d_xt_ajs(:,:,:)=d_xt_ajs(:,:,:)+d_xt_the(:,:,:)
+#endif
 
 !  incrementation des variables meteo
@@ -320,13 +384,31 @@
             qmemoire(:,:)=q_seri(:,:)
             q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
+#ifdef ISO
+            xtmemoire(:,:,:)=xt_seri(:,:,:)
+            xt_seri(:,:,:) = xt_seri(:,:,:) + d_xt_the(:,:,:)
+#ifdef ISOVERIF
+!      write(*,*) 'calltherm 350 tmp: ajout d_xt_the'
+      if (iso_HDO.gt.0) then
+!      i=479
+!      k=4
+!      write(*,*) 'xt_seri(iso_hdo,i,k),q_seri(i,k)=', &
+!     &   xt_seri(iso_hdo,i,k),q_seri(i,k)
+!      write(*,*) 'd_xt_the(iso_hdo,i,k),d_q_the(i,k)=', &
+!     &   d_xt_the(iso_hdo,i,k),d_q_the(i,k)
+      call iso_verif_aberrant_enc_vect2D( &
+     &        xt_seri,q_seri, &
+     &        'calltherm 353, apres ajout d_xt_the',ntraciso,klon,klev)
+      endif      
+#endif
+#endif 
            if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
 
        DO i=1,klon
             fm_therm(i,klev+1)=0.
-            Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
-!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
-            Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
-!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
-        if(prt_level.GE.10) print*,'calltherm i Alp_bl Alp Ale_bl Ale',i,Alp_bl(i),Alp(i),Ale_bl(i),Ale(i)
+            ale_bl(i)=ale_bl(i)+ale(i)/REAL(nsplit_thermals)
+!            write(22,*)'ALE CALLTHERM',ale_bl(i),ale(i)
+            alp_bl(i)=alp_bl(i)+alp(i)/REAL(nsplit_thermals)
+!            write(23,*)'ALP CALLTHERM',alp_bl(i),alp(i)
+        if(prt_level.GE.10) print*,'calltherm i alp_bl alp ale_bl ale',i,alp_bl(i),alp(i),ale_bl(i),ale(i)
        ENDDO
 
@@ -341,4 +423,9 @@
                 q_seri(i,k)=1.e-15
                 nbptspb=nbptspb+1
+#ifdef ISO
+                do ixt=1,ntraciso
+                  xt_seri(ixt,i,k)=1.e-15*(xtmemoire(ixt,i,k)/qmemoire(i,k))
+                enddo
+#endif
 !                if (prt_level.ge.10) then
 !                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
@@ -348,4 +435,14 @@
             ENDDO
             ENDDO
+#ifdef ISO
+#ifdef ISOVERIF
+      if (iso_HDO.gt.0) then
+      call iso_verif_aberrant_enc_vect2D( &
+     &        xt_seri,q_seri, &
+     &        'calltherm 393, apres bidouille q<0',ntraciso,klon,klev)
+      endif      
+#endif
+#endif
+
         IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
 ! tests sur les valeurs de la temperature
Index: /LMDZ6/trunk/libf/phylmd/calwake.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/calwake.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/calwake.F90	(revision 4089)
@@ -210,4 +210,16 @@
 
 
+! SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, &
+!                 tenv0, qe0, omgb, &
+!                 dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
+!                 sigd_con, Cin, &
+!                 deltatw, deltaqw, sigmaw, awdens, wdens, &                  ! state variables
+!                 dth, hw, wape, fip, gfl, &
+!                 dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
+!                 dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
+!                 d_deltat_gw, &
+!                 d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)     ! tendencies
+! 
+  ! retour a un Pupper fixe                                *
   CALL wake(klon,klev,znatsurf, p, ph, pi, dtime, &
     te, qe, omgbe, &
Index: /LMDZ6/trunk/libf/phylmd/clesphys.h
===================================================================
--- /LMDZ6/trunk/libf/phylmd/clesphys.h	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/clesphys.h	(revision 4089)
@@ -15,4 +15,7 @@
        INTEGER nbapp_rad, iflag_con, nbapp_cv, nbapp_wk, iflag_ener_conserv
        REAL co2_ppm, co2_ppm0, solaire
+       INTEGER iflag_thermals,nsplit_thermals
+       REAL tau_thermals
+
 !FC
        REAL Cd_frein
@@ -145,6 +148,6 @@
      &     , ok_chlorophyll,ok_conserv_q, adjust_tropopause             &
      &     , ok_daily_climoz, ok_all_xml, ok_lwoff                      &
-     &     , iflag_phytrac, ok_new_lscp
-     
+     &     , iflag_phytrac, ok_new_lscp                                 &
+     &     ,  iflag_thermals,nsplit_thermals, tau_thermals    
        save /clesphys/
 !$OMP THREADPRIVATE(/clesphys/)
Index: /LMDZ6/trunk/libf/phylmd/cloudth_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 4089)
@@ -24,5 +24,4 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-#include "thermcell.h"
 #include "nuage.h"
 
@@ -269,5 +268,4 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-#include "thermcell.h"
 #include "nuage.h"
       
@@ -609,5 +607,4 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-#include "thermcell.h"
 #include "nuage.h"
 
@@ -833,5 +830,4 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-#include "thermcell.h"
 #include "nuage.h"
       
@@ -1295,5 +1291,4 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-#include "thermcell.h"
 #include "nuage.h"
 
@@ -1562,5 +1557,4 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-#include "thermcell.h"
 #include "nuage.h"
       
Index: /LMDZ6/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 4089)
@@ -39,5 +39,5 @@
     INCLUDE "YOMCST.h"
     INCLUDE "YOMCST2.h"
-    INCLUDE "thermcell.h"
+    INCLUDE "alpale.h"
 
     !IM : on inclut/initialise les taux de CH4, N2O, CFC11 et CFC12
@@ -142,6 +142,4 @@
     REAL,SAVE :: seuil_inversion_omp
 
-    INTEGER,SAVE :: iflag_thermals_ed_omp,iflag_thermals_optflux_omp,iflag_thermals_closure_omp
-    REAL, SAVE :: fact_thermals_ed_dz_omp
     INTEGER,SAVE :: iflag_thermals_omp,nsplit_thermals_omp
     REAL,SAVE :: tau_thermals_omp,alp_bl_k_omp
@@ -1718,42 +1716,7 @@
     CALL getin('iflag_thermals',iflag_thermals_omp)
     !
-    !Config Key  = iflag_thermals_ed
-    !Config Desc =
-    !Config Def  = 0
-    !Config Help =
-    !
-    fact_thermals_ed_dz_omp = 0.1
-
-    CALL getin('fact_thermals_ed_dz',fact_thermals_ed_dz_omp)
-    !
-    !
-    !Config Key  = iflag_thermals_ed
-    !Config Desc =
-    !Config Def  = 0
-    !Config Help =
-    !
-    iflag_thermals_ed_omp = 0
-    CALL getin('iflag_thermals_ed',iflag_thermals_ed_omp)
-    !
-    !
-    !Config Key  = iflag_thermals_optflux
-    !Config Desc =
-    !Config Def  = 0
-    !Config Help =
-    !
-    iflag_thermals_optflux_omp = 0
-    CALL getin('iflag_thermals_optflux',iflag_thermals_optflux_omp)
-    !
-    !Config Key  = iflag_thermals_closure
-    !Config Desc =
-    !Config Def  = 0
-    !Config Help =
-    !
-    iflag_thermals_closure_omp = 1
-    CALL getin('iflag_thermals_closure',iflag_thermals_closure_omp)
-    !
     !Config Key  = nsplit_thermals
     !Config Desc =
-    !Config Def  = 1
+    !Config Def  = 0
     !Config Help =
     !
@@ -2633,8 +2596,4 @@
     ip_ebil_phy = ip_ebil_phy_omp
     iflag_thermals = iflag_thermals_omp
-    iflag_thermals_ed = iflag_thermals_ed_omp
-    fact_thermals_ed_dz = fact_thermals_ed_dz_omp
-    iflag_thermals_optflux = iflag_thermals_optflux_omp
-    iflag_thermals_closure = iflag_thermals_closure_omp
     nsplit_thermals = nsplit_thermals_omp
     tau_thermals = tau_thermals_omp
@@ -3002,8 +2961,4 @@
     WRITE(lunout,*) ' iflag_order2_sollw = ', iflag_order2_sollw
     WRITE(lunout,*) ' iflag_thermals = ', iflag_thermals
-    WRITE(lunout,*) ' iflag_thermals_ed = ', iflag_thermals_ed
-    WRITE(lunout,*) ' fact_thermals_ed_dz = ', fact_thermals_ed_dz
-    WRITE(lunout,*) ' iflag_thermals_optflux = ', iflag_thermals_optflux
-    WRITE(lunout,*) ' iflag_thermals_closure = ', iflag_thermals_closure
     WRITE(lunout,*) ' iflag_clos = ', iflag_clos
     WRITE(lunout,*) ' coef_clos_ls = ', coef_clos_ls
Index: /LMDZ6/trunk/libf/phylmd/phyetat0.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phyetat0.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/phyetat0.F90	(revision 4089)
@@ -43,5 +43,5 @@
   include "dimsoil.h"
   include "clesphys.h"
-  include "thermcell.h"
+  include "alpale.h"
   include "compbl.h"
   include "YOMCST.h"
Index: /LMDZ6/trunk/libf/phylmd/phyredem.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 4089)
@@ -47,5 +47,5 @@
   include "dimsoil.h"
   include "clesphys.h"
-  include "thermcell.h"
+  include "alpale.h"
   include "compbl.h"
   !======================================================================
Index: /LMDZ6/trunk/libf/phylmd/phys_output_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_mod.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_mod.F90	(revision 4089)
@@ -53,5 +53,4 @@
     IMPLICIT NONE
     include "clesphys.h"
-    include "thermcell.h"
     include "YOMCST.h"
 
Index: /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4089)
@@ -426,5 +426,5 @@
 
     INCLUDE "clesphys.h"
-    INCLUDE "thermcell.h"
+    INCLUDE "alpale.h"
     INCLUDE "compbl.h"
     INCLUDE "YOMCST.h"
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4089)
@@ -79,4 +79,5 @@
     USE lscp_mod, ONLY : lscp
     USE wake_ini_mod, ONLY : wake_ini
+    USE thermcell_ini_mod, ONLY : thermcell_ini
 
     !USE cmp_seri_mod
@@ -355,5 +356,5 @@
     include "dimsoil.h"
     include "clesphys.h"
-    include "thermcell.h"
+    include "alpale.h"
     include "dimpft.h"
     !======================================================================
@@ -1732,4 +1733,7 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        CALL wake_ini(rg,rd,rv,prt_level)
+       CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
+   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
Index: /LMDZ6/trunk/libf/phylmd/phytrac_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phytrac_mod.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/phytrac_mod.F90	(revision 4089)
@@ -153,5 +153,4 @@
     INCLUDE "YOMCST.h"
     INCLUDE "clesphys.h"
-    INCLUDE "thermcell.h"
     !==========================================================================
     !                   -- ARGUMENT DESCRIPTION --
Index: /LMDZ6/trunk/libf/phylmd/thermcell_alim.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/thermcell_alim.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/thermcell_alim.F90	(revision 4089)
@@ -11,10 +11,4 @@
 !--------------------------------------------------------------------------
 
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"
-#include "thermcell.h"
-
-!      fort(10) ptimestep,ztv,zthl,po,zl,rhobarz,zlev,pplev,pphi,zpspsk,f0
       INTEGER, INTENT(IN) :: ngrid,klev
       REAL, INTENT(IN) :: ztv(ngrid,klev)
@@ -41,5 +35,5 @@
 
 !-------------------------------------------------------------------------
-! Definition de l'alimentation a l'origine dans thermcell_init
+! Definition de l'alimentation
 !-------------------------------------------------------------------------
    IF (flag==0) THEN ! CMIP5 version
Index: /LMDZ6/trunk/libf/phylmd/thermcell_alp.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/thermcell_alp.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/thermcell_alp.F90	(revision 4089)
@@ -1,53 +1,31 @@
 ! $Id: thermcell_main.F90 2351 2015-08-25 15:14:59Z emillour $
 !
-      SUBROUTINE thermcell_alp(ngrid,nlay,ptimestep  &
-     &                  ,pplay,pplev  &
-     &                  ,fm0,entr0,lmax  &
-     &                  ,ale_bl,alp_bl,lalim_conv,wght_th &
-     &                  ,zw2,fraca &
-!!! ncessaire en plus
-     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
-!!! nrlmd le 10/04/2012
-     &                  ,pbl_tke,pctsrf,omega,airephy &
-     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
-     &                  ,n2,s2,ale_bl_stat &
-     &                  ,therm_tke_max,env_tke_max &
-     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
-     &                  ,alp_bl_conv,alp_bl_stat &
-!!! fin nrlmd le 10/04/2012
+      SUBROUTINE thermcell_alp(ngrid,nlay,ptimestep  &                         ! in
+     &                  ,pplay,pplev  &                                        ! in
+     &                  ,fm0,entr0,lmax  &                                     ! in
+     &                  ,pbl_tke,pctsrf,omega,airephy &                        ! in
+     &                  ,zw2,fraca &                                           ! in
+     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
+!
+     &                  ,ale_bl,alp_bl,lalim_conv,wght_th &                    ! out
+     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &   ! out
+     &                  ,n2,s2,ale_bl_stat &                                   ! out
+     &                  ,therm_tke_max,env_tke_max &                           ! out
+     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
+     &                  ,alp_bl_conv,alp_bl_stat &                             ! out
      &)
 
-      USE dimphy
       USE indice_sol_mod
       IMPLICIT NONE
 
 !=======================================================================
-!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
-!   Version du 09.02.07
-!   Calcul du transport vertical dans la couche limite en presence
-!   de "thermiques" explicitement representes avec processus nuageux
 !
-!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
-!
-!   le thermique est suppose homogene et dissipe par melange avec
-!   son environnement. la longueur l_mix controle l'efficacite du
-!   melange
-!
-!   Le calcul du transport des differentes especes se fait en prenant
-!   en compte:
-!     1. un flux de masse montant
-!     2. un flux de masse descendant
-!     3. un entrainement
-!     4. un detrainement
-!
-! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
-!    Introduction of an implicit computation of vertical advection in
-!    the environment of thermal plumes in thermcell_dq
-!    impl =     0 : explicit, 1 : implicit, -1 : old version
-!    controled by iflag_thermals =
-!       15, 16 run with impl=-1 : numerical convergence with NPv3
-!       17, 18 run with impl=1  : more stable
-!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
-!
+!   Auteurs: Catherine Rio
+!   Modifications :
+!   Nicolas Rochetin et Jean-Yves Grandpeix
+!         pour la fermeture stochastique. 2012
+!   Frédéric Hourdin :
+!         netoyage informatique. 2022
+!   
 !=======================================================================
 !-----------------------------------------------------------------------
@@ -58,106 +36,74 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-#include "thermcell.h"
+#include "alpale.h"
 
 !   arguments:
 !   ----------
 
-!IM 140508
-
-      INTEGER ngrid,nlay
-      real ptimestep
-      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
-
-!   local:
-!   ------
-
+!------Entrees
+      integer, intent(in) :: ngrid,nlay
+      real, intent(in) :: ptimestep
+      real, intent(in) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
+      integer, intent(in), dimension(ngrid) ::lmax,lalim
+      real, intent(in), dimension(ngrid) :: zmax
+      real, intent(in), dimension(ngrid,nlay+1) :: zw2
+      real, intent(in), dimension(ngrid,nlay+1) :: fraca
+      real, intent(in), dimension(ngrid,nlay) :: wth3
+      real, intent(in), dimension(ngrid,nlay) :: rhobarz
+      real, intent(in), dimension(ngrid) :: wmax_sec
+      real, intent(in), dimension(ngrid,nlay) :: entr0
+      real, intent(in), dimension(ngrid,nlay+1) :: fm0,fm
+      real, intent(in), dimension(ngrid) :: pcon
+      real, intent(in), dimension(ngrid,nlay) :: alim_star
+      real, intent(in), dimension(ngrid,nlay+1,nbsrf) :: pbl_tke
+      real, intent(in), dimension(ngrid,nbsrf) :: pctsrf
+      real, intent(in), dimension(ngrid,nlay) :: omega
+      real, intent(in), dimension(ngrid) :: airephy
+!------Sorties
+      real, intent(out), dimension(ngrid) :: ale_bl,alp_bl
+      real, intent(out), dimension(ngrid,nlay) :: wght_th
+      integer, intent(out), dimension(ngrid) :: lalim_conv
+      real, intent(out), dimension(ngrid) :: zlcl,fraca0,w0,w_conv
+      real, intent(out), dimension(ngrid) :: therm_tke_max0,env_tke_max0,n2,s2,ale_bl_stat
+      real, intent(out), dimension(ngrid,nlay) :: therm_tke_max,env_tke_max
+      real, intent(out), dimension(ngrid) :: alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke
+      real, intent(out), dimension(ngrid) :: alp_bl_conv,alp_bl_stat
+
+!=============================================================================================
+!------Local
+!=============================================================================================
 
       REAL susqr2pi, reuler
-
       INTEGER ig,k,l
-      INTEGER lmax(klon),lalim(klon)
-      real zmax(klon),zw2(klon,klev+1)
-
-!on garde le zmax du pas de temps precedent
-
-
-      real fraca(klon,klev+1)
-      real wth3(klon,klev)
-! FH probleme de dimensionnement avec l'allocation dynamique
-!     common/comtherm/thetath2,wth2
-      real rhobarz(klon,klev)
-
-      real wmax_sec(klon)
-      real fm0(klon,klev+1),entr0(klon,klev)
-      real fm(klon,klev+1)
-
-!niveau de condensation
-      real pcon(klon)
-
-      real alim_star(klon,klev)
-
-!!! nrlmd le 10/04/2012
-
-!------Entrées
-      real pbl_tke(klon,klev+1,nbsrf)
-      real pctsrf(klon,nbsrf)
-      real omega(klon,klev)
-      real airephy(klon)
-!------Sorties
-      real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
-      real therm_tke_max0(klon),env_tke_max0(klon)
-      real n2(klon),s2(klon)
-      real ale_bl_stat(klon)
-      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
-      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
-!------Local
       integer nsrf
-      real rhobarz0(klon)                    ! Densité au LCL
-      logical ok_lcl(klon)                   ! Existence du LCL des thermiques
-      integer klcl(klon)                     ! Niveau du LCL
-      real interp(klon)                      ! Coef d'interpolation pour le LCL
+      real rhobarz0(ngrid)                    ! Densité au LCL
+      logical ok_lcl(ngrid)                   ! Existence du LCL des thermiques
+      integer klcl(ngrid)                     ! Niveau du LCL
+      real interp(ngrid)                      ! Coef d'interpolation pour le LCL
 !--Triggering
-      real Su                                ! Surface unité: celle d'un updraft élémentaire
-      parameter(Su=4e4)
-      real hcoef                             ! Coefficient directeur pour le calcul de s2
-      parameter(hcoef=1)
-      real hmincoef                          ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 
-      parameter(hmincoef=0.3)
-      real eps1                              ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
-      parameter(eps1=0.3)
-      real hmin(ngrid)                       ! Ordonnée à l'origine pour le calcul de s2
-      real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
-      real zmax_moy_coef
-      parameter(zmax_moy_coef=0.33)
-      real depth(klon)                       ! Epaisseur moyenne du cumulus
-      real w_max(klon)                       ! Vitesse max statistique 
-      real s_max(klon)
+      real, parameter :: su_cst=4e4              ! Surface unite: celle d'un updraft élémentaire
+      real, parameter :: hcoef=1             ! Coefficient directeur pour le calcul de s2
+      real, parameter :: hmincoef=0.3        ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 
+      real, parameter :: eps1=0.3            ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
+      real, dimension(ngrid) :: hmin         ! Ordonnée à l'origine pour le calcul de s2
+      real, dimension(ngrid) :: zmax_moy     ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
+      real, parameter :: zmax_moy_coef=0.33
+      real, dimension(ngrid) :: depth        ! Epaisseur moyenne du cumulus
+      real, dimension(ngrid) ::  w_max                 ! Vitesse max statistique 
+      real, dimension(ngrid) ::  s_max(ngrid)
 !--Closure
-      real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne 
-      real pbl_tke_max0(klon)                ! TKE moyenne au LCL
-      real w_ls(klon,klev)                   ! Vitesse verticale grande échelle (m/s)
-      real coef_m                            ! On considère un rendement pour alp_bl_fluct_m
-      parameter(coef_m=1.)
-      real coef_tke                          ! On considère un rendement pour alp_bl_fluct_tke
-      parameter(coef_tke=1.)
-
-!!! fin nrlmd le 10/04/2012
-
-!
-      !nouvelles variables pour la convection
-      real ale_bl(klon)
-      real alp_bl(klon)
-      real alp_int(klon),dp_int(klon),zdp
-      real fm_tot(klon)
-      real wght_th(klon,klev)
-      integer lalim_conv(klon)
-!v1d     logical therm
-!v1d     save therm
-
+      real, dimension(ngrid,nlay) :: pbl_tke_max       ! Profil de TKE moyenne 
+      real, dimension(ngrid) :: pbl_tke_max0           ! TKE moyenne au LCL
+      real, dimension(ngrid,nlay) :: w_ls              ! Vitesse verticale grande échelle (m/s)
+      real, parameter :: coef_m=1.            ! On considère un rendement pour alp_bl_fluct_m
+      real, parameter :: coef_tke=1.          ! On considère un rendement pour alp_bl_fluct_tke
+      real :: zdp
+      real, dimension(ngrid) :: alp_int,dp_int
+      real, dimension(ngrid) :: fm_tot
 
 !------------------------------------------------------------
 !  Initialize output arrays related to stochastic triggering
 !------------------------------------------------------------
-  DO ig = 1,klon
+  DO ig = 1,ngrid
      zlcl(ig) = 0.
      fraca0(ig) = 0.
@@ -175,17 +121,15 @@
      alp_bl_stat(ig) = 0.
   ENDDO
-  DO l = 1,klev
-    DO ig = 1,klon
+  DO l = 1,nlay
+    DO ig = 1,ngrid
      therm_tke_max(ig,l) = 0.
      env_tke_max(ig,l) = 0.
     ENDDO
   ENDDO
-!------------------------------------------------------------
-
 
 !------------Test sur le LCL des thermiques
     do ig=1,ngrid
       ok_lcl(ig)=.false.
-      if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
+      if ( (pcon(ig) .gt. pplay(ig,nlay-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
     enddo
 
@@ -207,12 +151,4 @@
     enddo
 
-!------------Hauteur des thermiques
-!!jyg le 27/04/2012
-!!    do ig =1,ngrid
-!!    rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
-!! &               -rhobarz(ig,klcl(ig)))*interp(ig)
-!!    zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
-!!      if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
-!!    enddo
     do ig =1,ngrid
 !CR:REHABILITATION ZMAX CONTINU
@@ -257,6 +193,6 @@
 
 !-----Calcul de la TKE transportée par les thermiques : therm_tke_max
-   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
-  &           rg,pplev,therm_tke_max)
+   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &  ! in
+  &           rg,pplev,therm_tke_max)                               ! out
 !   print *,' thermcell_tke_transport -> '   !!jyg
 
@@ -330,16 +266,5 @@
 !   print *,'avant Calcul de Wmax '    !!jyg
 
-!-----Calcul de Wmax et ALE_BL_STAT associée
-!!jyg le 30/04/2012
-!!   do ig=1,ngrid
-!!     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
-!!     w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/su)-log(2.*3.14)-log(2.*log(s_max(ig)/su)-log(2.*3.14))))
-!!     ale_bl_stat(ig)=0.5*w_max(ig)**2
-!!     else
-!!     w_max(ig)=0.
-!!     ale_bl_stat(ig)=0.
-!!     endif
-!!   enddo
-   susqr2pi=su*sqrt(2.*Rpi)
+   susqr2pi=su_cst*sqrt(2.*Rpi)
    reuler=exp(1.)
    do ig=1,ngrid
@@ -409,5 +334,5 @@
       lalim_conv(:)=lalim(:)
 
-      do k=1,klev
+      do k=1,nlay
          do ig=1,ngrid
             if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
@@ -417,5 +342,5 @@
 ! assez bizarre car, si on est dans la couche d'alim et que alim_star et
 ! plus petit que 1.e-10, on prend wght_th=1.
-      do k=1,klev
+      do k=1,nlay
          do ig=1,ngrid
             if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
Index: /LMDZ6/trunk/libf/phylmd/thermcell_closure.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/thermcell_closure.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/thermcell_closure.F90	(revision 4089)
@@ -17,5 +17,4 @@
       IMPLICIT NONE
 
-#include "thermcell.h"
 INTEGER ngrid,nlay
 INTEGER ig,k       
Index: /LMDZ6/trunk/libf/phylmd/thermcell_env.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/thermcell_env.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/thermcell_env.F90	(revision 4089)
@@ -1,3 +1,3 @@
-      SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
+   SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
      &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out)
 
@@ -7,36 +7,21 @@
 !--------------------------------------------------------------
 
-      USE print_control_mod, ONLY: prt_level
-      IMPLICIT NONE
 
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"      
+   USE thermcell_ini_mod, ONLY : prt_level,RLvCp,RKAPPA,RETV
+   IMPLICIT NONE
 
-      INTEGER ngrid,nlay
-      REAL po(ngrid,nlay)
-      REAL pt(ngrid,nlay)
-      REAL pu(ngrid,nlay)
-      REAL pv(ngrid,nlay)
-      REAL pplay(ngrid,nlay)
-      REAL pplev(ngrid,nlay+1)
-      integer lev_out                           ! niveau pour les print
+! arguments
 
-      REAL zo(ngrid,nlay)
-      REAL zl(ngrid,nlay)
-      REAL zh(ngrid,nlay)
-      REAL ztv(ngrid,nlay)
-      REAL zthl(ngrid,nlay)
-      REAL zpspsk(ngrid,nlay)
-      REAL zu(ngrid,nlay)
-      REAL zv(ngrid,nlay)
-      REAL pqsat(ngrid,nlay)
+   integer,intent(in) :: ngrid,nlay,lev_out
+   real,intent(in), dimension(ngrid,nlay) :: po,pt,pu,pv,pplay
+   real,intent(in), dimension(ngrid,nlay+1) :: pplev
+   real,intent(out), dimension(ngrid,nlay) :: zo,zl,zh,ztv,zthl
+   real,intent(out), dimension(ngrid,nlay) :: zpspsk,zu,zv,pqsat
+   
+! Local
 
-      INTEGER ig,ll
-
-      real dqsat_dT
-      real RLvCp
-
-logical mask(ngrid,nlay)
+   integer ig,ll
+   real dqsat_dT
+   logical mask(ngrid,nlay)
 
 
@@ -45,44 +30,39 @@
 !------------------
 
-mask(:,:)=.true.
-RLvCp = RLVTT/RCPD
+   mask(:,:)=.true.
 
 !
 ! calcul des caracteristiques de l environnement
-       DO  ll=1,nlay
-         DO ig=1,ngrid
-            zo(ig,ll)=po(ig,ll)
-            zl(ig,ll)=0.
-            zh(ig,ll)=pt(ig,ll)
-         EndDO
-       EndDO
-!
-!
+   DO  ll=1,nlay
+     DO ig=1,ngrid
+        zo(ig,ll)=po(ig,ll)
+        zl(ig,ll)=0.
+        zh(ig,ll)=pt(ig,ll)
+     enddo
+   enddo
+
 ! Condensation :
 !---------------
 ! Calcul de l'humidite a saturation et de la condensation
 
-call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
-DO ll=1,nlay
-   DO ig=1,ngrid
-      zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll))
-      zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)         !   T = Tl + Lv/Cp ql
-      zo(ig,ll) = po(ig,ll)-zl(ig,ll)
-   ENDDO
-ENDDO
-!
-!
+   call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
+   do ll=1,nlay
+      do ig=1,ngrid
+         zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll))
+         zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)         !   T = Tl + Lv/Cp ql
+         zo(ig,ll) = po(ig,ll)-zl(ig,ll)
+      enddo
+   enddo
+
 !-----------------------------------------------------------------------
+   if (prt_level.ge.1) print*,'0 OK convect8'
 
-      if (prt_level.ge.1) print*,'0 OK convect8'
-
-      DO ll=1,nlay
-         DO ig=1,ngrid
-             zpspsk(ig,ll)=(pplay(ig,ll)/100000.)**RKAPPA
-             zu(ig,ll)=pu(ig,ll)
-             zv(ig,ll)=pv(ig,ll)
+   do ll=1,nlay
+      do ig=1,ngrid
+          zpspsk(ig,ll)=(pplay(ig,ll)/100000.)**RKAPPA
+          zu(ig,ll)=pu(ig,ll)
+          zv(ig,ll)=pv(ig,ll)
 !attention zh est maintenant le profil de T et plus le profil de theta !
 ! Quelle horreur ! A eviter.
-!
 !   T-> Theta
             ztv(ig,ll)=zh(ig,ll)/zpspsk(ig,ll)
@@ -92,7 +72,7 @@
             zthl(ig,ll)=pt(ig,ll)/zpspsk(ig,ll)
 !            
-         ENDDO
-      ENDDO
+      enddo
+   enddo
  
-      RETURN
-      END
+   RETURN
+   END
Index: /LMDZ6/trunk/libf/phylmd/thermcell_flux2.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/thermcell_flux2.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/thermcell_flux2.F90	(revision 4089)
@@ -13,7 +13,6 @@
 !---------------------------------------------------------------------------
 
-      USE print_control_mod, ONLY: prt_level
+      USE thermcell_ini_mod, ONLY : prt_level,iflag_thermals_optflux
       IMPLICIT NONE
-#include "thermcell.h"
       
       INTEGER ig,l
Index: /LMDZ6/trunk/libf/phylmd/thermcell_height.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/thermcell_height.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/thermcell_height.F90	(revision 4089)
@@ -6,5 +6,4 @@
 !-----------------------------------------------------------------------------
       IMPLICIT NONE
-#include "thermcell.h"
 
       INTEGER ig,l
Index: /LMDZ6/trunk/libf/phylmd/thermcell_main.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/thermcell_main.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/thermcell_main.F90	(revision 4089)
@@ -1,6 +1,6 @@
-!
+
 ! $Id$
 !
-      SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep  &
+      subroutine thermcell_main(itap,ngrid,nlay,ptimestep  &
      &                  ,pplay,pplev,pphi,debut  &
      &                  ,pu,pv,pt,po  &
@@ -8,21 +8,27 @@
      &                  ,fm0,entr0,detr0,zqta,zqla,lmax  &
      &                  ,ratqscth,ratqsdiff,zqsatth  &
-     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
      &                  ,zmax0, f0,zw2,fraca,ztv &
-     &                  ,zpspsk,ztla,zthl &
-!!! nrlmd le 10/04/2012
-     &                  ,pbl_tke,pctsrf,omega,airephy &
-     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
-     &                  ,n2,s2,ale_bl_stat &
-     &                  ,therm_tke_max,env_tke_max &
-     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
-     &                  ,alp_bl_conv,alp_bl_stat &
-!!! fin nrlmd le 10/04/2012
-     &                  ,ztva  )
-
-      USE dimphy
-      USE ioipsl
-      USE indice_sol_mod
-      USE print_control_mod, ONLY: lunout,prt_level
+     &                  ,zpspsk,ztla,zthl,ztva &
+     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
+#ifdef ISO         
+     &      ,xtpo,xtpdoadj &
+#endif         
+     &   )
+
+
+      USE thermcell_ini_mod, ONLY: thermcell_ini,dqimpl,dvdq,prt_level,lunout,prt_level
+      USE thermcell_ini_mod, ONLY: iflag_thermals_closure,iflag_thermals_ed,tau_thermals,r_aspect_thermals
+      USE thermcell_ini_mod, ONLY: RD,RG
+
+#ifdef ISO
+  USE infotrac_phy, ONLY : ntraciso
+#ifdef ISOVERIF
+  USE isotopes_mod, ONLY : iso_eau,iso_HDO
+  USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
+        iso_verif_aberrant_encadre
+#endif
+#endif
+
+
       IMPLICIT NONE
 
@@ -62,35 +68,29 @@
 !   -------------
 
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"
-#include "thermcell.h"
 
 !   arguments:
 !   ----------
-
-!IM 140508
-      INTEGER itap
-
-      INTEGER ngrid,nlay
-      real ptimestep
-      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)
-      LOGICAL debut
+      integer, intent(in) :: itap,ngrid,nlay
+      real, intent(in) ::  ptimestep
+      real, intent(in), dimension(ngrid,nlay)    :: pt,pu,pv,po,pplay,pphi,zpspsk
+      real, intent(in), dimension(ngrid,nlay+1)  :: pplev
+      real, intent(out), dimension(ngrid,nlay)   :: pdtadj,pduadj,pdvadj,pdoadj,entr0,detr0
+      real, intent(out), dimension(ngrid,nlay)   :: ztla,zqla,zqta,zqsatth,zthl
+      real, intent(out), dimension(ngrid,nlay+1) :: fm0,zw2,fraca
+      real, intent(out), dimension(ngrid) :: zmax0,f0
+      real, intent(out), dimension(ngrid,nlay) :: ztva,ztv
+      logical, intent(in) :: debut
+
+      real, intent(out), dimension(ngrid) :: pcon
+      real, intent(out), dimension(ngrid,nlay) :: rhobarz,wth3
+      real, intent(out), dimension(ngrid) :: wmax_sec
+      integer,intent(out), dimension(ngrid) :: lalim
+      real, intent(out), dimension(ngrid,nlay+1) :: fm
+      real, intent(out), dimension(ngrid,nlay) :: alim_star
+      real, intent(out), dimension(ngrid) :: zmax
 
 !   local:
 !   ------
 
-      integer icount
-
-      integer, save :: dvdq=1,dqimpl=-1
-!$OMP THREADPRIVATE(dvdq,dqimpl)
-      data icount/0/
-      save icount
-!$OMP THREADPRIVATE(icount)
 
       integer,save :: igout=1
@@ -101,150 +101,33 @@
 !$OMP THREADPRIVATE(lev_out)
 
-      REAL susqr2pi, Reuler
-
-      INTEGER ig,k,l,ll,ierr
-      real zsortie1d(klon)
-      INTEGER lmax(klon),lmin(klon),lalim(klon)
-      INTEGER lmix(klon)
-      INTEGER lmix_bis(klon)
-      real linter(klon)
-      real zmix(klon)
-      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
-!      real fraca(klon,klev)
-
-      real zmax_sec(klon)
-!on garde le zmax du pas de temps precedent
-      real zmax0(klon)
-!FH/IM     save zmax0
-
-      real lambda
-
-      real zlev(klon,klev+1),zlay(klon,klev)
-      real deltaz(klon,klev)
-      REAL zh(klon,klev)
-      real zthl(klon,klev),zdthladj(klon,klev)
-      REAL ztv(klon,klev)
-      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
-      real zl(klon,klev)
-      real zsortie(klon,klev)
-      real zva(klon,klev)
-      real zua(klon,klev)
-      real zoa(klon,klev)
-
-      real zta(klon,klev)
-      real zha(klon,klev)
-      real fraca(klon,klev+1)
-      real zf,zf2
-      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
-      real q2(klon,klev)
-! FH probleme de dimensionnement avec l'allocation dynamique
-!     common/comtherm/thetath2,wth2
-      real wq(klon,klev)
-      real wthl(klon,klev)
-      real wthv(klon,klev)
-    
-      real ratqscth(klon,klev)
-      real var
-      real vardiff
-      real ratqsdiff(klon,klev)
-
+      real lambda, zf,zf2,var,vardiff,CHI
+      integer ig,k,l,ierr,ll
       logical sorties
-      real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
-      real zpspsk(klon,klev)
-
-      real wmax(klon)
-      real wmax_tmp(klon)
-      real wmax_sec(klon)
-      real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
-      real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
-
-      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
-!niveau de condensation
-      integer nivcon(klon)
-      real zcon(klon)
-      REAL CHI
-      real zcon2(klon)
-      real pcon(klon)
-      real zqsat(klon,klev)
-      real zqsatth(klon,klev) 
-
-      real f_star(klon,klev+1),entr_star(klon,klev)
-      real detr_star(klon,klev)
-      real alim_star_tot(klon)
-      real alim_star(klon,klev)
-      real alim_star_clos(klon,klev)
-      real f(klon), f0(klon)
-!FH/IM     save f0
-      real zlevinter(klon)
-       real seuil
-      real csc(klon,klev)
-
-!!! nrlmd le 10/04/2012
-
-!------Entrées
-      real pbl_tke(klon,klev+1,nbsrf)
-      real pctsrf(klon,nbsrf)
-      real omega(klon,klev)
-      real airephy(klon)
-!------Sorties
-      real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
-      real therm_tke_max0(klon),env_tke_max0(klon)
-      real n2(klon),s2(klon)
-      real ale_bl_stat(klon)
-      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
-      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
-!------Local
-      integer nsrf
-      real rhobarz0(klon)                    ! Densité au LCL
-      logical ok_lcl(klon)                   ! Existence du LCL des thermiques
-      integer klcl(klon)                     ! Niveau du LCL
-      real interp(klon)                      ! Coef d'interpolation pour le LCL
-!--Triggering
-      real Su                                ! Surface unité: celle d'un updraft élémentaire
-      parameter(Su=4e4)
-      real hcoef                             ! Coefficient directeur pour le calcul de s2
-      parameter(hcoef=1)
-      real hmincoef                          ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 
-      parameter(hmincoef=0.3)
-      real eps1                              ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
-      parameter(eps1=0.3)
-      real hmin(ngrid)                       ! Ordonnée à l'origine pour le calcul de s2
-      real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
-      real zmax_moy_coef
-      parameter(zmax_moy_coef=0.33)
-      real depth(klon)                       ! Epaisseur moyenne du cumulus
-      real w_max(klon)                       ! Vitesse max statistique 
-      real s_max(klon)
-!--Closure
-      real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne 
-      real pbl_tke_max0(klon)                ! TKE moyenne au LCL
-      real w_ls(klon,klev)                   ! Vitesse verticale grande échelle (m/s)
-      real coef_m                            ! On considère un rendement pour alp_bl_fluct_m
-      parameter(coef_m=1.)
-      real coef_tke                          ! On considère un rendement pour alp_bl_fluct_tke
-      parameter(coef_tke=1.)
-
-!!! fin nrlmd le 10/04/2012
-
-!
-      !nouvelles variables pour la convection
-      real Ale_bl(klon)
-      real Alp_bl(klon)
-      real alp_int(klon),dp_int(klon),zdp
-      real ale_int(klon)
-      integer n_int(klon)
-      real fm_tot(klon)
-      real wght_th(klon,klev)
-      integer lalim_conv(klon)
-!v1d     logical therm
-!v1d     save therm
-
-      character*2 str2
-      character*10 str10
+      real, dimension(ngrid) :: linter,zmix, zmax_sec
+      integer,dimension(ngrid) :: lmax,lmin,lmix,lmix_bis,nivcon
+      real, dimension(ngrid,nlay) :: ztva_est
+      real, dimension(ngrid,nlay) :: deltaz,zlay,zh,zdthladj,zu,zv,zo,zl,zva,zua,zoa
+      real, dimension(ngrid,nlay) :: zta,zha,q2,wq,wthl,wthv,thetath2,wth2
+      real, dimension(ngrid,nlay) :: ratqscth,ratqsdiff,rho,masse
+      real, dimension(ngrid,nlay+1) :: zw_est,zlev
+      real, dimension(ngrid) :: wmax,wmax_tmp
+      real, dimension(ngrid,nlay+1) :: f_star
+      real, dimension(ngrid,nlay) :: entr,detr,entr_star,detr_star,alim_star_clos
+      real, dimension(ngrid,nlay) :: zqsat,csc
+      real, dimension(ngrid) :: zcon,zcon2,alim_star_tot,f
 
       character (len=20) :: modname='thermcell_main'
       character (len=80) :: abort_message
 
-      EXTERNAL SCOPY
+
+#ifdef ISO
+      REAL xtpo(ntraciso,ngrid,nlay),xtpdoadj(ntraciso,ngrid,nlay)
+      REAL xtzo(ntraciso,ngrid,nlay)
+      REAL xtpdoadj_tmp(ngrid,nlay)
+      REAL xtpo_tmp(ngrid,nlay)
+      REAL xtzo_tmp(ngrid,nlay)
+      integer ixt
+#endif
+
 !
 
@@ -253,43 +136,21 @@
 !   ---------------
 !
-
-   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
+    print*,'NEW THERMCELL cool'
+
+
    fm=0. ; entr=0. ; detr=0.
-   icount=icount+1
-
-!IM 090508 beg
-!print*,'====================================================================='
-!print*,'====================================================================='
-!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
-!print*,'====================================================================='
-!print*,'====================================================================='
-!IM 090508 end
 
       if (prt_level.ge.1) print*,'thermcell_main V4'
 
        sorties=.true.
-      IF(ngrid.NE.klon) THEN
+      IF(ngrid.NE.ngrid) THEN
          PRINT*
          PRINT*,'STOP dans convadj'
          PRINT*,'ngrid    =',ngrid
-         PRINT*,'klon  =',klon
+         PRINT*,'ngrid  =',ngrid
       ENDIF
 !
 !     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
-     do ig=1,klon
+     do ig=1,ngrid
          f0(ig)=max(f0(ig),1.e-2)
          zmax0(ig)=max(zmax0(ig),40.)
@@ -336,34 +197,26 @@
          zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
       enddo
-         zlev(:,1)=0.
-         zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
+      zlev(:,1)=0.
+      zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
       do l=1,nlay
          zlay(:,l)=pphi(:,l)/RG
       enddo
-!calcul de l epaisseur des couches
       do l=1,nlay
          deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
       enddo
 
-!     print*,'2 OK convect8'
 !-----------------------------------------------------------------------
-!   Calcul des densites
+!   Calcul des densites et masses
 !-----------------------------------------------------------------------
 
-     rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))
-
-     if (prt_level.ge.10)write(lunout,*)                                &
-    &    'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
+      rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))
+      if (prt_level.ge.10) write(lunout,*) 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
       rhobarz(:,1)=rho(:,1)
-
       do l=2,nlay
          rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
       enddo
-
-!calcul de la masse
       do l=1,nlay
          masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
       enddo
-
       if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
 
@@ -480,6 +333,6 @@
       if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
 
-      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
+      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
+      call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
 
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
@@ -509,8 +362,8 @@
 
 
-      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
+      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
+      call test_ltherm(ngrid,nlay,pplay,lmin ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
+      call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
+      call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
 
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
@@ -526,6 +379,6 @@
 
  
-call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
-call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
+call test_ltherm(ngrid,nlay,pplay,lmin,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
+call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
 
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
@@ -582,5 +435,4 @@
 !-------------------------------------------------------------------------------
 !deduction des flux
-!-------------------------------------------------------------------------------
 
       CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
@@ -591,6 +443,6 @@
 
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
-      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
+      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
+      call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
 
 !------------------------------------------------------------------
@@ -620,13 +472,54 @@
      &                   po,pdoadj,zoa,lev_out)
 
+#ifdef ISO
+        ! C Risi: on utilise directement la même routine
+        do ixt=1,ntraciso
+          do ll=1,nlay
+            DO ig=1,ngrid
+                xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll)
+                xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll)
+            enddo
+          enddo
+          call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
+     &                   xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out)
+          do ll=1,nlay
+            DO ig=1,ngrid
+                xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll)
+            enddo
+          enddo
+        enddo !do ixt=1,ntraciso
+#endif
+
+#ifdef ISO      
+#ifdef ISOVERIF
+      DO  ll=1,nlay
+        DO ig=1,ngrid
+          if (iso_eau.gt.0) then
+              call iso_verif_egalite(xtpo(iso_eau,ig,ll), &
+     &          po(ig,ll),'thermcell_main 594')
+              call iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), &
+     &          pdoadj(ig,ll),'thermcell_main 596')
+          endif
+          if (iso_HDO.gt.0) then
+              call iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) &
+     &           /po(ig,ll),'thermcell_main 610')
+          endif
+        enddo
+      enddo !DO  ll=1,nlay 
+      write(*,*) 'thermcell_main 600 tmp: apres thermcell_dq'
+#endif      
+#endif
+
+
+
 !------------------------------------------------------------------
 ! Calcul de la fraction de l'ascendance
 !------------------------------------------------------------------
-      do ig=1,klon
+      do ig=1,ngrid
          fraca(ig,1)=0.
          fraca(ig,nlay+1)=0.
       enddo
       do l=2,nlay
-         do ig=1,klon
+         do ig=1,ngrid
             if (zw2(ig,l).gt.1.e-10) then
             fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
@@ -760,25 +653,4 @@
          enddo
       enddo
-!
-! $Id$
-!
-      CALL thermcell_alp(ngrid,nlay,ptimestep  &
-     &                  ,pplay,pplev  &
-     &                  ,fm0,entr0,lmax  &
-     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
-     &                  ,zw2,fraca &
-!!! necessire en plus
-     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
-!!! nrlmd le 10/04/2012
-     &                  ,pbl_tke,pctsrf,omega,airephy &
-     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
-     &                  ,n2,s2,ale_bl_stat &
-     &                  ,therm_tke_max,env_tke_max &
-     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
-     &                  ,alp_bl_conv,alp_bl_stat &
-!!! fin nrlmd le 10/04/2012
-     &                   )
-
-
 
 !calcul du ratqscdiff
@@ -788,5 +660,5 @@
       ratqsdiff(:,:)=0.
 
-      do l=1,klev
+      do l=1,nlay
          do ig=1,ngrid
             if (l<=lalim(ig)) then
@@ -798,5 +670,5 @@
       if (prt_level.ge.1) print*,'14f OK convect8'
 
-      do l=1,klev
+      do l=1,nlay
          do ig=1,ngrid
             if (l<=lalim(ig)) then
@@ -809,15 +681,9 @@
 
       if (prt_level.ge.1) print*,'14g OK convect8'
-      do l=1,nlay
-         do ig=1,ngrid
-            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
-!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
-         enddo
-      enddo 
-!--------------------------------------------------------------------    
-!
-!ecriture des fichiers sortie
-!     print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
-
+         do l=1,nlay
+            do ig=1,ngrid
+               ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
+            enddo
+         enddo 
       endif
 
@@ -825,23 +691,22 @@
 
       return
-      end
-
-!-----------------------------------------------------------------------------
-
-      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
-      USE print_control_mod, ONLY: prt_level
+      end subroutine thermcell_main
+
+!=============================================================================
+!/////////////////////////////////////////////////////////////////////////////
+!=============================================================================
+      subroutine test_ltherm(ngrid,nlay,pplay,long,ztv,po,ztva, &  ! in
+    &            zqla,f_star,zw2,comment)                          ! in
+!=============================================================================
+      USE thermcell_ini_mod, ONLY: prt_level
       IMPLICIT NONE
 
-      integer i, k, klon,klev
-      real pplev(klon,klev+1),pplay(klon,klev)
-      real ztv(klon,klev)
-      real po(klon,klev)
-      real ztva(klon,klev)
-      real zqla(klon,klev)
-      real f_star(klon,klev)
-      real zw2(klon,klev)
-      integer long(klon)
+      integer i, k, ngrid,nlay
+      real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,po,ztva,zqla
+      real, intent(in), dimension(ngrid,nlay) :: f_star,zw2
+      integer, intent(in), dimension(ngrid) :: long
       real seuil
       character*21 comment
+      seuil=0.25
 
       if (prt_level.ge.1) THEN
@@ -851,10 +716,10 @@
 
 !  test sur la hauteur des thermiques ...
-         do i=1,klon
+         do i=1,ngrid
 !IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
            if (prt_level.ge.10) then
                print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
                print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
-               do k=1,klev
+               do k=1,nlay
                   write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
                enddo
@@ -866,10 +731,16 @@
       end
 
-!!! nrlmd le 10/04/2012                          Transport de la TKE par le thermique moyen pour la fermeture en ALP 
-!                                                         On transporte pbl_tke pour donner therm_tke
-!                                          Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
-      subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
-     &           rg,pplev,therm_tke_max)
-      USE print_control_mod, ONLY: prt_level
+! nrlmd le 10/04/2012   Transport de la TKE par le thermique moyen pour la fermeture en ALP 
+!                       On transporte pbl_tke pour donner therm_tke
+!                       Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
+
+!=======================================================================
+!///////////////////////////////////////////////////////////////////////
+!=======================================================================
+
+      subroutine thermcell_tke_transport( &
+     &     ngrid,nlay,ptimestep,fm0,entr0,rg,pplev,  &   ! in
+     &     therm_tke_max)                                ! out
+      USE thermcell_ini_mod, ONLY: prt_level
       implicit none
 
@@ -882,13 +753,14 @@
 !=======================================================================
 
-      integer ngrid,nlay,nsrf
-
-      real ptimestep
-      real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
-      real entr0(ngrid,nlay),rg
-      real therm_tke_max(ngrid,nlay)
+      integer ngrid,nlay
+
+      real, intent(in) :: ptimestep
+      real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev
+      real, intent(in), dimension(ngrid,nlay) :: entr0
+      real, intent(in) :: rg
+      real, intent(out), dimension(ngrid,nlay) :: therm_tke_max
+
       real detr0(ngrid,nlay)
-
-
+      real masse0(ngrid,nlay)
       real masse(ngrid,nlay),fm(ngrid,nlay+1)
       real entr(ngrid,nlay)
@@ -897,9 +769,5 @@
 
       real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
-
-      real zzm
-
       integer ig,k
-      integer isrf
 
 
@@ -929,14 +797,5 @@
       fm(:,nlay+1)=0.
 
-!!! nrlmd le 16/09/2010
-!   calcul de la valeur dans les ascendances
-!       do ig=1,ngrid
-!          qa(ig,1)=q(ig,1)
-!       enddo
-!!! 
-
-!do isrf=1,nsrf
-
-!   q(:,:)=therm_tke(:,:,isrf)
+
    q(:,:)=therm_tke_max(:,:)
 !!! nrlmd le 16/09/2010
Index: /LMDZ6/trunk/libf/phylmd/thermcell_plume.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/thermcell_plume.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/thermcell_plume.F90	(revision 4089)
@@ -21,14 +21,11 @@
 !   = 29 : an other way to compute the modified buoyancy (to be tested)
 !--------------------------------------------------------------------------
-USE IOIPSL, ONLY : getin
-USE ioipsl_getin_p_mod, ONLY : getin_p
-
-       USE print_control_mod, ONLY: prt_level
+       USE thermcell_ini_mod, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
+       USE thermcell_ini_mod, ONLY: fact_epsilon, betalpha, afact, fact_shell
+       USE thermcell_ini_mod, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
+       USE thermcell_ini_mod, ONLY: mix0, thermals_flag_alim
+
+
        IMPLICIT NONE
-
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"
-#include "thermcell.h"
 
       INTEGER itap
@@ -101,28 +98,8 @@
       real zbetalpha, coefzlmel
       real eps
-      REAL REPS,RLvCp,DDT0
-      PARAMETER (DDT0=.01)
       logical Zsat
       LOGICAL active(ngrid),activetmp(ngrid)
       REAL fact_gamma,fact_gamma2,fact_epsilon2
 
-      REAL, SAVE :: fact_epsilon=0.002
-      REAL, SAVE :: betalpha=0.9
-      REAL, SAVE :: afact=2./3.
-      REAL, SAVE :: fact_shell=1.
-      REAL,SAVE :: detr_min=1.e-5
-      REAL,SAVE :: entr_min=1.e-5
-      REAL,SAVE :: detr_q_coef=0.012
-      REAL,SAVE :: detr_q_power=0.5
-      REAL,SAVE :: mix0=0.
-      INTEGER,SAVE :: thermals_flag_alim=0
-
-!$OMP THREADPRIVATE(fact_epsilon, betalpha, afact, fact_shell)
-!$OMP THREADPRIVATE(detr_min, entr_min, detr_q_coef, detr_q_power)
-!$OMP THREADPRIVATE( mix0, thermals_flag_alim)
-
-      LOGICAL, SAVE :: first=.true.
-  !$OMP THREADPRIVATE(first)
-
 
       REAL c2(ngrid,klev)
@@ -132,21 +109,4 @@
 ! Initialisation
 
-      RLvCp = RLVTT/RCPD
-      IF (first) THEN
-
-     CALL getin_p('thermals_fact_epsilon',fact_epsilon)
-     CALL getin_p('thermals_betalpha',betalpha)
-     CALL getin_p('thermals_afact',afact)
-     CALL getin_p('thermals_fact_shell',fact_shell)
-     CALL getin_p('thermals_detr_min',detr_min)
-     CALL getin_p('thermals_entr_min',entr_min)
-     CALL getin_p('thermals_detr_q_coef',detr_q_coef)
-     CALL getin_p('thermals_detr_q_power',detr_q_power)
-     CALL getin_p('thermals_mix0',mix0)
-     CALL getin_p('thermals_flag_alim',thermals_flag_alim)
-
-
-      first=.false.
-      ENDIF
 
       zbetalpha=betalpha/(1.+betalpha)
Index: /LMDZ6/trunk/libf/phylmd/thermcell_plume_6A.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/thermcell_plume_6A.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/thermcell_plume_6A.F90	(revision 4089)
@@ -11,61 +11,45 @@
 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
 !--------------------------------------------------------------------------
-USE IOIPSL, ONLY : getin
-USE ioipsl_getin_p_mod, ONLY : getin_p
-
-       USE print_control_mod, ONLY: prt_level
+
+       USE thermcell_ini_mod, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
+       USE thermcell_ini_mod, ONLY: fact_epsilon, betalpha, afact, fact_shell
+       USE thermcell_ini_mod, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
+       USE thermcell_ini_mod, ONLY: mix0, thermals_flag_alim
+
        IMPLICIT NONE
 
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"
-#include "thermcell.h"
-
-      INTEGER itap
-      INTEGER lunout1,igout
-      INTEGER ngrid,klev
-      REAL ptimestep
-      REAL ztv(ngrid,klev)
-      REAL zthl(ngrid,klev)
-      REAL po(ngrid,klev)
-      REAL zl(ngrid,klev)
-      REAL rhobarz(ngrid,klev)
-      REAL zlev(ngrid,klev+1)
-      REAL pplev(ngrid,klev+1)
-      REAL pphi(ngrid,klev)
-      REAL zpspsk(ngrid,klev)
-      REAL alim_star(ngrid,klev)
-      REAL f0(ngrid)
-      INTEGER lalim(ngrid)
-      integer lev_out                           ! niveau pour les print
-      integer nbpb
-    
-      real alim_star_tot(ngrid)
-
-      REAL ztva(ngrid,klev)
-      REAL ztla(ngrid,klev)
-      REAL zqla(ngrid,klev)
-      REAL zqta(ngrid,klev)
-      REAL zha(ngrid,klev)
-
-      REAL detr_star(ngrid,klev)
-      REAL coefc
-      REAL entr_star(ngrid,klev)
-      REAL detr(ngrid,klev)
-      REAL entr(ngrid,klev)
-
-      REAL csc(ngrid,klev)
-
-      REAL zw2(ngrid,klev+1)
-      REAL w_est(ngrid,klev+1)
-      REAL f_star(ngrid,klev+1)
-      REAL wa_moy(ngrid,klev+1)
-
-      REAL ztva_est(ngrid,klev)
-      REAL ztv_est(ngrid,klev)
-      REAL zqla_est(ngrid,klev)
-      REAL zqsatth(ngrid,klev)
-      REAL zta_est(ngrid,klev)
-      REAL ztemp(ngrid),zqsat(ngrid)
+      integer,intent(in) :: itap,lev_out,lunout1,igout,ngrid,klev
+      real,intent(in) :: ptimestep
+      real,intent(in),dimension(ngrid,klev) :: ztv
+      real,intent(in),dimension(ngrid,klev) :: zthl
+      real,intent(in),dimension(ngrid,klev) :: po
+      real,intent(in),dimension(ngrid,klev) :: zl
+      real,intent(in),dimension(ngrid,klev) :: rhobarz
+      real,intent(in),dimension(ngrid,klev+1) :: zlev
+      real,intent(in),dimension(ngrid,klev+1) :: pplev
+      real,intent(in),dimension(ngrid,klev) :: pphi
+      real,intent(in),dimension(ngrid,klev) :: zpspsk
+      real,intent(in),dimension(ngrid) :: f0
+
+      integer,intent(out) :: lalim(ngrid)
+      real,intent(out),dimension(ngrid,klev) :: alim_star
+      real,intent(out),dimension(ngrid) :: alim_star_tot
+      real,intent(out),dimension(ngrid,klev) :: detr_star
+      real,intent(out),dimension(ngrid,klev) :: entr_star
+      real,intent(out),dimension(ngrid,klev+1) :: f_star
+      real,intent(out),dimension(ngrid,klev) :: csc
+      real,intent(out),dimension(ngrid,klev) :: ztva
+      real,intent(out),dimension(ngrid,klev) :: ztla
+      real,intent(out),dimension(ngrid,klev) :: zqla
+      real,intent(out),dimension(ngrid,klev) :: zqta
+      real,intent(out),dimension(ngrid,klev) :: zha
+      real,intent(out),dimension(ngrid,klev+1) :: zw2
+      real,intent(out),dimension(ngrid,klev+1) :: w_est
+      real,intent(out),dimension(ngrid,klev) :: ztva_est
+      real,intent(out),dimension(ngrid,klev) :: zqsatth
+      integer,intent(out),dimension(ngrid) :: lmix(ngrid)
+      integer,intent(out),dimension(ngrid) :: lmix_bis(ngrid)
+      real,intent(out),dimension(ngrid) :: linter(ngrid)
+
       REAL zdw2,zdw2bis
       REAL zw2modif
@@ -73,10 +57,16 @@
       REAL zeps(ngrid,klev)
 
-      REAL linter(ngrid)
-      INTEGER lmix(ngrid)
-      INTEGER lmix_bis(ngrid)
       REAL    wmaxa(ngrid)
 
       INTEGER ig,l,k,lt,it,lm
+      integer nbpb
+
+      real,dimension(ngrid,klev) :: detr
+      real,dimension(ngrid,klev) :: entr
+      real,dimension(ngrid,klev+1) :: wa_moy
+      real,dimension(ngrid,klev) :: ztv_est
+      real,dimension(ngrid) :: ztemp,zqsat
+      real,dimension(ngrid,klev) :: zqla_est
+      real,dimension(ngrid,klev) :: zta_est
 
       real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
@@ -91,29 +81,8 @@
       real zbetalpha, coefzlmel
       real eps
-      REAL REPS,RLvCp,DDT0
-      PARAMETER (DDT0=.01)
       logical Zsat
       LOGICAL active(ngrid),activetmp(ngrid)
       REAL fact_gamma,fact_gamma2,fact_epsilon2
-
-      REAL, SAVE :: fact_epsilon=0.002
-      REAL, SAVE :: betalpha=0.9
-      REAL, SAVE :: afact=2./3.
-      REAL, SAVE :: fact_shell=1.
-      REAL,SAVE :: detr_min=1.e-5
-      REAL,SAVE :: entr_min=1.e-5
-      REAL,SAVE :: detr_q_coef=0.012
-      REAL,SAVE :: detr_q_power=0.5
-      REAL,SAVE :: mix0=0.
-      INTEGER,SAVE :: thermals_flag_alim=0
-
-!$OMP THREADPRIVATE(fact_epsilon, betalpha, afact, fact_shell)
-!$OMP THREADPRIVATE(detr_min, entr_min, detr_q_coef, detr_q_power)
-!$OMP THREADPRIVATE( mix0, thermals_flag_alim)
-
-      LOGICAL, SAVE :: first=.true.
-  !$OMP THREADPRIVATE(first)
-
-
+      REAL coefc
       REAL c2(ngrid,klev)
 
@@ -122,21 +91,4 @@
 ! Initialisation
 
-      RLvCp = RLVTT/RCPD
-      IF (first) THEN
-
-     CALL getin_p('thermals_fact_epsilon',fact_epsilon)
-     CALL getin_p('thermals_betalpha',betalpha)
-     CALL getin_p('thermals_afact',afact)
-     CALL getin_p('thermals_fact_shell',fact_shell)
-     CALL getin_p('thermals_detr_min',detr_min)
-     CALL getin_p('thermals_entr_min',entr_min)
-     CALL getin_p('thermals_detr_q_coef',detr_q_coef)
-     CALL getin_p('thermals_detr_q_power',detr_q_power)
-     CALL getin_p('thermals_mix0',mix0)
-     CALL getin_p('thermals_flag_alim',thermals_flag_alim)
-
-
-      first=.false.
-      ENDIF
 
       zbetalpha=betalpha/(1.+betalpha)
@@ -786,11 +738,6 @@
 !--------------------------------------------------------------------------
 
-      USE print_control_mod, ONLY: prt_level
+      USE thermcell_ini_mod, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
       IMPLICIT NONE
-
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"
-#include "thermcell.h"
 
       INTEGER itap
@@ -857,6 +804,4 @@
       real betalpha,zbetalpha
       real eps, afact
-      REAL REPS,RLvCp,DDT0
-      PARAMETER (DDT0=.01)
       logical Zsat
       LOGICAL active(ngrid),activetmp(ngrid)
@@ -866,5 +811,4 @@
 ! Initialisation
 
-      RLvCp = RLVTT/RCPD
       fact_epsilon=0.002
       betalpha=0.9 
@@ -923,5 +867,5 @@
 
 !-------------------------------------------------------------------------
-! Definition de l'alimentation a l'origine dans thermcell_init
+! Definition de l'alimentation
 !-------------------------------------------------------------------------
       do l=1,klev-1
Index: /LMDZ6/trunk/libf/phylmd/wake.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/wake.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmd/wake.F90	(revision 4089)
@@ -280,4 +280,7 @@
   ! Initialisations
   ! -------------------------------------------------------------------------
+  ! ALON = 3.e5
+  ! alon = 1.E6
+
   !  Provisionnal; to be suppressed when f_shear is parameterized
   f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
Index: /LMDZ6/trunk/libf/phylmdiso/calltherm.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/calltherm.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmdiso/calltherm.F90	(revision 4089)
@@ -7,5 +7,5 @@
      &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs  &
      &      ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,  &
-     &       ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, &
+     &       ratqsdiff,zqsatth,ale_bl,alp_bl,lalim_conv,wght_th, &
      &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl &
 !!! nrlmd le 10/04/2012
@@ -39,5 +39,6 @@
 
       implicit none
-      include "thermcell.h"
+      include "clesphys.h"
+      include "thermcell_old.h"
 
 
@@ -94,8 +95,8 @@
       real zqsatth(klon,klev)  
 !nouvelles variables pour la convection
-      real Ale_bl(klon)
-      real Alp_bl(klon)
-      real Ale(klon)
-      real Alp(klon)
+      real ale_bl(klon)
+      real alp_bl(klon)
+      real ale(klon)
+      real alp(klon)
 !RC
       !on garde le zmax du pas de temps precedent
@@ -117,4 +118,13 @@
 !********************************************************
 
+      real, dimension(klon) :: pcon
+      real, dimension(klon,klev) :: rhobarz,wth3
+      integer,dimension(klon) :: lalim
+      real, dimension(klon,klev+1) :: fm
+      real, dimension(klon,klev) :: alim_star
+      real, dimension(klon) :: zmax
+
+
+
 
 ! variables locales
@@ -130,6 +140,7 @@
       character (len=80) :: abort_message
 
-      integer i,k
+      integer i,k,isplit
       logical, save :: first=.true.
+      logical :: new_thermcell
 
 #ifdef ISO
@@ -173,6 +184,6 @@
          detr_therm(:,:)=0.
 
-         Ale_bl(:)=0.
-         Alp_bl(:)=0.
+         ale_bl(:)=0.
+         alp_bl(:)=0.
          if (prt_level.ge.10) then
           print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
@@ -207,5 +218,10 @@
          if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
 
-#ifdef ISO
+
+         new_thermcell=iflag_thermals>=15.and.iflag_thermals<=18
+#ifdef ISO
+      if (.not.new_thermcell) then
+           CALL abort_gcm('calltherm 234','isos pas prevus ici',1)
+      endif
 #ifdef ISOVERIF
       if (iso_eau.gt.0) then
@@ -217,10 +233,9 @@
 #endif
          zdt=dtime/REAL(nsplit_thermals)
+
+
          do isplit=1,nsplit_thermals
 
           if (iflag_thermals>=1000) then
-#ifdef ISO
-              CALL abort_gcm('calltherm 173','isos pas prevus ici',1)
-#endif
             CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
      &      ,pplay,paprs,pphi  &
@@ -231,7 +246,4 @@
      &      ,tau_thermals)
           else if (iflag_thermals.eq.2) then
-#ifdef ISO
-              CALL abort_gcm('calltherm 186','isos pas prevus ici',1)  
-#endif
             CALL thermcell_sec(klon,klev,zdt  &
      &      ,pplay,paprs,pphi,zlev  &
@@ -242,8 +254,4 @@
      &      ,tau_thermals)
           else if (iflag_thermals.eq.3) then
-#ifdef ISO
-              write(*,*) 'calltherm 199: isos pas prévus ici'
-              stop
-#endif
             CALL thermcell(klon,klev,zdt  &
      &      ,pplay,paprs,pphi  &
@@ -254,7 +262,4 @@
      &      ,tau_thermals)
           else if (iflag_thermals.eq.10) then
-#ifdef ISO
-              CALL abort_gcm('calltherm 212','isos pas prevus ici',1)  
-#endif
             CALL thermcell_eau(klon,klev,zdt  &
      &      ,pplay,paprs,pphi  &
@@ -264,22 +269,8 @@
      &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
      &      ,tau_thermals)
-#ifdef ISO
-              CALL abort_gcm('calltherm 267','isos pas prevus ici',1)  
-#endif
           else if (iflag_thermals.eq.11) then
               abort_message = 'cas non prevu dans calltherm'
               CALL abort_physic (modname,abort_message,1)
-
-!           CALL thermcell_pluie(klon,klev,zdt  &
-!   &      ,pplay,paprs,pphi,zlev  &
-!    &      ,u_seri,v_seri,t_seri,q_seri  &
-!    &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
-!    &      ,zfm_therm,zentr_therm,zqla  &
-!    &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
-!    &      ,tau_thermals,3)
           else if (iflag_thermals.eq.12) then
-#ifdef ISO
-              CALL abort_gcm('calltherm 282','isos pas prevus ici',1)  
-#endif
             CALL calcul_sec(klon,klev,zdt  &
      &      ,pplay,paprs,pphi,zlev  &
@@ -289,19 +280,7 @@
      &      ,tau_thermals)
           else if (iflag_thermals==13.or.iflag_thermals==14) then
-#ifdef ISO
-              CALL abort_gcm('calltherm 292','isos pas prevus ici',1)  
-#endif
-            CALL thermcellV0_main(itap,klon,klev,zdt  &
-     &      ,pplay,paprs,pphi,debut  &
-     &      ,u_seri,v_seri,t_seri,q_seri  &
-     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
-     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
-     &      ,ratqscth,ratqsdiff,zqsatth  &
-     &      ,r_aspect_thermals,l_mix_thermals  &
-     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
-     &      ,zmax0,f0,zw2,fraca)
-          else if (iflag_thermals>=15.and.iflag_thermals<=18) then
-
-!            print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
+              abort_message = 'thermcellV0_main enleve svn>2084'
+              CALL abort_physic (modname,abort_message,1)
+          else if (new_thermcell) then
             CALL thermcell_main(itap,klon,klev,zdt  &
      &      ,pplay,paprs,pphi,debut  &
@@ -310,22 +289,26 @@
      &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
      &      ,ratqscth,ratqsdiff,zqsatth  &
-!    &      ,r_aspect_thermals,l_mix_thermals &
-!    &      ,tau_thermals,iflag_thermals_ed,iflag_coupl &
-     &      ,Ale,Alp,lalim_conv,wght_th &
      &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
-     &      ,ztla,zthl &
-!!! nrlmd le 10/04/2012
-     &      ,pbl_tke,pctsrf,omega,airephy &
-     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
-     &      ,n2,s2,ale_bl_stat &
-     &      ,therm_tke_max,env_tke_max &
-     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
-     &      ,alp_bl_conv,alp_bl_stat &
-!!! fin nrlmd le 10/04/2012
-     &      ,ztva &
+     &      ,ztla,zthl,ztva &
+     &      ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
 #ifdef ISO         
      &      ,xt_seri,d_xt_the &
 #endif         
      &   ) 
+
+            CALL thermcell_alp(klon,klev,zdt  &                      ! in
+     &        ,pplay,paprs  &                                        ! in
+     &        ,zfm_therm,zentr_therm,lmax  &                         ! in
+     &        ,pbl_tke,pctsrf,omega,airephy &                        ! in
+     &        ,zw2,fraca &                                           ! in
+     &        ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
+     &        ,ale,alp,lalim_conv,wght_th &                          ! out
+     &        ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out
+     &        ,n2,s2,ale_bl_stat &                                   ! out
+     &        ,therm_tke_max,env_tke_max &                           ! out
+     &        ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
+     &        ,alp_bl_conv,alp_bl_stat &                             ! out
+     &        )
+
            if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
          else
@@ -423,9 +406,9 @@
        DO i=1,klon
             fm_therm(i,klev+1)=0.
-            Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
-!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
-            Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
-!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
-        if(prt_level.GE.10) print*,'calltherm i Alp_bl Alp Ale_bl Ale',i,Alp_bl(i),Alp(i),Ale_bl(i),Ale(i)
+            ale_bl(i)=ale_bl(i)+ale(i)/REAL(nsplit_thermals)
+!            write(22,*)'ALE CALLTHERM',ale_bl(i),ale(i)
+            alp_bl(i)=alp_bl(i)+alp(i)/REAL(nsplit_thermals)
+!            write(23,*)'ALP CALLTHERM',alp_bl(i),alp(i)
+        if(prt_level.GE.10) print*,'calltherm i alp_bl alp ale_bl ale',i,alp_bl(i),alp(i),ale_bl(i),ale(i)
        ENDDO
 
Index: /LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90	(revision 4089)
@@ -16016,5 +16016,4 @@
 #include "dimsoil.h"
 #include "clesphys.h"
-#include "thermcell.h"
 #include "compbl.h"      
 
@@ -16201,5 +16200,5 @@
 #include "dimsoil.h"
 #include "clesphys.h"
-#include "thermcell.h"
+!  #include "thermcell.h"
 #include "compbl.h"
 
@@ -16589,5 +16588,5 @@
 #include "dimsoil.h"
 #include "clesphys.h"
-#include "thermcell.h"
+! #include "thermcell.h"
 #include "compbl.h"    
 
Index: /LMDZ6/trunk/libf/phylmdiso/phyetat0.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phyetat0.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmdiso/phyetat0.F90	(revision 4089)
@@ -59,5 +59,5 @@
   include "dimsoil.h"
   include "clesphys.h"
-  include "thermcell.h"
+  include "alpale.h"
   include "compbl.h"
   include "YOMCST.h"
Index: /LMDZ6/trunk/libf/phylmdiso/phyredem.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phyredem.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmdiso/phyredem.F90	(revision 4089)
@@ -56,5 +56,5 @@
   include "dimsoil.h"
   include "clesphys.h"
-  include "thermcell.h"
+  include "alpale.h"
   include "compbl.h"
   !======================================================================
@@ -504,5 +504,5 @@
 #include "dimsoil.h"
 #include "clesphys.h"
-#include "thermcell.h"
+#include "alpale.h"
 #include "compbl.h"      
       ! inputs
Index: /LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90	(revision 4089)
@@ -59,5 +59,5 @@
     IMPLICIT NONE
     include "clesphys.h"
-    include "thermcell.h"
+    include "alpale.h"
     include "YOMCST.h"
 
Index: /LMDZ6/trunk/libf/phylmdiso/phys_state_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_state_var_mod.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_state_var_mod.F90	(revision 4089)
@@ -1,4 +1,4 @@
 !
-! $Id: phys_state_var_mod.F90 3888 2021-05-05 10:50:37Z jyg $
+! $Id: phys_state_var_mod.F90 4088 2022-03-10 07:03:20Z fhourdin $
 !
       MODULE phys_state_var_mod
Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4089)
@@ -78,4 +78,5 @@
     USE write_field_phy
     USE lscp_mod, ONLY : lscp
+    USE thermcell_ini_mod, ONLY : thermcell_ini
 
     !USE cmp_seri_mod
@@ -421,5 +422,5 @@
     include "dimsoil.h"
     include "clesphys.h"
-    include "thermcell.h"
+    include "alpale.h"
     include "dimpft.h"
     !======================================================================
@@ -1868,4 +1869,9 @@
 
        CALL iniradia(klon,klev,paprs(1,1:klev+1))
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
+   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
        !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Index: /LMDZ6/trunk/libf/phylmdiso/thermcell_main.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/thermcell_main.F90	(revision 4088)
+++ /LMDZ6/trunk/libf/phylmdiso/thermcell_main.F90	(revision 4089)
@@ -1,6 +1,6 @@
-!
+
 ! $Id: thermcell_main.F90 3451 2019-01-27 11:07:30Z fhourdin $
 !
-      SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep  &
+      subroutine thermcell_main(itap,ngrid,nlay,ptimestep  &
      &                  ,pplay,pplev,pphi,debut  &
      &                  ,pu,pv,pt,po  &
@@ -8,16 +8,7 @@
      &                  ,fm0,entr0,detr0,zqta,zqla,lmax  &
      &                  ,ratqscth,ratqsdiff,zqsatth  &
-     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
      &                  ,zmax0, f0,zw2,fraca,ztv &
-     &                  ,zpspsk,ztla,zthl &
-!!! nrlmd le 10/04/2012
-     &                  ,pbl_tke,pctsrf,omega,airephy &
-     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
-     &                  ,n2,s2,ale_bl_stat &
-     &                  ,therm_tke_max,env_tke_max &
-     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
-     &                  ,alp_bl_conv,alp_bl_stat &
-!!! fin nrlmd le 10/04/2012
-     &                  ,ztva & 
+     &                  ,zpspsk,ztla,zthl,ztva &
+     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
 #ifdef ISO         
      &      ,xtpo,xtpdoadj &
@@ -25,8 +16,9 @@
      &   )
 
-      USE dimphy
-      USE ioipsl
-      USE indice_sol_mod
-      USE print_control_mod, ONLY: lunout,prt_level
+
+      USE thermcell_ini_mod, ONLY: thermcell_ini,dqimpl,dvdq,prt_level,lunout,prt_level
+      USE thermcell_ini_mod, ONLY: iflag_thermals_closure,iflag_thermals_ed,tau_thermals,r_aspect_thermals
+      USE thermcell_ini_mod, ONLY: RD,RG
+
 #ifdef ISO
   USE infotrac_phy, ONLY : ntraciso
@@ -37,4 +29,6 @@
 #endif
 #endif
+
+
       IMPLICIT NONE
 
@@ -74,35 +68,29 @@
 !   -------------
 
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"
-#include "thermcell.h"
 
 !   arguments:
 !   ----------
-
-!IM 140508
-      INTEGER itap
-
-      INTEGER ngrid,nlay
-      real ptimestep
-      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)
-      LOGICAL debut
+      integer, intent(in) :: itap,ngrid,nlay
+      real, intent(in) ::  ptimestep
+      real, intent(in), dimension(ngrid,nlay)    :: pt,pu,pv,po,pplay,pphi,zpspsk
+      real, intent(in), dimension(ngrid,nlay+1)  :: pplev
+      real, intent(out), dimension(ngrid,nlay)   :: pdtadj,pduadj,pdvadj,pdoadj,entr0,detr0
+      real, intent(out), dimension(ngrid,nlay)   :: ztla,zqla,zqta,zqsatth,zthl
+      real, intent(out), dimension(ngrid,nlay+1) :: fm0,zw2,fraca
+      real, intent(out), dimension(ngrid) :: zmax0,f0
+      real, intent(out), dimension(ngrid,nlay) :: ztva,ztv
+      logical, intent(in) :: debut
+
+      real, intent(out), dimension(ngrid) :: pcon
+      real, intent(out), dimension(ngrid,nlay) :: rhobarz,wth3
+      real, intent(out), dimension(ngrid) :: wmax_sec
+      integer,intent(out), dimension(ngrid) :: lalim
+      real, intent(out), dimension(ngrid,nlay+1) :: fm
+      real, intent(out), dimension(ngrid,nlay) :: alim_star
+      real, intent(out), dimension(ngrid) :: zmax
 
 !   local:
 !   ------
 
-      integer icount
-
-      integer, save :: dvdq=1,dqimpl=-1
-!$OMP THREADPRIVATE(dvdq,dqimpl)
-      data icount/0/
-      save icount
-!$OMP THREADPRIVATE(icount)
 
       integer,save :: igout=1
@@ -113,159 +101,33 @@
 !$OMP THREADPRIVATE(lev_out)
 
-      REAL susqr2pi, Reuler
-
-      INTEGER ig,k,l,ll,ierr
-      real zsortie1d(klon)
-      INTEGER lmax(klon),lmin(klon),lalim(klon)
-      INTEGER lmix(klon)
-      INTEGER lmix_bis(klon)
-      real linter(klon)
-      real zmix(klon)
-      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
-!      real fraca(klon,klev)
-
-      real zmax_sec(klon)
-!on garde le zmax du pas de temps precedent
-      real zmax0(klon)
-!FH/IM     save zmax0
-
-      real lambda
-
-      real zlev(klon,klev+1),zlay(klon,klev)
-      real deltaz(klon,klev)
-      REAL zh(klon,klev)
-      real zthl(klon,klev),zdthladj(klon,klev)
-      REAL ztv(klon,klev)
-      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
-      real zl(klon,klev)
-      real zsortie(klon,klev)
-      real zva(klon,klev)
-      real zua(klon,klev)
-      real zoa(klon,klev)
-
-      real zta(klon,klev)
-      real zha(klon,klev)
-      real fraca(klon,klev+1)
-      real zf,zf2
-      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
-      real q2(klon,klev)
-! FH probleme de dimensionnement avec l'allocation dynamique
-!     common/comtherm/thetath2,wth2
-      real wq(klon,klev)
-      real wthl(klon,klev)
-      real wthv(klon,klev)
-    
-      real ratqscth(klon,klev)
-      real var
-      real vardiff
-      real ratqsdiff(klon,klev)
-
+      real lambda, zf,zf2,var,vardiff,CHI
+      integer ig,k,l,ierr,ll
       logical sorties
-      real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
-      real zpspsk(klon,klev)
-
-      real wmax(klon)
-      real wmax_tmp(klon)
-      real wmax_sec(klon)
-      real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
-      real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
-
-      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
-!niveau de condensation
-      integer nivcon(klon)
-      real zcon(klon)
-      REAL CHI
-      real zcon2(klon)
-      real pcon(klon)
-      real zqsat(klon,klev)
-      real zqsatth(klon,klev) 
-
-      real f_star(klon,klev+1),entr_star(klon,klev)
-      real detr_star(klon,klev)
-      real alim_star_tot(klon)
-      real alim_star(klon,klev)
-      real alim_star_clos(klon,klev)
-      real f(klon), f0(klon)
-!FH/IM     save f0
-      real zlevinter(klon)
-       real seuil
-      real csc(klon,klev)
-
-!!! nrlmd le 10/04/2012
-
-!------Entrées
-      real pbl_tke(klon,klev+1,nbsrf)
-      real pctsrf(klon,nbsrf)
-      real omega(klon,klev)
-      real airephy(klon)
-!------Sorties
-      real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
-      real therm_tke_max0(klon),env_tke_max0(klon)
-      real n2(klon),s2(klon)
-      real ale_bl_stat(klon)
-      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
-      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
-!------Local
-      integer nsrf
-      real rhobarz0(klon)                    ! Densité au LCL
-      logical ok_lcl(klon)                   ! Existence du LCL des thermiques
-      integer klcl(klon)                     ! Niveau du LCL
-      real interp(klon)                      ! Coef d'interpolation pour le LCL
-!--Triggering
-      real Su                                ! Surface unité: celle d'un updraft élémentaire
-      parameter(Su=4e4)
-      real hcoef                             ! Coefficient directeur pour le calcul de s2
-      parameter(hcoef=1)
-      real hmincoef                          ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 
-      parameter(hmincoef=0.3)
-      real eps1                              ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
-      parameter(eps1=0.3)
-      real hmin(ngrid)                       ! Ordonnée à l'origine pour le calcul de s2
-      real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
-      real zmax_moy_coef
-      parameter(zmax_moy_coef=0.33)
-      real depth(klon)                       ! Epaisseur moyenne du cumulus
-      real w_max(klon)                       ! Vitesse max statistique 
-      real s_max(klon)
-!--Closure
-      real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne 
-      real pbl_tke_max0(klon)                ! TKE moyenne au LCL
-      real w_ls(klon,klev)                   ! Vitesse verticale grande échelle (m/s)
-      real coef_m                            ! On considère un rendement pour alp_bl_fluct_m
-      parameter(coef_m=1.)
-      real coef_tke                          ! On considère un rendement pour alp_bl_fluct_tke
-      parameter(coef_tke=1.)
-
-!!! fin nrlmd le 10/04/2012
-
-!
-      !nouvelles variables pour la convection
-      real Ale_bl(klon)
-      real Alp_bl(klon)
-      real alp_int(klon),dp_int(klon),zdp
-      real ale_int(klon)
-      integer n_int(klon)
-      real fm_tot(klon)
-      real wght_th(klon,klev)
-      integer lalim_conv(klon)
-!v1d     logical therm
-!v1d     save therm
-
-      character*2 str2
-      character*10 str10
+      real, dimension(ngrid) :: linter,zmix, zmax_sec
+      integer,dimension(ngrid) :: lmax,lmin,lmix,lmix_bis,nivcon
+      real, dimension(ngrid,nlay) :: ztva_est
+      real, dimension(ngrid,nlay) :: deltaz,zlay,zh,zdthladj,zu,zv,zo,zl,zva,zua,zoa
+      real, dimension(ngrid,nlay) :: zta,zha,q2,wq,wthl,wthv,thetath2,wth2
+      real, dimension(ngrid,nlay) :: ratqscth,ratqsdiff,rho,masse
+      real, dimension(ngrid,nlay+1) :: zw_est,zlev
+      real, dimension(ngrid) :: wmax,wmax_tmp
+      real, dimension(ngrid,nlay+1) :: f_star
+      real, dimension(ngrid,nlay) :: entr,detr,entr_star,detr_star,alim_star_clos
+      real, dimension(ngrid,nlay) :: zqsat,csc
+      real, dimension(ngrid) :: zcon,zcon2,alim_star_tot,f
 
       character (len=20) :: modname='thermcell_main'
       character (len=80) :: abort_message
 
-      EXTERNAL SCOPY
 
 #ifdef ISO
       REAL xtpo(ntraciso,ngrid,nlay),xtpdoadj(ntraciso,ngrid,nlay)
-      REAL xtzo(ntraciso,klon,klev)
+      REAL xtzo(ntraciso,ngrid,nlay)
       REAL xtpdoadj_tmp(ngrid,nlay)
-      REAL xtpo_tmp(klon,klev)
-      REAL xtzo_tmp(klon,klev)
+      REAL xtpo_tmp(ngrid,nlay)
+      REAL xtzo_tmp(ngrid,nlay)
       integer ixt
 #endif
+
 !
 
@@ -274,43 +136,21 @@
 !   ---------------
 !
-
-   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
+    print*,'NEW THERMCELL cool'
+
+
    fm=0. ; entr=0. ; detr=0.
-   icount=icount+1
-
-!IM 090508 beg
-!print*,'====================================================================='
-!print*,'====================================================================='
-!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
-!print*,'====================================================================='
-!print*,'====================================================================='
-!IM 090508 end
 
       if (prt_level.ge.1) print*,'thermcell_main V4'
 
        sorties=.true.
-      IF(ngrid.NE.klon) THEN
+      IF(ngrid.NE.ngrid) THEN
          PRINT*
          PRINT*,'STOP dans convadj'
          PRINT*,'ngrid    =',ngrid
-         PRINT*,'klon  =',klon
+         PRINT*,'ngrid  =',ngrid
       ENDIF
 !
 !     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
-     do ig=1,klon
+     do ig=1,ngrid
          f0(ig)=max(f0(ig),1.e-2)
          zmax0(ig)=max(zmax0(ig),40.)
@@ -357,34 +197,26 @@
          zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
       enddo
-         zlev(:,1)=0.
-         zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
+      zlev(:,1)=0.
+      zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
       do l=1,nlay
          zlay(:,l)=pphi(:,l)/RG
       enddo
-!calcul de l epaisseur des couches
       do l=1,nlay
          deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
       enddo
 
-!     print*,'2 OK convect8'
 !-----------------------------------------------------------------------
-!   Calcul des densites
+!   Calcul des densites et masses
 !-----------------------------------------------------------------------
 
-     rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))
-
-     if (prt_level.ge.10)write(lunout,*)                                &
-    &    'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
+      rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))
+      if (prt_level.ge.10) write(lunout,*) 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
       rhobarz(:,1)=rho(:,1)
-
       do l=2,nlay
          rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
       enddo
-
-!calcul de la masse
       do l=1,nlay
          masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
       enddo
-
       if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
 
@@ -501,6 +333,6 @@
       if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
 
-      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
+      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
+      call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
 
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
@@ -530,8 +362,8 @@
 
 
-      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
+      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
+      call test_ltherm(ngrid,nlay,pplay,lmin ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
+      call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
+      call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
 
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
@@ -547,6 +379,6 @@
 
  
-call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
-call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
+call test_ltherm(ngrid,nlay,pplay,lmin,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
+call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
 
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
@@ -603,5 +435,4 @@
 !-------------------------------------------------------------------------------
 !deduction des flux
-!-------------------------------------------------------------------------------
 
       CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
@@ -612,6 +443,6 @@
 
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
-      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
+      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
+      call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
 
 !------------------------------------------------------------------
@@ -640,4 +471,5 @@
       call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
      &                   po,pdoadj,zoa,lev_out)
+
 #ifdef ISO
         ! C Risi: on utilise directement la même routine
@@ -675,17 +507,19 @@
         enddo
       enddo !DO  ll=1,nlay 
-      write(*,*) 'thermcell_main 600 tmp: apres thermcell_dq' 
+      write(*,*) 'thermcell_main 600 tmp: apres thermcell_dq'
 #endif      
 #endif
 
+
+
 !------------------------------------------------------------------
 ! Calcul de la fraction de l'ascendance
 !------------------------------------------------------------------
-      do ig=1,klon
+      do ig=1,ngrid
          fraca(ig,1)=0.
          fraca(ig,nlay+1)=0.
       enddo
       do l=2,nlay
-         do ig=1,klon
+         do ig=1,ngrid
             if (zw2(ig,l).gt.1.e-10) then
             fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
@@ -819,25 +653,4 @@
          enddo
       enddo
-!
-! $Id: thermcell_main.F90 3451 2019-01-27 11:07:30Z fhourdin $
-!
-      CALL thermcell_alp(ngrid,nlay,ptimestep  &
-     &                  ,pplay,pplev  &
-     &                  ,fm0,entr0,lmax  &
-     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
-     &                  ,zw2,fraca &
-!!! necessire en plus
-     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
-!!! nrlmd le 10/04/2012
-     &                  ,pbl_tke,pctsrf,omega,airephy &
-     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
-     &                  ,n2,s2,ale_bl_stat &
-     &                  ,therm_tke_max,env_tke_max &
-     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
-     &                  ,alp_bl_conv,alp_bl_stat &
-!!! fin nrlmd le 10/04/2012
-     &                   )
-
-
 
 !calcul du ratqscdiff
@@ -847,5 +660,5 @@
       ratqsdiff(:,:)=0.
 
-      do l=1,klev
+      do l=1,nlay
          do ig=1,ngrid
             if (l<=lalim(ig)) then
@@ -857,5 +670,5 @@
       if (prt_level.ge.1) print*,'14f OK convect8'
 
-      do l=1,klev
+      do l=1,nlay
          do ig=1,ngrid
             if (l<=lalim(ig)) then
@@ -868,15 +681,9 @@
 
       if (prt_level.ge.1) print*,'14g OK convect8'
-      do l=1,nlay
-         do ig=1,ngrid
-            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
-!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
-         enddo
-      enddo 
-!--------------------------------------------------------------------    
-!
-!ecriture des fichiers sortie
-!     print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
-
+         do l=1,nlay
+            do ig=1,ngrid
+               ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
+            enddo
+         enddo 
       endif
 
@@ -884,23 +691,22 @@
 
       return
-      end
-
-!-----------------------------------------------------------------------------
-
-      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
-      USE print_control_mod, ONLY: prt_level
+      end subroutine thermcell_main
+
+!=============================================================================
+!/////////////////////////////////////////////////////////////////////////////
+!=============================================================================
+      subroutine test_ltherm(ngrid,nlay,pplay,long,ztv,po,ztva, &  ! in
+    &            zqla,f_star,zw2,comment)                          ! in
+!=============================================================================
+      USE thermcell_ini_mod, ONLY: prt_level
       IMPLICIT NONE
 
-      integer i, k, klon,klev
-      real pplev(klon,klev+1),pplay(klon,klev)
-      real ztv(klon,klev)
-      real po(klon,klev)
-      real ztva(klon,klev)
-      real zqla(klon,klev)
-      real f_star(klon,klev)
-      real zw2(klon,klev)
-      integer long(klon)
+      integer i, k, ngrid,nlay
+      real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,po,ztva,zqla
+      real, intent(in), dimension(ngrid,nlay) :: f_star,zw2
+      integer, intent(in), dimension(ngrid) :: long
       real seuil
       character*21 comment
+      seuil=0.25
 
       if (prt_level.ge.1) THEN
@@ -910,10 +716,10 @@
 
 !  test sur la hauteur des thermiques ...
-         do i=1,klon
+         do i=1,ngrid
 !IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
            if (prt_level.ge.10) then
                print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
                print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
-               do k=1,klev
+               do k=1,nlay
                   write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
                enddo
@@ -925,10 +731,16 @@
       end
 
-!!! nrlmd le 10/04/2012                          Transport de la TKE par le thermique moyen pour la fermeture en ALP 
-!                                                         On transporte pbl_tke pour donner therm_tke
-!                                          Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
-      subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
-     &           rg,pplev,therm_tke_max)
-      USE print_control_mod, ONLY: prt_level
+! nrlmd le 10/04/2012   Transport de la TKE par le thermique moyen pour la fermeture en ALP 
+!                       On transporte pbl_tke pour donner therm_tke
+!                       Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
+
+!=======================================================================
+!///////////////////////////////////////////////////////////////////////
+!=======================================================================
+
+      subroutine thermcell_tke_transport( &
+     &     ngrid,nlay,ptimestep,fm0,entr0,rg,pplev,  &   ! in
+     &     therm_tke_max)                                ! out
+      USE thermcell_ini_mod, ONLY: prt_level
       implicit none
 
@@ -941,13 +753,14 @@
 !=======================================================================
 
-      integer ngrid,nlay,nsrf
-
-      real ptimestep
-      real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
-      real entr0(ngrid,nlay),rg
-      real therm_tke_max(ngrid,nlay)
+      integer ngrid,nlay
+
+      real, intent(in) :: ptimestep
+      real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev
+      real, intent(in), dimension(ngrid,nlay) :: entr0
+      real, intent(in) :: rg
+      real, intent(out), dimension(ngrid,nlay) :: therm_tke_max
+
       real detr0(ngrid,nlay)
-
-
+      real masse0(ngrid,nlay)
       real masse(ngrid,nlay),fm(ngrid,nlay+1)
       real entr(ngrid,nlay)
@@ -956,9 +769,5 @@
 
       real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
-
-      real zzm
-
       integer ig,k
-      integer isrf
 
 
@@ -988,14 +797,5 @@
       fm(:,nlay+1)=0.
 
-!!! nrlmd le 16/09/2010
-!   calcul de la valeur dans les ascendances
-!       do ig=1,ngrid
-!          qa(ig,1)=q(ig,1)
-!       enddo
-!!! 
-
-!do isrf=1,nsrf
-
-!   q(:,:)=therm_tke(:,:,isrf)
+
    q(:,:)=therm_tke_max(:,:)
 !!! nrlmd le 16/09/2010
