Index: LMDZ6/branches/Amaury_dev/libf/misc/lmdz_cppkeys_wrapper.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/misc/lmdz_cppkeys_wrapper.F90	(revision 5101)
+++ LMDZ6/branches/Amaury_dev/libf/misc/lmdz_cppkeys_wrapper.F90	(revision 5102)
@@ -5,9 +5,10 @@
 ! and imported through USE ..., ONLY: ... elsewhere
 ! CPP keys supported (key -> fortran variables associated):
-!      NC_DOUBLE -> nf90_format
-!      CPP_PHYS  -> CPPKEY_PHYS
-!      INCA      -> CPPKEY_INCA
-!      CPP_StratAer-> CPPKEY_STRATAER
-!      CPP_DUST  -> CPPKEY_DUST
+!      NC_DOUBLE      -> nf90_format
+!      CPP_PHYS       -> CPPKEY_PHYS
+!      INCA           -> CPPKEY_INCA
+!      CPP_StratAer   -> CPPKEY_STRATAER
+!      CPP_DUST       -> CPPKEY_DUST
+!      CPP_INLANDSIS  -> CPPKEY_INLANDSIS
 ! ---------------------------------------------
 
@@ -17,5 +18,5 @@
   IMPLICIT NONE; PRIVATE
   PUBLIC nf90_format, CPPKEY_PHYS, CPPKEY_INCA, CPPKEY_STRATAER, CPPKEY_DUST, &
-      CPPKEY_DEBUGIO
+      CPPKEY_DEBUGIO, CPPKEY_INLANDSIS
 
 #ifdef NC_DOUBLE
@@ -55,3 +56,9 @@
 #endif
 
