Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90	(revision 1292)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90	(revision 1294)
@@ -1,6 +1,9 @@
+!
+! $Header$
+!
       SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
-     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
-     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
-     &           ztla,zqla,zqta,zha,zw2,w_est,zqsatth,lmix,lmix_bis,linter &
+     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
+     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
+     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
      &           ,lev_out,lunout1,igout)
 
@@ -31,8 +34,5 @@
       REAL zpspsk(ngrid,klev)
       REAL alim_star(ngrid,klev)
-      REAL zmax_sec(ngrid)
       REAL f0(ngrid)
-      REAL l_mix
-      REAL r_aspect
       INTEGER lalim(ngrid)
       integer lev_out                           ! niveau pour les print
@@ -44,5 +44,4 @@
       REAL ztla(ngrid,klev)
       REAL zqla(ngrid,klev)
-      REAL zqla0(ngrid,klev)
       REAL zqta(ngrid,klev)
       REAL zha(ngrid,klev)
@@ -50,16 +49,9 @@
       REAL detr_star(ngrid,klev)
       REAL coefc
-      REAL detr_stara(ngrid,klev)
-      REAL detr_starb(ngrid,klev)
-      REAL detr_starc(ngrid,klev)
-      REAL detr_star0(ngrid,klev)
-      REAL detr_star1(ngrid,klev)
-      REAL detr_star2(ngrid,klev)
-
       REAL entr_star(ngrid,klev)
-      REAL entr_star1(ngrid,klev)
-      REAL entr_star2(ngrid,klev)
       REAL detr(ngrid,klev)
       REAL entr(ngrid,klev)
+
+      REAL csc(ngrid,klev)
 
       REAL zw2(ngrid,klev+1)
@@ -72,4 +64,7 @@
       REAL zqsatth(ngrid,klev)
       REAL zta_est(ngrid,klev)
+      REAL zdw2
+      REAL zw2modif
+      REAL zeps
 
       REAL linter(ngrid)
@@ -80,4 +75,5 @@
       INTEGER ig,l,k
 
+      real zdz,zfact,zbuoy,zalpha
       real zcor,zdelta,zcvm5,qlbef
       real Tbef,qsatbef
@@ -86,7 +82,9 @@
       PARAMETER (DDT0=.01)
       logical Zsat
-      REAL fact_gamma,fact_epsilon
+      LOGICAL active(ngrid),activetmp(ngrid)
+      REAL fact_gamma,fact_epsilon,fact_gamma2
       REAL c2(ngrid,klev)
 
+      REAL zw2fact,expa
       Zsat=.false.
 ! Initialisation
@@ -97,124 +95,130 @@
          fact_epsilon=1.
       else if (iflag_thermals_ed==1)  then
+! Valeurs au moment de la premiere soumission des papiers
          fact_gamma=1.
-         fact_epsilon=1.
+         fact_epsilon=0.002
+         fact_gamma2=0.6
+! Suggestions des Fleurs, Septembre 2009
+         fact_epsilon=0.015
+!test cr
+!         fact_epsilon=0.002
+         fact_gamma=0.9
+         fact_gamma2=0.7
+
       else if (iflag_thermals_ed==2)  then
          fact_gamma=1.
          fact_epsilon=2.
+      else if (iflag_thermals_ed==3)  then
+         fact_gamma=3./4.
+         fact_epsilon=3.
       endif
 
