!
!
!
SUBROUTINE thermcell_alp(ngrid,nlay,ptimestep,pplay,pplev,                    &
                         fm0,entr0,lmax,ale_bl,alp_bl,lalim_conv,             &
                         wght_th,zw2,fraca,                                   &
                         pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,  &
                         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)
      
      
!==============================================================================
!   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"
!
!==============================================================================
      
      USE thermcell_mod
      
      IMPLICIT NONE
      
      
!==============================================================================
! Declarations
!==============================================================================
      
!      Inputs:
!      -------
      
      INTEGER,INTENT(in) :: ngrid
      INTEGER,INTENT(in) :: nlay
      INTEGER lmax(ngrid)
      INTEGER lalim(ngrid)
      
      REAL,INTENT(in) :: ptimestep
      REAL,INTENT(in) :: pplay(ngrid,nlay)
      REAL,INTENT(in) :: pplev(ngrid,nlay+1)
! On garde le zmax du pas de temps precedent
      REAL zmax(ngrid)
      REAL zw2(ngrid,nlay+1)
      REAL fraca(ngrid,nlay+1)
      real fm0(ngrid,nlay+1)
      real entr0(ngrid,nlay)
      real pbl_tke(ngrid,nlay+1,nbsrf)
      real pctsrf(ngrid,nbsrf)
      real omega(ngrid,nlay)
      real airephy(ngrid)
      real alim_star(ngrid,nlay)
      
!      Outputs:
!     ---------
      
      REAL susqr2pi
      REAL reuler
      real wth3(ngrid,nlay)
! FH probleme de dimensionnement avec l'allocation dynamique
!     common/comtherm/thetath2,wth2
      real rhobarz(ngrid,nlay)
      real wmax_sec(ngrid)
      real fm(ngrid,nlay+1)
      
! niveau de condensation
      real pcon(ngrid)
      
      real zlcl(ngrid)
      real fraca0(ngrid)
      real w0(ngrid)
      real w_conv(ngrid)
      real therm_tke_max0(ngrid)
      real env_tke_max0(ngrid)
      real n2(ngrid)
      real s2(ngrid)
      real ale_bl_stat(ngrid)
      real therm_tke_max(ngrid,nlay)
      real env_tke_max(ngrid,nlay)
      real alp_bl_det(ngrid)
      real alp_bl_fluct_m(ngrid)
      real alp_bl_fluct_tke(ngrid)
      real alp_bl_conv(ngrid)
      real alp_bl_stat(ngrid)
      
!      local:
!      ------
      
      INTEGER ig,k,l
      integer nsrf
      integer klcl(ngrid)                     ! Niveau du LCL
      
      real rhobarz0(ngrid)                    ! Densité au LCL
      logical ok_lcl(ngrid)                   ! Existence du LCL des thermiques
      real interp(ngrid)                      ! Coef d'interpolation pour le LCL
!--Triggering
      real,parameter :: Su = 4.e4             ! Surface unité: celle d'un updraft élémentaire
      real,parameter :: hcoef = 1.            ! Coefficient directeur pour le calcul de s2
      real,parameter :: hmincoef = 0.3        ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 
      real,parameter :: eps1 = 0.3            ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
      real,parameter :: zmax_moy_coef=0.33
      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 depth(ngrid)                       ! Epaisseur moyenne du cumulus
      real w_max(ngrid)                       ! Vitesse max statistique 
      real s_max(ngrid)
!--Closure
      real pbl_tke_max(ngrid,nlay)            ! Profil de TKE moyenne 
      real pbl_tke_max0(ngrid)                ! TKE moyenne au LCL
      real w_ls(ngrid,nlay)                   ! Vitesse verticale grande échelle (m/s)
      real,parameter :: coef_m=1.             ! On considère un rendement pour alp_bl_fluct_m
      real,parameter :: coef_tke = 1.         ! On considère un rendement pour alp_bl_fluct_tke
      
! nouvelles variables pour la convection
      real ale_bl(ngrid)
      real alp_bl(ngrid)
      real alp_int(ngrid),dp_int(ngrid),zdp
      real fm_tot(ngrid)
      real wght_th(ngrid,nlay)
      integer lalim_conv(ngrid)
!v1d     logical therm
!v1d     save therm
      
!==============================================================================
! Initialization
!==============================================================================
      
      DO ig = 1,ngrid
         zlcl(ig) = 0.
         fraca0(ig) = 0.
         w0(ig) = 0.
         w_conv(ig) = 0.
         therm_tke_max0(ig) = 0.
         env_tke_max0(ig) = 0.
         n2(ig) = 0.
         s2(ig) = 0.
         ale_bl_stat(ig) = 0.
         alp_bl_det(ig) = 0.
         alp_bl_fluct_m(ig) = 0.
         alp_bl_fluct_tke(ig) = 0.
         alp_bl_conv(ig) = 0.
         alp_bl_stat(ig) = 0.
      ENDDO
      
      DO l = 1,nlay
         DO ig = 1,ngrid
            therm_tke_max(ig,l) = 0.
            env_tke_max(ig,l) = 0.
         ENDDO
      ENDDO
      
!------------------------------------------------------------
!------------Test sur le LCL des thermiques
      do ig=1,ngrid
         ok_lcl(ig)=.false.
         if ( (pcon(ig) .gt. pplay(ig,nlay-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 
!ATTENTION,zw2 calcule en pplev
!               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
               if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
                  klcl(ig)=l
                  interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(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)
!         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
!CR:REHABILITATION ZMAX CONTINU
         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(:)=100. ! some low value, arbitrary
         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
               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.2) 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.2)
      
!!! fin nrlmd le 10/04/2012
      
!      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
      
! ale sec (max de wmax/2 sous la zone d'inhibition) dans le cas iflag_trig_bl=3
      
      IF (iflag_trig_bl==3) then
         ale_bl(:)=0.5*wmax_sec(:)**2
      ENDIF
      
!test:calcul de la ponderation des couches pour KE
!initialisations
      
      fm_tot(:)=0.
      wght_th(:,:)=1.
      lalim_conv(:)=lalim(:)
      
      do k=1,nlay
         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,nlay
         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(:)
      
      
return
end