source: LMDZ5/trunk/libf/phylmd/thermcell_main.F90 @ 1851

!
! $Id: thermcell_main.F90 1790 2013-07-17 09:24:04Z emillour $
!
      SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep  &
     &                  ,pplay,pplev,pphi,debut  &
     &                  ,pu,pv,pt,po  &
     &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
     &                  ,fm0,entr0,detr0,zqta,zqla,lmax  &
     &                  ,ratqscth,ratqsdiff,zqsatth  &
     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
     &                  ,zmax0, f0,zw2,fraca,ztv &
     &                  ,zpspsk,ztla,zthl &
!!! nrlmd le 10/04/2012
     &                  ,pbl_tke,pctsrf,omega,airephy &
     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
     &                  ,n2,s2,ale_bl_stat &
     &                  ,therm_tke_max,env_tke_max &
     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
     &                  ,alp_bl_conv,alp_bl_stat &
!!! fin nrlmd le 10/04/2012
     &                  ,ztva  )

      USE dimphy
      USE ioipsl
      USE comgeomphy , ONLY:rlond,rlatd
      USE indice_sol_mod
      IMPLICIT NONE

!=======================================================================
!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
!   Version du 09.02.07
!   Calcul du transport vertical dans la couche limite en presence
!   de "thermiques" explicitement representes avec processus nuageux
!
!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
!
!   le thermique est suppose homogene et dissipe par melange avec
!   son environnement. la longueur l_mix controle l'efficacite du
!   melange
!
!   Le calcul du transport des differentes especes se fait en prenant
!   en compte:
!     1. un flux de masse montant
!     2. un flux de masse descendant
!     3. un entrainement
!     4. un detrainement
!
! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
!    Introduction of an implicit computation of vertical advection in
!    the environment of thermal plumes in thermcell_dq
!    impl =     0 : explicit, 1 : implicit, -1 : old version
!    controled by iflag_thermals =
!       15, 16 run with impl=-1 : numerical convergence with NPv3
!       17, 18 run with impl=1  : more stable
!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
!
!=======================================================================


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

#include "dimensions.h"
#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"
#include "iniprint.h"
#include "thermcell.h"

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

!IM 140508
      INTEGER itap

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

!   local:
!   ------

      integer icount

      integer, save :: dvdq=1,dqimpl=-1
!$OMP THREADPRIVATE(dvdq,dqimpl)
      data icount/0/
      save icount
!$OMP THREADPRIVATE(icount)

      integer,save :: igout=1
!$OMP THREADPRIVATE(igout)
      integer,save :: lunout1=6
!$OMP THREADPRIVATE(lunout1)
      integer,save :: lev_out=10
!$OMP THREADPRIVATE(lev_out)

      REAL susqr2pi, Reuler

      INTEGER ig,k,l,ll,ierr
      real zsortie1d(klon)
      INTEGER lmax(klon),lmin(klon),lalim(klon)
      INTEGER lmix(klon)
      INTEGER lmix_bis(klon)
      real linter(klon)
      real zmix(klon)
      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
!      real fraca(klon,klev)

      real zmax_sec(klon)
!on garde le zmax du pas de temps precedent
      real zmax0(klon)
!FH/IM     save zmax0

      real lambda

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

      real zta(klon,klev)
      real zha(klon,klev)
      real fraca(klon,klev+1)
      real zf,zf2
      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
      real q2(klon,klev)
! FH probleme de dimensionnement avec l'allocation dynamique
!     common/comtherm/thetath2,wth2
      real wq(klon,klev)
      real wthl(klon,klev)
      real wthv(klon,klev)
    
      real ratqscth(klon,klev)
      real var
      real vardiff
      real ratqsdiff(klon,klev)

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

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

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

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

!!! nrlmd le 10/04/2012

!------Entrées
      real pbl_tke(klon,klev+1,nbsrf)
      real pctsrf(klon,nbsrf)
      real omega(klon,klev)
      real airephy(klon)
!------Sorties
      real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
      real therm_tke_max0(klon),env_tke_max0(klon)
      real n2(klon),s2(klon)
      real ale_bl_stat(klon)
      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