+#ifdef CPP_INLANDSIS
+  LOGICAL, PARAMETER :: CPPKEY_INLANDSIS = .TRUE.
+#else
+  LOGICAL, PARAMETER :: CPPKEY_INLANDSIS = .FALSE.
+#endif
+
 END MODULE lmdz_cppkeys_wrapper
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume.F90	(revision 5101)
+++ 	(revision )
@@ -1,453 +1,0 @@
-MODULE lmdz_thermcell_plume
-
-! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $
-
-CONTAINS
-
-      SUBROUTINE thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
-             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)
-!     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
-!--------------------------------------------------------------------------
-! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam
-
-!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
-!   This versions starts from a cleaning of thermcell_plume_6A (2019/01/20)
-!   thermcell_plume_6A is activate for flag_thermas_ed < 10
-!   thermcell_plume_5B for flag_thermas_ed < 20
-!   thermcell_plume for flag_thermals_ed>= 20
-!   Various options are controled by the flag_thermals_ed parameter
-!   = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8
-!   = 21 : the Jam strato-cumulus modif is not activated in detrainment
-!   = 29 : an other way to compute the modified buoyancy (to be tested)
-!--------------------------------------------------------------------------
-
-       USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
-       USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
-       USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
-       USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
-       USE lmdz_thermcell_alim, ONLY: thermcell_alim
-       USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
-
-
-       IMPLICIT NONE
-
-      integer,intent(in) :: itap,lev_out,lunout1,igout,ngrid,nlay
-      real,intent(in) :: ptimestep
-      real,intent(in),dimension(ngrid,nlay) :: ztv
-      real,intent(in),dimension(ngrid,nlay) :: zthl
-      real,intent(in),dimension(ngrid,nlay) :: po
-      real,intent(in),dimension(ngrid,nlay) :: zl
-      real,intent(in),dimension(ngrid,nlay) :: rhobarz
-      real,intent(in),dimension(ngrid,nlay+1) :: zlev
-      real,intent(in),dimension(ngrid,nlay+1) :: pplev
-      real,intent(in),dimension(ngrid,nlay) :: pphi
-      real,intent(in),dimension(ngrid,nlay) :: zpspsk
-      real,intent(in),dimension(ngrid) :: f0
-
-      integer,intent(out) :: lalim(ngrid)
-      real,intent(out),dimension(ngrid,nlay) :: alim_star
-      real,intent(out),dimension(ngrid) :: alim_star_tot
-      real,intent(out),dimension(ngrid,nlay) :: detr_star
-      real,intent(out),dimension(ngrid,nlay) :: entr_star
-      real,intent(out),dimension(ngrid,nlay+1) :: f_star
-      real,intent(out),dimension(ngrid,nlay) :: csc
-      real,intent(out),dimension(ngrid,nlay) :: ztva
-      real,intent(out),dimension(ngrid,nlay) :: ztla
-      real,intent(out),dimension(ngrid,nlay) :: zqla
-      real,intent(out),dimension(ngrid,nlay) :: zqta
-      real,intent(out),dimension(ngrid,nlay) :: zha
-      real,intent(out),dimension(ngrid,nlay+1) :: zw2
-      real,intent(out),dimension(ngrid,nlay+1) :: w_est
-      real,intent(out),dimension(ngrid,nlay) :: ztva_est
-      real,intent(out),dimension(ngrid,nlay) :: zqsatth
-      integer,intent(out),dimension(ngrid) :: lmix(ngrid)
-      integer,intent(out),dimension(ngrid) :: lmix_bis(ngrid)
-      real,intent(out),dimension(ngrid) :: linter(ngrid)
-
-
-      REAL,dimension(ngrid,nlay+1) :: wa_moy
-      REAL,dimension(ngrid,nlay) :: entr,detr
-      REAL,dimension(ngrid,nlay) :: ztv_est
-      REAL,dimension(ngrid,nlay) :: zqla_est
-      REAL,dimension(ngrid,nlay) :: zta_est
-      REAL,dimension(ngrid) :: ztemp,zqsat
-      REAL zdw2,zdw2bis
-      REAL zw2modif
-      REAL zw2fact,zw2factbis
-      REAL,dimension(ngrid,nlay) :: zeps
-
-      REAL,dimension(ngrid) :: wmaxa
-
-      INTEGER ig,l,k,lt,it,lm,nbpb
-
-      real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt
-      real zdz,zalpha,zw2m
-      real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam
-      real zdz2,zdz3,lmel,entrbis,zdzbis
-      real,dimension(ngrid) :: d_temp
-      real ztv1,ztv2,factinv,zinv,zlmel
-      real zlmelup,zlmeldwn,zlt,zltdwn,zltup
-      real atv1,atv2,btv1,btv2
-      real ztv_est1,ztv_est2
-      real zcor,zdelta,zcvm5,qlbef
-      real zbetalpha, coefzlmel
-      real eps
-      logical Zsat
-      LOGICAL,dimension(ngrid) :: active,activetmp
-      REAL fact_gamma,fact_gamma2,fact_epsilon2
-
-
-      REAL,dimension(ngrid,nlay) :: c2
-
-      if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
-      Zsat=.false.
-! Initialisation
-
-
-      zbetalpha=betalpha/(1.+betalpha)
-
-
-! Initialisations des variables r?elles
-if (1==1) then
-      ztva(:,:)=ztv(:,:)
-      ztva_est(:,:)=ztva(:,:)
-      ztv_est(:,:)=ztv(:,:)
-      ztla(:,:)=zthl(:,:)
-      zqta(:,:)=po(:,:)
-      zqla(:,:)=0.
-      zha(:,:) = ztva(:,:)
-else
-      ztva(:,:)=0.
-      ztv_est(:,:)=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.
-      zbuoy(:,:)=0.
-      zbuoyjam(:,:)=0.
-      gamma(:,:)=0.
-      zeps(:,:)=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.
-
-
-!-------------------------------------------------------------------------
-! On ne considere comme actif que les colonnes dont les deux premieres
-! couches sont instables.
-!-------------------------------------------------------------------------
-
-      active(:)=ztv(:,1)>ztv(:,2)
-      d_temp(:)=0. ! Pour activer un contraste de temperature a la base
-                   ! du panache
-!  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
-      CALL thermcell_alim(thermals_flag_alim,ngrid,nlay,ztv,d_temp,zlev,alim_star,lalim)
-
-!------------------------------------------------------------------------------
-! 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.
-    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
-
-!==============================================================================
-!boucle de calcul de la vitesse verticale dans le thermique
-!==============================================================================
-do l=2,nlay-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
-
-
-
-!---------------------------------------------------------------------------
-! calcul des proprietes thermodynamiques et de la vitesse de la couche l
-! sans tenir compte du detrainement et de l'entrainement dans cette
-! couche
-! C'est a dire qu'on suppose 
-! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
-! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
-! avant) a l'alimentation pour avoir un calcul plus propre
-!---------------------------------------------------------------------------
-
-   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
-   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
-    do ig=1,ngrid 
-!       print*,'active',active(ig),ig,l
-        if(active(ig)) then 
-        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
-        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)
-        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
-        -zqla_est(ig,l))-zqla_est(ig,l))
- 
-
-!Modif AJAM
-
-        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) 
-        zdz=zlev(ig,l+1)-zlev(ig,l)         
-        lmel=fact_thermals_ed_dz*zlev(ig,l)
-!        lmel=0.09*zlev(ig,l)
-        zlmel=zlev(ig,l)+lmel
-        zlmelup=zlmel+(zdz/2)
-        zlmeldwn=zlmel-(zdz/2)
-
-        lt=l+1
-        zlt=zlev(ig,lt)
-        zdz3=zlev(ig,lt+1)-zlt
-        zltdwn=zlt-zdz3/2
-        zltup=zlt+zdz3/2
-         
-!=========================================================================
-! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
-!=========================================================================
-
-!--------------------------------------------------
-        lt=l+1
-        zlt=zlev(ig,lt)
-        zdz2=zlev(ig,lt)-zlev(ig,l)
-
-        do while (lmel>zdz2)
-           lt=lt+1
-           zlt=zlev(ig,lt)
-           zdz2=zlt-zlev(ig,l)
-        enddo
-        zdz3=zlev(ig,lt+1)-zlt
-        zltdwn=zlev(ig,lt)-zdz3/2
-        zlmelup=zlmel+(zdz/2)
-        coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
-        zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- & 
-     ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
-     ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
-
-!------------------------------------------------
-!AJAM:nouveau calcul de w?  
-!------------------------------------------------
-        zdz=zlev(ig,l+1)-zlev(ig,l)
-        zdzbis=zlev(ig,l)-zlev(ig,l-1)
-        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-        zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
-        zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
-        zdw2=afact*zbuoy(ig,l)/fact_epsilon
-        zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
-        lm=Max(1,l-2)
-        w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
-       endif
-    enddo
-
-
-!-------------------------------------------------
-!calcul des taux d'entrainement et de detrainement
-!-------------------------------------------------
-
-     do ig=1,ngrid
-        if (active(ig)) then
-
-!          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
-          zw2m=w_est(ig,l+1)
-          zdz=zlev(ig,l+1)-zlev(ig,l)
-          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
-          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
-
-!=========================================================================
-! 4. Calcul de l'entrainement et du detrainement
-!=========================================================================
-
-          detr_star(ig,l)=f_star(ig,l)*zdz             &
-       *( mix0 * 0.1 / (zalpha+0.001)               &
-       + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
-       + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
-
-          if ( iflag_thermals_ed == 20 ) then
-             entr_star(ig,l)=f_star(ig,l)*zdz* (         &
-            mix0 * 0.1 / (zalpha+0.001)               &
-          + zbetalpha*MAX(entr_min,                   &
-          afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
-          else
-             entr_star(ig,l)=f_star(ig,l)*zdz* (         &
-            mix0 * 0.1 / (zalpha+0.001)               &
-          + zbetalpha*MAX(entr_min,                   &
-          afact*zbuoy(ig,l)/zw2m - fact_epsilon))
-          endif
-          
-! En dessous de lalim, on prend le max de alim_star et entr_star pour
-! alim_star et 0 sinon
-        if (l<lalim(ig)) then
-          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
-          entr_star(ig,l)=0.
-        endif
-        f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
-                -detr_star(ig,l)
-
-      endif
-   enddo
-
-
-!============================================================================
-! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
-!===========================================================================
-
-   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)+  &
-              (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
-              /(f_star(ig,l+1)+detr_star(ig,l))
-
-        endif
-    enddo
-
-   ztemp(:)=zpspsk(:,l)*ztla(:,l)
-   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
-   do ig=1,ngrid
-      if (activetmp(ig)) then
-! on ecrit de maniere conservative (sat ou non)
-!          T = Tl +Lv/Cp ql
-           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
-           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
-           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
-!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
-           zha(ig,l) = ztva(ig,l)
-           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
-                -zqla(ig,l))-zqla(ig,l))
-           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
-           zdz=zlev(ig,l+1)-zlev(ig,l)
-           zdzbis=zlev(ig,l)-zlev(ig,l-1)
-           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
-           zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
-           zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
-           zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
-           zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
-           zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
-      endif
-   enddo
-
-   if (prt_level>=20) print*,'coucou calcul detr 460: ig, l',ig, l
-
-!===========================================================================
-! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max 
-!===========================================================================
-
-   nbpb=0
-   do ig=1,ngrid
-            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1)<1.e-10) then
-!               stop'On tombe sur le cas particulier de thermcell_dry'
-!               print*,'On tombe sur le cas particulier de thermcell_plume'
-                nbpb=nbpb+1
-                zw2(ig,l+1)=0.
-                linter(ig)=l+1
-            endif
-
-        if (zw2(ig,l+1)<0.) then
-           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
-                 -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
-           zw2(ig,l+1)=0.
-!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
-        elseif (f_star(ig,l+1)<0.) then
-           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
-                 -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
-           zw2(ig,l+1)=0.
-!fin CR:04/05/12
-        endif
-
-           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 
-
-        if (wa_moy(ig,l+1)>wmaxa(ig)) then
-!   lmix est le niveau de la couche ou w (wa_moy) est maximum
-!on rajoute le calcul de lmix_bis
-            if (zqla(ig,l)<1.e-10) then
-               lmix_bis(ig)=l+1
-            endif
-            lmix(ig)=l+1
-            wmaxa(ig)=wa_moy(ig,l+1)
-        endif
-   enddo
-
-   if (nbpb>0) then
-   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
-   endif
-
-!=========================================================================
-! FIN DE LA BOUCLE VERTICALE
-      enddo
-!=========================================================================
-
-!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>=20) print*,'coucou calcul detr 470: ig, l', ig, l
-
-#undef wrgrads_thermcell
-#ifdef wrgrads_thermcell
-         call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
-         call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
-         call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
-         call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
-         call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
-         call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
-         call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
-#endif
-
-
- RETURN
-     end
-END MODULE lmdz_thermcell_plume
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume.f90	(revision 5102)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume.f90	(revision 5102)
@@ -0,0 +1,435 @@
+MODULE lmdz_thermcell_plume
+
+  ! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $
+
+CONTAINS
+
+  SUBROUTINE thermcell_plume(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
+          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)
+    !     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
+    !--------------------------------------------------------------------------
+    ! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam
+
+    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
+    !   This versions starts from a cleaning of thermcell_plume_6A (2019/01/20)
+    !   thermcell_plume_6A is activate for flag_thermas_ed < 10
+    !   thermcell_plume_5B for flag_thermas_ed < 20
+    !   thermcell_plume for flag_thermals_ed>= 20
+    !   Various options are controled by the flag_thermals_ed parameter
+    !   = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8
+    !   = 21 : the Jam strato-cumulus modif is not activated in detrainment
+    !   = 29 : an other way to compute the modified buoyancy (to be tested)
+    !--------------------------------------------------------------------------
+
+    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
+    USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
+    USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
+    USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
+    USE lmdz_thermcell_alim, ONLY: thermcell_alim
+    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
+
+    IMPLICIT NONE
+
+    integer, intent(in) :: itap, lev_out, lunout1, igout, ngrid, nlay
+    real, intent(in) :: ptimestep
+    real, intent(in), dimension(ngrid, nlay) :: ztv
+    real, intent(in), dimension(ngrid, nlay) :: zthl
+    real, intent(in), dimension(ngrid, nlay) :: po
+    real, intent(in), dimension(ngrid, nlay) :: zl
+    real, intent(in), dimension(ngrid, nlay) :: rhobarz
+    real, intent(in), dimension(ngrid, nlay + 1) :: zlev
+    real, intent(in), dimension(ngrid, nlay + 1) :: pplev
+    real, intent(in), dimension(ngrid, nlay) :: pphi
+    real, intent(in), dimension(ngrid, nlay) :: zpspsk
+    real, intent(in), dimension(ngrid) :: f0
+
+    integer, intent(out) :: lalim(ngrid)
+    real, intent(out), dimension(ngrid, nlay) :: alim_star
+    real, intent(out), dimension(ngrid) :: alim_star_tot
+    real, intent(out), dimension(ngrid, nlay) :: detr_star
+    real, intent(out), dimension(ngrid, nlay) :: entr_star
+    real, intent(out), dimension(ngrid, nlay + 1) :: f_star
+    real, intent(out), dimension(ngrid, nlay) :: csc
+    real, intent(out), dimension(ngrid, nlay) :: ztva
+    real, intent(out), dimension(ngrid, nlay) :: ztla
+    real, intent(out), dimension(ngrid, nlay) :: zqla
+    real, intent(out), dimension(ngrid, nlay) :: zqta
+    real, intent(out), dimension(ngrid, nlay) :: zha
+    real, intent(out), dimension(ngrid, nlay + 1) :: zw2
+    real, intent(out), dimension(ngrid, nlay + 1) :: w_est
+    real, intent(out), dimension(ngrid, nlay) :: ztva_est
+    real, intent(out), dimension(ngrid, nlay) :: zqsatth
+    integer, intent(out), dimension(ngrid) :: lmix(ngrid)
+    integer, intent(out), dimension(ngrid) :: lmix_bis(ngrid)
+    real, intent(out), dimension(ngrid) :: linter(ngrid)
+
+    REAL, dimension(ngrid, nlay + 1) :: wa_moy
+    REAL, dimension(ngrid, nlay) :: entr, detr
+    REAL, dimension(ngrid, nlay) :: ztv_est
+    REAL, dimension(ngrid, nlay) :: zqla_est
+    REAL, dimension(ngrid, nlay) :: zta_est
+    REAL, dimension(ngrid) :: ztemp, zqsat
+    REAL zdw2, zdw2bis
+    REAL zw2modif
+    REAL zw2fact, zw2factbis
+    REAL, dimension(ngrid, nlay) :: zeps
+
+    REAL, dimension(ngrid) :: wmaxa
+
+    INTEGER ig, l, k, lt, it, lm, nbpb
+
+    real, dimension(ngrid, nlay) :: zbuoy, gamma, zdqt
+    real zdz, zalpha, zw2m
+    real, dimension(ngrid, nlay) :: zbuoyjam, zdqtjam
+    real zdz2, zdz3, lmel, entrbis, zdzbis
+    real, dimension(ngrid) :: d_temp
+    real ztv1, ztv2, factinv, zinv, zlmel
+    real zlmelup, zlmeldwn, zlt, zltdwn, zltup
+    real atv1, atv2, btv1, btv2
+    real ztv_est1, ztv_est2
+    real zcor, zdelta, zcvm5, qlbef
+    real zbetalpha, coefzlmel
+    real eps
+    logical Zsat
+    LOGICAL, dimension(ngrid) :: active, activetmp
+    REAL fact_gamma, fact_gamma2, fact_epsilon2
+
+    REAL, dimension(ngrid, nlay) :: c2
+
+    if (ngrid==1) print*, 'THERMCELL PLUME MODIFIE 2014/07/11'
+    Zsat = .false.
+    ! Initialisation
+
+    zbetalpha = betalpha / (1. + betalpha)
+
+
+    ! Initialisations des variables r?elles
+    if (1==1) then
+      ztva(:, :) = ztv(:, :)
+      ztva_est(:, :) = ztva(:, :)
+      ztv_est(:, :) = ztv(:, :)
+      ztla(:, :) = zthl(:, :)
+      zqta(:, :) = po(:, :)
+      zqla(:, :) = 0.
+      zha(:, :) = ztva(:, :)
+    else
+      ztva(:, :) = 0.
+      ztv_est(:, :) = 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.
+    zbuoy(:, :) = 0.
+    zbuoyjam(:, :) = 0.
+    gamma(:, :) = 0.
+    zeps(:, :) = 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.
+
+
+    !-------------------------------------------------------------------------
+    ! On ne considere comme actif que les colonnes dont les deux premieres
+    ! couches sont instables.
+    !-------------------------------------------------------------------------
+
+    active(:) = ztv(:, 1)>ztv(:, 2)
+    d_temp(:) = 0. ! Pour activer un contraste de temperature a la base
+    ! du panache
+    !  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
+    CALL thermcell_alim(thermals_flag_alim, ngrid, nlay, ztv, d_temp, zlev, alim_star, lalim)
+
+    !------------------------------------------------------------------------------
+    ! 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.
+      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
+
+    !==============================================================================
+    !boucle de calcul de la vitesse verticale dans le thermique
+    !==============================================================================
+    do l = 2, nlay - 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
+
+
+
+      !---------------------------------------------------------------------------
+      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
+      ! sans tenir compte du detrainement et de l'entrainement dans cette
+      ! couche
+      ! C'est a dire qu'on suppose
+      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
+      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
+      ! avant) a l'alimentation pour avoir un calcul plus propre
+      !---------------------------------------------------------------------------
+
+      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
+      call thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
+      do ig = 1, ngrid
+        !       print*,'active',active(ig),ig,l
+        if(active(ig)) then
+          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
+          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)
+          ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1)  &
+                  - zqla_est(ig, l)) - zqla_est(ig, l))
+
+
+          !Modif AJAM
+
+          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          lmel = fact_thermals_ed_dz * zlev(ig, l)
+          !        lmel=0.09*zlev(ig,l)
+          zlmel = zlev(ig, l) + lmel
+          zlmelup = zlmel + (zdz / 2)
+          zlmeldwn = zlmel - (zdz / 2)
+
+          lt = l + 1
+          zlt = zlev(ig, lt)
+          zdz3 = zlev(ig, lt + 1) - zlt
+          zltdwn = zlt - zdz3 / 2
+          zltup = zlt + zdz3 / 2
+
+          !=========================================================================
+          ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
+          !=========================================================================
+
+          !--------------------------------------------------
+          lt = l + 1
+          zlt = zlev(ig, lt)
+          zdz2 = zlev(ig, lt) - zlev(ig, l)
+
+          do while (lmel>zdz2)
+            lt = lt + 1
+            zlt = zlev(ig, lt)
+            zdz2 = zlt - zlev(ig, l)
+          enddo
+          zdz3 = zlev(ig, lt + 1) - zlt
+          zltdwn = zlev(ig, lt) - zdz3 / 2
+          zlmelup = zlmel + (zdz / 2)
+          coefzlmel = Min(1., (zlmelup - zltdwn) / zdz)
+          zbuoyjam(ig, l) = 1. * RG * (coefzlmel * (ztva_est(ig, l) - &
+                  ztv(ig, lt)) / ztv(ig, lt) + (1. - coefzlmel) * (ztva_est(ig, l) - &
+                  ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + 0. * zbuoy(ig, l)
+
+          !------------------------------------------------
+          !AJAM:nouveau calcul de w?
+          !------------------------------------------------
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
+          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
+          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
+          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
+          zdw2 = afact * zbuoy(ig, l) / fact_epsilon
+          zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
+          lm = Max(1, l - 2)
+          w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
+        endif
+      enddo
+
+
+      !-------------------------------------------------
+      !calcul des taux d'entrainement et de detrainement
+      !-------------------------------------------------
+
+      do ig = 1, ngrid
+        if (active(ig)) then
+
+          !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
+          zw2m = w_est(ig, l + 1)
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
+          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
+          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)
+
+          !=========================================================================
+          ! 4. Calcul de l'entrainement et du detrainement
+          !=========================================================================
+
+          detr_star(ig, l) = f_star(ig, l) * zdz             &
+                  * (mix0 * 0.1 / (zalpha + 0.001)               &
+                          + MAX(detr_min, -afact * zbetalpha * zbuoyjam(ig, l) / zw2m   &
+                                  + detr_q_coef * (zdqt(ig, l) / zw2m)**detr_q_power))
+
+          if (iflag_thermals_ed == 20) then
+            entr_star(ig, l) = f_star(ig, l) * zdz * (&
+                    mix0 * 0.1 / (zalpha + 0.001)               &
+                            + zbetalpha * MAX(entr_min, &
+                            afact * zbuoyjam(ig, l) / zw2m - fact_epsilon))
+          else
+            entr_star(ig, l) = f_star(ig, l) * zdz * (&
+                    mix0 * 0.1 / (zalpha + 0.001)               &
+                            + zbetalpha * MAX(entr_min, &
+                            afact * zbuoy(ig, l) / zw2m - fact_epsilon))
+          endif
+
+          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
+          ! alim_star et 0 sinon
+          if (l<lalim(ig)) then
+            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
+            entr_star(ig, l) = 0.
+          endif
+          f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l)  &
+                  - detr_star(ig, l)
+
+        endif
+      enddo
+
+
+      !============================================================================
+      ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
+      !===========================================================================
+
+      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) + &
+                  (alim_star(ig, l) + entr_star(ig, l)) * po(ig, l))  &
+                  / (f_star(ig, l + 1) + detr_star(ig, l))
+
+        endif
+      enddo
+
+      ztemp(:) = zpspsk(:, l) * ztla(:, l)
+      call thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
+      do ig = 1, ngrid
+        if (activetmp(ig)) then
+          ! on ecrit de maniere conservative (sat ou non)
+          !          T = Tl +Lv/Cp ql
+          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
+          ztva(ig, l) = ztla(ig, l) * zpspsk(ig, l) + RLvCp * zqla(ig, l)
+          ztva(ig, l) = ztva(ig, l) / zpspsk(ig, l)
+          !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
+          zha(ig, l) = ztva(ig, l)
+          ztva(ig, l) = ztva(ig, l) * (1. + RETV * (zqta(ig, l)  &
+                  - zqla(ig, l)) - zqla(ig, l))
+          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
+          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
+          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
+          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
+          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
+          zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
+          zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
+        endif
+      enddo
+
+      if (prt_level>=20) print*, 'coucou calcul detr 460: ig, l', ig, l
+
+      !===========================================================================
+      ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
+      !===========================================================================
+
+      nbpb = 0
+      do ig = 1, ngrid
+        if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then
+          !               stop'On tombe sur le cas particulier de thermcell_dry'
+          !               print*,'On tombe sur le cas particulier de thermcell_plume'
+          nbpb = nbpb + 1
+          zw2(ig, l + 1) = 0.
+          linter(ig) = l + 1
+        endif
+
+        if (zw2(ig, l + 1)<0.) then
+          linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l))  &
+                  - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l))
+          zw2(ig, l + 1) = 0.
+          !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
+        elseif (f_star(ig, l + 1)<0.) then
+          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
+                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
+          zw2(ig, l + 1) = 0.
+          !fin CR:04/05/12
+        endif
+
+        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))
+
+        if (wa_moy(ig, l + 1)>wmaxa(ig)) then
+          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
+          !on rajoute le calcul de lmix_bis
+          if (zqla(ig, l)<1.e-10) then
+            lmix_bis(ig) = l + 1
+          endif
+          lmix(ig) = l + 1
+          wmaxa(ig) = wa_moy(ig, l + 1)
+        endif
+      enddo
+
+      if (nbpb>0) then
+        print*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
+      endif
+
+      !=========================================================================
+      ! FIN DE LA BOUCLE VERTICALE
+    enddo
+    !=========================================================================
+
+    !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>=20) print*, 'coucou calcul detr 470: ig, l', ig, l
+    RETURN
+  end
+END MODULE lmdz_thermcell_plume
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume_6A.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume_6A.F90	(revision 5101)
+++ 	(revision )
@@ -1,1116 +1,0 @@
-MODULE lmdz_thermcell_plume_6A
-
-! $Id$
-
-CONTAINS
-
-      SUBROUTINE thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
-             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)
-!     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
-!--------------------------------------------------------------------------
-!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
-!--------------------------------------------------------------------------
-
-       USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
-       USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
-       USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
-       USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
-       USE lmdz_thermcell_alim, ONLY: thermcell_alim
-       USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
-
-
-       IMPLICIT NONE
-
-      integer,intent(in) :: itap,lev_out,lunout1,igout,ngrid,nlay
-      real,intent(in) :: ptimestep
-      real,intent(in),dimension(ngrid,nlay) :: ztv
-      real,intent(in),dimension(ngrid,nlay) :: zthl
-      real,intent(in),dimension(ngrid,nlay) :: po
-      real,intent(in),dimension(ngrid,nlay) :: zl
-      real,intent(in),dimension(ngrid,nlay) :: rhobarz
-      real,intent(in),dimension(ngrid,nlay+1) :: zlev
-      real,intent(in),dimension(ngrid,nlay+1) :: pplev
-      real,intent(in),dimension(ngrid,nlay) :: pphi
-      real,intent(in),dimension(ngrid,nlay) :: zpspsk
-      real,intent(in),dimension(ngrid) :: f0
-
-      integer,intent(out) :: lalim(ngrid)
-      real,intent(out),dimension(ngrid,nlay) :: alim_star
-      real,intent(out),dimension(ngrid) :: alim_star_tot
-      real,intent(out),dimension(ngrid,nlay) :: detr_star
-      real,intent(out),dimension(ngrid,nlay) :: entr_star
-      real,intent(out),dimension(ngrid,nlay+1) :: f_star
-      real,intent(out),dimension(ngrid,nlay) :: csc
-      real,intent(out),dimension(ngrid,nlay) :: ztva
-      real,intent(out),dimension(ngrid,nlay) :: ztla
-      real,intent(out),dimension(ngrid,nlay) :: zqla
-      real,intent(out),dimension(ngrid,nlay) :: zqta
-      real,intent(out),dimension(ngrid,nlay) :: zha
-      real,intent(out),dimension(ngrid,nlay+1) :: zw2
-      real,intent(out),dimension(ngrid,nlay+1) :: w_est
-      real,intent(out),dimension(ngrid,nlay) :: ztva_est
-      real,intent(out),dimension(ngrid,nlay) :: zqsatth
-      integer,intent(out),dimension(ngrid) :: lmix
-      integer,intent(out),dimension(ngrid) :: lmix_bis
-      real,intent(out),dimension(ngrid) :: linter
-
-      REAL zdw2,zdw2bis
-      REAL zw2modif
-      REAL zw2fact,zw2factbis
-      REAL,dimension(ngrid,nlay) :: zeps
-
-      REAL, dimension(ngrid) ::    wmaxa
-
-      INTEGER ig,l,k,lt,it,lm
-      integer nbpb
-
-      real,dimension(ngrid,nlay) :: detr
-      real,dimension(ngrid,nlay) :: entr
-      real,dimension(ngrid,nlay+1) :: wa_moy
-      real,dimension(ngrid,nlay) :: ztv_est
-      real,dimension(ngrid) :: ztemp,zqsat
-      real,dimension(ngrid,nlay) :: zqla_est
-      real,dimension(ngrid,nlay) :: zta_est
-
-      real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt
-      real zdz,zalpha,zw2m
-      real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam
-      real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
-      real, dimension(ngrid) :: d_temp
-      real ztv1,ztv2,factinv,zinv,zlmel
-      real zlmelup,zlmeldwn,zlt,zltdwn,zltup
-      real atv1,atv2,btv1,btv2
-      real ztv_est1,ztv_est2
-      real zcor,zdelta,zcvm5,qlbef
-      real zbetalpha, coefzlmel
-      real eps
-      logical Zsat
-      LOGICAL,dimension(ngrid) :: active,activetmp
-      REAL fact_gamma,fact_gamma2,fact_epsilon2
-      REAL coefc
-      REAL,dimension(ngrid,nlay) :: c2
-
-      if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
-      Zsat=.false.
-! Initialisation
-
-
-      zbetalpha=betalpha/(1.+betalpha)
-
-
-! Initialisations des variables r?elles
-if (1==1) then
-      ztva(:,:)=ztv(:,:)
-      ztva_est(:,:)=ztva(:,:)
-      ztv_est(:,:)=ztv(:,:)
-      ztla(:,:)=zthl(:,:)
-      zqta(:,:)=po(:,:)
-      zqla(:,:)=0.
-      zha(:,:) = ztva(:,:)
-else
-      ztva(:,:)=0.
-      ztv_est(:,:)=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.
-      zbuoy(:,:)=0.
-      zbuoyjam(:,:)=0.
-      gamma(:,:)=0.
-      zeps(:,:)=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.
-
-! Initialisation a 0  en cas de sortie dans replay
-      zqsat(:)=0.
-      zta_est(:,:)=0.
-      zdqt(:,:)=0.
-      zdqtjam(:,:)=0.
-      c2(:,:)=0.
-
-
-!-------------------------------------------------------------------------
-! On ne considere comme actif que les colonnes dont les deux premieres
-! couches sont instables.
-!-------------------------------------------------------------------------
-
-      active(:)=ztv(:,1)>ztv(:,2)
-      d_temp(:)=0. ! Pour activer un contraste de temperature a la base
-                   ! du panache
-!  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
-      CALL thermcell_alim(thermals_flag_alim,ngrid,nlay,ztv,d_temp,zlev,alim_star,lalim)
-
-!------------------------------------------------------------------------------
-! 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.
-    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
-
-!==============================================================================
-!boucle de calcul de la vitesse verticale dans le thermique
-!==============================================================================
-do l=2,nlay-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
-
-
-
-!---------------------------------------------------------------------------
-! calcul des proprietes thermodynamiques et de la vitesse de la couche l
-! sans tenir compte du detrainement et de l'entrainement dans cette
-! couche
-! C'est a dire qu'on suppose 
-! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
-! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
-! avant) a l'alimentation pour avoir un calcul plus propre
-!---------------------------------------------------------------------------
-
-   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
-   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
-    do ig=1,ngrid 
-!       print*,'active',active(ig),ig,l
-        if(active(ig)) then 
-        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
-        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)
-        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
-        -zqla_est(ig,l))-zqla_est(ig,l))
- 
-
-!Modif AJAM
-
-        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) 
-        zdz=zlev(ig,l+1)-zlev(ig,l)         
-        lmel=fact_thermals_ed_dz*zlev(ig,l)
-!        lmel=0.09*zlev(ig,l)
-        zlmel=zlev(ig,l)+lmel
-        zlmelup=zlmel+(zdz/2)
-        zlmeldwn=zlmel-(zdz/2)
-
-        lt=l+1
-        zlt=zlev(ig,lt)
-        zdz3=zlev(ig,lt+1)-zlt
-        zltdwn=zlt-zdz3/2
-        zltup=zlt+zdz3/2
-         
-!=========================================================================
-! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
-!=========================================================================
-
-!--------------------------------------------------
-        if (iflag_thermals_ed<8) then
-!--------------------------------------------------
-!AJ052014: J'ai remplac?? la boucle do par un do while
-! afin de faire moins de calcul dans la boucle
-!--------------------------------------------------
-            do while (zlmelup>zltup)
-               lt=lt+1
-               zlt=zlev(ig,lt)
-               zdz3=zlev(ig,lt+1)-zlt
-               zltdwn=zlt-zdz3/2
-               zltup=zlt+zdz3/2        
-            enddo
-!--------------------------------------------------
-!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
-! on cherche o?? se trouve l'altitude d'inversion 
-! en calculant ztv1 (interpolation de la valeur de 
-! theta au niveau lt en utilisant les niveaux lt-1 et
-! lt-2) et ztv2 (interpolation avec les niveaux lt+1
-! et lt+2). Si theta r??ellement calcul??e au niveau lt
-! comprise entre ztv1 et ztv2, alors il y a inversion
-! et on calcule son altitude zinv en supposant que ztv(lt)
-! est une combinaison lineaire de ztv1 et ztv2.
-! Ensuite, on calcule la flottabilite en comparant 
-! la temperature de la couche l a celle de l'air situe 
-! l+lmel plus haut, ce qui necessite de savoir quel fraction 
-! de cet air est au-dessus ou en-dessous de l'inversion   
-!--------------------------------------------------
-            atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2))
-            btv1=(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
-            /(zlev(ig,lt-1)-zlev(ig,lt-2))
-            atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1))
-            btv2=(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
-            /(zlev(ig,lt+2)-zlev(ig,lt+1))
-
-             ztv1=atv1*zlt+btv1
-             ztv2=atv2*zlt+btv2
-
-             if (ztv(ig,lt)>ztv1.and.ztv(ig,lt)<ztv2) then
-
-!--------------------------------------------------
-!AJ052014: D??calage de zinv qui est entre le haut
-!          et le bas de la couche lt
-!--------------------------------------------------
-                factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
-                zinv=zltdwn+zdz3*factinv
-
-          
-                if (zlmeldwn>=zinv) then
-                   ztv_est(ig,l)=atv2*zlmel+btv2
-                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
-                      +(1.-fact_shell)*zbuoy(ig,l)
-                elseif (zlmelup>=zinv) then
-                 ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2
-                   ztv_est1=atv1*0.5*(zinv+zlmeldwn)+btv1
-                   ztv_est(ig,l)=((zlmelup-zinv)/zdz)*ztv_est2+((zinv-zlmeldwn)/zdz)*ztv_est1
-
-                   zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- & 
-              ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
-              ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
-
-                else 
-                   ztv_est(ig,l)=atv1*zlmel+btv1
-                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) & 
-                             +(1.-fact_shell)*zbuoy(ig,l)
-                endif
-
-             else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
-
-                if (zlmeldwn>zltdwn) then
-                   zbuoyjam(ig,l)=fact_shell*RG*((ztva_est(ig,l)- & 
-                  ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l)
-                else
-                   zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & 
-                  ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
-                  ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
-
-                endif
-
-!          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & 
-!    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
-!    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
-!         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- & 
-!    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
-!     &          po(ig,lt-1))/po(ig,lt-1))
-          endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
-
-        else  !   if (iflag_thermals_ed.lt.8) then
-           lt=l+1
-           zlt=zlev(ig,lt)
-           zdz2=zlev(ig,lt)-zlev(ig,l)
-
-           do while (lmel>zdz2)
-             lt=lt+1
-             zlt=zlev(ig,lt)
-             zdz2=zlt-zlev(ig,l)
-           enddo
-           zdz3=zlev(ig,lt+1)-zlt
-           zltdwn=zlev(ig,lt)-zdz3/2
-           zlmelup=zlmel+(zdz/2)
-           coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
-           zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- & 
-            ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
-            ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
-        endif !   if (iflag_thermals_ed.lt.8) then
-
-!------------------------------------------------
-!AJAM:nouveau calcul de w?  
-!------------------------------------------------
-              zdz=zlev(ig,l+1)-zlev(ig,l)
-              zdzbis=zlev(ig,l)-zlev(ig,l-1)
-              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-
-              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
-              zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
-              zdw2=afact*zbuoy(ig,l)/fact_epsilon
-              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
-!              zdw2bis=0.5*(zdw2+zdw2bis)
-              lm=Max(1,l-2)
-!              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
-!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) 
-!              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
-!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
-!             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
-!             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
-!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
-!    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
-!              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
-
-!--------------------------------------------------
-!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
-!--------------------------------------------------
-         if (iflag_thermals_ed==8) then
-! Ancienne version
-!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
-!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
-!    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
-
-            w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
-
-! Nouvelle version Arnaud
-         else
-!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
-!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
-!    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
-
-            w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
-
-!             w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
-!    &                     (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
-!    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
-
-
-
-!            w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))
-
-!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
-!    &                      (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
-
-!             w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
-
-         endif
-
-
-         if (iflag_thermals_ed<6) then
-             zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
-!              fact_epsilon=0.0005/(zalpha+0.025)**0.5
-!              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
-              fact_epsilon=0.0002/(zalpha+0.1)
-              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
-              zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
-              zdw2=afact*zbuoy(ig,l)/fact_epsilon
-              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
-!              w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
-
-!              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
-!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
-!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
-
-            w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
-
-
-         endif
-!--------------------------------------------------
-!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
-!on fait max(0.0001,.....)
-!--------------------------------------------------         
-
-!             if (w_est(ig,l+1).lt.0.) then
-!               w_est(ig,l+1)=zw2(ig,l)
-!                w_est(ig,l+1)=0.0001
-!             endif
-
-       endif
-    enddo
-
-
-!-------------------------------------------------
-!calcul des taux d'entrainement et de detrainement
-!-------------------------------------------------
-
-     do ig=1,ngrid
-        if (active(ig)) then
-
-!          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
-          zw2m=w_est(ig,l+1)
-!          zw2m=zw2(ig,l)
-          zdz=zlev(ig,l+1)-zlev(ig,l)
-          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
-          zbuoybis=zbuoy(ig,l)
-          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
-          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
-
-          
-!          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
-!    &     afact*zbuoybis/zw2m - fact_epsilon )
-
-!          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
-!    &     afact*zbuoybis/zw2m - fact_epsilon )
-
-
-
-!          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-
-!=========================================================================
-! 4. Calcul de l'entrainement et du detrainement
-!=========================================================================
-
-!          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
-!    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon ) 
-!          entrbis=entr_star(ig,l)
-
-          if (iflag_thermals_ed<6) then
-          fact_epsilon=0.0002/(zalpha+0.1)
-          endif
-          
-
-
-          detr_star(ig,l)=f_star(ig,l)*zdz             &
-       *( mix0 * 0.1 / (zalpha+0.001)               &
-       + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
-       + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
-
-!          detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
-!    &                          ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
-
-          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-
-          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
-         mix0 * 0.1 / (zalpha+0.001)               &
-       + zbetalpha*MAX(entr_min,                   &
-       afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
-
-
-!          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
-!    &       mix0 * 0.1 / (zalpha+0.001)               &
-!    &     + MAX(entr_min,                   &
-!    &     zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon +  & 
-!    &     detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
-
-
-!          entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
-!    &                          ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
-
-!          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &     
-!    &     afact*zbuoy(ig,l)/zw2m &
-!    &     - 1.*fact_epsilon)
-
-          
-! En dessous de lalim, on prend le max de alim_star et entr_star pour
-! alim_star et 0 sinon
-        if (l<lalim(ig)) then
-          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
-          entr_star(ig,l)=0.
-        endif
-!        if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
-!          alim_star(ig,l)=entrbis
-!        endif
-
-!        print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
-! 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)
-
-      endif
-   enddo
-
-
-!============================================================================
-! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
-!===========================================================================
-
-   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)+  &
-              (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
-              /(f_star(ig,l+1)+detr_star(ig,l))
-
-        endif
-    enddo
-
-   ztemp(:)=zpspsk(:,l)*ztla(:,l)
-   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
-   do ig=1,ngrid
-      if (activetmp(ig)) then
-! on ecrit de maniere conservative (sat ou non)
-!          T = Tl +Lv/Cp ql
-           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
-           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
-           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
-!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
-           zha(ig,l) = ztva(ig,l)
-           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
-                -zqla(ig,l))-zqla(ig,l))
-           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
-           zdz=zlev(ig,l+1)-zlev(ig,l)
-           zdzbis=zlev(ig,l)-zlev(ig,l-1)
-           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
-!!!!!!!          fact_epsilon=0.002
-            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
-            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
-            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
-            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
-!              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
-!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) 
-!              lm=Max(1,l-2)
-!              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
-!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
-!            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
-!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
-!     &                   (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
-!            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
-!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
-!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
-!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
-            if (iflag_thermals_ed==8) then
-            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
-            else
-            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
-            endif
-!            zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
-!    &                     (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
-!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
-
-
-           if (iflag_thermals_ed<6) then
-           zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
-!           fact_epsilon=0.0005/(zalpha+0.025)**0.5
-!           fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
-           fact_epsilon=0.0002/(zalpha+0.1)**1
-            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
-            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
-            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
-            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
-
-!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
-!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
-!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
-!            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
-            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
-
-           endif
-
-
-      endif
-   enddo
-
-        if (prt_level>=20) print*,'coucou calcul detr 460: ig, l',ig, l
-
-!===========================================================================
-! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max 
-!===========================================================================
-
-   nbpb=0
-   do ig=1,ngrid
-            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1)<1.e-10) then
-!               stop'On tombe sur le cas particulier de thermcell_dry'
-!               print*,'On tombe sur le cas particulier de thermcell_plume'
-                nbpb=nbpb+1
-                zw2(ig,l+1)=0.
-                linter(ig)=l+1
-            endif
-
-        if (zw2(ig,l+1)<0.) then
-           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
-                 -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
-           zw2(ig,l+1)=0.
-!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
-        elseif (f_star(ig,l+1)<0.) then
-           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
-                 -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
-           zw2(ig,l+1)=0.
-!fin CR:04/05/12
-        endif
-
-           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 
-
-        if (wa_moy(ig,l+1)>wmaxa(ig)) then
-!   lmix est le niveau de la couche ou w (wa_moy) est maximum
-!on rajoute le calcul de lmix_bis
-            if (zqla(ig,l)<1.e-10) then
-               lmix_bis(ig)=l+1
-            endif
-            lmix(ig)=l+1
-            wmaxa(ig)=wa_moy(ig,l+1)
-        endif
-   enddo
-
-   if (nbpb>0) then
-   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
-   endif
-
-!=========================================================================
-! FIN DE LA BOUCLE VERTICALE
-      enddo
-!=========================================================================
-
-!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>=20) print*,'coucou calcul detr 470: ig, l', ig, l
-
-#undef wrgrads_thermcell
-#ifdef wrgrads_thermcell
-         call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
-         call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
-         call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
-         call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
-         call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
-         call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
-         call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
-#endif
-
-
- RETURN
-     end
-
-
-
-
-
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- SUBROUTINE thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
-&           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)
-!&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
-
-!--------------------------------------------------------------------------
-!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
-! Version conforme a l'article de Rio et al. 2010.
-! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
-!--------------------------------------------------------------------------
-
-      USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
-       USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
-      IMPLICIT NONE
-
-      INTEGER itap
-      INTEGER lunout1,igout
-      INTEGER ngrid,nlay
-      REAL ptimestep
-      REAL ztv(ngrid,nlay)
-      REAL zthl(ngrid,nlay)
-      REAL, INTENT(IN) :: po(ngrid,nlay)
-      REAL zl(ngrid,nlay)
-      REAL rhobarz(ngrid,nlay)
-      REAL zlev(ngrid,nlay+1)
-      REAL pplev(ngrid,nlay+1)
-      REAL pphi(ngrid,nlay)
-      REAL zpspsk(ngrid,nlay)
-      REAL alim_star(ngrid,nlay)
-      REAL f0(ngrid)
-      INTEGER lalim(ngrid)
-      integer lev_out                           ! niveau pour les print
-      integer nbpb
-    
-      real alim_star_tot(ngrid)
-
-      REAL ztva(ngrid,nlay)
-      REAL ztla(ngrid,nlay)
-      REAL zqla(ngrid,nlay)
-      REAL zqta(ngrid,nlay)
-      REAL zha(ngrid,nlay)
-
-      REAL detr_star(ngrid,nlay)
-      REAL coefc
-      REAL entr_star(ngrid,nlay)
-      REAL detr(ngrid,nlay)
-      REAL entr(ngrid,nlay)
-
-      REAL csc(ngrid,nlay)
-
-      REAL zw2(ngrid,nlay+1)
-      REAL w_est(ngrid,nlay+1)
-      REAL f_star(ngrid,nlay+1)
-      REAL wa_moy(ngrid,nlay+1)
-
-      REAL ztva_est(ngrid,nlay)
-      REAL zqla_est(ngrid,nlay)
-      REAL zqsatth(ngrid,nlay)
-      REAL zta_est(ngrid,nlay)
-      REAL zbuoyjam(ngrid,nlay)
-      REAL ztemp(ngrid),zqsat(ngrid)
-      REAL zdw2
-      REAL zw2modif
-      REAL zw2fact
-      REAL zeps(ngrid,nlay)
-
-      REAL linter(ngrid)
-      INTEGER lmix(ngrid)
-      INTEGER lmix_bis(ngrid)
-      REAL    wmaxa(ngrid)
-
-      INTEGER ig,l,k
-
-      real zdz,zbuoy(ngrid,nlay),zalpha,gamma(ngrid,nlay),zdqt(ngrid,nlay),zw2m
-      real zbuoybis
-      real zcor,zdelta,zcvm5,qlbef,zdz2
-      real betalpha,zbetalpha
-      real eps, afact
-      logical Zsat
-      LOGICAL active(ngrid),activetmp(ngrid)
-      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
-      REAL c2(ngrid,nlay)
-      Zsat=.false.
-! Initialisation
-
-      fact_epsilon=0.002
-      betalpha=0.9 
-      afact=2./3.            
-
-      zbetalpha=betalpha/(1.+betalpha)
-
-
-! Initialisations des variables reeles
-if (1==1) 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.
-      zbuoy(:,:)=0.
-      zbuoyjam(:,:)=0.
-      gamma(:,:)=0.
-      zeps(:,:)=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
-!-------------------------------------------------------------------------
-      do l=1,nlay-1
-         do ig=1,ngrid
-            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(ig)=l+1
-               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
-            endif
-         enddo
-      enddo
-      do l=1,nlay
-         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
-      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.
-    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
-
-!==============================================================================
-!boucle de calcul de la vitesse verticale dans le thermique
-!==============================================================================
-do l=2,nlay-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
-
-
-
-!---------------------------------------------------------------------------
-! calcul des proprietes thermodynamiques et de la vitesse de la couche l
-! sans tenir compte du detrainement et de l'entrainement dans cette
-! couche
-! C'est a dire qu'on suppose 
-! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
-! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
-! avant) a l'alimentation pour avoir un calcul plus propre
-!---------------------------------------------------------------------------
-
-   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
-   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
-
-    do ig=1,ngrid 
-!       print*,'active',active(ig),ig,l
-        if(active(ig)) then 
-        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
-        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)
-        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
-        -zqla_est(ig,l))-zqla_est(ig,l))
-
-!------------------------------------------------
-!AJAM:nouveau calcul de w?  
-!------------------------------------------------
-              zdz=zlev(ig,l+1)-zlev(ig,l)
-              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-
-              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
-              zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
-              w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
- 
-
-             if (w_est(ig,l+1)<0.) then
-                w_est(ig,l+1)=zw2(ig,l)
-             endif
-       endif
-    enddo
-
-
-!-------------------------------------------------
-!calcul des taux d'entrainement et de detrainement
-!-------------------------------------------------
-
-     do ig=1,ngrid
-        if (active(ig)) then
-
-          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
-          zw2m=w_est(ig,l+1)
-          zdz=zlev(ig,l+1)-zlev(ig,l)
-          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
-          zbuoybis=zbuoy(ig,l)
-          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
-          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
-
-          
-          entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
-       afact*zbuoybis/zw2m - fact_epsilon )
-
-
-          detr_star(ig,l)=f_star(ig,l)*zdz                        &
-       *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
-       + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
-          
-! En dessous de lalim, on prend le max de alim_star et entr_star pour
-! alim_star et 0 sinon
-        if (l<lalim(ig)) then
-          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
-          entr_star(ig,l)=0.
-        endif
-
-! 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)
-
-      endif
-   enddo
-
-
-!----------------------------------------------------------------------------
-!calcul de la vitesse verticale en melangeant Tl et qt du thermique
-!---------------------------------------------------------------------------
-   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)+  &
-              (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
-              /(f_star(ig,l+1)+detr_star(ig,l))
-
-        endif
-    enddo
-
-   ztemp(:)=zpspsk(:,l)*ztla(:,l)
-   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
-
-   do ig=1,ngrid
-      if (activetmp(ig)) then
-! on ecrit de maniere conservative (sat ou non)
-!          T = Tl +Lv/Cp ql
-           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
-           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
-           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
-!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
-           zha(ig,l) = ztva(ig,l)
-           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
-                -zqla(ig,l))-zqla(ig,l))
-           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
-           zdz=zlev(ig,l+1)-zlev(ig,l)
-           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
-
-            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
-            zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
-            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) 
-      endif
-   enddo
-
-        if (prt_level>=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 
-!---------------------------------------------------------------------------
-
-   nbpb=0
-   do ig=1,ngrid
-            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1)<1.e-10) then
-!               stop'On tombe sur le cas particulier de thermcell_dry'
-!               print*,'On tombe sur le cas particulier de thermcell_plume'
-                nbpb=nbpb+1
-                zw2(ig,l+1)=0.
-                linter(ig)=l+1
-            endif
-
-        if (zw2(ig,l+1)<0.) then
-           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
-                 -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
-           zw2(ig,l+1)=0.
-        elseif (f_star(ig,l+1)<0.) then
-           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
-                 -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
-!           print*,"linter plume", linter(ig)
-           zw2(ig,l+1)=0.
-        endif
-
-           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 
-
-        if (wa_moy(ig,l+1)>wmaxa(ig)) then
-!   lmix est le niveau de la couche ou w (wa_moy) est maximum
-!on rajoute le calcul de lmix_bis
-            if (zqla(ig,l)<1.e-10) then
-               lmix_bis(ig)=l+1
-            endif
-            lmix(ig)=l+1
-            wmaxa(ig)=wa_moy(ig,l+1)
-        endif
-   enddo
-
-   if (nbpb>0) then
-   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
-   endif
-
-!=========================================================================
-! FIN DE LA BOUCLE VERTICALE
-      enddo
-!=========================================================================
-
-!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>=20) print*,'coucou calcul detr 470: ig, l', ig, l
-
-#undef wrgrads_thermcell
-#ifdef wrgrads_thermcell
-         call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
-         call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
-         call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
-         call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
-         call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
-         call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
-         call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
-#endif
-
-
-     return 
-     end
-END MODULE lmdz_thermcell_plume_6A
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume_6A.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume_6A.f90	(revision 5102)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume_6A.f90	(revision 5102)
@@ -0,0 +1,1065 @@
+! $Id$
+
+
+MODULE lmdz_thermcell_plume_6A
+
+CONTAINS
+
+  SUBROUTINE thermcell_plume_6A(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
+          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)
+    !     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
+    !--------------------------------------------------------------------------
+    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
+    !--------------------------------------------------------------------------
+
+    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
+    USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
+    USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
+    USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
+    USE lmdz_thermcell_alim, ONLY: thermcell_alim
+    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
+
+    IMPLICIT NONE
+
+    integer, intent(in) :: itap, lev_out, lunout1, igout, ngrid, nlay
+    real, intent(in) :: ptimestep
+    real, intent(in), dimension(ngrid, nlay) :: ztv
+    real, intent(in), dimension(ngrid, nlay) :: zthl
+    real, intent(in), dimension(ngrid, nlay) :: po
+    real, intent(in), dimension(ngrid, nlay) :: zl
+    real, intent(in), dimension(ngrid, nlay) :: rhobarz
+    real, intent(in), dimension(ngrid, nlay + 1) :: zlev
+    real, intent(in), dimension(ngrid, nlay + 1) :: pplev
+    real, intent(in), dimension(ngrid, nlay) :: pphi
+    real, intent(in), dimension(ngrid, nlay) :: zpspsk
+    real, intent(in), dimension(ngrid) :: f0
+
+    integer, intent(out) :: lalim(ngrid)
+    real, intent(out), dimension(ngrid, nlay) :: alim_star
+    real, intent(out), dimension(ngrid) :: alim_star_tot
+    real, intent(out), dimension(ngrid, nlay) :: detr_star
+    real, intent(out), dimension(ngrid, nlay) :: entr_star
+    real, intent(out), dimension(ngrid, nlay + 1) :: f_star
+    real, intent(out), dimension(ngrid, nlay) :: csc
+    real, intent(out), dimension(ngrid, nlay) :: ztva
+    real, intent(out), dimension(ngrid, nlay) :: ztla
+    real, intent(out), dimension(ngrid, nlay) :: zqla
+    real, intent(out), dimension(ngrid, nlay) :: zqta
+    real, intent(out), dimension(ngrid, nlay) :: zha
+    real, intent(out), dimension(ngrid, nlay + 1) :: zw2
+    real, intent(out), dimension(ngrid, nlay + 1) :: w_est
+    real, intent(out), dimension(ngrid, nlay) :: ztva_est
+    real, intent(out), dimension(ngrid, nlay) :: zqsatth
+    integer, intent(out), dimension(ngrid) :: lmix
+    integer, intent(out), dimension(ngrid) :: lmix_bis
+    real, intent(out), dimension(ngrid) :: linter
+
+    REAL zdw2, zdw2bis
+    REAL zw2modif
+    REAL zw2fact, zw2factbis
+    REAL, dimension(ngrid, nlay) :: zeps
+
+    REAL, dimension(ngrid) :: wmaxa
+
+    INTEGER ig, l, k, lt, it, lm
+    integer nbpb
+
+    real, dimension(ngrid, nlay) :: detr
+    real, dimension(ngrid, nlay) :: entr
+    real, dimension(ngrid, nlay + 1) :: wa_moy
+    real, dimension(ngrid, nlay) :: ztv_est
+    real, dimension(ngrid) :: ztemp, zqsat
+    real, dimension(ngrid, nlay) :: zqla_est
+    real, dimension(ngrid, nlay) :: zta_est
+
+    real, dimension(ngrid, nlay) :: zbuoy, gamma, zdqt
+    real zdz, zalpha, zw2m
+    real, dimension(ngrid, nlay) :: zbuoyjam, zdqtjam
+    real zbuoybis, zdz2, zdz3, lmel, entrbis, zdzbis
+    real, dimension(ngrid) :: d_temp
+    real ztv1, ztv2, factinv, zinv, zlmel
+    real zlmelup, zlmeldwn, zlt, zltdwn, zltup
+    real atv1, atv2, btv1, btv2
+    real ztv_est1, ztv_est2
+    real zcor, zdelta, zcvm5, qlbef
+    real zbetalpha, coefzlmel
+    real eps
+    logical Zsat
+    LOGICAL, dimension(ngrid) :: active, activetmp
+    REAL fact_gamma, fact_gamma2, fact_epsilon2
+    REAL coefc
+    REAL, dimension(ngrid, nlay) :: c2
+
+    if (ngrid==1) print*, 'THERMCELL PLUME MODIFIE 2014/07/11'
+    Zsat = .false.
+    ! Initialisation
+
+    zbetalpha = betalpha / (1. + betalpha)
+
+
+    ! Initialisations des variables r?elles
+    if (1==1) then
+      ztva(:, :) = ztv(:, :)
+      ztva_est(:, :) = ztva(:, :)
+      ztv_est(:, :) = ztv(:, :)
+      ztla(:, :) = zthl(:, :)
+      zqta(:, :) = po(:, :)
+      zqla(:, :) = 0.
+      zha(:, :) = ztva(:, :)
+    else
+      ztva(:, :) = 0.
+      ztv_est(:, :) = 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.
+    zbuoy(:, :) = 0.
+    zbuoyjam(:, :) = 0.
+    gamma(:, :) = 0.
+    zeps(:, :) = 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.
+
+    ! Initialisation a 0  en cas de sortie dans replay
+    zqsat(:) = 0.
+    zta_est(:, :) = 0.
+    zdqt(:, :) = 0.
+    zdqtjam(:, :) = 0.
+    c2(:, :) = 0.
+
+
+    !-------------------------------------------------------------------------
+    ! On ne considere comme actif que les colonnes dont les deux premieres
+    ! couches sont instables.
+    !-------------------------------------------------------------------------
+
+    active(:) = ztv(:, 1)>ztv(:, 2)
+    d_temp(:) = 0. ! Pour activer un contraste de temperature a la base
+    ! du panache
+    !  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
+    CALL thermcell_alim(thermals_flag_alim, ngrid, nlay, ztv, d_temp, zlev, alim_star, lalim)
+
+    !------------------------------------------------------------------------------
+    ! 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.
+      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
+
+    !==============================================================================
+    !boucle de calcul de la vitesse verticale dans le thermique
+    !==============================================================================
+    do l = 2, nlay - 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
+
+
+
+      !---------------------------------------------------------------------------
+      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
+      ! sans tenir compte du detrainement et de l'entrainement dans cette
+      ! couche
+      ! C'est a dire qu'on suppose
+      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
+      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
+      ! avant) a l'alimentation pour avoir un calcul plus propre
+      !---------------------------------------------------------------------------
+
+      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
+      call thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
+      do ig = 1, ngrid
+        !       print*,'active',active(ig),ig,l
+        if(active(ig)) then
+          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
+          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)
+          ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1)  &
+                  - zqla_est(ig, l)) - zqla_est(ig, l))
+
+
+          !Modif AJAM
+
+          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          lmel = fact_thermals_ed_dz * zlev(ig, l)
+          !        lmel=0.09*zlev(ig,l)
+          zlmel = zlev(ig, l) + lmel
+          zlmelup = zlmel + (zdz / 2)
+          zlmeldwn = zlmel - (zdz / 2)
+
+          lt = l + 1
+          zlt = zlev(ig, lt)
+          zdz3 = zlev(ig, lt + 1) - zlt
+          zltdwn = zlt - zdz3 / 2
+          zltup = zlt + zdz3 / 2
+
+          !=========================================================================
+          ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
+          !=========================================================================
+
+          !--------------------------------------------------
+          if (iflag_thermals_ed<8) then
+            !--------------------------------------------------
+            !AJ052014: J'ai remplac?? la boucle do par un do while
+            ! afin de faire moins de calcul dans la boucle
+            !--------------------------------------------------
+            do while (zlmelup>zltup)
+              lt = lt + 1
+              zlt = zlev(ig, lt)
+              zdz3 = zlev(ig, lt + 1) - zlt
+              zltdwn = zlt - zdz3 / 2
+              zltup = zlt + zdz3 / 2
+            enddo
+            !--------------------------------------------------
+            !AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
+            ! on cherche o?? se trouve l'altitude d'inversion
+            ! en calculant ztv1 (interpolation de la valeur de
+            ! theta au niveau lt en utilisant les niveaux lt-1 et
+            ! lt-2) et ztv2 (interpolation avec les niveaux lt+1
+            ! et lt+2). Si theta r??ellement calcul??e au niveau lt
+            ! comprise entre ztv1 et ztv2, alors il y a inversion
+            ! et on calcule son altitude zinv en supposant que ztv(lt)
+            ! est une combinaison lineaire de ztv1 et ztv2.
+            ! Ensuite, on calcule la flottabilite en comparant
+            ! la temperature de la couche l a celle de l'air situe
+            ! l+lmel plus haut, ce qui necessite de savoir quel fraction
+            ! de cet air est au-dessus ou en-dessous de l'inversion
+            !--------------------------------------------------
+            atv1 = (ztv(ig, lt - 1) - ztv(ig, lt - 2)) / (zlev(ig, lt - 1) - zlev(ig, lt - 2))
+            btv1 = (ztv(ig, lt - 2) * zlev(ig, lt - 1) - ztv(ig, lt - 1) * zlev(ig, lt - 2)) &
+                    / (zlev(ig, lt - 1) - zlev(ig, lt - 2))
+            atv2 = (ztv(ig, lt + 2) - ztv(ig, lt + 1)) / (zlev(ig, lt + 2) - zlev(ig, lt + 1))
+            btv2 = (ztv(ig, lt + 1) * zlev(ig, lt + 2) - ztv(ig, lt + 2) * zlev(ig, lt + 1)) &
+                    / (zlev(ig, lt + 2) - zlev(ig, lt + 1))
+
+            ztv1 = atv1 * zlt + btv1
+            ztv2 = atv2 * zlt + btv2
+
+            if (ztv(ig, lt)>ztv1.and.ztv(ig, lt)<ztv2) then
+
+              !--------------------------------------------------
+              !AJ052014: D??calage de zinv qui est entre le haut
+              !          et le bas de la couche lt
+              !--------------------------------------------------
+              factinv = (ztv2 - ztv(ig, lt)) / (ztv2 - ztv1)
+              zinv = zltdwn + zdz3 * factinv
+
+              if (zlmeldwn>=zinv) then
+                ztv_est(ig, l) = atv2 * zlmel + btv2
+                zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) &
+                        + (1. - fact_shell) * zbuoy(ig, l)
+              elseif (zlmelup>=zinv) then
+                ztv_est2 = atv2 * 0.5 * (zlmelup + zinv) + btv2
+                ztv_est1 = atv1 * 0.5 * (zinv + zlmeldwn) + btv1
+                ztv_est(ig, l) = ((zlmelup - zinv) / zdz) * ztv_est2 + ((zinv - zlmeldwn) / zdz) * ztv_est1
+
+                zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zinv) / zdz) * (ztva_est(ig, l) - &
+                        ztv_est2) / ztv_est2 + ((zinv - zlmeldwn) / zdz) * (ztva_est(ig, l) - &
+                        ztv_est1) / ztv_est1) + (1. - fact_shell) * zbuoy(ig, l)
+
+              else
+                ztv_est(ig, l) = atv1 * zlmel + btv1
+                zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) &
+                        + (1. - fact_shell) * zbuoy(ig, l)
+              endif
+
+            else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
+
+              if (zlmeldwn>zltdwn) then
+                zbuoyjam(ig, l) = fact_shell * RG * ((ztva_est(ig, l) - &
+                        ztv(ig, lt)) / ztv(ig, lt)) + (1. - fact_shell) * zbuoy(ig, l)
+              else
+                zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zltdwn) / zdz) * (ztva_est(ig, l) - &
+                        ztv(ig, lt)) / ztv(ig, lt) + ((zltdwn - zlmeldwn) / zdz) * (ztva_est(ig, l) - &
+                        ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + (1. - fact_shell) * zbuoy(ig, l)
+
+              endif
+
+              !          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
+              !    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
+              !    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
+              !         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
+              !    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
+              !     &          po(ig,lt-1))/po(ig,lt-1))
+            endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
+
+          else  !   if (iflag_thermals_ed.lt.8) then
+            lt = l + 1
+            zlt = zlev(ig, lt)
+            zdz2 = zlev(ig, lt) - zlev(ig, l)
+
+            do while (lmel>zdz2)
+              lt = lt + 1
+              zlt = zlev(ig, lt)
+              zdz2 = zlt - zlev(ig, l)
+            enddo
+            zdz3 = zlev(ig, lt + 1) - zlt
+            zltdwn = zlev(ig, lt) - zdz3 / 2
+            zlmelup = zlmel + (zdz / 2)
+            coefzlmel = Min(1., (zlmelup - zltdwn) / zdz)
+            zbuoyjam(ig, l) = 1. * RG * (coefzlmel * (ztva_est(ig, l) - &
+                    ztv(ig, lt)) / ztv(ig, lt) + (1. - coefzlmel) * (ztva_est(ig, l) - &
+                    ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + 0. * zbuoy(ig, l)
+          endif !   if (iflag_thermals_ed.lt.8) then
+
+          !------------------------------------------------
+          !AJAM:nouveau calcul de w?
+          !------------------------------------------------
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
+          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
+
+          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
+          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
+          zdw2 = afact * zbuoy(ig, l) / fact_epsilon
+          zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
+          !              zdw2bis=0.5*(zdw2+zdw2bis)
+          lm = Max(1, l - 2)
+          !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
+          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
+          !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
+          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
+          !             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
+          !             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
+          !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
+          !    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
+          !              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
+
+          !--------------------------------------------------
+          !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
+          !--------------------------------------------------
+          if (iflag_thermals_ed==8) then
+            ! Ancienne version
+            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
+            !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
+            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
+
+            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
+
+            ! Nouvelle version Arnaud
+          else
+            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
+            !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
+            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
+
+            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2)
+
+            !             w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
+            !    &                     (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
+            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
+
+
+
+            !            w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))
+
+            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
+            !    &                      (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
+
+            !             w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
+
+          endif
+
+          if (iflag_thermals_ed<6) then
+            zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
+            !              fact_epsilon=0.0005/(zalpha+0.025)**0.5
+            !              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
+            fact_epsilon = 0.0002 / (zalpha + 0.1)
+            zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
+            zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
+            zdw2 = afact * zbuoy(ig, l) / fact_epsilon
+            zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
+            !              w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
+
+            !              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
+            !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
+            !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
+
+            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2)
+
+          endif
+          !--------------------------------------------------
+          !AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
+          !on fait max(0.0001,.....)
+          !--------------------------------------------------
+
+          !             if (w_est(ig,l+1).lt.0.) then
+          !               w_est(ig,l+1)=zw2(ig,l)
+          !                w_est(ig,l+1)=0.0001
+          !             endif
+
+        endif
+      enddo
+
+
+      !-------------------------------------------------
+      !calcul des taux d'entrainement et de detrainement
+      !-------------------------------------------------
+
+      do ig = 1, ngrid
+        if (active(ig)) then
+
+          !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
+          zw2m = w_est(ig, l + 1)
+          !          zw2m=zw2(ig,l)
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
+          !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
+          zbuoybis = zbuoy(ig, l)
+          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
+          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)
+
+
+          !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
+          !    &     afact*zbuoybis/zw2m - fact_epsilon )
+
+          !          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
+          !    &     afact*zbuoybis/zw2m - fact_epsilon )
+
+
+
+          !          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
+
+          !=========================================================================
+          ! 4. Calcul de l'entrainement et du detrainement
+          !=========================================================================
+
+          !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
+          !    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
+          !          entrbis=entr_star(ig,l)
+
+          if (iflag_thermals_ed<6) then
+            fact_epsilon = 0.0002 / (zalpha + 0.1)
+          endif
+
+          detr_star(ig, l) = f_star(ig, l) * zdz             &
+                  * (mix0 * 0.1 / (zalpha + 0.001)               &
+                          + MAX(detr_min, -afact * zbetalpha * zbuoyjam(ig, l) / zw2m   &
+                                  + detr_q_coef * (zdqt(ig, l) / zw2m)**detr_q_power))
+
+          !          detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
+          !    &                          ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
+
+          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
+
+          entr_star(ig, l) = f_star(ig, l) * zdz * (&
+                  mix0 * 0.1 / (zalpha + 0.001)               &
+                          + zbetalpha * MAX(entr_min, &
+                          afact * zbuoyjam(ig, l) / zw2m - fact_epsilon))
+
+
+          !          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
+          !    &       mix0 * 0.1 / (zalpha+0.001)               &
+          !    &     + MAX(entr_min,                   &
+          !    &     zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon +  &
+          !    &     detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
+
+
+          !          entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
+          !    &                          ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
+
+          !          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &
+          !    &     afact*zbuoy(ig,l)/zw2m &
+          !    &     - 1.*fact_epsilon)
+
+
+          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
+          ! alim_star et 0 sinon
+          if (l<lalim(ig)) then
+            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
+            entr_star(ig, l) = 0.
+          endif
+          !        if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
+          !          alim_star(ig,l)=entrbis
+          !        endif
+
+          !        print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
+          ! 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)
+
+        endif
+      enddo
+
+
+      !============================================================================
+      ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
+      !===========================================================================
+
+      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) + &
+                  (alim_star(ig, l) + entr_star(ig, l)) * po(ig, l))  &
+                  / (f_star(ig, l + 1) + detr_star(ig, l))
+
+        endif
+      enddo
+
+      ztemp(:) = zpspsk(:, l) * ztla(:, l)
+      call thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
+      do ig = 1, ngrid
+        if (activetmp(ig)) then
+          ! on ecrit de maniere conservative (sat ou non)
+          !          T = Tl +Lv/Cp ql
+          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
+          ztva(ig, l) = ztla(ig, l) * zpspsk(ig, l) + RLvCp * zqla(ig, l)
+          ztva(ig, l) = ztva(ig, l) / zpspsk(ig, l)
+          !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
+          zha(ig, l) = ztva(ig, l)
+          ztva(ig, l) = ztva(ig, l) * (1. + RETV * (zqta(ig, l)  &
+                  - zqla(ig, l)) - zqla(ig, l))
+          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
+          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
+          !!!!!!!          fact_epsilon=0.002
+          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
+          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
+          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
+          zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
+          !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
+          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
+          !              lm=Max(1,l-2)
+          !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
+          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
+          !            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
+          !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
+          !     &                   (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
+          !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
+          !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
+          !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
+          !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
+          if (iflag_thermals_ed==8) then
+            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
+          else
+            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2)
+          endif
+          !            zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
+          !    &                     (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
+          !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
+
+          if (iflag_thermals_ed<6) then
+            zalpha = f0(ig) * f_star(ig, l) / sqrt(zw2(ig, l + 1)) / rhobarz(ig, l)
+            !           fact_epsilon=0.0005/(zalpha+0.025)**0.5
+            !           fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
+            fact_epsilon = 0.0002 / (zalpha + 0.1)**1
+            zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
+            zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
+            zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
+            zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
+
+            !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
+            !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
+            !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
+            !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
+            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2)
+
+          endif
+
+        endif
+      enddo
+
+      if (prt_level>=20) print*, 'coucou calcul detr 460: ig, l', ig, l
+
+      !===========================================================================
+      ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
+      !===========================================================================
+
+      nbpb = 0
+      do ig = 1, ngrid
+        if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then
+          !               stop'On tombe sur le cas particulier de thermcell_dry'
+          !               print*,'On tombe sur le cas particulier de thermcell_plume'
+          nbpb = nbpb + 1
+          zw2(ig, l + 1) = 0.
+          linter(ig) = l + 1
+        endif
+
+        if (zw2(ig, l + 1)<0.) then
+          linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l))  &
+                  - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l))
+          zw2(ig, l + 1) = 0.
+          !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
+        elseif (f_star(ig, l + 1)<0.) then
+          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
+                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
+          zw2(ig, l + 1) = 0.
+          !fin CR:04/05/12
+        endif
+
+        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))
+
+        if (wa_moy(ig, l + 1)>wmaxa(ig)) then
+          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
+          !on rajoute le calcul de lmix_bis
+          if (zqla(ig, l)<1.e-10) then
+            lmix_bis(ig) = l + 1
+          endif
+          lmix(ig) = l + 1
+          wmaxa(ig) = wa_moy(ig, l + 1)
+        endif
+      enddo
+
+      if (nbpb>0) then
+        print*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
+      endif
+
+      !=========================================================================
+      ! FIN DE LA BOUCLE VERTICALE
+    enddo
+    !=========================================================================
+
+    !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>=20) print*, 'coucou calcul detr 470: ig, l', ig, l
+    RETURN
+  end
+
+
+  SUBROUTINE thermcell_plume_5B(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
+          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)
+    !&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
+
+    !--------------------------------------------------------------------------
+    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
+    ! Version conforme a l'article de Rio et al. 2010.
+    ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
+    !--------------------------------------------------------------------------
+
+    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
+    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
+    IMPLICIT NONE
+
+    INTEGER itap
+    INTEGER lunout1, igout
+    INTEGER ngrid, nlay
+    REAL ptimestep
+    REAL ztv(ngrid, nlay)
+    REAL zthl(ngrid, nlay)
+    REAL, INTENT(IN) :: po(ngrid, nlay)
+    REAL zl(ngrid, nlay)
+    REAL rhobarz(ngrid, nlay)
+    REAL zlev(ngrid, nlay + 1)
+    REAL pplev(ngrid, nlay + 1)
+    REAL pphi(ngrid, nlay)
+    REAL zpspsk(ngrid, nlay)
+    REAL alim_star(ngrid, nlay)
+    REAL f0(ngrid)
+    INTEGER lalim(ngrid)
+    integer lev_out                           ! niveau pour les print
+    integer nbpb
+
+    real alim_star_tot(ngrid)
+
+    REAL ztva(ngrid, nlay)
+    REAL ztla(ngrid, nlay)
+    REAL zqla(ngrid, nlay)
+    REAL zqta(ngrid, nlay)
+    REAL zha(ngrid, nlay)
+
+    REAL detr_star(ngrid, nlay)
+    REAL coefc
+    REAL entr_star(ngrid, nlay)
+    REAL detr(ngrid, nlay)
+    REAL entr(ngrid, nlay)
+
+    REAL csc(ngrid, nlay)
+
+    REAL zw2(ngrid, nlay + 1)
+    REAL w_est(ngrid, nlay + 1)
+    REAL f_star(ngrid, nlay + 1)
+    REAL wa_moy(ngrid, nlay + 1)
+
+    REAL ztva_est(ngrid, nlay)
+    REAL zqla_est(ngrid, nlay)
+    REAL zqsatth(ngrid, nlay)
+    REAL zta_est(ngrid, nlay)
+    REAL zbuoyjam(ngrid, nlay)
+    REAL ztemp(ngrid), zqsat(ngrid)
+    REAL zdw2
+    REAL zw2modif
+    REAL zw2fact
+    REAL zeps(ngrid, nlay)
+
+    REAL linter(ngrid)
+    INTEGER lmix(ngrid)
+    INTEGER lmix_bis(ngrid)
+    REAL    wmaxa(ngrid)
+
+    INTEGER ig, l, k
+
+    real zdz, zbuoy(ngrid, nlay), zalpha, gamma(ngrid, nlay), zdqt(ngrid, nlay), zw2m
+    real zbuoybis
+    real zcor, zdelta, zcvm5, qlbef, zdz2
+    real betalpha, zbetalpha
+    real eps, afact
+    logical Zsat
+    LOGICAL active(ngrid), activetmp(ngrid)
+    REAL fact_gamma, fact_epsilon, fact_gamma2, fact_epsilon2
+    REAL c2(ngrid, nlay)
+    Zsat = .false.
+    ! Initialisation
+
+    fact_epsilon = 0.002
+    betalpha = 0.9
+    afact = 2. / 3.
+
+    zbetalpha = betalpha / (1. + betalpha)
+
+
+    ! Initialisations des variables reeles
+    if (1==1) 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.
+    zbuoy(:, :) = 0.
+    zbuoyjam(:, :) = 0.
+    gamma(:, :) = 0.
+    zeps(:, :) = 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
+    !-------------------------------------------------------------------------
+    do l = 1, nlay - 1
+      do ig = 1, ngrid
+        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(ig) = l + 1
+          alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
+        endif
+      enddo
+    enddo
+    do l = 1, nlay
+      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
+    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.
+      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
+
+    !==============================================================================
+    !boucle de calcul de la vitesse verticale dans le thermique
+    !==============================================================================
+    do l = 2, nlay - 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
+
+
+
+      !---------------------------------------------------------------------------
+      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
+      ! sans tenir compte du detrainement et de l'entrainement dans cette
+      ! couche
+      ! C'est a dire qu'on suppose
+      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
+      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
+      ! avant) a l'alimentation pour avoir un calcul plus propre
+      !---------------------------------------------------------------------------
+
+      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
+      call thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
+
+      do ig = 1, ngrid
+        !       print*,'active',active(ig),ig,l
+        if(active(ig)) then
+          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
+          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)
+          ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1)  &
+                  - zqla_est(ig, l)) - zqla_est(ig, l))
+
+          !------------------------------------------------
+          !AJAM:nouveau calcul de w?
+          !------------------------------------------------
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
+
+          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
+          zdw2 = (afact) * zbuoy(ig, l) / (fact_epsilon)
+          w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
+
+          if (w_est(ig, l + 1)<0.) then
+            w_est(ig, l + 1) = zw2(ig, l)
+          endif
+        endif
+      enddo
+
+
+      !-------------------------------------------------
+      !calcul des taux d'entrainement et de detrainement
+      !-------------------------------------------------
+
+      do ig = 1, ngrid
+        if (active(ig)) then
+
+          zw2m = max(0.5 * (w_est(ig, l) + w_est(ig, l + 1)), 0.1)
+          zw2m = w_est(ig, l + 1)
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
+          !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
+          zbuoybis = zbuoy(ig, l)
+          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
+          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)
+
+          entr_star(ig, l) = f_star(ig, l) * zdz * zbetalpha * MAX(0., &
+                  afact * zbuoybis / zw2m - fact_epsilon)
+
+          detr_star(ig, l) = f_star(ig, l) * zdz                        &
+                  * MAX(1.e-3, -afact * zbetalpha * zbuoy(ig, l) / zw2m          &
+                          + 0.012 * (zdqt(ig, l) / zw2m)**0.5)
+
+          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
+          ! alim_star et 0 sinon
+          if (l<lalim(ig)) then
+            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
+            entr_star(ig, l) = 0.
+          endif
+
+          ! 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)
+
+        endif
+      enddo
+
+
+      !----------------------------------------------------------------------------
+      !calcul de la vitesse verticale en melangeant Tl et qt du thermique
+      !---------------------------------------------------------------------------
+      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) + &
+                  (alim_star(ig, l) + entr_star(ig, l)) * po(ig, l))  &
+                  / (f_star(ig, l + 1) + detr_star(ig, l))
+
+        endif
+      enddo
+
+      ztemp(:) = zpspsk(:, l) * ztla(:, l)
+      call thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
+
+      do ig = 1, ngrid
+        if (activetmp(ig)) then
+          ! on ecrit de maniere conservative (sat ou non)
+          !          T = Tl +Lv/Cp ql
+          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
+          ztva(ig, l) = ztla(ig, l) * zpspsk(ig, l) + RLvCp * zqla(ig, l)
+          ztva(ig, l) = ztva(ig, l) / zpspsk(ig, l)
+          !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
+          zha(ig, l) = ztva(ig, l)
+          ztva(ig, l) = ztva(ig, l) * (1. + RETV * (zqta(ig, l)  &
+                  - zqla(ig, l)) - zqla(ig, l))
+          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
+          zdz = zlev(ig, l + 1) - zlev(ig, l)
+          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
+
+          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
+          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
+          zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
+        endif
+      enddo
+
+      if (prt_level>=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
+      !---------------------------------------------------------------------------
+
+      nbpb = 0
+      do ig = 1, ngrid
+        if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then
+          !               stop'On tombe sur le cas particulier de thermcell_dry'
+          !               print*,'On tombe sur le cas particulier de thermcell_plume'
+          nbpb = nbpb + 1
+          zw2(ig, l + 1) = 0.
+          linter(ig) = l + 1
+        endif
+
+        if (zw2(ig, l + 1)<0.) then
+          linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l))  &
+                  - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l))
+          zw2(ig, l + 1) = 0.
+        elseif (f_star(ig, l + 1)<0.) then
+          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
+                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
+          !           print*,"linter plume", linter(ig)
+          zw2(ig, l + 1) = 0.
+        endif
+
+        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))
+
+        if (wa_moy(ig, l + 1)>wmaxa(ig)) then
+          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
+          !on rajoute le calcul de lmix_bis
+          if (zqla(ig, l)<1.e-10) then
+            lmix_bis(ig) = l + 1
+          endif
+          lmix(ig) = l + 1
+          wmaxa(ig) = wa_moy(ig, l + 1)
+        endif
+      enddo
+
+      if (nbpb>0) then
+        print*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
+      endif
+
+      !=========================================================================
+      ! FIN DE LA BOUCLE VERTICALE
+    enddo
+    !=========================================================================
+
+    !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>=20) print*, 'coucou calcul detr 470: ig, l', ig, l
+  end
+END MODULE lmdz_thermcell_plume_6A
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/surf_landice_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/surf_landice_mod.F90	(revision 5101)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/surf_landice_mod.F90	(revision 5102)
@@ -1,42 +1,41 @@
-
 MODULE surf_landice_mod
