Sources/phylmd/hbtm.f

module HBTM_m

  IMPLICIT none

contains

  SUBROUTINE HBTM(knon, paprs, pplay, t2m, q2m, ustar, flux_t, &
       flux_q, u, v, t, q, pblh, cape, EauLiq, ctei, pblT, therm, trmb1, &
       trmb2, trmb3, plcl)

    ! D'apr\'es Holstag et Boville et Troen et Mahrt
    ! JAS 47 BLM
    ! Algorithme th\'ese Anne Mathieu
    ! Crit\'ere d'entra\^inement Peter Duynkerke (JAS 50)
    ! written by: Anne MATHIEU and Alain LAHELLEC, 22nd November 1999
    ! features : implem. exces Mathieu

    ! modifications : decembre 99 passage th a niveau plus bas. voir fixer
    ! la prise du th a z/Lambda = -.2 (max Ray)
    ! Autre algo : entrainement ~ Theta+v =cste mais comment=>The?
    ! on peut fixer q a .7 qsat (cf. non adiabatique) => T2 et The2
    ! voir aussi //KE pblh = niveau The_e ou l = env.

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

    ! Adaptation a LMDZ version couplee Pour le moment on fait passer
    ! en argument les grandeurs de surface : flux, t, q2m, t, on va
    ! utiliser systematiquement les grandeurs a 2m mais on garde la
    ! possibilit\'e de changer si besoin est (jusqu'\`a pr\'esent la
    ! forme de HB avec le 1er niveau modele etait conservee)

    USE dimphy, ONLY: klev, klon
    USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rtt
    USE yoethf_m, ONLY: r2es, rvtmp2
    USE fcttre, ONLY: foeew

    ! Arguments:

    ! nombre de points a calculer
    INTEGER, intent(in):: knon

    REAL, intent(in):: t2m(klon) ! temperature a 2 m
    ! q a 2 et 10m
    REAL q2m(klon)
    REAL ustar(klon)
    ! pression a inter-couche (Pa)
    REAL paprs(klon, klev+1)
    ! pression au milieu de couche (Pa)
    REAL pplay(klon, klev)
    ! Flux
    REAL flux_t(klon, klev), flux_q(klon, klev)
    ! vitesse U (m/s)
    REAL u(klon, klev)
    ! vitesse V (m/s)
    REAL v(klon, klev)
    ! temperature (K)
    REAL t(klon, klev)
    ! vapeur d'eau (kg/kg)
    REAL q(klon, klev)

    INTEGER isommet
    ! limite max sommet pbl
    PARAMETER (isommet=klev)
    REAL vk
    ! Von Karman => passer a .41 ! cf U.Olgstrom
    PARAMETER (vk=0.35)
    REAL ricr
    PARAMETER (ricr=0.4)
    REAL fak
    ! b calcul du Prandtl et de dTetas
    PARAMETER (fak=8.5)
    REAL fakn
    ! a
    PARAMETER (fakn=7.2)
    REAL onet
    PARAMETER (onet=1.0/3.0)
    REAL t_coup
    PARAMETER(t_coup=273.15)
    REAL zkmin
    PARAMETER (zkmin=0.01)
    REAL betam
    ! pour Phim / h dans la S.L stable
    PARAMETER (betam=15.0)
    REAL betah
    PARAMETER (betah=15.0)
    REAL betas
    ! Phit dans la S.L. stable (mais 2 formes /
    PARAMETER (betas=5.0)
    ! z/OBL<>1
    REAL sffrac
    ! S.L. = z/h < .1
    PARAMETER (sffrac=0.1)
    REAL binm
    PARAMETER (binm=betam*sffrac)
    REAL binh
    PARAMETER (binh=betah*sffrac)
    REAL ccon
    PARAMETER (ccon=fak*sffrac*vk)

    REAL q_star, t_star
    ! Lambert correlations T' q' avec T* q*
    REAL b1, b2, b212, b2sr
    PARAMETER (b1=70., b2=20.)

    REAL z(klon, klev)

    REAL zref
    ! Niveau de ref a 2m peut eventuellement
    PARAMETER (zref=2.)
    ! etre choisi a 10m

    INTEGER i, k
    REAL zxt
    ! surface kinematic heat flux [mK/s]
    REAL khfs(klon)
    ! sfc kinematic constituent flux [m/s]
    REAL kqfs(klon)
    ! surface virtual heat flux
    REAL heatv(klon)
    ! bulk Richardon no. mais en Theta_v
    REAL rhino(klon, klev)
    ! pts w/unstbl pbl (positive virtual ht flx)
    LOGICAL unstbl(klon)
    LOGICAL check(klon) ! Richardson number > critical
    ! flag de prolongerment cape pour pt Omega
    LOGICAL omegafl(klon)
    REAL pblh(klon)
    REAL pblT(klon)
    REAL plcl(klon)

    ! Monin-Obukhov lengh
    REAL obklen(klon)

    REAL zdu2
    ! thermal virtual temperature excess
    REAL therm(klon)
    REAL trmb1(klon), trmb2(klon), trmb3(klon)
    ! Algorithme thermique
    REAL s(klon, klev) ! [P/Po]^Kappa milieux couches
    ! total water of thermal
    REAL qT_th(klon)
    ! T thermique niveau precedent
    REAL qsatbef(klon)
    ! le thermique est sature
    LOGICAL Zsat(klon)
    ! Cape du thermique
    REAL Cape(klon)
    ! Eau liqu integr du thermique
    REAL EauLiq(klon)
    ! Critere d'instab d'entrainmt des nuages de
    REAL ctei(klon)
    REAL zthvd, zthvu, qqsat
    REAL t2

    ! inverse phi function for momentum
    REAL phiminv(klon)
    ! turbulent velocity scale for momentum
    REAL wm(klon)
    ! current level height + one level up
    REAL zp(klon)
    REAL zcor

    REAL pblmin

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

    ! initialisations
    q_star = 0
    t_star = 0

    b212=sqrt(b1*b2)
    b2sr=sqrt(b2)

    ! Calculer les hauteurs de chaque couche
    ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
    ! pourquoi ne pas utiliser Phi/RG ?
    DO i = 1, knon
       z(i, 1) = RD * t(i, 1) / (0.5*(paprs(i, 1)+pplay(i, 1))) &
            * (paprs(i, 1)-pplay(i, 1)) / RG
       s(i, 1) = (pplay(i, 1)/paprs(i, 1))**RKappa
    ENDDO
    ! 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)

    DO k = 2, klev
       DO i = 1, knon
          z(i, k) = z(i, k-1) &
               + RD * 0.5*(t(i, k-1)+t(i, k)) / paprs(i, k) &
               * (pplay(i, k-1)-pplay(i, k)) / RG
          s(i, k) = (pplay(i, k) / paprs(i, 1))**RKappa
       ENDDO
    ENDDO

    ! Determination des grandeurs de surface
    DO i = 1, knon
       ! Niveau de ref choisi a 2m
       zxt = t2m(i)

       ! convention >0 vers le bas ds lmdz
       khfs(i) = - flux_t(i, 1)*zxt*Rd / (RCPD*paprs(i, 1))
       kqfs(i) = - flux_q(i, 1)*zxt*Rd / (paprs(i, 1))
       ! verifier que khfs et kqfs sont bien de la forme w'l'
       heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
       ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
       ! Theta et qT du thermique sans exces (interpolin vers surf)
       ! chgt de niveau du thermique (jeudi 30/12/1999)
       ! (interpolation lineaire avant integration phi_h)
       qT_th(i) = q2m(i)
    ENDDO

    DO i = 1, knon
       ! Global Richardson
       rhino(i, 1) = 0.0
       check(i) = .TRUE.
       ! on initialise pblh a l'altitude du 1er niv
       pblh(i) = z(i, 1)
       plcl(i) = 6000.
       ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
       obklen(i) = -t(i, 1)*ustar(i)**3/(RG*vk*heatv(i))
       trmb1(i) = 0.
       trmb2(i) = 0.
       trmb3(i) = 0.
    ENDDO

    ! 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)
    DO k = 2, isommet
       DO i = 1, knon
          IF (check(i)) THEN
             ! pourquoi / niveau 1 (au lieu du sol) et le terme en u*^2 ?
             zdu2 = u(i, k)**2+v(i, k)**2
             zdu2 = max(zdu2, 1.0e-20)
             ! Theta_v environnement
             zthvd=t(i, k)/s(i, k)*(1.+RETV*q(i, k))

             ! therm Theta_v sans exces (avec hypothese fausse de H&B, sinon,
             ! passer par Theta_e et virpot)
             zthvu = T2m(i)*(1.+RETV*qT_th(i))
             ! Le Ri par Theta_v
             ! On a nveau de ref a 2m ???
             rhino(i, k) = (z(i, k)-zref)*RG*(zthvd-zthvu) &
                  /(zdu2*0.5*(zthvd+zthvu))

             IF (rhino(i, k).GE.ricr) THEN
                pblh(i) = z(i, k-1) + (z(i, k-1)-z(i, k)) * &
                     (ricr-rhino(i, k-1))/(rhino(i, k-1)-rhino(i, k))
                ! test04
                pblh(i) = pblh(i) + 100.
                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))
                check(i) = .FALSE.
             ENDIF
          ENDIF
       ENDDO
    ENDDO

    ! Set pbl height to maximum value where computation exceeds number of
    ! layers allowed
    DO i = 1, knon
       if (check(i)) pblh(i) = z(i, isommet)
    ENDDO

    ! Improve estimate of pbl height for the unstable points.
    ! Find unstable points (sensible heat flux is upward):
    DO i = 1, knon
       IF (heatv(i) > 0.) THEN
          unstbl(i) = .TRUE.
          check(i) = .TRUE.
       ELSE
          unstbl(i) = .FALSE.
          check(i) = .FALSE.
       ENDIF
    ENDDO

    ! For the unstable case, compute velocity scale and the
    ! convective temperature excess:
    DO i = 1, knon
       IF (check(i)) THEN
          phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet

          ! CALCUL DE wm
          ! Ici on considerera que l'on est dans la couche de surf jusqu'a 100
          ! On prend svt couche de surface=0.1*h mais on ne connait pas h
          ! Dans la couche de surface
          wm(i)= ustar(i)*phiminv(i)

          ! forme Mathieu :
          q_star = kqfs(i)/wm(i)
          t_star = khfs(i)/wm(i)

          therm(i) = sqrt( b1*(1.+2.*RETV*qT_th(i))*t_star**2 &
               + (RETV*T2m(i))**2*b2*q_star**2 &
               + max(0., 2.*RETV*T2m(i)*b212*q_star*t_star))

          ! Theta et qT du thermique (forme H&B) avec exces
          ! (attention, on ajoute therm(i) qui est virtuelle ...)
          ! pourquoi pas sqrt(b1)*t_star ?
          qT_th(i) = qT_th(i) + b2sr*q_star
          ! new on differre le calcul de Theta_e
          rhino(i, 1) = 0.
       ENDIF
    ENDDO

    ! Improve pblh estimate for unstable conditions using the
    ! convective temperature excess :
    DO k = 2, isommet
       DO i = 1, knon
          IF (check(i)) THEN
             zdu2 = u(i, k)**2 + v(i, k)**2
             zdu2 = max(zdu2, 1e-20)
             ! Theta_v environnement
             zthvd=t(i, k)/s(i, k)*(1.+RETV*q(i, k))

             ! et therm Theta_v (avec hypothese de constance de H&B,
             zthvu = T2m(i)*(1.+RETV*qT_th(i)) + therm(i)

             ! Le Ri par Theta_v
             ! Niveau de ref 2m
             rhino(i, k) = (z(i, k)-zref)*RG*(zthvd-zthvu) &
                  /(zdu2*0.5*(zthvd+zthvu))

             IF (rhino(i, k).GE.ricr) THEN
                pblh(i) = z(i, k-1) + (z(i, k-1)-z(i, k)) * &
                     (ricr-rhino(i, k-1))/(rhino(i, k-1)-rhino(i, k))
                ! test04
                pblh(i) = pblh(i) + 100.
                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))
                check(i) = .FALSE.
             ENDIF
          ENDIF
       ENDDO
    ENDDO

    ! Set pbl height to maximum value where computation exceeds number of
    ! layers allowed
    DO i = 1, knon
       if (check(i)) pblh(i) = z(i, isommet)
    ENDDO

    ! 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.
    DO i = 1, knon
       pblmin = 700. * ustar(i)
       pblh(i) = MAX(pblh(i), pblmin)
       ! par exemple :
       pblT(i) = t(i, 2) + (t(i, 3)-t(i, 2)) * &
            (pblh(i)-z(i, 2))/(z(i, 3)-z(i, 2))
    ENDDO

    ! pblh is now available; do preparation for diffusivity calculation:
    DO i = 1, knon
       check(i) = .TRUE.
       Zsat(i) = .FALSE.
       ! omegafl utilise pour prolongement CAPE
       omegafl(i) = .FALSE.
       Cape(i) = 0.
       EauLiq(i) = 0.
       CTEI(i) = 0.

       ! Do additional preparation for unstable cases only, set temperature
       ! and moisture perturbations depending on stability.
       ! Remarque : les formule sont prises dans leur forme CS
       IF (unstbl(i)) THEN
          ! Niveau de ref du thermique
          zxt=(T2m(i)-zref*0.5*RG/RCPD/(1.+RVTMP2*qT_th(i))) &
               *(1.+RETV*qT_th(i))
          phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
          wm(i) = ustar(i)*phiminv(i)
       ENDIF
    ENDDO

    ! Main level loop to compute the diffusivities and
    ! counter-gradient terms:
    loop_level: DO k = 2, isommet
       ! Find levels within boundary layer:
       DO i = 1, knon
          zp(i) = z(i, k)
          IF (zkmin == 0. .AND. zp(i) > pblh(i)) zp(i) = pblh(i)
       ENDDO

       ! For all layers, compute integral info and CTEI
       DO i = 1, knon
          if (check(i) .or. omegafl(i)) then
             if (.not. Zsat(i)) then
                T2 = T2m(i) * s(i, k)
                ! thermodyn functions
                qqsat= r2es * FOEEW(T2, RTT >= T2) / pplay(i, k)
                qqsat=MIN(0.5, qqsat)
                zcor=1./(1.-retv*qqsat)
                qqsat=qqsat*zcor

                if (qqsat < qT_th(i)) then
                   ! on calcule lcl
                   if (k == 2) then
                      plcl(i) = z(i, k)
                   else
                      plcl(i) = z(i, k-1) + (z(i, k-1)-z(i, k)) &
                           * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)
                   endif
                   Zsat(i) = .true.
                endif
             endif
             qsatbef(i) = qqsat
             ! cette ligne a deja ete faite normalement ?
          endif
       ENDDO
    end DO loop_level

  END SUBROUTINE HBTM

end module HBTM_m
1	guez	37	module HBTM_m
2	guez	3
3	guez	37	IMPLICIT none
4	guez	3
5	guez	37	contains
6	guez	3
7	guez	149	SUBROUTINE HBTM(knon, paprs, pplay, t2m, q2m, ustar, flux_t, &
8	guez	37	flux_q, u, v, t, q, pblh, cape, EauLiq, ctei, pblT, therm, trmb1, &
9			trmb2, trmb3, plcl)
10	guez	3
11	guez	149	! D'apr\'es Holstag et Boville et Troen et Mahrt
12	guez	37	! JAS 47 BLM
13	guez	149	! Algorithme th\'ese Anne Mathieu
14			! Crit\'ere d'entra\^inement Peter Duynkerke (JAS 50)
15	guez	37	! written by: Anne MATHIEU and Alain LAHELLEC, 22nd November 1999
16			! features : implem. exces Mathieu
17	guez	3
18	guez	37	! modifications : decembre 99 passage th a niveau plus bas. voir fixer
19			! la prise du th a z/Lambda = -.2 (max Ray)
20			! Autre algo : entrainement ~ Theta+v =cste mais comment=>The?
21			! on peut fixer q a .7 qsat (cf. non adiabatique) => T2 et The2
22			! voir aussi //KE pblh = niveau The_e ou l = env.
23	guez	3
24	guez	37	! fin therm a la HBTM passage a forme Mathieu 12/09/2001
25	guez	3
26	guez	149	! Adaptation a LMDZ version couplee Pour le moment on fait passer
27			! en argument les grandeurs de surface : flux, t, q2m, t, on va
28			! utiliser systematiquement les grandeurs a 2m mais on garde la
29	guez	178	! possibilit\'e de changer si besoin est (jusqu'\`a pr\'esent la
30	guez	149	! forme de HB avec le 1er niveau modele etait conservee)
31	guez	3
32	guez	71	USE dimphy, ONLY: klev, klon
33	guez	178	USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rtt
34	guez	49	USE yoethf_m, ONLY: r2es, rvtmp2
35			USE fcttre, ONLY: foeew
36
37	guez	37	! Arguments:
38	guez	3
39	guez	37	! nombre de points a calculer
40			INTEGER, intent(in):: knon
41	guez	3
42	guez	37	REAL, intent(in):: t2m(klon) ! temperature a 2 m
43			! q a 2 et 10m
44	guez	149	REAL q2m(klon)
45	guez	37	REAL ustar(klon)
46			! pression a inter-couche (Pa)
47			REAL paprs(klon, klev+1)
48			! pression au milieu de couche (Pa)
49			REAL pplay(klon, klev)
50			! Flux
51			REAL flux_t(klon, klev), flux_q(klon, klev)
52			! vitesse U (m/s)
53			REAL u(klon, klev)
54			! vitesse V (m/s)
55			REAL v(klon, klev)
56			! temperature (K)
57			REAL t(klon, klev)
58			! vapeur d'eau (kg/kg)
59			REAL q(klon, klev)
60	guez	3
61	guez	37	INTEGER isommet
62			! limite max sommet pbl
63			PARAMETER (isommet=klev)
64			REAL vk
65			! Von Karman => passer a .41 ! cf U.Olgstrom
66			PARAMETER (vk=0.35)
67			REAL ricr
68			PARAMETER (ricr=0.4)
69			REAL fak
70			! b calcul du Prandtl et de dTetas
71			PARAMETER (fak=8.5)
72			REAL fakn
73			! a
74			PARAMETER (fakn=7.2)
75			REAL onet
76			PARAMETER (onet=1.0/3.0)
77			REAL t_coup
78			PARAMETER(t_coup=273.15)
79			REAL zkmin
80			PARAMETER (zkmin=0.01)
81			REAL betam
82			! pour Phim / h dans la S.L stable
83			PARAMETER (betam=15.0)
84			REAL betah
85			PARAMETER (betah=15.0)
86			REAL betas
87			! Phit dans la S.L. stable (mais 2 formes /
88			PARAMETER (betas=5.0)
89			! z/OBL<>1
90			REAL sffrac
91			! S.L. = z/h < .1
92			PARAMETER (sffrac=0.1)
93			REAL binm
94			PARAMETER (binm=betam*sffrac)
95			REAL binh
96			PARAMETER (binh=betah*sffrac)
97			REAL ccon
98			PARAMETER (ccon=faksffracvk)
99
100			REAL q_star, t_star
101			! Lambert correlations T' q' avec T* q*
102			REAL b1, b2, b212, b2sr
103			PARAMETER (b1=70., b2=20.)
104
105			REAL z(klon, klev)
106
107			REAL zref
108			! Niveau de ref a 2m peut eventuellement
109			PARAMETER (zref=2.)
110			! etre choisi a 10m
111
112			INTEGER i, k
113			REAL zxt
114			! surface kinematic heat flux [mK/s]
115			REAL khfs(klon)
116			! sfc kinematic constituent flux [m/s]
117			REAL kqfs(klon)
118			! surface virtual heat flux
119			REAL heatv(klon)
120			! bulk Richardon no. mais en Theta_v
121			REAL rhino(klon, klev)
122			! pts w/unstbl pbl (positive virtual ht flx)
123			LOGICAL unstbl(klon)
124			LOGICAL check(klon) ! Richardson number > critical
125			! flag de prolongerment cape pour pt Omega
126			LOGICAL omegafl(klon)
127			REAL pblh(klon)
128			REAL pblT(klon)
129			REAL plcl(klon)
130
131			! Monin-Obukhov lengh
132			REAL obklen(klon)
133
134			REAL zdu2
135			! thermal virtual temperature excess
136			REAL therm(klon)
137			REAL trmb1(klon), trmb2(klon), trmb3(klon)
138			! Algorithme thermique
139			REAL s(klon, klev) ! [P/Po]^Kappa milieux couches
140			! total water of thermal
141			REAL qT_th(klon)
142			! T thermique niveau precedent
143			REAL qsatbef(klon)
144			! le thermique est sature
145			LOGICAL Zsat(klon)
146			! Cape du thermique
147			REAL Cape(klon)
148			! Eau liqu integr du thermique
149			REAL EauLiq(klon)
150			! Critere d'instab d'entrainmt des nuages de
151			REAL ctei(klon)
152	guez	178	REAL zthvd, zthvu, qqsat
153	guez	149	REAL t2
154	guez	37
155			! inverse phi function for momentum
156			REAL phiminv(klon)
157			! turbulent velocity scale for momentum
158			REAL wm(klon)
159			! current level height + one level up
160			REAL zp(klon)
161	guez	149	REAL zcor
162	guez	37
163	guez	178	REAL pblmin
164	guez	37
165			!-----------------------------------------------------------------
166
167			! initialisations
168			q_star = 0
169			t_star = 0
170
171			b212=sqrt(b1*b2)
172			b2sr=sqrt(b2)
173
174			! Calculer les hauteurs de chaque couche
175			! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
176			! pourquoi ne pas utiliser Phi/RG ?
177			DO i = 1, knon
178			z(i, 1) = RD * t(i, 1) / (0.5*(paprs(i, 1)+pplay(i, 1))) &
179			* (paprs(i, 1)-pplay(i, 1)) / RG
180			s(i, 1) = (pplay(i, 1)/paprs(i, 1))**RKappa
181			ENDDO
182			! s(k) = [pplay(k)/ps]^kappa
183			! + + + + + + + + + pplay <-> s(k) t dp=pplay(k-1)-pplay(k)
184			! ----------------- paprs <-> sig(k)
185			! + + + + + + + + + pplay <-> s(k-1)
186			! + + + + + + + + + pplay <-> s(1) t dp=paprs-pplay z(1)
187			! ----------------- paprs <-> sig(1)
188
189			DO k = 2, klev
190			DO i = 1, knon
191			z(i, k) = z(i, k-1) &
192			+ RD * 0.5*(t(i, k-1)+t(i, k)) / paprs(i, k) &
193			* (pplay(i, k-1)-pplay(i, k)) / RG
194			s(i, k) = (pplay(i, k) / paprs(i, 1))**RKappa
195			ENDDO
196			ENDDO
197
198			! Determination des grandeurs de surface
199			DO i = 1, knon
200			! Niveau de ref choisi a 2m
201			zxt = t2m(i)
202
203			! convention >0 vers le bas ds lmdz
204			khfs(i) = - flux_t(i, 1)zxtRd / (RCPD*paprs(i, 1))
205			kqfs(i) = - flux_q(i, 1)zxtRd / (paprs(i, 1))
206			! verifier que khfs et kqfs sont bien de la forme w'l'
207			heatv(i) = khfs(i) + 0.608zxtkqfs(i)
208			! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
209			! Theta et qT du thermique sans exces (interpolin vers surf)
210			! chgt de niveau du thermique (jeudi 30/12/1999)
211			! (interpolation lineaire avant integration phi_h)
212			qT_th(i) = q2m(i)
213			ENDDO
214
215			DO i = 1, knon
216			! Global Richardson
217			rhino(i, 1) = 0.0
218			check(i) = .TRUE.
219			! on initialise pblh a l'altitude du 1er niv
220			pblh(i) = z(i, 1)
221			plcl(i) = 6000.
222			! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
223			obklen(i) = -t(i, 1)ustar(i)3/(RGvk*heatv(i))
224			trmb1(i) = 0.
225			trmb2(i) = 0.
226			trmb3(i) = 0.
227			ENDDO
228
229			! PBL height calculation: Search for level of pbl. Scan upward
230			! until the Richardson number between the first level and the
231			! current level exceeds the "critical" value. (bonne idee Nu de
232			! separer le Ric et l'exces de temp du thermique)
233			DO k = 2, isommet
234			DO i = 1, knon
235			IF (check(i)) THEN
236			! pourquoi / niveau 1 (au lieu du sol) et le terme en u*^2 ?
237			zdu2 = u(i, k)2+v(i, k)2
238			zdu2 = max(zdu2, 1.0e-20)
239			! Theta_v environnement
240			zthvd=t(i, k)/s(i, k)(1.+RETVq(i, k))
241
242			! therm Theta_v sans exces (avec hypothese fausse de H&B, sinon,
243			! passer par Theta_e et virpot)
244			zthvu = T2m(i)(1.+RETVqT_th(i))
245			! Le Ri par Theta_v
246			! On a nveau de ref a 2m ???
247			rhino(i, k) = (z(i, k)-zref)RG(zthvd-zthvu) &
248			/(zdu20.5(zthvd+zthvu))
249
250			IF (rhino(i, k).GE.ricr) THEN
251			pblh(i) = z(i, k-1) + (z(i, k-1)-z(i, k)) * &
252			(ricr-rhino(i, k-1))/(rhino(i, k-1)-rhino(i, k))
253			! test04
254			pblh(i) = pblh(i) + 100.
255			pblT(i) = t(i, k-1) + (t(i, k)-t(i, k-1)) * &
256			(pblh(i)-z(i, k-1))/(z(i, k)-z(i, k-1))
257			check(i) = .FALSE.
258			ENDIF
259			ENDIF
260			ENDDO
261			ENDDO
262
263			! Set pbl height to maximum value where computation exceeds number of
264			! layers allowed
265			DO i = 1, knon
266			if (check(i)) pblh(i) = z(i, isommet)
267			ENDDO
268
269			! Improve estimate of pbl height for the unstable points.
270			! Find unstable points (sensible heat flux is upward):
271			DO i = 1, knon
272			IF (heatv(i) > 0.) THEN
273			unstbl(i) = .TRUE.
274			check(i) = .TRUE.
275			ELSE
276			unstbl(i) = .FALSE.
277			check(i) = .FALSE.
278			ENDIF
279			ENDDO
280
281			! For the unstable case, compute velocity scale and the
282			! convective temperature excess:
283			DO i = 1, knon
284			IF (check(i)) THEN
285			phiminv(i) = (1.-binmpblh(i)/obklen(i))*onet
286
287			! CALCUL DE wm
288			! Ici on considerera que l'on est dans la couche de surf jusqu'a 100
289			! On prend svt couche de surface=0.1*h mais on ne connait pas h
290			! Dans la couche de surface
291			wm(i)= ustar(i)*phiminv(i)
292
293			! forme Mathieu :
294			q_star = kqfs(i)/wm(i)
295			t_star = khfs(i)/wm(i)
296
297			therm(i) = sqrt( b1(1.+2.RETVqT_th(i))t_star**2 &
298			+ (RETVT2m(i))2b2q_star*2 &
299			+ max(0., 2.RETVT2m(i)b212q_star*t_star))
300
301			! Theta et qT du thermique (forme H&B) avec exces
302			! (attention, on ajoute therm(i) qui est virtuelle ...)
303			! pourquoi pas sqrt(b1)*t_star ?
304			qT_th(i) = qT_th(i) + b2sr*q_star
305			! new on differre le calcul de Theta_e
306			rhino(i, 1) = 0.
307			ENDIF
308			ENDDO
309
310			! Improve pblh estimate for unstable conditions using the
311			! convective temperature excess :
312			DO k = 2, isommet
313			DO i = 1, knon
314			IF (check(i)) THEN
315			zdu2 = u(i, k)2 + v(i, k)2
316			zdu2 = max(zdu2, 1e-20)
317			! Theta_v environnement
318			zthvd=t(i, k)/s(i, k)(1.+RETVq(i, k))
319
320			! et therm Theta_v (avec hypothese de constance de H&B,
321			zthvu = T2m(i)(1.+RETVqT_th(i)) + therm(i)
322
323			! Le Ri par Theta_v
324			! Niveau de ref 2m
325			rhino(i, k) = (z(i, k)-zref)RG(zthvd-zthvu) &
326			/(zdu20.5(zthvd+zthvu))
327
328			IF (rhino(i, k).GE.ricr) THEN
329			pblh(i) = z(i, k-1) + (z(i, k-1)-z(i, k)) * &
330			(ricr-rhino(i, k-1))/(rhino(i, k-1)-rhino(i, k))
331			! test04
332			pblh(i) = pblh(i) + 100.
333			pblT(i) = t(i, k-1) + (t(i, k)-t(i, k-1)) * &
334			(pblh(i)-z(i, k-1))/(z(i, k)-z(i, k-1))
335			check(i) = .FALSE.
336			ENDIF
337			ENDIF
338			ENDDO
339			ENDDO
340
341			! Set pbl height to maximum value where computation exceeds number of
342			! layers allowed
343			DO i = 1, knon
344			if (check(i)) pblh(i) = z(i, isommet)
345			ENDDO
346
347			! PBL height must be greater than some minimum mechanical mixing depth
348			! Several investigators have proposed minimum mechanical mixing depth
349			! relationships as a function of the local friction velocity, u*. We
350			! make use of a linear relationship of the form h = c u* where c=700.
351			! The scaling arguments that give rise to this relationship most often
352			! represent the coefficient c as some constant over the local coriolis
353			! parameter. Here we make use of the experimental results of Koracin
354			! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
355			! where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
356			! latitude value for f so that c = 0.07/f = 700.
357			DO i = 1, knon
358			pblmin = 700. * ustar(i)
359			pblh(i) = MAX(pblh(i), pblmin)
360			! par exemple :
361			pblT(i) = t(i, 2) + (t(i, 3)-t(i, 2)) * &
362			(pblh(i)-z(i, 2))/(z(i, 3)-z(i, 2))
363			ENDDO
364
365			! pblh is now available; do preparation for diffusivity calculation:
366			DO i = 1, knon
367			check(i) = .TRUE.
368			Zsat(i) = .FALSE.
369			! omegafl utilise pour prolongement CAPE
370			omegafl(i) = .FALSE.
371			Cape(i) = 0.
372			EauLiq(i) = 0.
373			CTEI(i) = 0.
374
375			! Do additional preparation for unstable cases only, set temperature
376			! and moisture perturbations depending on stability.
377			! Remarque : les formule sont prises dans leur forme CS
378			IF (unstbl(i)) THEN
379			! Niveau de ref du thermique
380			zxt=(T2m(i)-zref0.5RG/RCPD/(1.+RVTMP2*qT_th(i))) &
381			(1.+RETVqT_th(i))
382	guez	3	phiminv(i) = (1. - binmpblh(i)/obklen(i))*onet
383	guez	37	wm(i) = ustar(i)*phiminv(i)
384			ENDIF
385			ENDDO
386	guez	3
387	guez	37	! Main level loop to compute the diffusivities and
388			! counter-gradient terms:
389	guez	49	loop_level: DO k = 2, isommet
390	guez	37	! Find levels within boundary layer:
391			DO i = 1, knon
392			zp(i) = z(i, k)
393			IF (zkmin == 0. .AND. zp(i) > pblh(i)) zp(i) = pblh(i)
394			ENDDO
395
396			! For all layers, compute integral info and CTEI
397			DO i = 1, knon
398	guez	49	if (check(i) .or. omegafl(i)) then
399			if (.not. Zsat(i)) then
400	guez	37	T2 = T2m(i) * s(i, k)
401			! thermodyn functions
402	guez	103	qqsat= r2es * FOEEW(T2, RTT >= T2) / pplay(i, k)
403	guez	37	qqsat=MIN(0.5, qqsat)
404			zcor=1./(1.-retv*qqsat)
405			qqsat=qqsat*zcor
406
407			if (qqsat < qT_th(i)) then
408			! on calcule lcl
409			if (k == 2) then
410			plcl(i) = z(i, k)
411			else
412			plcl(i) = z(i, k-1) + (z(i, k-1)-z(i, k)) &
413			* (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)
414			endif
415			Zsat(i) = .true.
416			endif
417			endif
418			qsatbef(i) = qqsat
419			! cette ligne a deja ete faite normalement ?
420	guez	3	endif
421	guez	37	ENDDO
422	guez	49	end DO loop_level
423	guez	37
424			END SUBROUTINE HBTM
425
426			end module HBTM_m