!------Local
      integer nsrf
      real rhobarz0(klon)                    ! Densité au LCL
      logical ok_lcl(klon)                   ! Existence du LCL des thermiques
      integer klcl(klon)                     ! Niveau du LCL
      real interp(klon)                      ! Coef d'interpolation pour le LCL
!--Triggering
      real Su                                ! Surface unité: celle d'un updraft élémentaire
      parameter(Su=4e4)
      real hcoef                             ! Coefficient directeur pour le calcul de s2
      parameter(hcoef=1)
      real hmincoef                          ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 
      parameter(hmincoef=0.3)
      real eps1                              ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
      parameter(eps1=0.3)
      real hmin(ngrid)                       ! Ordonnée à l'origine pour le calcul de s2
      real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
      real zmax_moy_coef
      parameter(zmax_moy_coef=0.33)
      real depth(klon)                       ! Epaisseur moyenne du cumulus
      real w_max(klon)                       ! Vitesse max statistique 
      real s_max(klon)
!--Closure
      real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne 
      real pbl_tke_max0(klon)                ! TKE moyenne au LCL
      real w_ls(klon,klev)                   ! Vitesse verticale grande échelle (m/s)
      real coef_m                            ! On considère un rendement pour alp_bl_fluct_m
      parameter(coef_m=1.)
      real coef_tke                          ! On considère un rendement pour alp_bl_fluct_tke
      parameter(coef_tke=1.)

!!! fin nrlmd le 10/04/2012

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

      character*2 str2
      character*10 str10

      character (len=20) :: modname='thermcell_main'
      character (len=80) :: abort_message

      EXTERNAL SCOPY
!

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

      seuil=0.25

      if (debut)  then
!        call getin('dvdq',dvdq)
!        call getin('dqimpl',dqimpl)

         if (iflag_thermals==15.or.iflag_thermals==16) then
            dvdq=0
            dqimpl=-1
         else
            dvdq=1
            dqimpl=1
         endif

         fm0=0.
         entr0=0.
         detr0=0.


#undef wrgrads_thermcell
#ifdef wrgrads_thermcell
! Initialisation des sorties grads pour les thermiques.
! Pour l'instant en 1D sur le point igout.
! Utilise par thermcell_out3d.h
         str10='therm'
         call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
     &   rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1.,   &
     &   ptimestep,str10,'therm ')
#endif



      endif

      fm=0. ; entr=0. ; detr=0.


      icount=icount+1

!IM 090508 beg
!print*,'====================================================================='
!print*,'====================================================================='
!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
!print*,'====================================================================='
!print*,'====================================================================='
!IM 090508 end

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

       sorties=.true.
      IF(ngrid.NE.klon) THEN
         PRINT*
         PRINT*,'STOP dans convadj'
         PRINT*,'ngrid    =',ngrid
         PRINT*,'klon  =',klon
      ENDIF
!
!     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
     do ig=1,klon
         f0(ig)=max(f0(ig),1.e-2)
         zmax0(ig)=max(zmax0(ig),40.)
!IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
     enddo

      if (prt_level.ge.20) then
       do ig=1,ngrid
          print*,'th_main ig f0',ig,f0(ig)
       enddo
      endif
!-----------------------------------------------------------------------
! 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 (prt_level.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
!-----------------------------------------------------------------------

     rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))

     if (prt_level.ge.10)write(lunout,*)                                &
    &    'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
      rhobarz(:,1)=rho(:,1)

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

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

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

!------------------------------------------------------------------
!
!             /|\
!    --------  |  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 a la base du
!     panache thermique, calcule a partir de la flotabilite 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.
      lmin=1

!-----------------------------------------------------------------------------
!  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 utilise pour determiner la geometrie du thermique.
!------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
!calcul du melange et des variables dans le thermique
!--------------------------------------------------------------------------------
!
      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
!IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &

! Gestion temporaire de plusieurs appels à thermcell_plume au travers
! de la variable iflag_thermals

!      print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
      if (iflag_thermals_ed<=9) then
!         print*,'THERM NOUVELLE/NOUVELLE Arnaud'
         CALL 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,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     &    ,lev_out,lunout1,igout)

      elseif (iflag_thermals_ed>9) then