-      do l=1,klev
+      write(lunout,*)'THERM 31H '
+
+! Initialisations des variables reeles
+if (1==0) then
+      ztva(:,:)=ztv(:,:)
+      ztva_est(:,:)=ztva(:,:)
+      ztla(:,:)=zthl(:,:)
+      zqta(:,:)=po(:,:)
+      zha(:,:) = ztva(:,:)
+else
+      ztva(:,:)=0.
+      ztva_est(:,:)=0.
+      ztla(:,:)=0.
+      zqta(:,:)=0.
+      zha(:,:) =0.
+endif
+
+      zqla_est(:,:)=0.
+      zqsatth(:,:)=0.
+      zqla(:,:)=0.
+      detr_star(:,:)=0.
+      entr_star(:,:)=0.
+      alim_star(:,:)=0.
+      alim_star_tot(:)=0.
+      csc(:,:)=0.
+      detr(:,:)=0.
+      entr(:,:)=0.
+      zw2(:,:)=0.
+      w_est(:,:)=0.
+      f_star(:,:)=0.
+      wa_moy(:,:)=0.
+      linter(:)=1.
+      linter(:)=1.
+
+! Initialisation des variables entieres
+      lmix(:)=1
+      lmix_bis(:)=2
+      wmaxa(:)=0.
+      lalim(:)=1
+
+!-------------------------------------------------------------------------
+! On ne considere comme actif que les colonnes dont les deux premieres
+! couches sont instables.
+!-------------------------------------------------------------------------
+      active(:)=ztv(:,1)>ztv(:,2)
+
+!-------------------------------------------------------------------------
+! Definition de l'alimentation a l'origine dans thermcell_init
+!-------------------------------------------------------------------------
+      do l=1,klev-1
          do ig=1,ngrid
-            zqla_est(ig,l)=0.
-            ztva_est(ig,l)=ztva(ig,l)
-            zqsatth(ig,l)=0.
+            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
+               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
+     &                       *sqrt(zlev(ig,l+1)) 
+               lalim(:)=l+1
+               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
+            endif
          enddo
       enddo
-
-!CR: attention test couche alim
-!     do l=2,klev
-!     do ig=1,ngrid
-!        alim_star(ig,l)=0.
-!     enddo
-!     enddo
-!AM:initialisations du thermique
-      do k=1,klev
-         do ig=1,ngrid
-            ztva(ig,k)=ztv(ig,k)
-            ztla(ig,k)=zthl(ig,k)
-            zqla(ig,k)=0.
-            zqta(ig,k)=po(ig,k)
-!
-            ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+RLvCp*zqla(ig,k)
-            ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
-            zha(ig,k) = ztva(ig,k)
-!
-         enddo
-      enddo 
-      do k=1,klev
-        do ig=1,ngrid
-           detr_star(ig,k)=0.
-           entr_star(ig,k)=0.
-
-           detr_stara(ig,k)=0.
-           detr_starb(ig,k)=0.
-           detr_starc(ig,k)=0.
-           detr_star0(ig,k)=0.
-           zqla0(ig,k)=0.
-           detr_star1(ig,k)=0.
-           detr_star2(ig,k)=0.
-           entr_star1(ig,k)=0.
-           entr_star2(ig,k)=0.
-
-           detr(ig,k)=0.
-           entr(ig,k)=0.
-        enddo
-      enddo
-      if (prt_level.ge.1) print*,'7 OK convect8'
-      do k=1,klev+1
-         do ig=1,ngrid
-            zw2(ig,k)=0.
-            w_est(ig,k)=0.
-            f_star(ig,k)=0.
-            wa_moy(ig,k)=0.
+      do l=1,klev
+         do ig=1,ngrid 
+            if (alim_star_tot(ig) > 1.e-10 ) then
+               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
+            endif
          enddo
       enddo