-  
+
   IMPLICIT NONE
 
 CONTAINS
 
-!****************************************************************************************
+  !****************************************************************************************
 
   SUBROUTINE surf_landice(itime, dtime, knon, knindex, &
-       rlon, rlat, debut, lafin, &
-       rmu0, lwdownm, albedo, pphi1, &
-       swnet, lwnet, tsurf, p1lay, &
-       cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
-       AcoefH, AcoefQ, BcoefH, BcoefQ, &
-       AcoefU, AcoefV, BcoefU, BcoefV, &
-       AcoefQBS, BcoefQBS, &
-       ps, u1, v1, gustiness, rugoro, pctsrf, &
-       snow, qsurf, qsol, qbs1, agesno, &
-       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
-       tsurf_new, dflux_s, dflux_l, &
-       alt, slope, cloudf, &
-       snowhgt, qsnow, to_ice, sissnow, &
-       alb3, runoff, &
-       flux_u1, flux_v1 &
+            rlon, rlat, debut, lafin, &
+            rmu0, lwdownm, albedo, pphi1, &
+            swnet, lwnet, tsurf, p1lay, &
+            cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
+            AcoefH, AcoefQ, BcoefH, BcoefQ, &
+            AcoefU, AcoefV, BcoefU, BcoefV, &
+            AcoefQBS, BcoefQBS, &
+            ps, u1, v1, gustiness, rugoro, pctsrf, &
+            snow, qsurf, qsol, qbs1, agesno, &
+            tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
+            tsurf_new, dflux_s, dflux_l, &
+            alt, slope, cloudf, &
+            snowhgt, qsnow, to_ice, sissnow, &
+            alb3, runoff, &
+            flux_u1, flux_v1 &
 #ifdef ISO
         ,xtprecip_rain, xtprecip_snow,xtspechum,Rland_ice &
         ,xtsnow,xtsol,xtevap &
 #endif               
-      )
+            )
 
     USE dimphy