!        print*,'THERM RIO et al 2010, version d Arnaud'
         CALL thermcellV1_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,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     &    ,lev_out,lunout1,igout)

      endif

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

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

      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
      if (prt_level.ge.10) then
         write(lunout1,*) 'Dans thermcell_main 2'
         write(lunout1,*) 'lmin ',lmin(igout)
         write(lunout1,*) 'lalim ',lalim(igout)
         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
         write(lunout1,'(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)
! Attention, w2 est transforme en sa racine carree dans cette routine
! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
      wmax_tmp=0.
      do  l=1,nlay
         wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
      enddo
!     print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax



      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 (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'

!-------------------------------------------------------------------------------
! Fermeture,determination de f
!-------------------------------------------------------------------------------
!
!
!!      write(lunout,*)'THERM NOUVEAU XXXXX'
      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 (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
      if (prt_level.ge.10) then
         write(lunout1,*) 'Dans thermcell_main 1b'
         write(lunout1,*) 'lmin ',lmin(igout)
         write(lunout1,*) 'lalim ',lalim(igout)
         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
         write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
     &    ,l=1,lalim(igout)+4)
      endif




! Choix de la fonction d'alimentation utilisee pour la fermeture.
! Apparemment sans importance
      alim_star_clos(:,:)=alim_star(:,:)
      alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)

! Appel avec la version seche
      CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &
     &   zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Appel avec les zmax et wmax tenant compte de la condensation
! Semble moins bien marcher
!     CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &
!    &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

      if (tau_thermals>1.) then
         lambda=exp(-ptimestep/tau_thermals)
         f0=(1.-lambda)*f+lambda*f0
      else
         f0=f
      endif

! Test valable seulement en 1D mais pas genant
      if (.not. (f0(1).ge.0.) ) then
              abort_message = '.not. (f0(1).ge.0.)'
              CALL abort_gcm (modname,abort_message,1)
      endif

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

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

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

!------------------------------------------------------------------
!   On ne prend pas directement les profils issus des calculs precedents
!   mais on s'autorise genereusement une relaxation vers ceci avec
!   une constante de temps tau_thermals (typiquement 1800s).
!------------------------------------------------------------------

      if (tau_thermals>1.) then
         lambda=exp(-ptimestep/tau_thermals)
         fm0=(1.-lambda)*fm+lambda*fm0
         entr0=(1.-lambda)*entr+lambda*entr0
         detr0=(1.-lambda)*detr+lambda*detr0
      else
         fm0=fm
         entr0=entr
         detr0=detr
      endif

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

      call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
     &                    zthl,zdthladj,zta,lev_out)
      call thermcell_dq(ngrid,nlay,dqimpl,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
!------------------------------------------------------------------

!IM 090508  
      if (dvdq == 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*dvdq,zmax &
     &    ,fraca,zmax &
     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)

      else

! calcul purement conservatif pour le transport de V
         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &
     &    ,zu,pduadj,zua,lev_out)
         call thermcell_dq(ngrid,nlay,dqimpl,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

      if (prt_level.ge.1) print*,'14 OK convect8'
!------------------------------------------------------------------
!   Calculs de diagnostiques pour les sorties
!------------------------------------------------------------------
!calcul de fraca pour les sorties
      
      if (sorties) then
      if (prt_level.ge.1) 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
!IM   do k=1,nlay
      do k=1,nlay-1
         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
!IM
      ierr=0
      do ig=1,ngrid
        if (pcon(ig).le.pplay(ig,nlay)) then 
           zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
           ierr=1
        endif
      enddo
      if (ierr==1) then
           abort_message = 'thermcellV0_main: les thermiques vont trop haut '
           CALL abort_gcm (modname,abort_message,1)
      endif

      if (prt_level.ge.1) 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
      if (prt_level.ge.1) 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      
      if (prt_level.ge.1) print*,'14d OK convect8'
      if (prt_level.ge.10)write(lunout,*)                                &
    &     'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
      do l=1,nlay
         do ig=1,ngrid
            zf=fraca(ig,l)
            zf2=zf/(1.-zf)
!
            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
            if(zw2(ig,l).gt.1.e-10) then
             wth2(ig,l)=zf2*(zw2(ig,l))**2
            else
             wth2(ig,l)=0.
            endif
            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 des flux: q, thetal et thetav
      do l=1,nlay
         do ig=1,ngrid
      wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
      wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
      wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
         enddo
      enddo
!

!!! nrlmd le 10/04/2012

!------------Test sur le LCL des thermiques
    do ig=1,ngrid
      ok_lcl(ig)=.false.
      if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
    enddo

!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
    do l=1,nlay-1
      do ig=1,ngrid
        if (ok_lcl(ig)) then 
          if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
          klcl(ig)=l
          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
          endif
        endif
      enddo
    enddo

!------------Hauteur des thermiques
!!jyg le 27/04/2012
!!    do ig =1,ngrid
!!    rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
!! &               -rhobarz(ig,klcl(ig)))*interp(ig)
!!    zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
!!    zmax(ig)=pphi(ig,lmax(ig))/rg
!!      if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
!!    enddo
    do ig =1,ngrid
     zmax(ig)=pphi(ig,lmax(ig))/rg
     if (ok_lcl(ig)) then 
      rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
 &               -rhobarz(ig,klcl(ig)))*interp(ig)
      zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
      zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
     else
      rhobarz0(ig)=0.
      zlcl(ig)=zmax(ig)
     endif
    enddo
