My Project
 All Classes Files Functions Variables Macros
hbtm.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4 
5  SUBROUTINE hbtm(knon, paprs, pplay,
6  . t2m,t10m,q2m,q10m,ustar,
7  . flux_t,flux_q,u,v,t,q,
8  . pblh,cape,eauliq,ctei,pblt,
9  . therm,trmb1,trmb2,trmb3,plcl)
10  USE dimphy
11  IMPLICIT none
12 
13 c***************************************************************
14 c* *
15 c* HBTM2 D'apres Holstag&Boville et Troen&Mahrt *
16 c* JAS 47 BLM *
17 c* Algorithme These Anne Mathieu *
18 c* Critere d'Entrainement Peter Duynkerke (JAS 50) *
19 c* written by : Anne MATHIEU & Alain LAHELLEC, 22/11/99 *
20 c* features : implem. exces Mathieu *
21 c***************************************************************
22 c* mods : decembre 99 passage th a niveau plus bas. voir fixer *
23 c* la prise du th a z/Lambda = -.2 (max Ray) *
24 c* Autre algo : entrainement ~ Theta+v =cste mais comment=>The?*
25 c* on peut fixer q a .7qsat(cf non adiab)=>T2 et The2 *
26 c* voir aussi //KE pblh = niveau The_e ou l = env. *
27 c***************************************************************
28 c* fin therm a la HBTM passage a forme Mathieu 12/09/2001 *
29 c***************************************************************
30 c*
31 c
32 c
33 cAM Fev 2003
34 c Adaptation a LMDZ version couplee
35 c
36 c Pour le moment on fait passer en argument les grdeurs de surface :
37 c flux, t,q2m, t,q10m, on va utiliser systematiquement les grdeurs a 2m ms
38 c on garde la possibilite de changer si besoin est (jusqu'a present la
39 c forme de HB avec le 1er niveau modele etait conservee)
40 c
41 c
42 c
43 c
44 c
45 cym#include "dimensions.h"
46 cym#include "dimphy.h"
47 #include "YOMCST.h"
48  REAL rlvcp, reps
49 c Arguments:
50 c
51  INTEGER knon ! nombre de points a calculer
52 cAM
53  REAL t2m(klon), t10m(klon) ! temperature a 2 et 10m
54  REAL q2m(klon), q10m(klon) ! q a 2 et 10m
55  REAL ustar(klon)
56  REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
57  REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
58  REAL flux_t(klon,klev), flux_q(klon,klev) ! Flux
59  REAL u(klon,klev) ! vitesse U (m/s)
60  REAL v(klon,klev) ! vitesse V (m/s)
61  REAL t(klon,klev) ! temperature (K)
62  REAL q(klon,klev) ! vapeur d'eau (kg/kg)
63 cAM REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
64 cAM REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
65 c
66  INTEGER isommet
67 cum PARAMETER (isommet=klev) ! limite max sommet pbl
68  REAL vk
69  parameter(vk=0.35) ! Von Karman => passer a .41 ! cf U.Olgstrom
70  REAL ricr
71  parameter(ricr=0.4)
72  REAL fak
73  parameter(fak=8.5) ! b calcul du Prandtl et de dTetas
74  REAL fakn
75  parameter(fakn=7.2) ! a
76  REAL onet
77  parameter(onet=1.0/3.0)
78  REAL t_coup
79  parameter(t_coup=273.15)
80  REAL zkmin
81  parameter(zkmin=0.01)
82  REAL betam
83  parameter(betam=15.0) ! pour Phim / h dans la S.L stable
84  REAL betah
85  parameter(betah=15.0)
86  REAL betas
87  parameter(betas=5.0) ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
88  REAL sffrac
89  parameter(sffrac=0.1) ! S.L. = z/h < .1
90  REAL binm
91  parameter(binm=betam*sffrac)
92  REAL binh
93  parameter(binh=betah*sffrac)
94  REAL ccon
95  parameter(ccon=fak*sffrac*vk)
96 c
97  REAL q_star,t_star
98  REAL b1,b2,b212,b2sr ! Lambert correlations T' q' avec T* q*
99  parameter(b1=70.,b2=20.)
100 c
101  REAL z(klon,klev)
102 cAM REAL pcfm(klon,klev), pcfh(klon,klev)
103 cAM
104  REAL zref
105  parameter(zref=2.) ! Niveau de ref a 2m peut eventuellement
106 c etre choisi a 10m
107 cMA
108 c
109  INTEGER i, k, j
110  REAL zxt
111 cAM REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
112 cAM REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
113  REAL khfs(klon) ! surface kinematic heat flux [mK/s]
114  REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
115  REAL heatv(klon) ! surface virtual heat flux
116  REAL rhino(klon,klev) ! bulk Richardon no. mais en Theta_v
117  LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
118  LOGICAL stblev(klon) ! stable pbl with levels within pbl
119  LOGICAL unslev(klon) ! unstbl pbl with levels within pbl
120  LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
121  LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
122  LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
123  LOGICAL omegafl(klon) ! flag de prolongerment cape pour pt Omega
124  REAL pblh(klon)
125  REAL pblt(klon)
126  REAL plcl(klon)
127 cAM REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
128 cAM REAL cgq(klon,2:klev) ! counter-gradient term for constituents
129 cAM REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
130  REAL obklen(klon) ! Monin-Obukhov lengh
131 cAM REAL ztvd, ztvu,
132  REAL zdu2
133  REAL therm(klon) ! thermal virtual temperature excess
134  REAL trmb1(klon),trmb2(klon),trmb3(klon)
135 C Algorithme thermique
136  REAL s(klon,klev) ! [P/Po]^Kappa milieux couches
137  REAL th_th(klon) ! potential temperature of thermal
138  REAL the_th(klon) ! equivalent potential temperature of thermal
139  REAL qt_th(klon) ! total water of thermal
140  REAL tbef(klon) ! T thermique niveau precedent
141  REAL qsatbef(klon)
142  LOGICAL zsat(klon) ! le thermique est sature
143  REAL cape(klon) ! Cape du thermique
144  REAL kape(klon) ! Cape locale
145  REAL eauliq(klon) ! Eau liqu integr du thermique
146  REAL ctei(klon) ! Critere d'instab d'entrainmt des nuages de CL
147  REAL the1,the2,aa,bb,zthvd,zthvu,xintpos,qqsat
148 cIM 091204 BEG
149  REAL a1,a2,a3
150 cIM 091204 END
151  REAL xhis,rnum,denom,th1,th2,thv1,thv2,ql2
152  REAL dqsat_dt,qsat2,qt1,q2,t1,t2,xnull,delt_the
153  REAL delt_qt,delt_2,quadsat,spblh,reduc
154 c
155  REAL phiminv(klon) ! inverse phi function for momentum
156  REAL phihinv(klon) ! inverse phi function for heat
157  REAL wm(klon) ! turbulent velocity scale for momentum
158  REAL fak1(klon) ! k*ustar*pblh
159  REAL fak2(klon) ! k*wm*pblh
160  REAL fak3(klon) ! fakn*wstr/wm
161  REAL pblk(klon) ! level eddy diffusivity for momentum
162  REAL pr(klon) ! Prandtl number for eddy diffusivities
163  REAL zl(klon) ! zmzp / Obukhov length
164  REAL zh(klon) ! zmzp / pblh
165  REAL zzh(klon) ! (1-(zmzp/pblh))**2
166  REAL wstr(klon) ! w*, convective velocity scale
167  REAL zm(klon) ! current level height
168  REAL zp(klon) ! current level height + one level up
169  REAL zcor, zdelta, zcvm5
170 cAM REAL zxqs
171  REAL fac, pblmin, zmzp, term
172 c
173 #include "YOETHF.h"
174 #include "FCTTRE.h"
175 
176 
177 
178 ! initialisations (Anne)
179  isommet=klev
180  th_th(:) = 0.
181  q_star = 0
182  t_star = 0
183 
184 
185  b212=sqrt(b1*b2)
186  b2sr=sqrt(b2)
187 c
188 C ============================================================
189 C Fonctions thermo implicites
190 C ============================================================
191 c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
192 c Tetens : pression partielle de vap d'eau e_sat(T)
193 c =================================================
194 C++ e_sat(T) = r2*exp( r3*(T-Tf)/(T-r4) ) id a r2*FOEWE
195 C++ avec :
196 C++ Tf = 273.16 K (Temp de fusion de la glace)
197 C++ r2 = 611.14 Pa
198 C++ r3 = 17.269 (liquide) 21.875 (solide) adim
199 C++ r4 = 35.86 7.66 Kelvin
200 C++ q_sat = eps*e_sat/(p-(1-eps)*e_sat)
201 C++ deriv� :
202 C++ =========
203 C++ r3*(Tf-r4)*q_sat(T,p)
204 C++ d_qsat_dT = --------------------------------
205 C++ (T-r4)^2*( 1-(1-eps)*e_sat(T)/p )
206 c++ pour zcvm5=Lv, c'est FOEDE
207 c++ Rq :(1.-REPS)*esarg/Parg id a RETV*Qsat
208 C ------------------------------------------------------------------
209 c
210 c Initialisation
211  rlvcp = rlvtt/rcpd
212  reps = rd/rv
213 
214 c
215 c DO i = 1, klon
216 c pcfh(i,1) = cd_h(i)
217 c pcfm(i,1) = cd_m(i)
218 c ENDDO
219 c DO k = 2, klev
220 c DO i = 1, klon
221 c pcfh(i,k) = zkmin
222 c pcfm(i,k) = zkmin
223 c cgs(i,k) = 0.0
224 c cgh(i,k) = 0.0
225 c cgq(i,k) = 0.0
226 c ENDDO
227 c ENDDO
228 c
229 c Calculer les hauteurs de chaque couche
230 c (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
231 c pourquoi ne pas utiliser Phi/RG ?
232  DO i = 1, knon
233  z(i,1) = rd * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
234  . * (paprs(i,1)-pplay(i,1)) / rg
235  s(i,1) = (pplay(i,1)/paprs(i,1))**rkappa
236  ENDDO
237 c s(k) = [pplay(k)/ps]^kappa
238 c + + + + + + + + + pplay <-> s(k) t dp=pplay(k-1)-pplay(k)
239 c
240 c ----------------- paprs <-> sig(k)
241 c
242 c + + + + + + + + + pplay <-> s(k-1)
243 c
244 c
245 c + + + + + + + + + pplay <-> s(1) t dp=paprs-pplay z(1)
246 c
247 c ----------------- paprs <-> sig(1)
248 c
249  DO k = 2, klev
250  DO i = 1, knon
251  z(i,k) = z(i,k-1)
252  . + rd * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
253  . * (pplay(i,k-1)-pplay(i,k)) / rg
254  s(i,k) = (pplay(i,k)/paprs(i,1))**rkappa
255  ENDDO
256  ENDDO
257 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
258 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
259 c +++ Determination des grandeurs de surface +++++++++++++++++++++
260 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
261 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
262  DO i = 1, knon
263 cAM IF (thermcep) THEN
264 cAM zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
265 c zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
266 c zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
267 cAM zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
268 cAM zxqs=MIN(0.5,zxqs)
269 cAM zcor=1./(1.-retv*zxqs)
270 cAM zxqs=zxqs*zcor
271 cAM ELSE
272 cAM IF (tsol(i).LT.t_coup) THEN
273 cAM zxqs = qsats(tsol(i)) / paprs(i,1)
274 cAM ELSE
275 cAM zxqs = qsatl(tsol(i)) / paprs(i,1)
276 cAM ENDIF
277 cAM ENDIF
278 c niveau de reference bulk; mais ici, c,a pourrait etre le niveau de ref du thermique
279 cAM zx_alf1 = 1.0
280 cAM zx_alf2 = 1.0 - zx_alf1
281 cAM zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
282 cAM . *(1.+RETV*q(i,1))*zx_alf1
283 cAM . + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
284 cAM . *(1.+RETV*q(i,2))*zx_alf2
285 cAM zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
286 cAM zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
287 cAM zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
288 cAM
289 cAMAM zxu = u10m(i)
290 cAMAM zxv = v10m(i)
291 cAMAM zxmod = 1.0+SQRT(zxu**2+zxv**2)
292 cAM Niveau de ref choisi a 2m
293  zxt = t2m(i)
294 
295 c ***************************************************
296 c attention, il doit s'agir de <w'theta'>
297 c ;Calcul de tcls virtuel et de w'theta'virtuel
298 c ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
299 c tcls=tcls*(1+.608*qcls)
300 c
301 c ;Pour avoir w'theta',
302 c ; il faut diviser par ro.Cp
303 c Cp=Cpd*(1+0.84*qcls)
304 c fcs=fcs/(ro_surf*Cp)
305 c ;On transforme w'theta' en w'thetav'
306 c Lv=(2.501-0.00237*(tcls-273.15))*1.E6
307 c xle=xle/(ro_surf*Lv)
308 c fcsv=fcs+.608*xle*tcls
309 c ***************************************************
310 cAM khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
311 cAM kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
312 cAM
313 cdif khfs est deja w't'_v / heatv(i) = khfs(i) + RETV*zxt*kqfs(i)
314 cAM calcule de Ro = paprs(i,1)/Rd zxt
315 cAM convention >0 vers le bas ds lmdz
316  khfs(i) = - flux_t(i,1)*zxt*rd / (rcpd*paprs(i,1))
317  kqfs(i) = - flux_q(i,1)*zxt*rd / (paprs(i,1))
318 cAM verifier que khfs et kqfs sont bien de la forme w'l'
319  heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
320 c a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
321 cAM heatv(i) = khfs(i)
322 cAM ustar est en entree
323 cAM taux = zxu *zxmod*cd_m(i)
324 cAM tauy = zxv *zxmod*cd_m(i)
325 cAM ustar(i) = SQRT(taux**2+tauy**2)
326 cAM ustar(i) = MAX(SQRT(ustar(i)),0.01)
327 c Theta et qT du thermique sans exces (interpolin vers surf)
328 c chgt de niveau du thermique (jeudi 30/12/1999)
329 c (interpolation lineaire avant integration phi_h)
330 cAM qT_th(i) = zxqs*beta(i) + 4./z(i,1)*(q(i,1)-zxqs*beta(i))
331 cAM qT_th(i) = max(qT_th(i),q(i,1))
332  qt_th(i) = q2m(i)
333 cn The_th restera la Theta du thermique sans exces jusqu'a 2eme calcul
334 cn reste a regler convention P) pour Theta
335 c The_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
336 c - + RLvCp*qT_th(i)
337 cAM Th_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
338  th_th(i) = t2m(i)
339  ENDDO
340 c
341  DO i = 1, knon
342  rhino(i,1) = 0.0 ! Global Richardson
343  check(i) = .true.
344  pblh(i) = z(i,1) ! on initialise pblh a l'altitude du 1er niveau
345  plcl(i) = 6000.
346 c Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
347  obklen(i) = -t(i,1)*ustar(i)**3/(rg*vk*heatv(i))
348  trmb1(i) = 0.
349  trmb2(i) = 0.
350  trmb3(i) = 0.
351  ENDDO
352 
353 C
354 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
355 C PBL height calculation:
356 C Search for level of pbl. Scan upward until the Richardson number between
357 C the first level and the current level exceeds the "critical" value.
358 C (bonne idee Nu de separer le Ric et l'exces de temp du thermique)
359 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
360  fac = 100.0
361  DO k = 2, isommet
362  DO i = 1, knon
363  IF (check(i)) THEN
364 ! pourquoi / niveau 1 (au lieu du sol) et le terme en u*^2 ?
365 ctest zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
366  zdu2 = u(i,k)**2+v(i,k)**2
367  zdu2 = max(zdu2,1.0e-20)
368 c Theta_v environnement
369  zthvd=t(i,k)/s(i,k)*(1.+retv*q(i,k))
370 c
371 c therm Theta_v sans exces (avec hypothese fausse de H&B, sinon,
372 c passer par Theta_e et virpot)
373 c zthvu=t(i,1)/s(i,1)*(1.+RETV*q(i,1))
374 cAM zthvu = Th_th(i)*(1.+RETV*q(i,1))
375  zthvu = th_th(i)*(1.+retv*qt_th(i))
376 c Le Ri par Theta_v
377 cAM rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
378 cAM . /(zdu2*0.5*(zthvd+zthvu))
379 cAM On a nveau de ref a 2m ???
380  rhino(i,k) = (z(i,k)-zref)*rg*(zthvd-zthvu)
381  . /(zdu2*0.5*(zthvd+zthvu))
382 c
383  IF (rhino(i,k).GE.ricr) THEN
384  pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
385  . (ricr-rhino(i,k-1))/(rhino(i,k-1)-rhino(i,k))
386 c test04
387  pblh(i) = pblh(i) + 100.
388  pblt(i) = t(i,k-1) + (t(i,k)-t(i,k-1)) *
389  . (pblh(i)-z(i,k-1))/(z(i,k)-z(i,k-1))
390  check(i) = .false.
391  ENDIF
392  ENDIF
393  ENDDO
394  ENDDO
395 
396 C
397 C Set pbl height to maximum value where computation exceeds number of
398 C layers allowed
399 C
400  DO i = 1, knon
401  if (check(i)) pblh(i) = z(i,isommet)
402  ENDDO
403 C
404 C Improve estimate of pbl height for the unstable points.
405 C Find unstable points (sensible heat flux is upward):
406 C
407  DO i = 1, knon
408  IF (heatv(i) .GT. 0.) THEN
409  unstbl(i) = .true.
410  check(i) = .true.
411  ELSE
412  unstbl(i) = .false.
413  check(i) = .false.
414  ENDIF
415  ENDDO
416 C
417 C For the unstable case, compute velocity scale and the
418 C convective temperature excess:
419 C
420  DO i = 1, knon
421  IF (check(i)) THEN
422  phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
423 c ***************************************************
424 c Wm ? et W* ? c'est la formule pour z/h < .1
425 c ;Calcul de w* ;;
426 c ;;;;;;;;;;;;;;;;
427 c w_star=((g/tcls)*fcsv*z(ind))^(1/3.) [ou prendre la premiere approx de h)
428 c ;; CALCUL DE wm ;;
429 c ;;;;;;;;;;;;;;;;;;
430 c ; Ici on considerera que l'on est dans la couche de surf jusqu'a 100m
431 c ; On prend svt couche de surface=0.1*h mais on ne connait pas h
432 c ;;;;;;;;;;;Dans la couche de surface
433 c if (z(ind) le 20) then begin
434 c Phim=(1.-15.*(z(ind)/L))^(-1/3.)
435 c wm=u_star/Phim
436 c ;;;;;;;;;;;En dehors de la couche de surface
437 c endif else if (z(ind) gt 20) then begin
438 c wm=(u_star^3+c1*w_star^3)^(1/3.)
439 c endif
440 c ***************************************************
441  wm(i)= ustar(i)*phiminv(i)
442 c======================================================================
443 cvaleurs de Dominique Lambert de la campagne SEMAPHORE :
444 c <T'^2> = 100.T*^2; <q'^2> = 20.q*^2 a 10m
445 c <Tv'^2> = (1+1.2q).100.T* + 1.2Tv.sqrt(20*100).T*.q* + (.608*Tv)^2*20.q*^2;
446 c et dTetavS = sqrt(<Tv'^2>) ainsi calculee.
447 c avec : T*=<w'T'>_s/w* et q*=<w'q'>/w*
448 c !!! on peut donc utiliser w* pour les fluctuations <-> Lambert
449 c(leur corellation pourrait dependre de beta par ex)
450 c if fcsv(i,j) gt 0 then begin
451 c dTetavs=b1*(1.+2.*.608*q_10(i,j))*(fcs(i,j)/wm(i,j))^2+$
452 c (.608*Thetav_10(i,j))^2*b2*(xle(i,j)/wm(i,j))^2+$
453 c 2.*.608*thetav_10(i,j)*sqrt(b1*b2)*(xle(i,j)/wm(i,j))*(fcs(i,j)/wm(i,j))
454 c dqs=b2*(xle(i,j)/wm(i,j))^2
455 c theta_s(i,j)=thetav_10(i,j)+sqrt(dTetavs)
456 c q_s(i,j)=q_10(i,j)+sqrt(dqs)
457 c endif else begin
458 c Theta_s(i,j)=thetav_10(i,j)
459 c q_s(i,j)=q_10(i,j)
460 c endelse
461 c======================================================================
462 c
463 cHBTM therm(i) = heatv(i)*fak/wm(i)
464 c forme Mathieu :
465  q_star = kqfs(i)/wm(i)
466  t_star = khfs(i)/wm(i)
467 cIM 091204 BEG
468  IF(1.EQ.0) THEN
469  IF(t_star.LT.0..OR.q_star.LT.0.) THEN
470  print*,'i t_star q_star khfs kqfs wm',i,t_star,q_star,
471  $ khfs(i),kqfs(i),wm(i)
472  ENDIF
473  ENDIF
474 cIM 091204 END
475 cAM Nveau cde ref 2m =>
476 cAM therm(i) = sqrt( b1*(1.+2.*RETV*q(i,1))*t_star**2
477 cAM + + (RETV*T(i,1))**2*b2*q_star**2
478 cAM + + 2.*RETV*T(i,1)*b212*q_star*t_star
479 cAM + )
480 cIM 091204 BEG
481  a1=b1*(1.+2.*retv*qt_th(i))*t_star**2
482  a2=(retv*th_th(i))**2*b2*q_star*q_star
483  a3=2.*retv*th_th(i)*b212*q_star*t_star
484  aa=a1+a2+a3
485  IF(1.EQ.0) THEN
486  IF (aa.LT.0.) THEN
487  print*,'i a1 a2 a3 aa',i,a1,a2,a3,aa
488  print*,'i qT_th Th_th t_star q_star RETV b1 b2 b212',
489  $ i,qt_th(i),th_th(i),t_star,q_star,retv,b1,b2,b212
490  ENDIF
491  ENDIF
492 cIM 091204 END
493  therm(i) = sqrt( b1*(1.+2.*retv*qt_th(i))*t_star**2
494  + + (retv*th_th(i))**2*b2*q_star*q_star
495 cIM 101204 + + 2.*RETV*Th_th(i)*b212*q_star*t_star
496  + + max(0.,2.*retv*th_th(i)*b212*q_star*t_star)
497  + )
498 c
499 c Theta et qT du thermique (forme H&B) avec exces
500 c (attention, on ajoute therm(i) qui est virtuelle ...)
501 c pourquoi pas sqrt(b1)*t_star ?
502 c dqs = b2sr*kqfs(i)/wm(i)
503  qt_th(i) = qt_th(i) + b2sr*q_star
504 cnew on differre le calcul de Theta_e
505 c The_th(i) = The_th(i) + therm(i) + RLvCp*qT_th(i)
506 c ou: The_th(i) = The_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i)
507  rhino(i,1) = 0.0
508  ENDIF
509  ENDDO
510 C
511 c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
512 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
513 C ++ Improve pblh estimate for unstable conditions using the +++++++
514 C ++ convective temperature excess : +++++++
515 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
516 c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
517 C
518  DO k = 2, isommet
519  DO i = 1, knon
520  IF (check(i)) THEN
521 ctest zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
522  zdu2 = u(i,k)**2+v(i,k)**2
523  zdu2 = max(zdu2,1.0e-20)
524 c Theta_v environnement
525  zthvd=t(i,k)/s(i,k)*(1.+retv*q(i,k))
526 c
527 c et therm Theta_v (avec hypothese de constance de H&B,
528 c zthvu=(t(i,1)+therm(i))/s(i,1)*(1.+RETV*q(i,1))
529  zthvu = th_th(i)*(1.+retv*qt_th(i)) + therm(i)
530 
531 c
532 c Le Ri par Theta_v
533 cAM Niveau de ref 2m
534 cAM rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
535 cAM . /(zdu2*0.5*(zthvd+zthvu))
536  rhino(i,k) = (z(i,k)-zref)*rg*(zthvd-zthvu)
537  . /(zdu2*0.5*(zthvd+zthvu))
538 c
539 c
540  IF (rhino(i,k).GE.ricr) THEN
541  pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
542  . (ricr-rhino(i,k-1))/(rhino(i,k-1)-rhino(i,k))
543 c test04
544  pblh(i) = pblh(i) + 100.
545  pblt(i) = t(i,k-1) + (t(i,k)-t(i,k-1)) *
546  . (pblh(i)-z(i,k-1))/(z(i,k)-z(i,k-1))
547  check(i) = .false.
548 cIM 170305 BEG
549  IF(1.EQ.0) THEN
550 c debug print -120;34 -34- 58 et 0;26 wamp
551  if (i.eq.950.or.i.eq.192.or.i.eq.624.or.i.eq.118) then
552  print*,' i,Th_th,Therm,qT :',i,th_th(i),therm(i),qt_th(i)
553  q_star = kqfs(i)/wm(i)
554  t_star = khfs(i)/wm(i)
555  print*,'q* t*, b1,b2,b212 ',q_star,t_star
556  - , b1*(1.+2.*retv*qt_th(i))*t_star**2
557  - , (retv*th_th(i))**2*b2*q_star**2
558  - , 2.*retv*th_th(i)*b212*q_star*t_star
559  print*,'zdu2 ,100.*ustar(i)**2',zdu2 ,fac*ustar(i)**2
560  endif
561  ENDIF !(1.EQ.0) THEN
562 cIM 170305 END
563 c q_star = kqfs(i)/wm(i)
564 c t_star = khfs(i)/wm(i)
565 c trmb1(i) = b1*(1.+2.*RETV*q(i,1))*t_star**2
566 c trmb2(i) = (RETV*T(i,1))**2*b2*q_star**2
567 c Omega now trmb3(i) = 2.*RETV*T(i,1)*b212*q_star*t_star
568  ENDIF
569  ENDIF
570  ENDDO
571  ENDDO
572 C
573 C Set pbl height to maximum value where computation exceeds number of
574 C layers allowed
575 C
576  DO i = 1, knon
577  if (check(i)) pblh(i) = z(i,isommet)
578  ENDDO
579 C
580 C PBL height must be greater than some minimum mechanical mixing depth
581 C Several investigators have proposed minimum mechanical mixing depth
582 C relationships as a function of the local friction velocity, u*. We
583 C make use of a linear relationship of the form h = c u* where c=700.
584 C The scaling arguments that give rise to this relationship most often
585 C represent the coefficient c as some constant over the local coriolis
586 C parameter. Here we make use of the experimental results of Koracin
587 C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
588 C where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
589 C latitude value for f so that c = 0.07/f = 700.
590 C
591  DO i = 1, knon
592  pblmin = 700.0*ustar(i)
593  pblh(i) = max(pblh(i),pblmin)
594 c par exemple :
595  pblt(i) = t(i,2) + (t(i,3)-t(i,2)) *
596  . (pblh(i)-z(i,2))/(z(i,3)-z(i,2))
597  ENDDO
598 
599 C ********************************************************************
600 C pblh is now available; do preparation for diffusivity calculation :
601 C ********************************************************************
602  DO i = 1, knon
603  check(i) = .true.
604  zsat(i) = .false.
605 c omegafl utilise pour prolongement CAPE
606  omegafl(i) = .false.
607  cape(i) = 0.
608  kape(i) = 0.
609  eauliq(i) = 0.
610  ctei(i) = 0.
611  pblk(i) = 0.0
612  fak1(i) = ustar(i)*pblh(i)*vk
613 C
614 C Do additional preparation for unstable cases only, set temperature
615 C and moisture perturbations depending on stability.
616 C *** Rq: les formule sont prises dans leur forme CS ***
617  IF (unstbl(i)) THEN
618 cAM Niveau de ref du thermique
619 cAM zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
620 cAM . *(1.+RETV*q(i,1))
621  zxt=(th_th(i)-zref*0.5*rg/rcpd/(1.+rvtmp2*qt_th(i)))
622  . *(1.+retv*qt_th(i))
623  phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
624  phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
625  wm(i) = ustar(i)*phiminv(i)
626  fak2(i) = wm(i)*pblh(i)*vk
627  wstr(i) = (heatv(i)*rg*pblh(i)/zxt)**onet
628  fak3(i) = fakn*wstr(i)/wm(i)
629  ENDIF
630 c Computes Theta_e for thermal (all cases : to be modified)
631 c attention ajout therm(i) = virtuelle
632  the_th(i) = th_th(i) + therm(i) + rlvcp*qt_th(i)
633 c ou: The_th(i) = Th_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i)
634  ENDDO
635 
636 C Main level loop to compute the diffusivities and
637 C counter-gradient terms:
638 C
639  DO 1000 k = 2, isommet
640 C
641 C Find levels within boundary layer:
642 C
643  DO i = 1, knon
644  unslev(i) = .false.
645  stblev(i) = .false.
646  zm(i) = z(i,k-1)
647  zp(i) = z(i,k)
648  IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
649  IF (zm(i) .LT. pblh(i)) THEN
650  zmzp = 0.5*(zm(i) + zp(i))
651 C debug
652 c if (i.EQ.1864) then
653 c print*,'i,pblh(1864),obklen(1864)',i,pblh(i),obklen(i)
654 c endif
655 
656  zh(i) = zmzp/pblh(i)
657  zl(i) = zmzp/obklen(i)
658  zzh(i) = 0.
659  IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
660 C
661 C stblev for points zm < plbh and stable and neutral
662 C unslev for points zm < plbh and unstable
663 C
664  IF (unstbl(i)) THEN
665  unslev(i) = .true.
666  ELSE
667  stblev(i) = .true.
668  ENDIF
669  ENDIF
670  ENDDO
671 c print*,'fin calcul niveaux'
672 C
673 C Stable and neutral points; set diffusivities; counter-gradient
674 C terms zero for stable case:
675 C
676  DO i = 1, knon
677  IF (stblev(i)) THEN
678  IF (zl(i).LE.1.) THEN
679  pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i))
680  ELSE
681  pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))
682  ENDIF
683 c pcfm(i,k) = pblk(i)
684 c pcfh(i,k) = pcfm(i,k)
685  ENDIF
686  ENDDO
687 C
688 C unssrf, unstable within surface layer of pbl
689 C unsout, unstable within outer layer of pbl
690 C
691  DO i = 1, knon
692  unssrf(i) = .false.
693  unsout(i) = .false.
694  IF (unslev(i)) THEN
695  IF (zh(i).lt.sffrac) THEN
696  unssrf(i) = .true.
697  ELSE
698  unsout(i) = .true.
699  ENDIF
700  ENDIF
701  ENDDO
702 C
703 C Unstable for surface layer; counter-gradient terms zero
704 C
705  DO i = 1, knon
706  IF (unssrf(i)) THEN
707  term = (1. - betam*zl(i))**onet
708  pblk(i) = fak1(i)*zh(i)*zzh(i)*term
709  pr(i) = term/sqrt(1. - betah*zl(i))
710  ENDIF
711  ENDDO
712 c print*,'fin counter-gradient terms zero'
713 C
714 C Unstable for outer layer; counter-gradient terms non-zero:
715 C
716  DO i = 1, knon
717  IF (unsout(i)) THEN
718  pblk(i) = fak2(i)*zh(i)*zzh(i)
719 c cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
720 c cgh(i,k) = khfs(i)*cgs(i,k)
721  pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
722 c cgq(i,k) = kqfs(i)*cgs(i,k)
723  ENDIF
724  ENDDO
725 c print*,'fin counter-gradient terms non zero'
726 C
727 C For all unstable layers, compute diffusivities and ctrgrad ter m
728 C
729 c DO i = 1, knon
730 c IF (unslev(i)) THEN
731 c pcfm(i,k) = pblk(i)
732 c pcfh(i,k) = pblk(i)/pr(i)
733 c etc cf original
734 c ENDIF
735 c ENDDO
736 C
737 C For all layers, compute integral info and CTEI
738 C
739  DO i = 1, knon
740  if (check(i).or.omegafl(i)) then
741  if (.not.zsat(i)) then
742 c Th2 = The_th(i) - RLvCp*qT_th(i)
743  th2 = th_th(i)
744  t2 = th2*s(i,k)
745 c thermodyn functions
746  zdelta=max(0.,sign(1.,rtt-t2))
747  qqsat= r2es * foeew(t2,zdelta)/pplay(i,k)
748  qqsat=min(0.5,qqsat)
749  zcor=1./(1.-retv*qqsat)
750  qqsat=qqsat*zcor
751 c
752  if (qqsat.lt.qt_th(i)) then
753 c on calcule lcl
754  if (k.eq.2) then
755  plcl(i) = z(i,k)
756  else
757  plcl(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
758  . (qt_th(i)-qsatbef(i))/(qsatbef(i)-qqsat)
759  endif
760  zsat(i) = .true.
761  tbef(i) = t2
762  endif
763 c
764  qsatbef(i) = qqsat ! bug dans la version orig ???
765  endif
766 camn ???? cette ligne a deja ete faite normalement ?
767  endif
768 c print*,'hbtm2 i,k=',i,k
769  ENDDO
770  1000 continue ! end of level loop
771 cIM 170305 BEG
772  IF(1.EQ.0) THEN
773  print*,'hbtm2 ok'
774  ENDIF !(1.EQ.0) THEN
775 cIM 170305 END
776  RETURN
777  END