-    USE geometry_mod,     ONLY: longitude,latitude
-    USE surface_data,     ONLY: type_ocean, calice, calsno, landice_opt, iflag_albcalc
-    USE fonte_neige_mod,  ONLY: fonte_neige,run_off_lic,fqcalving_global,ffonte_global,fqfonte_global,runofflic_global
-    USE cpl_mod,          ONLY: cpl_send_landice_fields
+    USE geometry_mod, ONLY: longitude, latitude
+    USE surface_data, ONLY: type_ocean, calice, calsno, landice_opt, iflag_albcalc
+    USE fonte_neige_mod, ONLY: fonte_neige, run_off_lic, fqcalving_global, ffonte_global, fqfonte_global, runofflic_global
+    USE cpl_mod, ONLY: cpl_send_landice_fields
     USE calcul_fluxs_mod
     USE phys_local_var_mod, ONLY: zxrhoslic, zxustartlic, zxqsaltlic
-    USE phys_output_var_mod, ONLY: snow_o,zfra_o
-#ifdef ISO   
+    USE phys_output_var_mod, ONLY: snow_o, zfra_o
+#ifdef ISO
     USE fonte_neige_mod,  ONLY: xtrun_off_lic
     USE infotrac_phy,     ONLY: ntiso,niso
@@ -47,40 +46,39 @@
 #endif
 #endif