!!jyg fin

!------------Calcul des propriétés du thermique au LCL 
  IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN 

  !-----Initialisation de la TKE moyenne 
   do l=1,nlay
    do ig=1,ngrid
     pbl_tke_max(ig,l)=0.
    enddo
   enddo

!-----Calcul de la TKE moyenne 
   do nsrf=1,nbsrf
    do l=1,nlay
     do ig=1,ngrid
     pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
     enddo
    enddo
   enddo

!-----Initialisations des TKE dans et hors des thermiques 
   do l=1,nlay
    do ig=1,ngrid
    therm_tke_max(ig,l)=pbl_tke_max(ig,l)
    env_tke_max(ig,l)=pbl_tke_max(ig,l)
    enddo
   enddo

!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
  &           rg,pplev,therm_tke_max)
!   print *,' thermcell_tke_transport -> '   !!jyg

!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
   do l=1,nlay
    do ig=1,ngrid
     pbl_tke_max(ig,l)=fraca(ig,l)*therm_tke_max(ig,l)+(1.-fraca(ig,l))*env_tke_max(ig,l)         !  Recalcul de TKE moyenne aprés transport de TKE_TH
     env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l))       !  Recalcul de TKE dans  l'environnement aprés transport de TKE_TH
     w_ls(ig,l)=-1.*omega(ig,l)/(RG*rhobarz(ig,l))                                                !  Vitesse verticale de grande échelle
    enddo
   enddo
!    print *,' apres w_ls = '   !!jyg

  do ig=1,ngrid
   if (ok_lcl(ig)) then
     fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
 &             -fraca(ig,klcl(ig)))*interp(ig)
     w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
 &         -zw2(ig,klcl(ig)))*interp(ig)
     w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
 &             -w_ls(ig,klcl(ig)))*interp(ig)
     therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
 &                     +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
     env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
 &                   -env_tke_max(ig,klcl(ig)))*interp(ig)
     pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
 &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
     if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
     if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
     if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
   else 
     fraca0(ig)=0.
     w0(ig)=0.
!!jyg le 27/04/2012
!!     zlcl(ig)=0.
!!
   endif
  enddo

  ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
!  print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg

!------------Triggering------------------
  IF (iflag_trig_bl.ge.1) THEN 

!-----Initialisations
   depth(:)=0.
   n2(:)=0.
   s2(:)=0.
   s_max(:)=0.

!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
   do ig=1,ngrid
     zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
     depth(ig)=zmax_moy(ig)-zlcl(ig)
     hmin(ig)=hmincoef*zlcl(ig)
     if (depth(ig).ge.10.) then 
       s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
       n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
!!
!!jyg le 27/04/2012
!!       s_max(ig)=s2(ig)*log(n2(ig))
!!       if (n2(ig) .lt. 1) s_max(ig)=0.
       s_max(ig)=s2(ig)*log(max(n2(ig),1.))
