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)

      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=871
      integer,save :: lunout=6
      integer,save :: lev_out=0

      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)
      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
      real 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)
      save f0
      real zlevinter(klon)
      logical debut

!

      character*2 str2
      character*10 str10

      EXTERNAL SCOPY
!

!-----------------------------------------------------------------------
!   initialisation:
!   ---------------
!
      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      
      if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
            f0(ig)=1.e-5
            zmax0(ig)=40.
      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'

!------------------------------------------------------------------
!   Calcul de w2, carre de w a partir de la cape
!
!   Indicages:
!   l'ascendance provenant du niveau k traverse l'interface l avec
!   une vitesse wa(k,l).
!
!                       --------------------
!
!                       + + + + + + + + + + 
!
!  wa(k,l)   ----       --------------------    l
!             /\
!            /||\       + + + + + + + + + + 
!             ||
!             ||        --------------------
!             ||
!             ||        + + + + + + + + + + 
!             ||
!             ||        --------------------
!             ||__
!             |___      + + + + + + + + + +     k
!
!                       --------------------
!
!
!
!------------------------------------------------------------------
!definition du profil d alimentation a partir de la flottabilite: 
!alim_star, alim_star_tot, lalim et lmin
!------------------------------------------------------------------
!
      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)

      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)
      endif

!-----------------------------------------------------------------------------------
!calcul des caracteristiques du thermique sec pour le calcul de detr et la fermeture
!-----------------------------------------------------------------------------------
!
      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
     &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)

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



!---------------------------------------------------------------------------------
!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)

      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,lmax(igout))
      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)

      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'

!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)

      if (1.eq.0) 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,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
      do ig=1,klon
         fraca(ig,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 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
      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
!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      
      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(q2(ig,l))/(po(ig,l)*1000.)
         enddo
      enddo
!calcul du ratqscdiff
      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
      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
      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

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