-
-      if (prt_level.ge.1) print*,'8 OK convect8'
-      do ig=1,ngrid
-         linter(ig)=1.
-         lmix(ig)=1
-         lmix_bis(ig)=2
-         wmaxa(ig)=0.
-      enddo
-
-!-----------------------------------------------------------------------------------
-!boucle de calcul de la vitesse verticale dans le thermique
-!-----------------------------------------------------------------------------------
-      do l=1,klev-1
-         do ig=1,ngrid
-
-
-
-! Calcul dans la premiere couche active du thermique (ce qu'on teste
-! en disant que la couche est instable et que w2 en bas de la couche
-! est nulle.
-
-            if (ztv(ig,l).gt.ztv(ig,l+1)  &
-     &         .and.alim_star(ig,l).gt.1.e-10  &
-     &         .and.zw2(ig,l).lt.1e-10) then
-
-
+      alim_star_tot(:)=1.
+
+
+!------------------------------------------------------------------------------
+! Calcul dans la premiere couche
+! On decide dans cette version que le thermique n'est actif que si la premiere
+! couche est instable.
+! Pourrait etre change si on veut que le thermiques puisse se déclencher
+! dans une couche l>1
+!------------------------------------------------------------------------------
+do ig=1,ngrid
 ! Le panache va prendre au debut les caracteristiques de l'air contenu
 ! dans cette couche.
-               ztla(ig,l)=zthl(ig,l) 
-               zqta(ig,l)=po(ig,l)
-               zqla(ig,l)=zl(ig,l)
-               f_star(ig,l+1)=alim_star(ig,l)
-
-               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
-     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
-     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
-               w_est(ig,l+1)=zw2(ig,l+1)
+    if (active(ig)) then
+    ztla(ig,1)=zthl(ig,1) 
+    zqta(ig,1)=po(ig,1)
+    zqla(ig,1)=zl(ig,1)
+!cr: attention, prise en compte de f*(1)=1
+    f_star(ig,2)=alim_star(ig,1)
+    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
+&                     *(zlev(ig,2)-zlev(ig,1))  &
+&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
+    w_est(ig,2)=zw2(ig,2)
+    endif
+enddo
 !
 
-
-            else if ((zw2(ig,l).ge.1e-10).and.  &
-     &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
-!estimation du detrainement a partir de la geometrie du pas precedent
-!tests sur la definition du detr
-!calcul de detr_star et entr_star
-
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! FH le test miraculeux de Catherine ? Le bout du tunel ?
-!               w_est(ig,3)=zw2(ig,2)*  &
-!    &                   ((f_star(ig,2))**2)  &
-!    &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
-!    &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
-!    &                   *(zlev(ig,3)-zlev(ig,2))
-!     if (l.gt.2) then
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!==============================================================================
+!boucle de calcul de la vitesse verticale dans le thermique
+!==============================================================================
+do l=2,klev-1
+!==============================================================================
+
+
+! On decide si le thermique est encore actif ou non
+! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
+    do ig=1,ngrid
+       active(ig)=active(ig) &
+&                 .and. zw2(ig,l)>1.e-10 &
+&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
+    enddo
 
 
@@ -222,62 +226,29 @@
 ! Premier calcul de la vitesse verticale a partir de la temperature
 ! potentielle virtuelle
-
-! FH CESTQUOI CA ????
-#define int1d2
-!#undef int1d2
-#ifdef int1d2
-      if (l.ge.2) then
-#else
-      if (l.gt.2) then
-#endif
-
-      if (1.eq.1) then
-          w_est(ig,3)=zw2(ig,2)* &
-     &      ((f_star(ig,2))**2) &
-     &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
-     &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
-!     &      *1./3. &
-     &      *(zlev(ig,3)-zlev(ig,2))
-       endif
+!     if (1.eq.1) then
+!         w_est(ig,3)=zw2(ig,2)* &
+!    &      ((f_star(ig,2))**2) &
+!    &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
+!    &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
+!    &      *(zlev(ig,3)-zlev(ig,2))
+!      endif
 
 
 !---------------------------------------------------------------------------
-!calcul de l entrainement et du detrainement lateral
+! calcul des proprietes thermodynamiques et de la vitesse de la couche l
+! sans tenir compte du detrainement et de l'entrainement dans cette
+! couche
+! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
+! avant) a l'alimentation pour avoir un calcul plus propre
 !---------------------------------------------------------------------------
-!
-!test:estimation de ztva_new_est sans entrainement
-
-               Tbef=ztla(ig,l-1)*zpspsk(ig,l)
-               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
-               qsatbef=MIN(0.5,qsatbef)
-               zcor=1./(1.-retv*qsatbef)
-               qsatbef=qsatbef*zcor
-               Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
-               if (Zsat) then
-               qlbef=max(0.,zqta(ig,l-1)-qsatbef)
-               DT = 0.5*RLvCp*qlbef
-               do while (abs(DT).gt.DDT0)
-                 Tbef=Tbef+DT
-                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
-                 qsatbef=MIN(0.5,qsatbef)
-                 zcor=1./(1.-retv*qsatbef)
-                 qsatbef=qsatbef*zcor
-                 qlbef=zqta(ig,l-1)-qsatbef
-
-                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
-                 zcor=1./(1.-retv*qsatbef)
-                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
-                 num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef
-                 denom=1.+RLvCp*dqsat_dT
-                 DT=num/denom
-               enddo
-                 zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef) 
-               endif
+
+     call thermcell_condens(ngrid,active, &
+&          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
+
+    do ig=1,ngrid
+        if(active(ig)) then
         ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
+        zta_est(ig,l)=ztva_est(ig,l)
         ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
-        zta_est(ig,l)=ztva_est(ig,l)
         ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
      &      -zqla_est(ig,l))-zqla_est(ig,l))
