!
! $Header$
!
      SUBROUTINE thermcell_main(ngrid,nlay,ptimestep  &
     &                  ,pplay,pplev,pphi,debut  &
     &                  ,pu,pv,pt,po  &
     &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
     &                  ,fm0,entr0,zqla,lmax  &
     &                  ,ratqscth,ratqsdiff,zqsatth  &
     &                  ,r_aspect,l_mix,w2di,tho &
     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
     &                  ,zmax0, f0)

      IMPLICIT NONE

!=======================================================================
!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
!   Version du 09.02.07
!   Calcul du transport vertical dans la couche limite en presence
!   de "thermiques" explicitement representes avec processus nuageux
!
!   Réécriture à partir d'un listing papier à Habas, le 14/02/00
!
!   le thermique est supposé homogène et dissipé par mélange avec
!   son environnement. la longueur l_mix contrôle l'efficacité du
!   mélange
!
!   Le calcul du transport des différentes espèces se fait en prenant
!   en compte:
!     1. un flux de masse montant
!     2. un flux de masse descendant
!     3. un entrainement
!     4. un detrainement
!
!=======================================================================

!-----------------------------------------------------------------------
!   declarations:
!   -------------

#include "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"

!   arguments:
!   ----------

      INTEGER ngrid,nlay,w2di,tho
      real ptimestep,l_mix,r_aspect
      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
      real pphi(ngrid,nlay)

!   local:
!   ------

      integer,save :: igout=1
      integer,save :: lunout=6
      integer,save :: lev_out=10

      INTEGER ig,k,l,ll
      real zsortie1d(klon)
      INTEGER lmax(klon),lmin(klon),lalim(klon)
      INTEGER lmix(klon)
      real linter(klon)
      real zmix(klon)
      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev)
      real zmax_sec(klon)
      real w_est(klon,klev+1)
!on garde le zmax du pas de temps precedent
      real zmax0(klon)
!FH/IM     save zmax0

      real zlev(klon,klev+1),zlay(klon,klev)
      real deltaz(klon,klev)
      REAL zh(klon,klev),zdhadj(klon,klev)
      real zthl(klon,klev),zdthladj(klon,klev)
      REAL ztv(klon,klev)
      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
      real zl(klon,klev)
      real zsortie(klon,klev)
      real zva(klon,klev)
      real zua(klon,klev)
      real zoa(klon,klev)

      real zta(klon,klev)
      real zha(klon,klev)
      real fraca(klon,klev+1)
      real zf,zf2
      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
      real q2(klon,klev)
      common/comtherm/thetath2,wth2
    
      real ratqscth(klon,klev)
      real var
      real vardiff
      real ratqsdiff(klon,klev)
      integer isplit,nsplit
      parameter (nsplit=10)
      data isplit/0/
      save isplit

      logical sorties
      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
      real zpspsk(klon,klev)

      real wmax(klon)
      real wmax_sec(klon)
      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
      real detr0(klon,klev)
      real fm(klon,klev+1),entr(klon,klev)

      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
!niveau de condensation
      integer nivcon(klon)
      real zcon(klon)
      REAL CHI
      real zcon2(klon)
      real pcon(klon)
      real zqsat(klon,klev)
      real zqsatth(klon,klev) 

      real f_star(klon,klev+1),entr_star(klon,klev)
      real detr_star(klon,klev)
      real alim_star_tot(klon),alim_star2(klon)
      real alim_star(klon,klev)
      real f(klon), f0(klon)
!FH/IM     save f0
      real zlevinter(klon)
      logical debut
       real seuil

!
      !nouvelles variables pour la convection
      real Ale_bl(klon)
      real Alp_bl(klon)
      real alp_int(klon)
      real ale_int(klon)
      integer n_int(klon)
      real fm_tot(klon)
      real wght_th(klon,klev)
      integer lalim_conv(klon)
!v1d     logical therm
!v1d     save therm

      character*2 str2
      character*10 str10

      EXTERNAL SCOPY
!

!-----------------------------------------------------------------------
!   initialisation:
!   ---------------
!

      seuil=0.25

      if (lev_out.ge.1) print*,'thermcell_main V4'

       sorties=.true.
      IF(ngrid.NE.klon) THEN
         PRINT*
         PRINT*,'STOP dans convadj'
         PRINT*,'ngrid    =',ngrid
         PRINT*,'klon  =',klon
      ENDIF
!
!Initialisation
!
      do ig=1,klon      
!FH/IM 130308     if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
      if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
            f0(ig)=1.e-5
            zmax0(ig)=40.