- 
-!FC
+
+    !FC
     USE ioipsl_getin_p_mod, ONLY: getin_p
     USE lmdz_blowing_snow_ini, ONLY: c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs
     USE lmdz_blowing_snow_ini, ONLY: rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs
-#ifdef CPP_INLANDSIS
-    USE surf_inlandsis_mod,  ONLY: surf_inlandsis
-#endif
+    USE surf_inlandsis_mod, ONLY: surf_inlandsis
+    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INLANDSIS
 
     USE indice_sol_mod
 
-!    INCLUDE "indicesol.h"
+    !    INCLUDE "indicesol.h"
     INCLUDE "dimsoil.h"
     INCLUDE "YOMCST.h"
     INCLUDE "clesphys.h"
 
-! Input variables 
-!****************************************************************************************
-    INTEGER, INTENT(IN)                           :: itime, knon
-    INTEGER, DIMENSION(klon), INTENT(in)          :: knindex
-    REAL, INTENT(in)                              :: dtime
-    REAL, DIMENSION(klon), INTENT(IN)             :: swnet ! net shortwave radiance
-    REAL, DIMENSION(klon), INTENT(IN)             :: lwnet ! net longwave radiance
-    REAL, DIMENSION(klon), INTENT(IN)             :: tsurf
-    REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
-    REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
-    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
-    REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
-    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
-    REAL, DIMENSION(klon), INTENT(IN)             :: BcoefH, BcoefQ
-    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
-    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefQBS, BcoefQBS
-    REAL, DIMENSION(klon), INTENT(IN)             :: ps
-    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
-    REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
-    REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
+    ! Input variables
+    !****************************************************************************************
+    INTEGER, INTENT(IN) :: itime, knon
+    INTEGER, DIMENSION(klon), INTENT(in) :: knindex
+    REAL, INTENT(in) :: dtime
+    REAL, DIMENSION(klon), INTENT(IN) :: swnet ! net shortwave radiance
+    REAL, DIMENSION(klon), INTENT(IN) :: lwnet ! net longwave radiance
+    REAL, DIMENSION(klon), INTENT(IN) :: tsurf
+    REAL, DIMENSION(klon), INTENT(IN) :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
+    REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow, precip_bs
+    REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
+    REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
+    REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
+    REAL, DIMENSION(klon), INTENT(IN) :: AcoefQBS, BcoefQBS
+    REAL, DIMENSION(klon), INTENT(IN) :: ps
+    REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness, qbs1
+    REAL, DIMENSION(klon), INTENT(IN) :: rugoro
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf
 #ifdef ISO
     REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtprecip_rain, xtprecip_snow 