@@ -287,412 +258,117 @@
      &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
      &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
-!     &                   *1./3. &
      &                   *(zlev(ig,l+1)-zlev(ig,l))
              if (w_est(ig,l+1).lt.0.) then
                 w_est(ig,l+1)=zw2(ig,l)
              endif
-!
-!calcul du detrainement
-!=======================
-
-!CR:on vire les modifs
-         if (iflag_thermals_ed==0) then
-
-! Modifications du calcul du detrainement.
-! Dans la version de la these de Catherine, on passe brusquement
-! de la version seche a la version nuageuse pour le detrainement
-! ce qui peut occasioner des oscillations.
-! dans la nouvelle version, on commence par calculer un detrainement sec.
-! Puis un autre en cas de nuages.
-! Puis on combine les deux lineairement en fonction de la quantite d'eau.
-
-#define int1d3
-!#undef int1d3
-#define RIO_TH
-#ifdef RIO_TH
-!1. Cas non nuageux
-! 1.1 on est sous le zmax_sec et w croit
-          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
-     &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
-#ifdef int1d3
-     &       (zqla_est(ig,l).lt.1.e-10)) then 
-#else
-     &       (zqla(ig,l-1).lt.1.e-10)) then 
-#endif
-             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
-     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
-     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
-     &       /(r_aspect*zmax_sec(ig)))
-             detr_stara(ig,l)=detr_star(ig,l)
-
-       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
-
-! 1.2 on est sous le zmax_sec et w decroit
-          else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
-#ifdef int1d3
-     &            (zqla_est(ig,l).lt.1.e-10)) then
-#else
-     &            (zqla(ig,l-1).lt.1.e-10)) then
-#endif
-             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
-     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
-     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
-     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
-     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
-     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
-     &       *((zmax_sec(ig)-zlev(ig,l))/  &
-     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
-             detr_starb(ig,l)=detr_star(ig,l)
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
-
+       endif
+    enddo
+
+!-------------------------------------------------
+!calcul des taux d'entrainement et de detrainement
+!-------------------------------------------------
+
+     do ig=1,ngrid
+        if (active(ig)) then
+          zdz=zlev(ig,l+1)-zlev(ig,l)
+          zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
+          zfact=fact_gamma/(1.+fact_gamma)
+ 
+! estimation de la fraction couverte par les thermiques
+          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
+
+!calcul de la soumission papier
+          if (1.eq.1) then
+             fact_epsilon=0.0007
+!             zfact=0.9/(1.+0.9)
+             zfact=0.3
+             fact_gamma=0.7
+             fact_gamma2=0.6
+             expa=0.25
+! Calcul  du taux d'entrainement entr_star (epsilon)
+           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
+     &     zbuoy/w_est(ig,l+1) )&
+!- fact_epsilon/zalpha**0.25  ) &
+     &     +0.000 )
+
+!           entr_star(ig,l)=f_star(ig,l)*zdz * (  1./3 * MAX(0.,  &     
+!     &     zbuoy/w_est(ig,l+1) - 1./zalpha**0.25  ) &
+!     &     +0.000 )
+! Calcul du taux de detrainement detr_star (delta)
+!           if (zqla_est(ig,l).lt.1.e-10) then
+!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
+!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
+!     &     +0.0006 )
+!           else
+!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
+!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
+!     &     +0.002 )
+!           endif
+           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
+     &     fact_gamma2 * MAX(0., &
+!fact_epsilon/zalpha**0.25
+     &      -zbuoy/w_est(ig,l+1) )       &
+!     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
+!test
+     &     +  0.006*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &    
+     &     +0.0000 )
           else
 