!v1d        therm=.false.
      endif
      enddo 


!-----------------------------------------------------------------------
! Calcul de T,q,ql a partir de Tl et qT dans l environnement
!   --------------------------------------------------------------------
!
      CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
       
      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_env'

!------------------------------------------------------------------------
!                       --------------------
!
!
!                       + + + + + + + + + + +
!
!
!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
!  wh,wt,wo ...
!
!                       + + + + + + + + + + +  zh,zu,zv,zo,rho
!
!
!                       --------------------   zlev(1)
!                       \\\\\\\\\\\\\\\\\\\\
!
!

!-----------------------------------------------------------------------
!   Calcul des altitudes des couches
!-----------------------------------------------------------------------

      do l=2,nlay
         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
      enddo
         zlev(:,1)=0.
         zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
      do l=1,nlay
         zlay(:,l)=pphi(:,l)/RG
      enddo
!calcul de l epaisseur des couches
      do l=1,nlay
         deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
      enddo

!     print*,'2 OK convect8'
!-----------------------------------------------------------------------
!   Calcul des densites
!-----------------------------------------------------------------------

      do l=1,nlay
         rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l))
      enddo

      do l=2,nlay
         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
      enddo

!calcul de la masse
      do l=1,nlay
         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
      enddo

      if (lev_out.ge.1) print*,'thermcell_main apres initialisation'

!------------------------------------------------------------------
!
!             /|\
!    --------  |  F_k+1 -------   
!                              ----> D_k
!             /|\              <---- E_k , A_k
!    --------  |  F_k --------- 
!                              ----> D_k-1
!                              <---- E_k-1 , A_k-1
!
!
!
!
!
!    ---------------------------
!
!    ----- F_lmax+1=0 ----------         \
!            lmax     (zmax)              |
!    ---------------------------          |
!                                         |
!    ---------------------------          |
!                                         |
!    ---------------------------          |
!                                         |
!    ---------------------------          |
!                                         |
!    ---------------------------          |
!                                         |  E
!    ---------------------------          |  D
!                                         |
!    ---------------------------          |
!                                         |
!    ---------------------------  \       |
!            lalim                 |      |
!    ---------------------------   |      |
!                                  |      |
!    ---------------------------   |      |
!                                  | A    |
!    ---------------------------   |      |
!                                  |      |
!    ---------------------------   |      |
!    lmin  (=1 pour le moment)     |      |
!    ----- F_lmin=0 ------------  /      /
!
!    ---------------------------
!    //////////////////////////
!
!
!=============================================================================
!  Calculs initiaux ne faisant pas intervenir les changements de phase
!=============================================================================

!------------------------------------------------------------------
!  1. alim_star est le profil vertical de l'alimentation à la base du
!     panache thermique, calculé à partir de la flotabilité de l'air sec
!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
!------------------------------------------------------------------
!
      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
      CALL thermcell_init(ngrid,nlay,ztv,zlev,  &
     &                    lalim,lmin,alim_star,alim_star_tot,lev_out)

call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin  ')
call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')


      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_init'
      if (lev_out.ge.10) then
         write(lunout,*) 'Dans thermcell_main 1'
         write(lunout,*) 'lmin ',lmin(igout)
         write(lunout,*) 'lalim ',lalim(igout)
         write(lunout,*) ' ig l alim_star thetav'
         write(lunout,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
     &   ,ztv(igout,l),l=1,lalim(igout)+4)
      endif

      do ig=1,klon
!v1d     if (alim_star(ig,1).gt.1.e-10) then
!v1d     therm=.true.
!v1d     endif
      enddo
!-----------------------------------------------------------------------------
!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
!     panache sec conservatif (e=d=0) alimente selon alim_star 
!     Il s'agit d'un calcul de type CAPE
!     zmax_sec est utilisé pour déterminer la géométrie du thermique.
!------------------------------------------------------------------------------
!
      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
     &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)

call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')

      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_dry'
      if (lev_out.ge.10) then
         write(lunout,*) 'Dans thermcell_main 1b'
         write(lunout,*) 'lmin ',lmin(igout)
         write(lunout,*) 'lalim ',lalim(igout)
         write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
         write(lunout,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
     &    ,l=1,lalim(igout)+4)
      endif



