! $Header$

    cin, eauliq, ctei, d_qt, d_thv, dlt_2, xhis, posint, omega, diagok)

  USE dimphy

  IMPLICIT NONE

  ! ***************************************************************

  ! *                                                             *

  ! * HBTM2L   D'apres Holstag&Boville et Troen&Mahrt             *

  ! *                 JAS 47              BLM                     *

  ! * Algorithmes These Anne Mathieu                              *

  ! * Critere d'Entrainement Peter Duynkerke (JAS 50)             *

  ! * written by  : Anne MATHIEU & Alain LAHELLEC, 22/11/99       *

  ! * features : implem. exces Mathieu                            *

  ! ***************************************************************

  ! * mods : decembre 99 passage th a niveau plus bas. voir fixer *

  ! ***************************************************************

  ! * fin therm a la HBTM passage a forme Mathieu 12/09/2001      *

  ! ***************************************************************

  ! AM Fev 2003 Adaptation a LMDZ version couplee                 *

  ! *

  ! Pour le moment on fait passer en argument les grdeurs de surface :

  ! flux, t,q2m, t,q10m, on va utiliser systematiquement les grdeurs a 2m ms

  ! on garde la possibilite de changer si besoin est (jusqu'a present la

  ! forme de HB avec le 1er niveau modele etait conservee)       *

  ! ***************************************************************

  ! * re-ecriture complete Alain Mars 2012 dans LMDZ5V5           *

  ! ***************************************************************

  include "YOMCST.h"

  REAL rlvcp, reps

  ! Arguments:

  INTEGER knon ! nombre de points a calculer

  ! AM

  REAL t2m(klon), t10m(klon) ! temperature a 2 et 10m

  REAL q2m(klon), q10m(klon) ! q a 2 et 10m

  REAL ustar(klon)

  REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)

  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)

  REAL flux_t(klon, klev), flux_q(klon, klev) ! Flux

  REAL u(klon, klev) ! vitesse U (m/s)

  REAL v(klon, klev) ! vitesse V (m/s)

  REAL t(klon, klev) ! temperature (K)

  REAL q(klon, klev) ! vapeur d'eau (kg/kg)

  INTEGER isommet

  REAL vk

  PARAMETER (vk=0.35) ! Von Karman => passer a .41 ! cf U.Olgstrom

  REAL ricr

  PARAMETER (ricr=0.4)

  REAL fak

  PARAMETER (fak=8.5) ! b calcul du Prandtl et de dTetas

  REAL fakn

  PARAMETER (fakn=7.2) ! a

  REAL onet

  PARAMETER (onet=1.0/3.0)

  REAL betam

  PARAMETER (betam=15.0) ! pour Phim / h dans la S.L stable

  REAL betah

  PARAMETER (betah=15.0)

  REAL betas

  PARAMETER (betas=5.0) ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1

  REAL sffrac

  PARAMETER (sffrac=0.1) ! S.L. = z/h < .1

  REAL binm

  PARAMETER (binm=betam*sffrac)

  REAL binh

  PARAMETER (binh=betah*sffrac)

  REAL q_star, t_star

  REAL b1, b2, b212, b2sr ! Lambert correlations T' q' avec T* q*

  PARAMETER (b1=70., b2=20.) ! b1 entre 70 et 100

  ! AM

  REAL zref, dt0

  PARAMETER (zref=2.) ! Niveau de ref a 2m

  PARAMETER (dt0=0.1) ! convergence do while

  INTEGER i, k, j

  REAL pblh(klon) ! PBL H               (m)

  REAL therm(klon) ! exces du thermique  (K)

  REAL plcl(klon) ! Lifted Cnd Level (Pa)

  REAL cape(klon) ! Cape

  REAL cin(klon) ! Inhibition

  REAL eauliq(klon) ! Eau Liqu integree

  REAL ctei(klon) ! Cld Top Entr. Instab.

  REAL d_qt(klon) ! Saut de qT a l'inversion

  REAL d_thv(klon) !         Theta_e

  REAL dlt_2(klon) ! Ordonnee a gauche de courbe de melange

  REAL xhis(klon) ! fraction de melange pour flottab nulle

  REAL posint(klon) ! partie positive de l'int. de Peter

  REAL omega(klon) ! point ultime de l'ascention du thermique

  REAL diagok(klon) ! pour traiter les sous-mailles sans info

  ! Algorithme thermique

  ! calcul de CTEI etc

  REAL the1, the2, aa, bb, zthvd, zthvu, qsat, chi, rh, zxt, zdu2

  REAL rnum, denom, th1, th2, tv1, tv2, thv1, thv2, ql1, ql2, dt

  REAL dqsat_dt, qsat2, qt1, q1, q2, t1, t2, tl1, te2, xnull, delt_the

  REAL delt_qt, quadsat, spblh, reduc

  ! diag      REAL dTv21(klon,klev)

  REAL zcor, zdelta, zcvm5

  REAL fac, pblmin

  REAL missing_val

  include "YOETHF.h"

  include "FCTTRE.h"

  ! c      missing_val=nf90_fill_real (avec include netcdf)

  missing_val = 0.

  ! initialisations (Anne)

  isommet = klev

  b212 = sqrt(b1*b2)

  b2sr = sqrt(b2)

  ! Initialisation thermo

  reps = rd/rv

  ! raz

  q_star = 0.

  t_star = 0.

  ! diag      dTv21(:,:)= missing_val

  ! Calculer les hauteurs de chaque couche

  END DO

  ! s(k) = [pplay(k)/ps]^kappa

  ! + + + + + + + + + pplay  <-> s(k)   t  dp=pplay(k-1)-pplay(k)

  ! -----------------  paprs <-> sig(k)

  ! + + + + + + + + + pplay  <-> s(k-1)

  ! + + + + + + + + + pplay  <-> s(1)   t  dp=paprs-pplay   z(1)

  ! -----------------  paprs <-> sig(1)

    END DO

  END DO

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  ! +++  Determination des grandeurs de surface  +++++++++++++++++++++

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    ! AM Niveau de ref choisi a 2m

    ! ***************************************************

    ! attention, il doit s'agir de <w'theta'>

    ! ;Calcul de tcls virtuel et de w'theta'virtuel

    ! ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    ! tcls=tcls*(1+.608*qcls)

    ! ;Pour avoir w'theta',

    ! ; il faut diviser par ro.Cp

    ! Cp=Cpd*(1+0.84*qcls)

    ! fcs=fcs/(ro_surf*Cp)

    ! ;On transforme w'theta' en w'thetav'

    ! Lv=(2.501-0.00237*(tcls-273.15))*1.E6

    ! xle=xle/(ro_surf*Lv)

    ! fcsv=fcs+.608*xle*tcls

    ! ***************************************************

    ! dif khfs est deja w't'_v / heatv(i) = khfs(i) + RETV*zxt*kqfs(i)

    ! AM calcul de Ro = paprs(i,1)/Rd zxt

    ! AM convention >0 vers le bas ds lmdz

    ! AM   verifier que khfs et kqfs sont bien de la forme w'l'

    ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv

    ! AM ustar est en entree (calcul dans stdlevvar avec t2m q2m)

    ! Theta et qT du thermique sans exces

    ! Al1 Th_th restera la Theta du thermique sans exces jusqu'au 3eme calcul

  END DO

    ! Attention Plcl est pression ou altitude ?

    ! plcl(i) = 6000. ! m

      ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>

    ELSE

      ! set pblh to the friction high (cf + bas)

    END IF

  END DO

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  ! PBL height calculation:

  ! Search for level of pbl. Scan upward until the Richardson number between

  ! the first level and the current level exceeds the "critical" value.

  ! (bonne idee Nu de separer le Ric et l'exces de temp du thermique)

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  fac = 100.0

        ! Theta_v environnement

        ! Le Ri bulk par Theta_v

          ! test04 (la pblh est encore ici sous-estime'e)

          ! pblT(i) = t(i,k-1) + (t(i,k)-t(i,k-1)) *

          ! .              (pblh(i)-z(i,k-1))/(z(i,k)-z(i,k-1))

        END IF

      END IF

    END DO

  END DO

  ! Set pbl height to maximum value where computation exceeds number of

  ! layers allowed

  END DO

  ! Improve estimate of pbl height for the unstable points.

  ! Find unstable points (sensible heat flux is upward):

    ELSE

    END IF

  END DO

  ! For the unstable case, compute velocity scale and the

  ! convective temperature excess:

      ! ***************************************************

      ! Wm ? et W* ? c'est la formule pour z/h < .1

      ! ;Calcul de w* ;;

      ! ;;;;;;;;;;;;;;;;

      ! w_star=((g/tcls)*fcsv*z(ind))^(1/3.) [ou prendre la premiere approx de h)

      ! ;; CALCUL DE wm ;;

      ! ;;;;;;;;;;;;;;;;;;

      ! ; Ici on considerera que l'on est dans la couche de surf jusqu'a 100m

      ! ; On prend svt couche de surface=0.1*h mais on ne connait pas h

      ! ;;;;;;;;;;;Dans la couche de surface

      ! if (z(ind) le 20) then begin

      ! Phim=(1.-15.*(z(ind)/L))^(-1/3.)

      ! wm=u_star/Phim

      ! ;;;;;;;;;;;En dehors de la couche de surface

      ! endif else if (z(ind) gt 20) then begin

      ! wm=(u_star^3+c1*w_star^3)^(1/3.)

      ! endif

      ! ***************************************************

      ! ======================================================================

      ! valeurs de Dominique Lambert de la campagne SEMAPHORE :

      ! <T'^2> = 100.T*^2; <q'^2> = 20.q*^2 a 10m

      ! <Tv'^2> = (1+1.2q).100.T* + 1.2Tv.sqrt(20*100).T*.q* + (.608*Tv)^2*20.q*^2;

      ! et dTetavS = sqrt(<Tv'^2>) ainsi calculee.

      ! avec : T*=<w'T'>_s/w* et q*=<w'q'>/w*

      ! !!! on peut donc utiliser w* pour les fluctuations <-> Lambert

      ! (leur corellation pourrait dependre de beta par ex)

      ! if fcsv(i,j) gt 0 then begin

      ! dTetavs=b1*(1.+2.*.608*q_10(i,j))*(fcs(i,j)/wm(i,j))^2+$

      ! (.608*Thetav_10(i,j))^2*b2*(xle(i,j)/wm(i,j))^2+$

      ! 2.*.608*thetav_10(i,j)*sqrt(b1*b2)*(xle(i,j)/wm(i,j))*(fcs(i,j)/wm(i,j))

      ! dqs=b2*(xle(i,j)/wm(i,j))^2

      ! theta_s(i,j)=thetav_10(i,j)+sqrt(dTetavs)

      ! q_s(i,j)=q_10(i,j)+sqrt(dqs)

      ! endif else begin

      ! Theta_s(i,j)=thetav_10(i,j)

      ! q_s(i,j)=q_10(i,j)

      ! endelse

      ! leur reference est le niveau a 10m, mais on prend 2m ici.

      ! ======================================================================

      ! Premier calcul de l'exces tu thermique

      ! ======================================================================

      ! HBTM        therm(i) = heatv(i)*fak/wm(i)

      ! forme Mathieu :

      ! Al1 Houston, we have a problem : il arrive en effet que heatv soit

      ! positif (=thermique instable) mais pas t_star : avec evaporation

      ! importante, il se peut qu'on refroidisse la 2m Que faire alors ?

      ! Garder le seul terme en q_star^2 ? ou rendre negatif le t_star^2 ?

      therm(i) = sqrt(b1*(1.+2.*retv*qt_th(i))*t_star**2+(retv*th_th(i))**2*b2*q_star*q_star+2.*retv*th_th(i)*b212* &

      ! Theta et qT du thermique (forme H&B) avec exces

      ! (attention, on ajoute therm(i) qui est virtuelle ...)

      ! pourquoi pas sqrt(b1)*t_star ?

      ! dqs = b2sr*kqfs(i)/wm(i)

    END IF

  END DO

  ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  ! ++ Improve pblh estimate for unstable conditions using the +++++++

  ! ++          convective temperature excess :                +++++++

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

        ! test     zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2

        ! Theta_v environnement

        ! et therm Theta_v (avec hypothese de constance de H&B,

        ! qui assimile qT a vapeur)

        ! Le Ri par Theta_v

        ! AM Niveau de ref 2m

        ! Niveau critique atteint

          ! test04

          ! pblT(i) = t(i,k-1) + (t(i,k)-t(i,k-1)) *

          ! .              (pblh(i)-z(i,k-1))/(z(i,k)-z(i,k-1))

        END IF

      END IF

    END DO

  END DO

  ! Set pbl height to maximum value where computation exceeds number of

  ! layers allowed (H&B)

  END DO

  ! PBL height must be greater than some minimum mechanical mixing depth

  ! Several investigators have proposed minimum mechanical mixing depth

  ! relationships as a function of the local friction velocity, u*.  We

  ! make use of a linear relationship of the form h = c u* where c=700.

  ! The scaling arguments that give rise to this relationship most often

  ! represent the coefficient c as some constant over the local coriolis

  ! parameter.  Here we make use of the experimental results of Koracin

  ! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f

  ! where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid

  ! latitude value for f so that c = 0.07/f = 700.  (H&B)

  ! Al1 calcul de pblT dans ce cas

  END DO

      ! et par exemple :

      ! pblT(i) = t(i,2) + (t(i,3)-t(i,2)) *

      ! .              (pblh(i)-z(i,2))/(z(i,3)-z(i,2))

    END IF

  END DO

  ! ********************************************************************

  ! pblh is now available; do preparation for final calculations :

  ! ********************************************************************

    ! omegafl utilise pour prolongement CAPE

    ! Do additional preparation for unstable cases only, set temperature

    ! and moisture perturbations depending on stability.

    ! Rq: les formules sont prises dans leur forme Couche de Surface

      ! Al pblh a change', on recalcule :

      zxt = (th_th(i)-zref*0.5*rg/rcpd/(1.+rvtmp2*qt_th(i)))*(1.+retv*qt_th(i))

    END IF

  END DO

  ! =======================================================

  ! last upward integration

  ! For all unstable layers, compute integral info and CTEI

  ! =======================================================

  ! 1/Recompute surface characteristics with the improved pblh

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

      ! from missing_value to zero

      therm(i) = sqrt(b1*(1.+2.*retv*qt_th(i))*t_star**2+(retv*th_th(i))**2*b2*q_star*q_star+2.*retv*th_th(i)*b212* &

      ! Al1diag

      ! trmb1(i) = b1*(1.+2.*RETV*qT_th(i))*t_star**2

      ! trmb2(i) = (RETV*Th_th(i))**2*b2*q_star*q_star

      ! trmb3(i) = 2.*RETV*Th_th(i)*b212*q_star*t_star

      ! Th_th will now be the thermal-theta (including exces)

      ! c         Th_th(i) = Th_th(i)+sqrt(b1)*max(0.,khfs(i)/wm(i))

      ! al1diag

      ! trmb2(i) = wm(i)

      ! trmb3(i) = phiminv(i)

      ! and computes Theta_e for thermal

    END IF ! unstbl

    ! Al1 compute a first guess of Plcl with the Bolton/Emanuel formula

    ! thermodyn functions

    ! relative humidity of thermal at 2m

    ! al1diag

    ! ctei(i) = Plcl(i)

    ! cape(i) = T2

    ! trmb1(i)= Chi

    ! select unstable columns (=thermals)

    ! diag

    ! dTv21(i,1) = T2*(1+RETV*qT_th(i))-t(i,1)*(1+RETV*q(i,1))

  END DO

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

  ! 2/ upward integration for thermals

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ k loop

        ! CC         if (pplay(i,k) .le. plcl(i)) then

        ! Environnement : calcul de Tv1 a partir de t(:,:)== T liquide

        ! ==============

        t1 = tl1

        ! thermodyn function (Tl2Tql)

          ! correction lineaire pour conserver Tl env

          ! << Tl = T1 + DT - RLvCp*(ql1 - dqsat/dT*DT >>

        END DO

        ! Thermique    : on atteint le seuil B/E de condensation

        ! ==============

          ! first guess from The_th(i) = Th_th(i) + RLvCp* [qv=qT_th(i)]

          ! a PBLH non sature

            ! lmdz : qT1 et Thv1

            ! on calcule pour le cas sans nuage un ctei en Delta Thv

            ! diag

            ! dTv21(i,k) = Tv2-Tv1

          END IF

        END IF

          ! thermodyn functions (Te2Tqsat)

          dt = 1.

            ! correction lineaire pour conserver Te_th

            ! << Te = T2 + DT + RLvCp*(qsatbef + dq/dT*DT >>

          END DO

          ! jusqu'a PBLH y compris

            ! mais a PBLH, interpolation et complements

              ! CAPE et EauLiq a pblH

              ! CTEI

              ! T1 est en realite la Tl env (on a donc strict The1)

              ! Calcul de la Cloud Top Entrainement Instability

              ! cf Mathieu Lahellec QJRMS (2005) Comments to DYCOMS-II

              ! saut a l'inversion :

              ! init ctei(i)

              ctei(i) = dlt_2(i)

                ! integrale de Peter :

                ! approx de Xhi_s et de l'integrale Xint=ctei(i)

                ! trmb1(i) = xhis

                ! trmb3(i) = dlt_2

                ELSE

                END IF

                ELSE

                END IF

              ELSE

              END IF

            END IF ! end a pblh

          END IF

        END IF ! Zsat

        ! KAPE : thermique / environnement

        ! diag

        ! dTv21(i,k) = Tv2-Tv1

        ! Kape courante

        ! Cin

        END IF

        END IF

          ! trmb3(i) = paprs(i,k)

          ! diag

          ! print*,'Tv2-Tv1 (k): ',i,(dTv21(i,j),j=1,k)

        END IF

        ! CC         EndIf !plcl

      END IF ! check(i)

    END DO

  END DO ! end of level loop

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ end k loop

END SUBROUTINE hbtm2l