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.

    ! Modifications : d\'ecembre 99 passage th \`a niveau plus bas. Voir fixer
    ! la prise du th \`a z/Lambda = -.2 (max Ray)
    ! Autre algorithme : entra\^inement ~ Theta + v =constante
    !  mais comment ? The ?
    ! On peut fixer q \`a 0.7 qsat (cf. non adiabatique) d'où T2 et The2.
    ! Voir aussi KE pblh = niveau The_e ou l = env.

    ! Adaptation \`a LMDZ version coupl\'ee. Pour le moment on fait
    ! passer en argument les grandeurs de surface : flux, t, q2m. On
    ! va utiliser syst\'ematiquement les grandeurs \`a 2 m mais on
    ! garde la possibilit\'e de changer si besoin (jusqu'\`a pr\'esent
    ! la forme de HB avec le premier niveau mod\`ele \'etait
    ! conserv\'ee).

    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

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

    REAL, intent(out):: pblh(:) ! (knon)
    ! 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 pblT(klon)
    ! thermal virtual temperature excess
    REAL therm(klon)
    REAL trmb1(klon), trmb2(klon), trmb3(klon)
    REAL plcl(klon)

    ! Local:

    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)

    ! Monin-Obukhov lengh
    REAL obklen(klon)

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