Sources/phylmd/hbtm.f

module HBTM_m

  IMPLICIT none

contains

  SUBROUTINE HBTM(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:

    ! 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)
    REAL, intent(in):: flux_t(:), flux_q(:) ! (knon) flux à la surface

    ! 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 knon ! nombre de points a calculer
    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)
    ! a
    REAL onet
    PARAMETER (onet=1.0/3.0)
    REAL zkmin
    PARAMETER (zkmin=0.01)
    REAL betam
    ! pour Phim / h dans la S.L stable
    PARAMETER (betam=15.0)
    ! z/OBL<>1
    REAL sffrac
    ! S.L. = z/h < .1
    PARAMETER (sffrac=0.1)
    REAL binm
    PARAMETER (binm=betam*sffrac)

    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

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

    knon = size(pblh)

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