@@ -89,19 +87,19 @@
 
 
-    LOGICAL,  INTENT(IN)                          :: debut   !true if first step
-    LOGICAL,  INTENT(IN)                          :: lafin   !true if last step
-    REAL, DIMENSION(klon), INTENT(IN)             :: rlon, rlat
-    REAL, DIMENSION(klon), INTENT(IN)             :: rmu0
-    REAL, DIMENSION(klon), INTENT(IN)             :: lwdownm !ylwdown
-    REAL, DIMENSION(klon), INTENT(IN)             :: albedo  !mean albedo
-    REAL, DIMENSION(klon), INTENT(IN)             :: pphi1    
-    REAL, DIMENSION(klon), INTENT(IN)             :: alt   !mean altitude of the grid box  
-    REAL, DIMENSION(klon), INTENT(IN)             :: slope   !mean slope in grid box  
-    REAL, DIMENSION(klon), INTENT(IN)             :: cloudf  !total cloud fraction
-
-! In/Output variables
-!****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
-    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    LOGICAL, INTENT(IN) :: debut   !true if first step
+    LOGICAL, INTENT(IN) :: lafin   !true if last step
+    REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
+    REAL, DIMENSION(klon), INTENT(IN) :: rmu0
+    REAL, DIMENSION(klon), INTENT(IN) :: lwdownm !ylwdown
+    REAL, DIMENSION(klon), INTENT(IN) :: albedo  !mean albedo
+    REAL, DIMENSION(klon), INTENT(IN) :: pphi1
+    REAL, DIMENSION(klon), INTENT(IN) :: alt   !mean altitude of the grid box
+    REAL, DIMENSION(klon), INTENT(IN) :: slope   !mean slope in grid box
+    REAL, DIMENSION(klon), INTENT(IN) :: cloudf  !total cloud fraction
+
+    ! In/Output variables
+    !****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
+    REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
     REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
 #ifdef ISO
@@ -111,26 +109,26 @@
 
 
-! Output variables
-!****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
-    REAL, DIMENSION(klon), INTENT(OUT)            :: z0m, z0h
-!albedo SB >>>
-!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1  ! new albedo in visible SW interval
-!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2  ! new albedo in near IR interval
-    REAL, DIMENSION(6), INTENT(IN)                :: SFRWL
-    REAL, DIMENSION(klon,nsw), INTENT(OUT)        :: alb_dir,alb_dif
-!albedo SB <<<
-    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
-    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
-    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
-
-    REAL, DIMENSION(klon), INTENT(OUT)           :: alb3
-    REAL, DIMENSION(klon), INTENT(OUT)           :: qsnow   !column water in snow [kg/m2]
-    REAL, DIMENSION(klon), INTENT(OUT)           :: snowhgt !Snow height (m)
-    REAL, DIMENSION(klon), INTENT(OUT)           :: to_ice
-    REAL, DIMENSION(klon), INTENT(OUT)           :: sissnow
-    REAL, DIMENSION(klon), INTENT(OUT)           :: runoff  !Land ice runoff
+    ! Output variables
+    !****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT) :: z0m, z0h
+    !albedo SB >>>
+    !    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1  ! new albedo in visible SW interval
+    !    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2  ! new albedo in near IR interval
+    REAL, DIMENSION(6), INTENT(IN) :: SFRWL
+    REAL, DIMENSION(klon, nsw), INTENT(OUT) :: alb_dir, alb_dif
+    !albedo SB <<<
+    REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT) :: fluxbs
+    REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
+    REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
+    REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
+
+    REAL, DIMENSION(klon), INTENT(OUT) :: alb3
+    REAL, DIMENSION(klon), INTENT(OUT) :: qsnow   !column water in snow [kg/m2]
+    REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt !Snow height (m)
+    REAL, DIMENSION(klon), INTENT(OUT) :: to_ice
+    REAL, DIMENSION(klon), INTENT(OUT) :: sissnow
+    REAL, DIMENSION(klon), INTENT(OUT) :: runoff  !Land ice runoff
 #ifdef ISO
     REAL, DIMENSION(ntiso,klon), INTENT(OUT)     :: xtevap     
@@ -138,16 +136,16 @@
 !    fonte_neige
 #endif
- 
-
-! Local variables
-!****************************************************************************************
-    REAL, DIMENSION(klon)    :: soilcap, soilflux
-    REAL, DIMENSION(klon)    :: cal, beta, dif_grnd
-    REAL, DIMENSION(klon)    :: zfra, alb_neig
-    REAL, DIMENSION(klon)    :: radsol
-    REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay, ustar
-    INTEGER                  :: i,j,nt
-    REAL, DIMENSION(klon)    :: fqfonte,ffonte
-    REAL, DIMENSION(klon)    :: run_off_lic_frac
+
+
+    ! Local variables
+    !****************************************************************************************
+    REAL, DIMENSION(klon) :: soilcap, soilflux
+    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
+    REAL, DIMENSION(klon) :: zfra, alb_neig
+    REAL, DIMENSION(klon) :: radsol
+    REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay, ustar
+    INTEGER :: i, j, nt
+    REAL, DIMENSION(klon) :: fqfonte, ffonte
+    REAL, DIMENSION(klon) :: run_off_lic_frac
 #ifdef ISO       
     REAL, PARAMETER          :: t_coup = 273.15
@@ -167,47 +165,46 @@
 
 
-    REAL, DIMENSION(klon)    :: emis_new                  !Emissivity
-    REAL, DIMENSION(klon)    :: swdown,lwdown
-    REAL, DIMENSION(klon)    :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
-    REAL, DIMENSION(klon)    :: erod                      !erosion of surface snow (flux, kg/m2/s like evap)
-    REAL, DIMENSION(klon)    :: zsl_height, wind_velo     !surface layer height, wind spd
-    REAL, DIMENSION(klon)    :: dens_air,  snow_cont_air  !air density; snow content air
-    REAL, DIMENSION(klon)    :: alb_soil                  !albedo of underlying ice
-    REAL, DIMENSION(klon)    :: pexner                    !Exner potential
-    REAL                     :: pref
-    REAL, DIMENSION(klon,nsoilmx) :: tsoil0               !modif
-    REAL                          :: dtis                ! subtimestep
-    LOGICAL                       :: debut_is, lafin_is  ! debut and lafin for inlandsis
-
-    CHARACTER (len = 20)                      :: modname = 'surf_landice'
-    CHARACTER (len = 80)                      :: abort_message
-
-
-    REAL,DIMENSION(klon) :: alb1,alb2
-    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
-    REAL, DIMENSION (klon,6) :: alb6
-    REAL                   :: esalt
-    REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
-    REAL                   :: tau_dens, maxerosion
-    REAL, DIMENSION(klon)  :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
-    REAL, DIMENSION(klon)  :: fluxbs_1, fluxbs_2, bsweight_fresh
+    REAL, DIMENSION(klon) :: emis_new                  !Emissivity
+    REAL, DIMENSION(klon) :: swdown, lwdown
+    REAL, DIMENSION(klon) :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
+    REAL, DIMENSION(klon) :: erod                      !erosion of surface snow (flux, kg/m2/s like evap)
+    REAL, DIMENSION(klon) :: zsl_height, wind_velo     !surface layer height, wind spd
+    REAL, DIMENSION(klon) :: dens_air, snow_cont_air  !air density; snow content air
+    REAL, DIMENSION(klon) :: alb_soil                  !albedo of underlying ice
+    REAL, DIMENSION(klon) :: pexner                    !Exner potential
+    REAL :: pref
+    REAL, DIMENSION(klon, nsoilmx) :: tsoil0               !modif
+    REAL :: dtis                ! subtimestep
+    LOGICAL :: debut_is, lafin_is  ! debut and lafin for inlandsis
+
+    CHARACTER (len = 20) :: modname = 'surf_landice'
+    CHARACTER (len = 80) :: abort_message
+
+    REAL, DIMENSION(klon) :: alb1, alb2
+    REAL, DIMENSION(klon) :: precip_totsnow, evap_totsnow
+    REAL, DIMENSION (klon, 6) :: alb6
+    REAL :: esalt
+    REAL :: lambdasalt, fluxsalt, csalt, nunu, aa, bb, cc
+    REAL :: tau_dens, maxerosion
+    REAL, DIMENSION(klon) :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
+    REAL, DIMENSION(klon) :: fluxbs_1, fluxbs_2, bsweight_fresh
     LOGICAL, DIMENSION(klon) :: ok_remaining_freshsnow