!!fin jyg
     else
       s2(ig)=0.
       n2(ig)=0.
       s_max(ig)=0.
     endif
   enddo
!   print *,'avant Calcul de Wmax '    !!jyg

!-----Calcul de Wmax et ALE_BL_STAT associée
!!jyg le 30/04/2012
!!   do ig=1,ngrid
!!     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
!!     w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/su)-log(2.*3.14)-log(2.*log(s_max(ig)/su)-log(2.*3.14))))
!!     ale_bl_stat(ig)=0.5*w_max(ig)**2
!!     else
!!     w_max(ig)=0.
!!     ale_bl_stat(ig)=0.
!!     endif
!!   enddo
   susqr2pi=su*sqrt(2.*Rpi)
   Reuler=exp(1.)
   do ig=1,ngrid
     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
      w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
      ale_bl_stat(ig)=0.5*w_max(ig)**2
     else
      w_max(ig)=0.
      ale_bl_stat(ig)=0.
     endif
   enddo

  ENDIF ! iflag_trig_bl
!  print *,'ENDIF  iflag_trig_bl'    !!jyg

!------------Closure------------------

  IF (iflag_clos_bl.ge.1) THEN 

!-----Calcul de ALP_BL_STAT
  do ig=1,ngrid
  alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
  alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
 &                   (w0(ig)**2)
  alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
 &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
    if (iflag_clos_bl.ge.2) then 
    alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
 &                   (w0(ig)**2)
    else
    alp_bl_conv(ig)=0.
    endif
  alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
  enddo

!-----Sécurité ALP infinie
  do ig=1,ngrid
   if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
  enddo

  ENDIF ! (iflag_clos_bl.ge.1)

!!! fin nrlmd le 10/04/2012

      if (prt_level.ge.10) then
         ig=igout
         do l=1,nlay
            print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
            print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
         enddo
      endif

!      print*,'avant calcul ale et alp' 
!calcul de ALE et ALP pour la convection
      Alp_bl(:)=0.
      Ale_bl(:)=0.
!          print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
      do l=1,nlay
      do ig=1,ngrid
           Alp_bl(ig)=max(Alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
           Ale_bl(ig)=max(Ale_bl(ig),0.5*zw2(ig,l)**2)
!          print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
      enddo
      enddo

!test:calcul de la ponderation des couches pour KE
!initialisations

      fm_tot(:)=0.
      wght_th(:,:)=1.
      lalim_conv(:)=lalim(:)

      do k=1,klev
         do ig=1,ngrid
            if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
         enddo
      enddo

! assez bizarre car, si on est dans la couche d'alim et que alim_star et
! plus petit que 1.e-10, on prend wght_th=1.
      do k=1,klev
         do ig=1,ngrid
            if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
               wght_th(ig,k)=alim_star(ig,k)
            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

!------------------------------------------------------------------------
! Modif CR/FH 20110310 : Alp integree sur la verticale.
! Integrale verticale de ALP.
! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
! couches
!------------------------------------------------------------------------

      alp_int(:)=0.
      dp_int(:)=0.
      do l=2,nlay
        do ig=1,ngrid
           if(l.LE.lmax(ig)) THEN
           zdp=pplay(ig,l-1)-pplay(ig,l)
           alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)*zdp
           dp_int(ig)=dp_int(ig)+zdp
           endif
        enddo
      enddo

      if (iflag_coupl>=3 .and. iflag_coupl<=5) then
      do ig=1,ngrid
!valeur integree de alp_bl * 0.5:
        if (dp_int(ig)>0.) then
        Alp_bl(ig)=alp_int(ig)/dp_int(ig)
        endif
      enddo!
      endif


! Facteur multiplicatif sur Alp_bl
      Alp_bl(:)=alp_bl_k*Alp_bl(:)

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