!---------------------------------------------------------------------------------
!calcul du melange et des variables dans le thermique
!--------------------------------------------------------------------------------
!
      CALL thermcell_plume(ngrid,nlay,ztv,zthl,po,zl,rhobarz,  &
     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
     &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')

      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_plume'
      if (lev_out.ge.10) then
         write(lunout,*) 'Dans thermcell_main 2'
         write(lunout,*) 'lmin ',lmin(igout)
         write(lunout,*) 'lalim ',lalim(igout)
         write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
         write(lunout,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
      endif

!-------------------------------------------------------------------------------
! Calcul des caracteristiques du thermique:zmax,zmix,wmax
!-------------------------------------------------------------------------------
!
      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)


      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')

      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_height'

!-------------------------------------------------------------------------------
! Fermeture,determination de f
!-------------------------------------------------------------------------------

      CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out)

      if(lev_out.ge.1)print*,'thermcell_closure apres thermcell_closure'

!-------------------------------------------------------------------------------
!deduction des flux
!-------------------------------------------------------------------------------

      CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
     &       lalim,lmax,alim_star,  &
     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
     &       detr,zqla,zmax,lev_out,lunout,igout)

      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_flux'
      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')

!c------------------------------------------------------------------
!   calcul du transport vertical
!------------------------------------------------------------------

      if (w2di.eq.1) then
         fm0=fm0+ptimestep*(fm-fm0)/float(tho)
         entr0=entr0+ptimestep*(entr-entr0)/float(tho)
      else
         fm0=fm
         entr0=entr
         detr0=detr
      endif

      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
     &                    zthl,zdthladj,zta,lev_out)
      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
     &                   po,pdoadj,zoa,lev_out)

!------------------------------------------------------------------
! Calcul de la fraction de l'ascendance
!------------------------------------------------------------------
      do ig=1,klon
         fraca(ig,1)=0.
         fraca(ig,nlay+1)=0.
      enddo
      do l=2,nlay
         do ig=1,klon
            if (zw2(ig,l).gt.1.e-10) then
            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
            else
            fraca(ig,l)=0.
            endif
         enddo
      enddo
     
!------------------------------------------------------------------
!  calcul du transport vertical du moment horizontal
!------------------------------------------------------------------

      if (1.eq.1) then


! Calcul du transport de V tenant compte d'echange par gradient
! de pression horizontal avec l'environnement

         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
     &    ,fraca,zmax  &
     &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
      else

! calcul purement conservatif pour le transport de V
         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
     &    ,zu,pduadj,zua,lev_out)
         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
     &    ,zv,pdvadj,zva,lev_out)
      endif

!     print*,'13 OK convect8'
      do l=1,nlay
         do ig=1,ngrid
           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)  
         enddo
      enddo

      print*,'14 OK convect8'
!------------------------------------------------------------------
!   Calculs de diagnostiques pour les sorties
!------------------------------------------------------------------
!calcul de fraca pour les sorties
      
      if (sorties) then
      print*,'14a OK convect8'
! calcul du niveau de condensation
! initialisation
      do ig=1,ngrid
         nivcon(ig)=0
         zcon(ig)=0.
      enddo 
!nouveau calcul
      do ig=1,ngrid
      CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
      pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
      enddo
      do k=1,nlay
         do ig=1,ngrid
         if ((pcon(ig).le.pplay(ig,k))  &
     &      .and.(pcon(ig).gt.pplay(ig,k+1))) then
            zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
         endif
         enddo
      enddo
      print*,'14b OK convect8'
      do k=nlay,1,-1
         do ig=1,ngrid
            if (zqla(ig,k).gt.1e-10) then
               nivcon(ig)=k
               zcon(ig)=zlev(ig,k)
            endif
         enddo
      enddo
      print*,'14c OK convect8'
!calcul des moments
!initialisation
      do l=1,nlay
         do ig=1,ngrid
            q2(ig,l)=0.
            wth2(ig,l)=0.
            wth3(ig,l)=0.
            ratqscth(ig,l)=0.
            ratqsdiff(ig,l)=0.
         enddo
      enddo      
      print*,'14d OK convect8'
      do l=1,nlay
         do ig=1,ngrid
            zf=fraca(ig,l)
            zf2=zf/(1.-zf)
            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
            wth2(ig,l)=zf2*(zw2(ig,l))**2
!           print*,'wth2=',wth2(ig,l)
            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
!test: on calcul q2/po=ratqsc
            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
         enddo
      enddo
!calcul de ale_bl et alp_bl
!pour le calcul d'une valeur intégrée entre la surface et lmax
      do ig=1,ngrid
      alp_int(ig)=0.
      ale_int(ig)=0.
      n_int(ig)=0
      enddo
      do ig=1,ngrid