-    REAL  :: ta1, ta2, ta3, z01, z02, z03, coefa, coefb, coefc, coefd
-
-
-! End definition
-!****************************************************************************************
-!FC 
-!FC
-   REAL,SAVE :: alb_vis_sno_lic
-  !$OMP THREADPRIVATE(alb_vis_sno_lic)
-   REAL,SAVE :: alb_nir_sno_lic
-  !$OMP THREADPRIVATE(alb_nir_sno_lic)
-  LOGICAL, SAVE :: firstcall = .TRUE.
-  !$OMP THREADPRIVATE(firstcall)
-
-
-!FC firtscall initializations
-!******************************************************************************************
+    REAL :: ta1, ta2, ta3, z01, z02, z03, coefa, coefb, coefc, coefd
+
+
+    ! End definition
+    !****************************************************************************************
+    !FC
+    !FC
+    REAL, SAVE :: alb_vis_sno_lic
+    !$OMP THREADPRIVATE(alb_vis_sno_lic)
+    REAL, SAVE :: alb_nir_sno_lic
+    !$OMP THREADPRIVATE(alb_nir_sno_lic)
+    LOGICAL, SAVE :: firstcall = .TRUE.
+    !$OMP THREADPRIVATE(firstcall)
+
+
+    !FC firtscall initializations
+    !******************************************************************************************
 #ifdef ISO
 #ifdef ISOVERIF
@@ -222,43 +219,42 @@
 #endif
 
-  IF (firstcall) THEN
-  alb_vis_sno_lic=0.77
-  CALL getin_p('alb_vis_sno_lic',alb_vis_sno_lic)
-           PRINT*, 'alb_vis_sno_lic',alb_vis_sno_lic
-  alb_nir_sno_lic=0.77
-  CALL getin_p('alb_nir_sno_lic',alb_nir_sno_lic)
-           PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
-  
-  firstcall=.false.
-  ENDIF
-!******************************************************************************************
-
-! Initialize output variables
+    IF (firstcall) THEN
+      alb_vis_sno_lic = 0.77
+      CALL getin_p('alb_vis_sno_lic', alb_vis_sno_lic)
+      PRINT*, 'alb_vis_sno_lic', alb_vis_sno_lic
+      alb_nir_sno_lic = 0.77
+      CALL getin_p('alb_nir_sno_lic', alb_nir_sno_lic)
+      PRINT*, 'alb_nir_sno_lic', alb_nir_sno_lic
+
+      firstcall = .false.
+    ENDIF
+    !******************************************************************************************
+
+    ! Initialize output variables
     alb3(:) = 999999.
     alb2(:) = 999999.
     alb1(:) = 999999.
-    fluxbs(:)=0.  
+    fluxbs(:) = 0.
     runoff(:) = 0.
-!****************************************************************************************
-! Calculate total absorbed radiance at surface
-
-!****************************************************************************************
+    !****************************************************************************************
+    ! Calculate total absorbed radiance at surface
+
+    !****************************************************************************************
     radsol(:) = 0.0
     radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
 
-!****************************************************************************************
-
-!****************************************************************************************
-!  landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ...  
-!  landice_opt = 1  : prepare and call INterace Lmdz SISvat (INLANDSIS)
-!****************************************************************************************
-
+    !****************************************************************************************
+
+    !****************************************************************************************
+    !  landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ...
+    !  landice_opt = 1  : prepare and call INterace Lmdz SISvat (INLANDSIS)
+    !****************************************************************************************
 
     IF (landice_opt == 1) THEN
 
-!****************************************************************************************    
-! CALL to INLANDSIS interface
-!****************************************************************************************
-#ifdef CPP_INLANDSIS
+      !****************************************************************************************
+      ! CALL to INLANDSIS interface
+      !****************************************************************************************
+      IF (CPPKEY_INLANDSIS) THEN
 
 #ifdef ISO
@@ -266,58 +262,55 @@
 #endif
 
-        debut_is=debut
-        lafin_is=.false.
+        debut_is = debut
+        lafin_is = .false.
         ! Suppose zero surface speed
-        u0(:)            = 0.0
-        v0(:)            = 0.0
-
+        u0(:) = 0.0
+        v0(:) = 0.0
 
         CALL calcul_flux_wind(knon, dtime, &
-         u0, v0, u1, v1, gustiness, cdragm, &
-         AcoefU, AcoefV, BcoefU, BcoefV, &
-         p1lay, temp_air, &
-         flux_u1, flux_v1)
-
-       
-       ! Set constants and compute some input for SISVAT
-       ! = 1000 hPa
-       ! and calculate incoming flux for SW and LW interval: swdown, lwdown
-       swdown(:)        = 0.0
-       lwdown(:)        = 0.0
-       snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model      
-       alb_soil(:)      = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
-       ustar(:)         = 0.
-       pref             = 100000.       
-       DO i = 1, knon
-          swdown(i)        = swnet(i)/(1-albedo(i))
-          lwdown(i)        = lwdownm(i)
-          wind_velo(i)     = u1(i)**2 + v1(i)**2
-          wind_velo(i)     = wind_velo(i)**0.5
-          pexner(i)        = (p1lay(i)/pref)**(RD/RCPD)
-          dens_air(i)      = p1lay(i)/RD/temp_air(i)  ! dry air density
-          zsl_height(i)    = pphi1(i)/RG      
-          tsoil0(i,:)      = tsoil(i,:)  
-          ustar(i)= (cdragm(i)*(wind_velo(i)**2))**0.5   
-       END DO
-       
-
-
-        dtis=dtime
-
-          IF (lafin) THEN
-            lafin_is=.true.
-          END IF
-
-          CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is,&
-            rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow,   &
-            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf,&
-            rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
-            radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow,agesno,   &
-            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
-            run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat,dflux_s, dflux_l, &
-            tsurf_new, alb1, alb2, alb3, alb6, &
-            emis_new, z0m, z0h, qsurf)
-
-          debut_is=.false.
+                u0, v0, u1, v1, gustiness, cdragm, &
+                AcoefU, AcoefV, BcoefU, BcoefV, &
+                p1lay, temp_air, &
+                flux_u1, flux_v1)
+
+
+        ! Set constants and compute some input for SISVAT
+        ! = 1000 hPa
+        ! and calculate incoming flux for SW and LW interval: swdown, lwdown
+        swdown(:) = 0.0
+        lwdown(:) = 0.0
+        snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model
+        alb_soil(:) = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
+        ustar(:) = 0.
+        pref = 100000.
+        DO i = 1, knon
+          swdown(i) = swnet(i) / (1 - albedo(i))
+          lwdown(i) = lwdownm(i)
+          wind_velo(i) = u1(i)**2 + v1(i)**2
+          wind_velo(i) = wind_velo(i)**0.5
+          pexner(i) = (p1lay(i) / pref)**(RD / RCPD)
+          dens_air(i) = p1lay(i) / RD / temp_air(i)  ! dry air density
+          zsl_height(i) = pphi1(i) / RG
+          tsoil0(i, :) = tsoil(i, :)
+          ustar(i) = (cdragm(i) * (wind_velo(i)**2))**0.5
+        END DO
+
+        dtis = dtime
+
+        IF (lafin) THEN
+          lafin_is = .true.
+        END IF
+
+        CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is, &
+                rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow, &
+                zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
+                rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
+                radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, &
+                AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
+                run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s, dflux_l, &
+                tsurf_new, alb1, alb2, alb3, alb6, &
+                emis_new, z0m, z0h, qsurf)
+
+        debut_is = .false.
 
 
@@ -325,64 +318,62 @@
 
         ! for consistency with standard LMDZ, add calving to run_off_lic
-        run_off_lic(:)=run_off_lic(:) + to_ice(:)
+        run_off_lic(:) = run_off_lic(:) + to_ice(:)
 
         DO i = 1, knon
-           ffonte_global(knindex(i),is_lic)    = ffonte(i)
-           fqfonte_global(knindex(i),is_lic)   = fqfonte(i)! net melting= melting - refreezing
-           fqcalving_global(knindex(i),is_lic) = to_ice(i) ! flux
-           runofflic_global(knindex(i)) = run_off_lic(i)
+          ffonte_global(knindex(i), is_lic) = ffonte(i)
+          fqfonte_global(knindex(i), is_lic) = fqfonte(i)! net melting= melting - refreezing
+          fqcalving_global(knindex(i), is_lic) = to_ice(i) ! flux
+          runofflic_global(knindex(i)) = run_off_lic(i)
         ENDDO
         ! Here, we assume that the calving term is equal to the to_ice term
         ! (no ice accumulation)
 
-
-#else
-       abort_message='Pb de coherence: landice_opt = 1  mais CPP_INLANDSIS = .false.'
-       CALL abort_physic(modname,abort_message,1)
-#endif
-
-
-    ELSE 
-
-!****************************************************************************************
-! Soil calculations
-
-!****************************************************************************************
-
-    ! EV: use calbeta
-    CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
-
-
-    ! use soil model and recalculate properly cal
-    IF (soil_model) THEN 
-       CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
-   longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
-       cal(1:knon) = RCPD / soilcap(1:knon)
-       radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
-    ELSE 
-       cal = RCPD * calice
-       WHERE (snow > 0.0) cal = RCPD * calsno
-    ENDIF
-
-
-!****************************************************************************************
-! Calulate fluxes
-
-!****************************************************************************************
-!    beta(:) = 1.0
-!    dif_grnd(:) = 0.0
-
-! Suppose zero surface speed
-    u0(:)=0.0
-    v0(:)=0.0
-    u1_lay(:) = u1(:) - u0(:)
-    v1_lay(:) = v1(:) - v0(:)
-
-    CALL calcul_fluxs(knon, is_lic, dtime, &
-         tsurf, p1lay, cal, beta, cdragh, cdragh, ps, &
-         precip_rain, precip_snow, snow, qsurf,  &
-         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
-         1.,AcoefH, AcoefQ, BcoefH, BcoefQ, &
-         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+      ELSE
+        abort_message = 'Pb de coherence: landice_opt = 1  mais CPP_INLANDSIS = .false.'
+        CALL abort_physic(modname, abort_message, 1)
+      END IF
+
+    ELSE
+
+      !****************************************************************************************
+      ! Soil calculations
+
+      !****************************************************************************************
+
+      ! EV: use calbeta
+      CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
+
+
+      ! use soil model and recalculate properly cal
+      IF (soil_model) THEN
+        CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
+                longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
+        cal(1:knon) = RCPD / soilcap(1:knon)
+        radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
+      ELSE
+        cal = RCPD * calice
+        WHERE (snow > 0.0) cal = RCPD * calsno
+      ENDIF
+
+
+      !****************************************************************************************
+      ! Calulate fluxes
+
+      !****************************************************************************************
+      !    beta(:) = 1.0
+      !    dif_grnd(:) = 0.0
+
+      ! Suppose zero surface speed
+      u0(:) = 0.0
+      v0(:) = 0.0
+      u1_lay(:) = u1(:) - u0(:)
+      v1_lay(:) = v1(:) - v0(:)
+
+      CALL calcul_fluxs(knon, is_lic, dtime, &
+              tsurf, p1lay, cal, beta, cdragh, cdragh, ps, &
+              precip_rain, precip_snow, snow, qsurf, &
+              radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
+              1., AcoefH, AcoefQ, BcoefH, BcoefQ, &
+              tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
 
 #ifdef ISO
@@ -411,210 +402,207 @@
 #endif         
 
-    CALL calcul_flux_wind(knon, dtime, &
-         u0, v0, u1, v1, gustiness, cdragm, &
-         AcoefU, AcoefV, BcoefU, BcoefV, &
-         p1lay, temp_air, &
-         flux_u1, flux_v1)
-
-
-!****************************************************************************************
-! Calculate albedo
-
-!****************************************************************************************
-
-!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
-!       alb1(1 : knon)  = 0.6 !IM cf FH/GK 
-!       alb1(1 : knon)  = 0.82
-!       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
-!       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
-!IM: KstaTER0.77 & LMD_ARMIP6    
-
-! Attantion: alb1 and alb2 are not the same!
-    alb1(1:knon)  = alb_vis_sno_lic
-    alb2(1:knon)  = alb_nir_sno_lic
-
-
-!****************************************************************************************
-! Rugosity
-
-!****************************************************************************************
-
-if (z0m_landice > 0.) then
-    z0m(1:knon) = z0m_landice
-    z0h(1:knon) = z0h_landice
-else
-    ! parameterization of z0=f(T) following measurements in Adelie Land by Amory et al 2018
-    coefa = 0.1658 !0.1862 !Ant
-    coefb = -50.3869 !-55.7718 !Ant
-    ta1 = 253.15 !255. Ant
-    ta2 = 273.15
-    ta3 = 273.15+3
-    z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm
-    z02 = exp(coefa*ta2 + coefb) !~6  !~7 mm
-    z03 = z01
-    coefc = log(z03/z02)/(ta3-ta2)
-    coefd = log(z03)-coefc*ta3
-    do j=1,knon 
-      if (temp_air(j) < ta1) then
-        z0m(j) = z01
-      else if (temp_air(j)>=ta1 .and. temp_air(j)<ta2) then
-        z0m(j) = exp(coefa*temp_air(j) + coefb)
-      else if (temp_air(j)>=ta2 .and. temp_air(j)<ta3) then
-        ! if st > 0, melting induce smooth surface
-        z0m(j) = exp(coefc*temp_air(j) + coefd)
+      CALL calcul_flux_wind(knon, dtime, &
+              u0, v0, u1, v1, gustiness, cdragm, &
+              AcoefU, AcoefV, BcoefU, BcoefV, &
+              p1lay, temp_air, &
+              flux_u1, flux_v1)
+
+
+      !****************************************************************************************
+      ! Calculate albedo
+
+      !****************************************************************************************
+
+      !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
+      !       alb1(1 : knon)  = 0.6 !IM cf FH/GK
+      !       alb1(1 : knon)  = 0.82
+      !       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
+      !       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
+      !IM: KstaTER0.77 & LMD_ARMIP6
+
+      ! Attantion: alb1 and alb2 are not the same!
+      alb1(1:knon) = alb_vis_sno_lic
+      alb2(1:knon) = alb_nir_sno_lic
+
+
+      !****************************************************************************************
+      ! Rugosity
+
+      !****************************************************************************************
+
+      if (z0m_landice > 0.) then
+        z0m(1:knon) = z0m_landice
+        z0h(1:knon) = z0h_landice
       else
-        z0m(j) = z03
+        ! parameterization of z0=f(T) following measurements in Adelie Land by Amory et al 2018
+        coefa = 0.1658 !0.1862 !Ant
+        coefb = -50.3869 !-55.7718 !Ant
+        ta1 = 253.15 !255. Ant
+        ta2 = 273.15
+        ta3 = 273.15 + 3
+        z01 = exp(coefa * ta1 + coefb) !~0.2 ! ~0.25 mm
+        z02 = exp(coefa * ta2 + coefb) !~6  !~7 mm
+        z03 = z01
+        coefc = log(z03 / z02) / (ta3 - ta2)
+        coefd = log(z03) - coefc * ta3
+        do j = 1, knon
+          if (temp_air(j) < ta1) then
+            z0m(j) = z01
+          else if (temp_air(j)>=ta1 .and. temp_air(j)<ta2) then
+            z0m(j) = exp(coefa * temp_air(j) + coefb)
+          else if (temp_air(j)>=ta2 .and. temp_air(j)<ta3) then
+            ! if st > 0, melting induce smooth surface
+            z0m(j) = exp(coefc * temp_air(j) + coefd)
+          else
+            z0m(j) = z03
+          endif
+          z0h(j) = z0m(j)
+        enddo
+
       endif
-      z0h(j)=z0m(j)
-    enddo
-
-endif    
- 
-
-!****************************************************************************************
-! Simple blowing snow param
-!****************************************************************************************
-! we proceed in 2 steps:
-! first we erode - if possible -the accumulated snow during the time step
-! then we update the density of the underlying layer and see if we can also erode
-! this layer
-
-
-   if (ok_bs) then
-       fluxbs(:)=0.
-       do j=1,knon
-          ws1(j)=(u1(j)**2+v1(j)**2)**0.5
-          ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
-          rhod(j)=p1lay(j)/RD/temp_air(j)
-          ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j))
-       enddo
-
-       ! 1st step: erosion of fresh snow accumulated during the time step
-       do j=1, knon
-       if (precip_snow(j) > 0.) then
-           rhos(j)=rhofresh_bs
-           ! blowing snow flux formula used in MAR
-           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
-           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 
-           ! computation of qbs at the top of the saltation layer
-           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
-           esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
-           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
-           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
-           ! calculation of erosion (flux positive towards the surface here) 
-           ! consistent with implicit resolution of turbulent mixing equation
-           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
-           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
-           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
-           ! (rho*qsalt*hsalt)
-           ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
-           maxerosion=min(precip_snow(j),hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs)
-
-           fluxbs_1(j)=rhod(j)*ws1(j)*cdragh(j)*zeta_bs*(AcoefQBS(j)-qsalt(j)) &
-                   / (1.-rhod(j)*ws1(j)*cdragh(j)*zeta_bs*BcoefQBS(j)*dtime)
-           fluxbs_1(j)=max(-maxerosion,fluxbs_1(j))
-
-           if (precip_snow(j) > abs(fluxbs_1(j))) then
-               ok_remaining_freshsnow(j)=.true.
-               bsweight_fresh(j)=1.
-           else
-               ok_remaining_freshsnow(j)=.false.
-               bsweight_fresh(j)=exp(-(abs(fluxbs_1(j))-precip_snow(j))/precip_snow(j))
-           endif
-       else
-           ok_remaining_freshsnow(j)=.false.
-           fluxbs_1(j)=0.
-           bsweight_fresh(j)=0.
-       endif
-       enddo
-
-
-       ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
-       ! this is done through the routine albsno
-       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs_1(:))
-
-       ! 2nd step:
-       ! computation of threshold friction velocity
-       ! which depends on surface snow density
-       do j = 1, knon
-        if (ok_remaining_freshsnow(j)) then
-           fluxbs_2(j)=0.
-        else
-           ! we start eroding the underlying layer
-           ! estimation of snow density
-           ! snow density increases with snow age and
-           ! increases even faster in case of sedimentation of blowing snow or rain
-           tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - & 
-                    abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.)))
-           rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)*(1.-exp(-agesno(j)*86400.0/tau_dens))
-           ! blowing snow flux formula used in MAR
-           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs)) 
-           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 
-           ! computation of qbs at the top of the saltation layer
-           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
-           esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
-           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
-           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
-           ! calculation of erosion (flux positive towards the surface here) 
-           ! consistent with implicit resolution of turbulent mixing equation
-           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
-           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
-           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
-           ! (rho*qsalt*hsalt)
-           maxerosion=hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs
-           fluxbs_2(j)=rhod(j)*ws1(j)*cdragh(j)*zeta_bs*(AcoefQBS(j)-qsalt(j)) &
-                   / (1.-rhod(j)*ws1(j)*cdragh(j)*zeta_bs*BcoefQBS(j)*dtime)
-           fluxbs_2(j)=max(-maxerosion,fluxbs_2(j))
-         endif
-       enddo
-
-
-
-
-       ! final flux and outputs       
-        do j=1, knon
-              ! total flux is the erosion of fresh snow + 
-              ! a fraction of the underlying snow (if all the fresh snow has been eroded)
-              ! the calculation of the fraction is quite delicate since we do not know
-              ! how much time was needed to erode the fresh snow. We assume that this time
-              ! is dt*exp(-(abs(fluxbs1)-precipsnow)/precipsnow)=dt*bsweight_fresh
-
-              fluxbs(j)=fluxbs_1(j)+fluxbs_2(j)*(1.-bsweight_fresh(j))
-              i = knindex(j)
-              zxustartlic(i) = ustart(j)
-              zxrhoslic(i) = rhos(j)
-              zxqsaltlic(i)=qsalt(j)
+
+
+      !****************************************************************************************
+      ! Simple blowing snow param
+      !****************************************************************************************
+      ! we proceed in 2 steps:
+      ! first we erode - if possible -the accumulated snow during the time step
+      ! then we update the density of the underlying layer and see if we can also erode
+      ! this layer
+
+      if (ok_bs) then
+        fluxbs(:) = 0.
+        do j = 1, knon
+          ws1(j) = (u1(j)**2 + v1(j)**2)**0.5
+          ustar(j) = (cdragm(j) * (u1(j)**2 + v1(j)**2))**0.5
+          rhod(j) = p1lay(j) / RD / temp_air(j)
+          ustart0(j) = (log(2.868) - log(1.625)) / 0.085 * sqrt(cdragm(j))
         enddo
 