-! 1.3 dans les autres cas
-             detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
-     &                      *(zlev(ig,l+1)-zlev(ig,l))
-             detr_starc(ig,l)=detr_star(ig,l)
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
-             
+! Calcul  du taux d'entrainement entr_star (epsilon)
+           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
+     &     zbuoy/w_est(ig,l+1) - fact_epsilon  ) &
+     &     +0.0000 )
+
+! Calcul du taux de detrainement detr_star (delta)
+           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
+     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
+     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
+     &     +0.0000 )
+
           endif
-
-#else
-
-! 1.1 on est sous le zmax_sec et w croit
-          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
-     &       (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
-             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
-     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
-     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
-     &       /(r_aspect*zmax_sec(ig)))
-
-       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
-
-! 1.2 on est sous le zmax_sec et w decroit
-          else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
-             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
-     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
-     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
-     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
-     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
-     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
-     &       *((zmax_sec(ig)-zlev(ig,l))/  &
-     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
-       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
-
-          else
-             detr_star=0.
-          endif
-
-! 1.3 dans les autres cas
-          detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
-     &                      *(zlev(ig,l+1)-zlev(ig,l))
-
-          coefc=min(zqla(ig,l-1)/1.e-3,1.)
-          if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
-          coefc=1.
-! il semble qu'il soit important de baser le calcul sur
-! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
-          detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc)
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
-
-#endif
-
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
-!IM 730508 beg
-!        if(itap.GE.7200) THEN
-!         print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
-!        endif
-!IM 730508 end
-         
-         zqla0(ig,l)=zqla_est(ig,l)
-         detr_star0(ig,l)=detr_star(ig,l)
-!IM 060508 beg
-!         if(detr_star(ig,l).GT.1.) THEN
-!          print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
-!   &      ,detr_starc(ig,l),coefc
-!         endif
-!IM 060508 end
-!IM 160508 beg
-!IM 160508       IF (f0(ig).NE.0.) THEN
-           detr_star(ig,l)=detr_star(ig,l)/f0(ig)
-!IM 160508       ELSE IF(detr_star(ig,l).EQ.0.) THEN
-!IM 160508        print*,'WARNING1  : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
-!IM 160508       ELSE
-!IM 160508        print*,'WARNING2  : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
-!IM 160508       ENDIF
-!IM 160508 end
-!IM 060508 beg
-!        if(detr_star(ig,l).GT.1.) THEN
-!         print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
-!   &     float(1)/f0(ig)
-!        endif
-!IM 060508 end
-        if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
-!
-!calcul de entr_star
-
-! #undef test2
-! #ifdef test2
-! La version test2 destabilise beaucoup le modele.
-! Il semble donc que ca aide d'avoir un entrainement important sous
-! le nuage.
-!         if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
-!          entr_star(ig,l)=0.4*detr_star(ig,l)
-!         else
-!          entr_star(ig,l)=0.
-!         endif
-! #else
-!
-! Deplacement du calcul de entr_star pour eviter d'avoir aussi
-! entr_star > fstar.
-! Redeplacer suite a la transformation du cas detr>f
-! FH
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
-#define int1d
-!FH 070508 #define int1d4
-!#undef int1d4
-! L'option int1d4 correspond au choix dans le cas ou le detrainement
-! devient trop grand.
-
-#ifdef int1d
-
-#ifdef int1d4
-#else
-       detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
-!FH 070508 plus
-       detr_star(ig,l)=min(detr_star(ig,l),1.)
-#endif
-
-       entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
-#ifdef int1d4
-! Si le detrainement excede le flux en bas + l'entrainement, le thermique
-! doit disparaitre.
-       if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
-          detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
-          f_star(ig,l+1)=0.
-          linter(ig)=l+1
-          zw2(ig,l+1)=-1.e-10
-       endif
-#endif
-
-
-#else
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
-        if(l.gt.lalim(ig)) then
-         entr_star(ig,l)=0.4*detr_star(ig,l)
-        else
-
-! FH :
-! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
-! en haut de la couche d'alimentation.
-! A remettre en questoin a la premiere occasion mais ca peut aider a 
-! ecrire un code robuste.
-! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
-! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
-! d'alimentation, ce qui n'est pas forcement heureux.
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
-#undef pre_int1c
-#ifdef pre_int1c
-         entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
-         detr_star(ig,l)=entr_star(ig,l)
-#else
-         entr_star(ig,l)=0.
-#endif
-
+!endif choix du calcul de E* et D*
+
+!cr test
+!           entr_star(ig,l)=entr_star(ig,l)+0.2*detr_star(ig,l)     
+
+! Prise en compte de la fraction
+!      detr_star(ig,l)=detr_star(ig,l)*sqrt(0.01/max(zalpha,1.e-5))
+
+! En dessous de lalim, on prend le max de alim_star et entr_star pour
+! alim_star et 0 sinon
+        if (l.lt.lalim(ig)) then
+          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
+          entr_star(ig,l)=0.
         endif
 