!calcul du ratqscdiff
      if (prt_level.ge.1) print*,'14e OK convect8'
      var=0.
      vardiff=0.
      ratqsdiff(:,:)=0.

      do l=1,klev
         do ig=1,ngrid
            if (l<=lalim(ig)) then
            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
            endif
         enddo
      enddo

      if (prt_level.ge.1) print*,'14f OK convect8'

      do l=1,klev
         do ig=1,ngrid
            if (l<=lalim(ig)) then
               zf=fraca(ig,l)
               zf2=zf/(1.-zf)
               vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2
            endif
         enddo
      enddo

      if (prt_level.ge.1) print*,'14g OK convect8'
      do l=1,nlay
         do ig=1,ngrid
            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
         enddo
      enddo 
!--------------------------------------------------------------------    
!
!ecriture des fichiers sortie
!     print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'

#ifdef wrgrads_thermcell
      if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
#include "thermcell_out3d.h"
#endif

      endif

      if (prt_level.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)
      IMPLICIT NONE
#include "iniprint.h"

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

      if (prt_level.ge.1) THEN
       print*,'WARNING !!! TEST ',comment
      endif
      return

!  test sur la hauteur des thermiques ...
         do i=1,klon
!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
           if (prt_level.ge.10) then
               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
               do k=1,klev
                  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
           endif
         enddo


      return
      end

!!! nrlmd le 10/04/2012                          Transport de la TKE par le thermique moyen pour la fermeture en ALP 
!                                                         On transporte pbl_tke pour donner therm_tke
!                                          Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
      subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
     &           rg,pplev,therm_tke_max)
      implicit none

#include "iniprint.h"
!=======================================================================
!
!   Calcul du transport verticale dans la couche limite en presence
!   de "thermiques" explicitement representes
!   calcul du dq/dt une fois qu'on connait les ascendances
!
!=======================================================================

      integer ngrid,nlay,nsrf

      real ptimestep
      real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
      real entr0(ngrid,nlay),rg
      real therm_tke_max(ngrid,nlay)
      real detr0(ngrid,nlay)


      real masse(ngrid,nlay),fm(ngrid,nlay+1)
      real entr(ngrid,nlay)
      real q(ngrid,nlay)
      integer lev_out                           ! niveau pour les print

      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)

      real zzm

      integer ig,k
      integer isrf


      lev_out=0


      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'

!   calcul du detrainement
      do k=1,nlay
         detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
         masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
      enddo


! Decalage vertical des entrainements et detrainements.
      masse(:,1)=0.5*masse0(:,1)
      entr(:,1)=0.5*entr0(:,1)
      detr(:,1)=0.5*detr0(:,1)
      fm(:,1)=0.
      do k=1,nlay-1
         masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
         entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
         detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
      enddo
      fm(:,nlay+1)=0.

!!! nrlmd le 16/09/2010
!   calcul de la valeur dans les ascendances
!       do ig=1,ngrid
!          qa(ig,1)=q(ig,1)
!       enddo
!!! 

!do isrf=1,nsrf

!   q(:,:)=therm_tke(:,:,isrf)
   q(:,:)=therm_tke_max(:,:)
!!! nrlmd le 16/09/2010
      do ig=1,ngrid
         qa(ig,1)=q(ig,1)
      enddo
!!!

    if (1==1) then
      do k=2,nlay
         do ig=1,ngrid
            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     &         1.e-5*masse(ig,k)) then
         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
     &         /(fm(ig,k+1)+detr(ig,k))
            else
               qa(ig,k)=q(ig,k)
            endif
            if (qa(ig,k).lt.0.) then
!               print*,'qa<0!!!'
            endif
            if (q(ig,k).lt.0.) then
!               print*,'q<0!!!'
            endif
         enddo
      enddo

! Calcul du flux subsident

      do k=2,nlay
         do ig=1,ngrid
            wqd(ig,k)=fm(ig,k)*q(ig,k)
            if (wqd(ig,k).lt.0.) then
!               print*,'wqd<0!!!'
            endif
         enddo
      enddo
      do ig=1,ngrid
         wqd(ig,1)=0.
         wqd(ig,nlay+1)=0.
      enddo

! Calcul des tendances
      do k=1,nlay
         do ig=1,ngrid
            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
     &               -wqd(ig,k)+wqd(ig,k+1))  &
     &               *ptimestep/masse(ig,k)
         enddo
      enddo

 endif

   therm_tke_max(:,:)=q(:,:)

      return
!!! fin nrlmd le 10/04/2012
     end