!      do l=nivcon(ig),lmax(ig)
      do l=1,lmax(ig)
      alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
      ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
      n_int(ig)=n_int(ig)+1
      enddo
      enddo
!      print*,'avant calcul ale et alp' 
!calcul de ALE et ALP pour la convection
      do ig=1,ngrid
!      Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig))
!          Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig))
!      Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig)) 
!     &           *0.1
!valeur integree de alp_bl * 0.5:
       if (n_int(ig).gt.0) then
       Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
!       if (Alp_bl(ig).lt.0.) then
!       Alp_bl(ig)=0.
       endif
!       endif
!         write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
!     s               wth3(ig,nivcon(ig)),Alp_bl(ig)
!            write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
!      Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
!      if (nivcon(ig).eq.1) then
!       Ale_bl(ig)=0.
!       else
!valeur max de ale_bl:
       Ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2 
!     & /2.
!     & *0.1
!        Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2 
!       if (n_int(ig).gt.0) then
!       Ale_bl(ig)=ale_int(ig)/n_int(ig)
!        Ale_bl(ig)=4.
!       endif
!       endif
!            Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
!          Ale_bl(ig)=wth2(ig,nivcon(ig))
!          write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
      enddo
!test:calcul de la ponderation des couches pour KE
!initialisations
!      print*,'ponderation'
      do ig=1,ngrid
           fm_tot(ig)=0.
      enddo
       do ig=1,ngrid
        do k=1,klev
           wght_th(ig,k)=1.
        enddo
       enddo
       do ig=1,ngrid
!         lalim_conv(ig)=lmix_bis(ig)
!la hauteur de la couche alim_conv = hauteur couche alim_therm
         lalim_conv(ig)=lalim(ig)
!         zentr(ig)=zlev(ig,lalim(ig))
      enddo
      do ig=1,ngrid
        do k=1,lalim_conv(ig)
           fm_tot(ig)=fm_tot(ig)+fm(ig,k)
        enddo
      enddo
      do ig=1,ngrid
        do k=1,lalim_conv(ig)
           if (fm_tot(ig).gt.1.e-10) then
!           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
           endif
!on pondere chaque couche par a*
             if (alim_star(ig,k).gt.1.e-10) then
             wght_th(ig,k)=alim_star(ig,k)
             else
             wght_th(ig,k)=1.
             endif
        enddo
      enddo
!      print*,'apres wght_th'
!test pour prolonger la convection
      do ig=1,ngrid
!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
      if ((alim_star(ig,1).lt.1.e-10)) then
      lalim_conv(ig)=1
      wght_th(ig,1)=1.
!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
      endif
      enddo

!calcul du ratqscdiff
      print*,'14e OK convect8'
      var=0.
      vardiff=0.
      ratqsdiff(:,:)=0.
      do ig=1,ngrid
         do l=1,lalim(ig)
            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
         enddo
      enddo
      print*,'14f OK convect8'
      do ig=1,ngrid
          do l=1,lalim(ig)
          zf=fraca(ig,l)
          zf2=zf/(1.-zf)
          vardiff=vardiff+alim_star(ig,l)  &
     &           *(zqta(ig,l)*1000.-var)**2
!         ratqsdiff=ratqsdiff+alim_star(ig,l)*
!     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
          enddo
      enddo
      print*,'14g OK convect8'
      do l=1,nlay
         do ig=1,ngrid
            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
         enddo
      enddo 
!--------------------------------------------------------------------    
!
!ecriture des fichiers sortie
!     print*,'15 OK convect8'

      isplit=isplit+1       


#ifdef und
      if (lev_out.ge.1) print*,'thermcell_main sorties 1D'
#include "thermcell_out1d.h"
#endif


! #define troisD
      if (lev_out.ge.1) print*,'thermcell_main sorties 3D'
#ifdef troisD
#include "thermcell_out3d.h"
#endif

      endif

      if (lev_out.ge.1) print*,'thermcell_main FIN  OK'

      return
      end

!-----------------------------------------------------------------------------

      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)

      integer klon,klev
      real pplev(klon,klev+1),pplay(klon,klev)
      real ztv(klon,klev)
      real po(klon,klev)
      real ztva(klon,klev)
      real zqla(klon,klev)
      real f_star(klon,klev)
      real zw2(klon,klev)
      integer long(klon)
      real seuil
      character*21 comment

      print*,'WARNING !!! TEST ',comment
      return

!  test sur la hauteur des thermiques ...
         do i=1,klon
            if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
               do k=1,klev
                  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
               enddo
!              stop
            endif
         enddo


      return
      end