-#endif
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
-        entr_star1(ig,l)=entr_star(ig,l)
-        detr_star1(ig,l)=detr_star(ig,l)
-!
-
-#ifdef int1d
-#else
-        if (detr_star(ig,l).gt.f_star(ig,l)) then
-
-!  Ce test est là entre autres parce qu'on passe par des valeurs
-!  delirantes de detr_star.
-!  ca vaut sans doute le coup de verifier pourquoi.
-
-           detr_star(ig,l)=f_star(ig,l)
-#ifdef pre_int1c
-           if (l.gt.lalim(ig)+1) then
-               entr_star(ig,l)=0.
-               alim_star(ig,l)=0.
-! FH ajout pour forcer a stoper le thermique juste sous le sommet
-! de la couche (voir calcul de finter)
-               zw2(ig,l+1)=-1.e-10
-               linter(ig)=l+1
-            else
-               entr_star(ig,l)=0.4*detr_star(ig,l)
-            endif
-#else
-           entr_star(ig,l)=0.4*detr_star(ig,l)
-#endif
-        endif
-#endif
-
-      else !l > 2
-         detr_star(ig,l)=0.
-         entr_star(ig,l)=0.
-      endif
-
-        entr_star2(ig,l)=entr_star(ig,l)
-        detr_star2(ig,l)=detr_star(ig,l)
-        if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
-
-       endif  ! iflag_thermals_ed==0
-
-!CR:nvlle def de entr_star et detr_star
-      if (iflag_thermals_ed>=1) then
-!      if (l.lt.lalim(ig)) then
-!      if (l.lt.2) then 
-!        entr_star(ig,l)=0.
-!        detr_star(ig,l)=0.
-!      else
-!      if (0.001.gt.(RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))/(2.*w_est(ig,l+1)))) then 
-!         entr_star(ig,l)=0.001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
-!      else
-!         entr_star(ig,l)=  &
-!     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
-!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
-!     &                *(zlev(ig,l+1)-zlev(ig,l))
-
- 
-         entr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &          
-     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
-     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)   &
-     &                *(zlev(ig,l+1)-zlev(ig,l))) &
-     &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
-
-        if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
-            alim_star_tot(ig)=alim_star_tot(ig)+entr_star(ig,l)
-            lalim(ig)=lmix_bis(ig)
-            if(prt_level.GE.10) print*,'alim_star_tot',alim_star_tot(ig),entr_star(ig,l)
-        endif
-
-        if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
-!        c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
-         c2(ig,l)=0.001
-         detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
-     &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
-     &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
-     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
-     &                *(zlev(ig,l+1)-zlev(ig,l)))                    &
-     &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
-
-       else
-!         c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
-          c2(ig,l)=0.003
-
-         detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
-     &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
-     &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
-     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
-     &                *(zlev(ig,l+1)-zlev(ig,l))) &
-     &                +0.0002*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
-       endif
-         
-           
-!        detr_star(ig,l)=detr_star(ig,l)*3.
-!        if (l.lt.lalim(ig)) then
-!          entr_star(ig,l)=0.
-!        endif
-!        if (l.lt.2) then
-!          entr_star(ig,l)=0.
-!          detr_star(ig,l)=0.
-!        endif
-
-
-!      endif 
-!      else if ((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10) then
-!      entr_star(ig,l)=MAX(0.,0.8*f_star(ig,l)/(2.*w_est(ig,l+1))        &
-!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
-!     &                *(zlev(ig,l+1)-zlev(ig,l))
-!      detr_star(ig,l)=0.002*f_star(ig,l)                         &
-!     &                *(zlev(ig,l+1)-zlev(ig,l))
-!      else
-!      entr_star(ig,l)=0.001*f_star(ig,l)                         &
-!     &                *(zlev(ig,l+1)-zlev(ig,l))
-!      detr_star(ig,l)=MAX(0.,-0.2*f_star(ig,l)/(2.*w_est(ig,l+1))       &
-!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))       &
-!     &                *(zlev(ig,l+1)-zlev(ig,l))                      &
-!     &                +0.002*f_star(ig,l)                             &
-!     &                *(zlev(ig,l+1)-zlev(ig,l))
-!      endif
-
-      endif   ! iflag_thermals_ed==1
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! FH inutile si on conserve comme on l'a fait plus haut entr=detr
-! dans la couche d'alimentation
-!pas d entrainement dans la couche alim
-!      if ((l.le.lalim(ig))) then
-!           entr_star(ig,l)=0.
-!      endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-!prise en compte du detrainement et de l entrainement dans le calcul du flux
-
+! Calcul du flux montant normalise
       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
      &              -detr_star(ig,l)
 
-!test sur le signe de f_star
-        if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
-       if (f_star(ig,l+1).gt.1.e-10) then 
+      endif
+   enddo
+
 !----------------------------------------------------------------------------
 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
 !---------------------------------------------------------------------------
-!
-       Zsat=.false.
-       ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
+   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
+   do ig=1,ngrid
+       if (activetmp(ig)) then 
+           Zsat=.false.
+           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
      &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
      &            /(f_star(ig,l+1)+detr_star(ig,l))
-!
-       zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
+           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
      &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
      &            /(f_star(ig,l+1)+detr_star(ig,l))
-!  
-               Tbef=ztla(ig,l)*zpspsk(ig,l)
-               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)               
-               qsatbef=MIN(0.5,qsatbef)
-               zcor=1./(1.-retv*qsatbef)
-               qsatbef=qsatbef*zcor
-               Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
-               if (Zsat) then
-               qlbef=max(0.,zqta(ig,l)-qsatbef)
-               DT = 0.5*RLvCp*qlbef
-               do while (abs(DT).gt.DDT0)
-                 Tbef=Tbef+DT
-                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
-                 qsatbef=MIN(0.5,qsatbef)
-                 zcor=1./(1.-retv*qsatbef)
-                 qsatbef=qsatbef*zcor
-                 qlbef=zqta(ig,l)-qsatbef
-
-                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
-                 zcor=1./(1.-retv*qsatbef)
-                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
-                 num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
-                 denom=1.+RLvCp*dqsat_dT
-                 DT=num/denom
-              enddo
-                 zqla(ig,l) = max(0.,qlbef) 
-              endif
-!    
+
+        endif
+    enddo
+
+   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
+
+
+   do ig=1,ngrid
+      if (activetmp(ig)) then
         if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
 ! on ecrit de maniere conservative (sat ou non)