-
-  else ! not ok_bs
-  ! those lines are useful to calculate the snow age
-       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
-
-  endif ! if ok_bs
-
-
-
-!****************************************************************************************
-! Calculate snow amount
-
-!****************************************************************************************
-    IF (ok_bs) THEN
-      precip_totsnow(:)=precip_snow(:)+precip_bs(:)
-      evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
-    ELSE
-      precip_totsnow(:)=precip_snow(:)
-      evap_totsnow(:)=evap(:)
-    ENDIF
-    
- 
-    CALL fonte_neige(knon, is_lic, knindex, dtime, &
-         tsurf, precip_rain, precip_totsnow, &
-         snow, qsol, tsurf_new, evap_totsnow &
+        ! 1st step: erosion of fresh snow accumulated during the time step
+        do j = 1, knon
+          if (precip_snow(j) > 0.) then
+            rhos(j) = rhofresh_bs
+            ! blowing snow flux formula used in MAR
+            ustart(j) = ustart0(j) * exp(max(rhoice_bs / rhofresh_bs - rhoice_bs / rhos(j), 0.)) * exp(max(0., rhos(j) - rhohard_bs))
+            ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
+            ! computation of qbs at the top of the saltation layer
+            ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
+            esalt = 1. / (c_esalt_bs * max(1.e-6, ustar(j)))
+            hsalt(j) = 0.08436 * (max(1.e-6, ustar(j))**1.27)
+            qsalt(j) = (max(ustar(j)**2 - ustart(j)**2, 0.)) / (RG * hsalt(j)) * esalt
+            ! calculation of erosion (flux positive towards the surface here)
+            ! consistent with implicit resolution of turbulent mixing equation
+            ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
+            ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
+            ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
+            ! (rho*qsalt*hsalt)
+            ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
+            maxerosion = min(precip_snow(j), hsalt(j) * qsalt(j) * rhod(j) / tau_eqsalt_bs)
+
+            fluxbs_1(j) = rhod(j) * ws1(j) * cdragh(j) * zeta_bs * (AcoefQBS(j) - qsalt(j)) &
+                    / (1. - rhod(j) * ws1(j) * cdragh(j) * zeta_bs * BcoefQBS(j) * dtime)
+            fluxbs_1(j) = max(-maxerosion, fluxbs_1(j))
+
+            if (precip_snow(j) > abs(fluxbs_1(j))) then
+              ok_remaining_freshsnow(j) = .true.
+              bsweight_fresh(j) = 1.
+            else
+              ok_remaining_freshsnow(j) = .false.
+              bsweight_fresh(j) = exp(-(abs(fluxbs_1(j)) - precip_snow(j)) / precip_snow(j))
+            endif
+          else
+            ok_remaining_freshsnow(j) = .false.
+            fluxbs_1(j) = 0.
+            bsweight_fresh(j) = 0.
+          endif
+        enddo
+
+
+        ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
+        ! this is done through the routine albsno
+        CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:) + fluxbs_1(:))
+
+        ! 2nd step:
+        ! computation of threshold friction velocity
+        ! which depends on surface snow density
+        do j = 1, knon
+          if (ok_remaining_freshsnow(j)) then
+            fluxbs_2(j) = 0.
+          else
+            ! we start eroding the underlying layer
+            ! estimation of snow density
+            ! snow density increases with snow age and
+            ! increases even faster in case of sedimentation of blowing snow or rain
+            tau_dens = max(tau_densmin_bs, tau_dens0_bs * exp(-abs(precip_bs(j)) / pbst_bs - &
+                    abs(precip_rain(j)) / prt_bs) * exp(-max(tsurf(j) - RTT, 0.)))
+            rhos(j) = rhofresh_bs + (rhohard_bs - rhofresh_bs) * (1. - exp(-agesno(j) * 86400.0 / tau_dens))
+            ! blowing snow flux formula used in MAR
+            ustart(j) = ustart0(j) * exp(max(rhoice_bs / rhofresh_bs - rhoice_bs / rhos(j), 0.)) * exp(max(0., rhos(j) - rhohard_bs))
+            ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
+            ! computation of qbs at the top of the saltation layer
+            ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
+            esalt = 1. / (c_esalt_bs * max(1.e-6, ustar(j)))
+            hsalt(j) = 0.08436 * (max(1.e-6, ustar(j))**1.27)
+            qsalt(j) = (max(ustar(j)**2 - ustart(j)**2, 0.)) / (RG * hsalt(j)) * esalt
+            ! calculation of erosion (flux positive towards the surface here)
+            ! consistent with implicit resolution of turbulent mixing equation
+            ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
+            ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
+            ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
+            ! (rho*qsalt*hsalt)
+            maxerosion = hsalt(j) * qsalt(j) * rhod(j) / tau_eqsalt_bs
+            fluxbs_2(j) = rhod(j) * ws1(j) * cdragh(j) * zeta_bs * (AcoefQBS(j) - qsalt(j)) &
+                    / (1. - rhod(j) * ws1(j) * cdragh(j) * zeta_bs * BcoefQBS(j) * dtime)
+            fluxbs_2(j) = max(-maxerosion, fluxbs_2(j))
+          endif
+        enddo
+
+
+
+
+        ! final flux and outputs
+        do j = 1, knon
+          ! total flux is the erosion of fresh snow +
+          ! a fraction of the underlying snow (if all the fresh snow has been eroded)
+          ! the calculation of the fraction is quite delicate since we do not know
+          ! how much time was needed to erode the fresh snow. We assume that this time
+          ! is dt*exp(-(abs(fluxbs1)-precipsnow)/precipsnow)=dt*bsweight_fresh
+
+          fluxbs(j) = fluxbs_1(j) + fluxbs_2(j) * (1. - bsweight_fresh(j))
+          i = knindex(j)
+          zxustartlic(i) = ustart(j)
+          zxrhoslic(i) = rhos(j)
+          zxqsaltlic(i) = qsalt(j)
+        enddo
+
+      else ! not ok_bs
+        ! those lines are useful to calculate the snow age
+        CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:))
+
+      endif ! if ok_bs
+
+
+
+      !****************************************************************************************
+      ! Calculate snow amount
+
+      !****************************************************************************************
+      IF (ok_bs) THEN
+        precip_totsnow(:) = precip_snow(:) + precip_bs(:)
+        evap_totsnow(:) = evap(:) - fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
+      ELSE
+        precip_totsnow(:) = precip_snow(:)
+        evap_totsnow(:) = evap(:)
+      ENDIF
+
+      CALL fonte_neige(knon, is_lic, knindex, dtime, &
+      tsurf, precip_rain, precip_totsnow, &
+      snow, qsol, tsurf_new, evap_totsnow &
 #ifdef ISO    
     ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag     &
     ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
 #endif
-     )
+      )
 
 
@@ -641,75 +629,73 @@
 
 #endif
-   
-    WHERE (snow(1 : knon) < 0.0001) agesno(1 : knon) = 0.
-    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))  
-
+
+      WHERE (snow(1:knon) < 0.0001) agesno(1:knon) = 0.
+      zfra(1:knon) = MAX(0.0, MIN(1.0, snow(1:knon) / (snow(1:knon) + 10.0)))
 
     END IF ! landice_opt
 
 
-!****************************************************************************************
-! Send run-off on land-ice to coupler if coupled ocean.
-! run_off_lic has been calculated in fonte_neige or surf_inlandsis
-! If landice_opt>=2, corresponding call is done from surf_land_orchidee
-!****************************************************************************************
+    !****************************************************************************************
+    ! Send run-off on land-ice to coupler if coupled ocean.
+    ! run_off_lic has been calculated in fonte_neige or surf_inlandsis
+    ! If landice_opt>=2, corresponding call is done from surf_land_orchidee
+    !****************************************************************************************
     IF (type_ocean=='couple' .AND. landice_opt < 2) THEN
-       ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
-       run_off_lic_frac(:)=0.0
-       DO j = 1, knon
-          i = knindex(j)
-          run_off_lic_frac(j) = pctsrf(i,is_lic)
-       ENDDO
-
-       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
+      ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
+      run_off_lic_frac(:) = 0.0
+      DO j = 1, knon
+        i = knindex(j)
+        run_off_lic_frac(j) = pctsrf(i, is_lic)
+      ENDDO
+
+      CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
     ENDIF
 
- ! transfer runoff rate [kg/m2/s](!) to physiq for output
-    runoff(1:knon)=run_off_lic(1:knon)/dtime
-
-       snow_o=0.
-       zfra_o = 0.
-       DO j = 1, knon
-           i = knindex(j)
-           snow_o(i) = snow(j)
-           zfra_o(i) = zfra(j)
-       ENDDO
-
-
-!albedo SB >>>
-     select case(NSW)
-     case(2)
-       alb_dir(1:knon,1)=alb1(1:knon)
-       alb_dir(1:knon,2)=alb2(1:knon)
-     case(4)
-       alb_dir(1:knon,1)=alb1(1:knon)
-       alb_dir(1:knon,2)=alb2(1:knon)
-       alb_dir(1:knon,3)=alb2(1:knon)
-       alb_dir(1:knon,4)=alb2(1:knon)
-     case(6)
-       alb_dir(1:knon,1)=alb1(1:knon)
-       alb_dir(1:knon,2)=alb1(1:knon)
-       alb_dir(1:knon,3)=alb1(1:knon)
-       alb_dir(1:knon,4)=alb2(1:knon)
-       alb_dir(1:knon,5)=alb2(1:knon)
-       alb_dir(1:knon,6)=alb2(1:knon)
-
-       IF ((landice_opt == 1) .AND. (iflag_albcalc == 2)) THEN
-       alb_dir(1:knon,1)=alb6(1:knon,1)
-       alb_dir(1:knon,2)=alb6(1:knon,2)
-       alb_dir(1:knon,3)=alb6(1:knon,3)
-       alb_dir(1:knon,4)=alb6(1:knon,4)
-       alb_dir(1:knon,5)=alb6(1:knon,5)
-       alb_dir(1:knon,6)=alb6(1:knon,6)
-       ENDIF
-
-     end select
-alb_dif=alb_dir
-!albedo SB <<<
-
+    ! transfer runoff rate [kg/m2/s](!) to physiq for output
+    runoff(1:knon) = run_off_lic(1:knon) / dtime
+
+    snow_o = 0.
+    zfra_o = 0.
+    DO j = 1, knon
+      i = knindex(j)
+      snow_o(i) = snow(j)
+      zfra_o(i) = zfra(j)
+    ENDDO
+
+
+    !albedo SB >>>
+    select case(NSW)
+    case(2)
+      alb_dir(1:knon, 1) = alb1(1:knon)
+      alb_dir(1:knon, 2) = alb2(1:knon)
+    case(4)
+      alb_dir(1:knon, 1) = alb1(1:knon)
+      alb_dir(1:knon, 2) = alb2(1:knon)
+      alb_dir(1:knon, 3) = alb2(1:knon)
+      alb_dir(1:knon, 4) = alb2(1:knon)
+    case(6)
+      alb_dir(1:knon, 1) = alb1(1:knon)
+      alb_dir(1:knon, 2) = alb1(1:knon)
+      alb_dir(1:knon, 3) = alb1(1:knon)
+      alb_dir(1:knon, 4) = alb2(1:knon)
+      alb_dir(1:knon, 5) = alb2(1:knon)
+      alb_dir(1:knon, 6) = alb2(1:knon)
+
+      IF ((landice_opt == 1) .AND. (iflag_albcalc == 2)) THEN
+        alb_dir(1:knon, 1) = alb6(1:knon, 1)
+        alb_dir(1:knon, 2) = alb6(1:knon, 2)
+        alb_dir(1:knon, 3) = alb6(1:knon, 3)
+        alb_dir(1:knon, 4) = alb6(1:knon, 4)
+        alb_dir(1:knon, 5) = alb6(1:knon, 5)
+        alb_dir(1:knon, 6) = alb6(1:knon, 6)
+      ENDIF
+
+    end select
+    alb_dif = alb_dir
+    !albedo SB <<<
 
   END SUBROUTINE surf_landice
 
-!****************************************************************************************
+  !****************************************************************************************
 
 END MODULE surf_landice_mod