@@ -707,43 +383,55 @@
 !on ecrit zqsat 
            zqsatth(ig,l)=qsatbef  
-!calcul de vitesse
-           zw2(ig,l+1)=zw2(ig,l)*  &
-     &                 ((f_star(ig,l))**2)  &
-!  Tests de Catherine
-!     &                 /(f_star(ig,l+1)+detr_star(ig,l))**2+             &
-     &      /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-fact_epsilon))**2+ &
-     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
-     &                 *fact_gamma &
-     &                 *(zlev(ig,l+1)-zlev(ig,l))
-!prise en compte des forces de pression que qd flottabilité<0
-!              zw2(ig,l+1)=zw2(ig,l)*  &
-!     &            1./(1.+2.*entr_star(ig,l)/f_star(ig,l)) + &        
-!     &                 (f_star(ig,l))**2 &
-!     &                 /(f_star(ig,l)+entr_star(ig,l))**2+ &
-!     &                 (f_star(ig,l)-2.*entr_star(ig,l))**2/(f_star(ig,l)+2.*entr_star(ig,l))**2+  &        
-!     &                 /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-2.))**2+ &
-!     &                 /(f_star(ig,l)**2+2.*2.*detr_star(ig,l)*f_star(ig,l)+2.*entr_star(ig,l)*f_star(ig,l))+ &
-!     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
-!     &                 *1./3. &
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!          zw2(ig,l+1)=&
+!     &                 zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) &
+!     &                 +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
+!     &                 *1./(1.+fact_gamma) &
 !     &                 *(zlev(ig,l+1)-zlev(ig,l))
-          
-!        write(30,*),l+1,zw2(ig,l+1)-zw2(ig,l), &
-!     &              -2.*entr_star(ig,l)/f_star(ig,l)*zw2(ig,l), &
-!     &               2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
-
- 
-!             zw2(ig,l+1)=zw2(ig,l)*  &
-!     &                 (2.-2.*entr_star(ig,l)/f_star(ig,l)) &  
-!     &                 -zw2(ig,l-1)+  &        
-!     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
-!     &                 *1./3. &
-!     &                 *(zlev(ig,l+1)-zlev(ig,l))             
-
-            endif
-        endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! La meme en plus modulaire :
+           zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
+           zdz=zlev(ig,l+1)-zlev(ig,l)
+
+
+           zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
+
+if (1==0) then
+           zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz)
+           zdw2=2.*zbuoy/(1.+fact_gamma)*zdz
+           zw2(ig,l+1)=zw2modif+zdw2
+else
+! Tentative de reecriture de l'equation de w2. A reprendre ...
+!           zdw2=2*zdz*zbuoy
+!           zw2modif=zw2(ig,l)*(1.-2.*zdz*(zeps+fact_epsilon))
+!!!!!      zw2(ig,l+1)=(zw2(ig,l)+zdw2)/(1.+2.*zfactw2(ig,l))
+!formulation Arnaud
+!           zw2fact=zbuoy*zalpha**expa/fact_epsilon
+!           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa*(1+fact_gamma))*zdz) &
+!      &    +zw2fact
+           if (zbuoy.gt.1.e-10) then
+           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma-zfact)
+           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
+      &    +zw2fact
+           else
+           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma)
+           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
+      &    +zw2fact
+
+           endif
+
+endif
+!           zw2(ig,l+1)=zw2modif+zdw2
+      endif
+   enddo
+
         if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
 !
+!---------------------------------------------------------------------------
 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max 
-
+!---------------------------------------------------------------------------
+
+   do ig=1,ngrid
             if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
 !               stop'On tombe sur le cas particulier de thermcell_dry'
@@ -753,5 +441,4 @@
             endif
 
-!        if ((zw2(ig,l).gt.0.).and. (zw2(ig,l+1).le.0.)) then
         if (zw2(ig,l+1).lt.0.) then 
            linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
@@ -771,28 +458,22 @@
             wmaxa(ig)=wa_moy(ig,l+1)
         endif
-        enddo
+   enddo
+
+!=========================================================================
+! FIN DE LA BOUCLE VERTICALE
       enddo
-
-!on remplace a* par e* ds premiere couche
-!      if (iflag_thermals_ed.ge.1) then
-!       do ig=1,ngrid
-!       do l=2,klev
-!          if (l.lt.lalim(ig)) then
-!             alim_star(ig,l)=entr_star(ig,l)
-!          endif
-!       enddo
-!       enddo
-!       do ig=1,ngrid
-!          lalim(ig)=lmix_bis(ig)
-!       enddo
-!      endif
-       if (iflag_thermals_ed.ge.1) then
-          do ig=1,ngrid
-             do l=2,lalim(ig)
-                alim_star(ig,l)=entr_star(ig,l)
-                entr_star(ig,l)=0.
-             enddo
-           enddo
-       endif
+!=========================================================================
+
+!on recalcule alim_star_tot
+       do ig=1,ngrid
+          alim_star_tot(ig)=0.
+       enddo
+       do ig=1,ngrid
+          do l=1,lalim(ig)-1
+          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
+          enddo
+       enddo
+       
+
         if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
 
