Context Navigation

hbtm2l.F90 @ 5220

Last change on this file since 5220 was 5153, checked in by abarral, 3 months ago
Revert FCTTRE to INCLUDE to assess impact of inlining
File size: 26.2 KB

Line
1	! $Header$
2
3	SUBROUTINE hbtm2l(knon, paprs, pplay, t2m, t10m, q2m, q10m, ustar, flux_t, flux_q, u, v, t, q, pblh, therm, plcl, cape, &
4	cin, eauliq, ctei, d_qt, d_thv, dlt_2, xhis, posint, omega, diagok)
5	USE dimphy
6	USE lmdz_yoethf
7
8	USE lmdz_yomcst
9
10	IMPLICIT NONE
11	INCLUDE "FCTTRE.h"
12
13	! ***************************************************************
14	! * *
15	! * HBTM2L D'apres Holstag&Boville et Troen&Mahrt *
16	! * JAS 47 BLM *
17	! * Algorithmes These Anne Mathieu *
18	! * Critere d'Entrainement Peter Duynkerke (JAS 50) *
19	! * written by : Anne MATHIEU & Alain LAHELLEC, 22/11/99 *
20	! * features : implem. exces Mathieu *
21	! ***************************************************************
22	! * mods : decembre 99 passage th a niveau plus bas. voir fixer *
23	! ***************************************************************
24	! * fin therm a la HBTM passage a forme Mathieu 12/09/2001 *
25	! ***************************************************************
26	! AM Fev 2003 Adaptation a LMDZ version couplee *
27	! *
28	! Pour le moment on fait passer en argument les grdeurs de surface :
29	! flux, t,q2m, t,q10m, on va utiliser systematiquement les grdeurs a 2m ms
30	! on garde la possibilite de changer si besoin est (jusqu'a present la
31	! forme de HB avec le 1er niveau modele etait conservee) *
32	! ***************************************************************
33	! * re-ecriture complete Alain Mars 2012 dans LMDZ5V5 *
34	! ***************************************************************
35	REAL rlvcp, reps
36	! Arguments:
37
38	INTEGER knon ! nombre de points a calculer
39	! AM
40	REAL t2m(klon), t10m(klon) ! temperature a 2 et 10m
41	REAL q2m(klon), q10m(klon) ! q a 2 et 10m
42	REAL ustar(klon)
43	REAL paprs(klon, klev + 1) ! pression a inter-couche (Pa)
44	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
45	REAL flux_t(klon, klev), flux_q(klon, klev) ! Flux
46	REAL u(klon, klev) ! vitesse U (m/s)
47	REAL v(klon, klev) ! vitesse V (m/s)
48	REAL t(klon, klev) ! temperature (K)
49	REAL q(klon, klev) ! vapeur d'eau (kg/kg)
50
51	INTEGER isommet
52	REAL vk
53	PARAMETER (vk = 0.35) ! Von Karman => passer a .41 ! cf U.Olgstrom
54	REAL ricr
55	PARAMETER (ricr = 0.4)
56	REAL fak
57	PARAMETER (fak = 8.5) ! b calcul du Prandtl et de dTetas
58	REAL fakn
59	PARAMETER (fakn = 7.2) ! a
60	REAL onet
61	PARAMETER (onet = 1.0 / 3.0)
62	REAL betam
63	PARAMETER (betam = 15.0) ! pour Phim / h dans la S.L stable
64	REAL betah
65	PARAMETER (betah = 15.0)
66	REAL betas
67	PARAMETER (betas = 5.0) ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
68	REAL sffrac
69	PARAMETER (sffrac = 0.1) ! S.L. = z/h < .1
70	REAL binm
71	PARAMETER (binm = betam * sffrac)
72	REAL binh
73	PARAMETER (binh = betah * sffrac)
74
75	REAL q_star, t_star
76	REAL b1, b2, b212, b2sr ! Lambert correlations T' q' avec T* q*
77	PARAMETER (b1 = 70., b2 = 20.) ! b1 entre 70 et 100
78
79	REAL z(klon, klev)
80	! AM
81	REAL zref, dt0
82	PARAMETER (zref = 2.) ! Niveau de ref a 2m
83	PARAMETER (dt0 = 0.1) ! convergence do while
84
85	INTEGER i, k, j
86	REAL khfs(klon) ! surface kinematic heat flux [mK/s]
87	REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
88	REAL heatv(klon) ! surface virtual heat flux
89	REAL rhino(klon, klev) ! bulk Richardon no. mais en Theta_v
90	LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
91	LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
92	LOGICAL omegafl(klon) ! flag de prolongement cape pour pt Omega
93	REAL obklen(klon) ! Monin-Obukhov lengh
94
95	REAL pblh(klon) ! PBL H (m)
96	REAL therm(klon) ! exces du thermique (K)
97	REAL plcl(klon) ! Lifted Cnd Level (Pa)
98	REAL cape(klon) ! Cape
99	REAL cin(klon) ! Inhibition
100	REAL eauliq(klon) ! Eau Liqu integree
101	REAL ctei(klon) ! Cld Top Entr. Instab.
102	REAL d_qt(klon) ! Saut de qT a l'inversion
103	REAL d_thv(klon) ! Theta_e
104	REAL dlt_2(klon) ! Ordonnee a gauche de courbe de melange
105	REAL xhis(klon) ! fraction de melange pour flottab nulle
106	REAL posint(klon) ! partie positive de l'int. de Peter
107	REAL omega(klon) ! point ultime de l'ascention du thermique
108	REAL diagok(klon) ! pour traiter les sous-mailles sans info
109	! Algorithme thermique
110	REAL s(klon, klev) ! [P/Po]^Kappa milieux couches
111	REAL th_th(klon) ! potential temperature of thermal
112	REAL the_th(klon) ! equivalent potential temperature of thermal
113	REAL qt_th(klon) ! total water of thermal
114	REAL tbef(klon) ! T thermique niveau ou calcul precedent
115	LOGICAL zsat(klon) ! le thermique est sature
116	LOGICAL zcin(klon) ! calcul d'inhibition
117	REAL kape(klon) ! Cape locale
118	REAL kin(klon) ! Cin locale
119	! calcul de CTEI etc
120	REAL the1, the2, aa, bb, zthvd, zthvu, qsat, chi, rh, zxt, zdu2
121	REAL rnum, denom, th1, th2, tv1, tv2, thv1, thv2, ql1, ql2, dt
122	REAL dqsat_dt, qsat2, qt1, q1, q2, t1, t2, tl1, te2, xnull, delt_the
123	REAL delt_qt, quadsat, spblh, reduc
124	! diag REAL dTv21(klon,klev)
125
126	REAL phiminv(klon) ! inverse phi function for momentum
127	REAL phihinv(klon) ! inverse phi function for heat
128	REAL wm(klon) ! turbulent velocity scale for momentum
129	REAL zm(klon) ! current level height
130	REAL zp(klon) ! current level height + one level up
131	REAL zcor, zdelta, zcvm5
132	REAL fac, pblmin
133	REAL missing_val
134
135	! c missing_val=nf90_fill_real (avec include netcdf)
136	missing_val = 0.
137
138	! initialisations (Anne)
139	isommet = klev
140	b212 = sqrt(b1 * b2)
141	b2sr = sqrt(b2)
142
143	! Initialisation thermo
144	rlvcp = rlvtt / rcpd
145	reps = rd / rv
146	! raz
147	q_star = 0.
148	t_star = 0.
149	cape(:) = missing_val
150	kape(:) = 0.
151	cin(:) = missing_val
152	eauliq(:) = missing_val
153	ctei(:) = missing_val
154	d_qt(:) = missing_val
155	d_thv(:) = missing_val
156	dlt_2(:) = missing_val
157	xhis(:) = missing_val
158	posint(:) = missing_val
159	kin(:) = missing_val
160	omega(:) = missing_val
161	diagok(:) = 0.
162	! diag dTv21(:,:)= missing_val
163
164	! Calculer les hauteurs de chaque couche
165	DO i = 1, knon
166	z(i, 1) = rd * t(i, 1) / (0.5 * (paprs(i, 1) + pplay(i, 1))) * (paprs(i, 1) - pplay(i, 1)) / rg
167	s(i, 1) = (pplay(i, 1) / paprs(i, 1))**rkappa
168	END DO
169	! s(k) = [pplay(k)/ps]^kappa
170	! + + + + + + + + + pplay <-> s(k) t dp=pplay(k-1)-pplay(k)
171
172	! ----------------- paprs <-> sig(k)
173
174	! + + + + + + + + + pplay <-> s(k-1)
175
176
177	! + + + + + + + + + pplay <-> s(1) t dp=paprs-pplay z(1)
178
179	! ----------------- paprs <-> sig(1)
180
181	DO k = 2, klev
182	DO i = 1, knon
183	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
184	s(i, k) = (pplay(i, k) / paprs(i, 1))**rkappa
185	END DO
186	END DO
187	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
188	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
189	! +++ Determination des grandeurs de surface +++++++++++++++++++++
190	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
191	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
192	DO i = 1, knon
193	! AM Niveau de ref choisi a 2m
194	zxt = t2m(i)
195
196	! ***************************************************
197	! attention, il doit s'agir de <w'theta'>
198	! ;Calcul de tcls virtuel et de w'theta'virtuel
199	! ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
200	! tcls=tcls(1+.608qcls)
201
202	! ;Pour avoir w'theta',
203	! ; il faut diviser par ro.Cp
204	! Cp=Cpd(1+0.84qcls)
205	! fcs=fcs/(ro_surf*Cp)
206	! ;On transforme w'theta' en w'thetav'
207	! Lv=(2.501-0.00237(tcls-273.15))1.E6
208	! xle=xle/(ro_surf*Lv)
209	! fcsv=fcs+.608xletcls
210	! ***************************************************
211	! dif khfs est deja w't'_v / heatv(i) = khfs(i) + RETVzxtkqfs(i)
212	! AM calcul de Ro = paprs(i,1)/Rd zxt
213	! AM convention >0 vers le bas ds lmdz
214	khfs(i) = -flux_t(i, 1) * zxt * rd / (rcpd * paprs(i, 1))
215	kqfs(i) = -flux_q(i, 1) * zxt * rd / (paprs(i, 1))
216	! AM verifier que khfs et kqfs sont bien de la forme w'l'
217	heatv(i) = khfs(i) + retv * zxt * kqfs(i)
218	! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
219	! AM ustar est en entree (calcul dans stdlevvar avec t2m q2m)
220	! Theta et qT du thermique sans exces
221	qt_th(i) = q2m(i)
222	! Al1 Th_th restera la Theta du thermique sans exces jusqu'au 3eme calcul
223	th_th(i) = t2m(i)
224	END DO
225
226	DO i = 1, knon
227	rhino(i, 1) = 0.0 ! Global Richardson
228	check(i) = .TRUE.
229	pblh(i) = z(i, 1) ! on initialise pblh a l'altitude du 1er niveau
230	! Attention Plcl est pression ou altitude ?
231	! plcl(i) = 6000. ! m
232	plcl(i) = 200. ! hPa
233	IF (heatv(i)>0.0001) THEN
234	! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
235	obklen(i) = -t(i, 1) * ustar(i)*3 / (rg vk * heatv(i))
236	ELSE
237	! set pblh to the friction high (cf + bas)
238	pblh(i) = 700.0 * ustar(i)
239	check(i) = .FALSE.
240	END IF
241	END DO
242
243
244	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
245	! PBL height calculation:
246	! Search for level of pbl. Scan upward until the Richardson number between
247	! the first level and the current level exceeds the "critical" value.
248	! (bonne idee Nu de separer le Ric et l'exces de temp du thermique)
249	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
250	fac = 100.0
251	DO k = 2, isommet
252	DO i = 1, knon
253	IF (check(i)) THEN
254	zdu2 = u(i, k)2 + v(i, k)2
255	zdu2 = max(zdu2, 1.0E-20)
256	! Theta_v environnement
257	zthvd = t(i, k) / s(i, k) * (1. + retv * q(i, k))
258	zthvu = th_th(i) * (1. + retv * qt_th(i))
259	! Le Ri bulk par Theta_v
260	rhino(i, k) = (z(i, k) - zref) * rg * (zthvd - zthvu) / (zdu2 * 0.5 * (zthvd + zthvu))
261
262	IF (rhino(i, k)>=ricr) THEN
263	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))
264	! test04 (la pblh est encore ici sous-estime'e)
265	pblh(i) = pblh(i) + 100.
266	! pblT(i) = t(i,k-1) + (t(i,k)-t(i,k-1)) *
267	! . (pblh(i)-z(i,k-1))/(z(i,k)-z(i,k-1))
268	check(i) = .FALSE.
269	END IF
270	END IF
271	END DO
272	END DO
273
274
275	! Set pbl height to maximum value where computation exceeds number of
276	! layers allowed
277
278	DO i = 1, knon
279	IF (check(i)) pblh(i) = z(i, isommet)
280	END DO
281
282	! Improve estimate of pbl height for the unstable points.
283	! Find unstable points (sensible heat flux is upward):
284
285	DO i = 1, knon
286	IF (heatv(i)>0.) THEN
287	unstbl(i) = .TRUE.
288	check(i) = .TRUE.
289	ELSE
290	unstbl(i) = .FALSE.
291	check(i) = .FALSE.
292	END IF
293	END DO
294
295	! For the unstable case, compute velocity scale and the
296	! convective temperature excess:
297
298	DO i = 1, knon
299	IF (check(i)) THEN
300	phiminv(i) = (1. - binm * pblh(i) / obklen(i))**onet
301	! ***************************************************
302	! Wm ? et W* ? c'est la formule pour z/h < .1
303	! ;Calcul de w* ;;
304	! ;;;;;;;;;;;;;;;;
305	! w_star=((g/tcls)fcsvz(ind))^(1/3.) [ou prendre la premiere approx de h)
306	! ;; CALCUL DE wm ;;
307	! ;;;;;;;;;;;;;;;;;;
308	! ; Ici on considerera que l'on est dans la couche de surf jusqu'a 100m
309	! ; On prend svt couche de surface=0.1*h mais on ne connait pas h
310	! ;;;;;;;;;;;Dans la couche de surface
311	! if (z(ind) le 20) then begin
312	! Phim=(1.-15.*(z(ind)/L))^(-1/3.)
313	! wm=u_star/Phim
314	! ;;;;;;;;;;;En dehors de la couche de surface
315	! END IF ELSE IF (z(ind) gt 20) then begin
316	! wm=(u_star^3+c1*w_star^3)^(1/3.)
317	! END IF
318	! ***************************************************
319	wm(i) = ustar(i) * phiminv(i)
320	! ======================================================================
321	! valeurs de Dominique Lambert de la campagne SEMAPHORE :
322	! <T'^2> = 100.T^2; <q'^2> = 20.q^2 a 10m
323	! <Tv'^2> = (1+1.2q).100.T* + 1.2Tv.sqrt(20100).T.q* + (.608Tv)^220.q*^2;
324	! et dTetavS = sqrt(<Tv'^2>) ainsi calculee.
325	! avec : T=<w'T'>_s/w et q=<w'q'>/w
326	! !!! on peut donc utiliser w* pour les fluctuations <-> Lambert
327	! (leur corellation pourrait dependre de beta par ex)
328	! if fcsv(i,j) gt 0 then begin
329	! dTetavs=b1(1.+2..608q_10(i,j))(fcs(i,j)/wm(i,j))^2+$
330	! (.608Thetav_10(i,j))^2b2*(xle(i,j)/wm(i,j))^2+$
331	! 2..608thetav_10(i,j)sqrt(b1b2)(xle(i,j)/wm(i,j))(fcs(i,j)/wm(i,j))
332	! dqs=b2*(xle(i,j)/wm(i,j))^2
333	! theta_s(i,j)=thetav_10(i,j)+sqrt(dTetavs)
334	! q_s(i,j)=q_10(i,j)+sqrt(dqs)
335	! END IF else begin
336	! Theta_s(i,j)=thetav_10(i,j)
337	! q_s(i,j)=q_10(i,j)
338	! endelse
339	! leur reference est le niveau a 10m, mais on prend 2m ici.
340	! ======================================================================
341	! Premier calcul de l'exces tu thermique
342	! ======================================================================
343	! HBTM therm(i) = heatv(i)*fak/wm(i)
344	! forme Mathieu :
345	q_star = max(0., kqfs(i) / wm(i))
346	t_star = max(0., khfs(i) / wm(i))
347	! Al1 Houston, we have a problem : il arrive en effet que heatv soit
348	! positif (=thermique instable) mais pas t_star : avec evaporation
349	! importante, il se peut qu'on refroidisse la 2m Que faire alors ?
350	! Garder le seul terme en q_star^2 ? ou rendre negatif le t_star^2 ?
351	therm(i) = sqrt(b1 * (1. + 2. * retv * qt_th(i)) * t_star*2 + (retv th_th(i))*2 b2 * q_star * q_star + 2. * retv * th_th(i) * b212 * &
352	q_star * t_star)
353
354	! Theta et qT du thermique (forme H&B) avec exces
355	! (attention, on ajoute therm(i) qui est virtuelle ...)
356	! pourquoi pas sqrt(b1)*t_star ?
357	! dqs = b2sr*kqfs(i)/wm(i)
358	qt_th(i) = qt_th(i) + b2sr * q_star
359	rhino(i, 1) = 0.0
360	END IF
361	END DO
362
363	! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
364	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
365	! ++ Improve pblh estimate for unstable conditions using the +++++++
366	! ++ convective temperature excess : +++++++
367	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
368	! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
369
370	DO k = 2, isommet
371	DO i = 1, knon
372	IF (check(i)) THEN
373	! test zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
374	zdu2 = u(i, k)2 + v(i, k)2
375	zdu2 = max(zdu2, 1.0E-20)
376	! Theta_v environnement
377	zthvd = t(i, k) / s(i, k) * (1. + retv * q(i, k))
378
379	! et therm Theta_v (avec hypothese de constance de H&B,
380	! qui assimile qT a vapeur)
381	zthvu = th_th(i) * (1. + retv * qt_th(i)) + therm(i)
382
383
384	! Le Ri par Theta_v
385	! AM Niveau de ref 2m
386	rhino(i, k) = (z(i, k) - zref) * rg * (zthvd - zthvu) / (zdu2 * 0.5 * (zthvd + zthvu))
387
388	! Niveau critique atteint
389	IF (rhino(i, k)>=ricr) THEN
390	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))
391	! test04
392	pblh(i) = pblh(i) + 100.
393	! pblT(i) = t(i,k-1) + (t(i,k)-t(i,k-1)) *
394	! . (pblh(i)-z(i,k-1))/(z(i,k)-z(i,k-1))
395	check(i) = .FALSE.
396	END IF
397	END IF
398	END DO
399	END DO
400
401	! Set pbl height to maximum value where computation exceeds number of
402	! layers allowed (H&B)
403
404	DO i = 1, knon
405	IF (check(i)) pblh(i) = z(i, isommet)
406	END DO
407
408	! PBL height must be greater than some minimum mechanical mixing depth
409	! Several investigators have proposed minimum mechanical mixing depth
410	! relationships as a function of the local friction velocity, u*. We
411	! make use of a linear relationship of the form h = c u* where c=700.
412	! The scaling arguments that give rise to this relationship most often
413	! represent the coefficient c as some constant over the local coriolis
414	! parameter. Here we make use of the experimental results of Koracin
415	! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
416	! where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
417	! latitude value for f so that c = 0.07/f = 700. (H&B)
418	! Al1 calcul de pblT dans ce cas
419	DO i = 1, knon
420	pblmin = 700.0 * ustar(i)
421	IF (pblh(i)<pblmin) check(i) = .TRUE.
422	END DO
423	DO i = 1, knon
424	IF (check(i)) THEN
425	pblh(i) = 700.0 * ustar(i)
426	! et par exemple :
427	! pblT(i) = t(i,2) + (t(i,3)-t(i,2)) *
428	! . (pblh(i)-z(i,2))/(z(i,3)-z(i,2))
429	END IF
430	END DO
431
432	! ********************************************************************
433	! pblh is now available; do preparation for final calculations :
434	! ********************************************************************
435	DO i = 1, knon
436	check(i) = .TRUE.
437	zsat(i) = .FALSE.
438	zcin(i) = .FALSE.
439	! omegafl utilise pour prolongement CAPE
440	omegafl(i) = .FALSE.
441
442	! Do additional preparation for unstable cases only, set temperature
443	! and moisture perturbations depending on stability.
444	! Rq: les formules sont prises dans leur forme Couche de Surface
445	IF (unstbl(i)) THEN
446	! Al pblh a change', on recalcule :
447	zxt = (th_th(i) - zref * 0.5 * rg / rcpd / (1. + rvtmp2 * qt_th(i))) * (1. + retv * qt_th(i))
448	phiminv(i) = (1. - binm * pblh(i) / obklen(i))**onet
449	phihinv(i) = sqrt(1. - binh * pblh(i) / obklen(i))
450	wm(i) = ustar(i) * phiminv(i)
451	END IF
452	END DO
453
454
455	! =======================================================
456	! last upward integration
457	! For all unstable layers, compute integral info and CTEI
458	! =======================================================
459
460	! 1/Recompute surface characteristics with the improved pblh
461	! ----------------------------------------------------------
462	DO i = 1, knon
463	IF (unstbl(i)) THEN
464	diagok(i) = 1.
465	! from missing_value to zero
466	cape(i) = 0.
467	cin(i) = 0.
468	eauliq(i) = 0.
469	ctei(i) = 0.
470	d_qt(i) = 0.
471	d_thv(i) = 0.
472	dlt_2(i) = 0.
473	xhis(i) = 0.
474	posint(i) = 0.
475	kin(i) = 0.
476	omega(i) = 0.
477
478	phiminv(i) = (1. - binm * pblh(i) / obklen(i))**onet
479	wm(i) = ustar(i) * phiminv(i)
480	q_star = max(0., kqfs(i) / wm(i))
481	t_star = max(0., khfs(i) / wm(i))
482	therm(i) = sqrt(b1 * (1. + 2. * retv * qt_th(i)) * t_star*2 + (retv th_th(i))*2 b2 * q_star * q_star + 2. * retv * th_th(i) * b212 * &
483	q_star * t_star)
484	! Al1diag
485	! trmb1(i) = b1(1.+2.RETVqT_th(i))t_star**2
486	! trmb2(i) = (RETVTh_th(i))2b2q_starq_star
487	! trmb3(i) = 2.RETVTh_th(i)b212q_star*t_star
488
489	! Th_th will now be the thermal-theta (including exces)
490	! c Th_th(i) = Th_th(i)+sqrt(b1)*max(0.,khfs(i)/wm(i))
491	th_th(i) = th_th(i) + therm(i)
492	! al1diag
493	! trmb2(i) = wm(i)
494	! trmb3(i) = phiminv(i)
495	! and computes Theta_e for thermal
496	the_th(i) = th_th(i) + rlvcp * qt_th(i)
497	END IF ! unstbl
498	! Al1 compute a first guess of Plcl with the Bolton/Emanuel formula
499	t2 = th_th(i)
500	! thermodyn functions
501	zdelta = max(0., sign(1., rtt - t2))
502	qsat = r2es * foeew(t2, zdelta) / paprs(i, 1)
503	qsat = min(0.5, qsat)
504	zcor = 1. / (1. - retv * qsat)
505	qsat = qsat * zcor
506	! relative humidity of thermal at 2m
507	rh = qt_th(i) / qsat
508	chi = t2 / (1669.0 - 122.0 * rh - t2)
509	plcl(i) = paprs(i, 1) * (rh**chi)
510	! al1diag
511	! ctei(i) = Plcl(i)
512	! cape(i) = T2
513	! trmb1(i)= Chi
514	! select unstable columns (=thermals)
515	check(i) = .FALSE.
516	IF (heatv(i)>0.) check(i) = .TRUE.
517	! diag
518	! dTv21(i,1) = T2(1+RETVqT_th(i))-t(i,1)(1+RETVq(i,1))
519	END DO
520	! ----------------------------------
521	! 2/ upward integration for thermals
522	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ k loop
523	DO k = 2, isommet
524	DO i = 1, knon
525	IF (check(i) .OR. omegafl(i)) THEN
526	! CC if (pplay(i,k) .le. plcl(i)) THEN
527	zm(i) = z(i, k - 1)
528	zp(i) = z(i, k)
529	! Environnement : calcul de Tv1 a partir de t(:,:)== T liquide
530	! ==============
531	tl1 = t(i, k)
532	t1 = tl1
533	zdelta = max(0., sign(1., rtt - t1))
534	qsat = r2es * foeew(t1, zdelta) / pplay(i, k)
535	qsat = min(0.5, qsat)
536	zcor = 1. / (1. - retv * qsat)
537	qsat = qsat * zcor
538	q1 = min(q(i, k), qsat)
539	ql1 = max(0., q(i, k) - q1)
540	! thermodyn function (Tl2Tql)
541	dt = rlvcp * ql1
542	DO WHILE (abs(dt)>=dt0)
543	t1 = t1 + dt
544	zdelta = max(0., sign(1., rtt - t1))
545	zcvm5 = r5les * (1. - zdelta) + r5ies * zdelta
546	qsat = r2es * foeew(t1, zdelta) / pplay(i, k)
547	qsat = min(0.5, qsat)
548	zcor = 1. / (1. - retv * qsat)
549	qsat = qsat * zcor
550	dqsat_dt = foede(t1, zdelta, zcvm5, qsat, zcor)
551	! correction lineaire pour conserver Tl env
552	! << Tl = T1 + DT - RLvCp(ql1 - dqsat/dTDT >>
553	denom = 1. + rlvcp * dqsat_dt
554	q1 = min(q(i, k), qsat)
555	ql1 = q(i, k) - q1 ! can be negative
556	rnum = tl1 - t1 + rlvcp * ql1
557	dt = rnum / denom
558	END DO
559	ql1 = max(0., ql1)
560	tv1 = t1 * (1. + retv * q1 - ql1)
561	! Thermique : on atteint le seuil B/E de condensation
562	! ==============
563
564	IF (.NOT. zsat(i)) THEN
565	! first guess from The_th(i) = Th_th(i) + RLvCp* [qv=qT_th(i)]
566	t2 = s(i, k) * the_th(i) - rlvcp * qt_th(i)
567	zdelta = max(0., sign(1., rtt - t2))
568	qsat = r2es * foeew(t2, zdelta) / pplay(i, k)
569	qsat = min(0.5, qsat)
570	zcor = 1. / (1. - retv * qsat)
571	qsat = qsat * zcor
572	q2 = min(qt_th(i), qsat)
573	ql2 = max(0., qt_th(i) - q2)
574	IF (ql2>0.0001) zsat(i) = .TRUE.
575	tbef(i) = t2
576	! a PBLH non sature
577	IF (zm(i)<pblh(i) .AND. zp(i)>=pblh(i)) THEN
578	reduc = (pblh(i) - zm(i)) / (zp(i) - zm(i))
579	spblh = s(i, k - 1) + reduc * (s(i, k) - s(i, k - 1))
580	! lmdz : qT1 et Thv1
581	t1 = (t(i, k - 1) + reduc * (t(i, k) - t(i, k - 1)))
582	thv1 = t1 * (1. + retv * q(i, k)) / spblh
583	! on calcule pour le cas sans nuage un ctei en Delta Thv
584	thv2 = t2 / spblh * (1. + retv * qt_th(i))
585	ctei(i) = thv1 - thv2
586	tv2 = t2 * (1. + retv * q2 - ql2)
587	! diag
588	! dTv21(i,k) = Tv2-Tv1
589	check(i) = .FALSE.
590	omegafl(i) = .TRUE.
591	END IF
592	END IF
593
594	IF (zsat(i)) THEN
595	! thermodyn functions (Te2Tqsat)
596	t2 = tbef(i)
597	dt = 1.
598	te2 = s(i, k) * the_th(i)
599	DO WHILE (abs(dt)>=dt0)
600	zdelta = max(0., sign(1., rtt - t2))
601	zcvm5 = r5les * (1. - zdelta) + r5ies * zdelta
602	qsat = r2es * foeew(t2, zdelta) / pplay(i, k)
603	qsat = min(0.5, qsat)
604	zcor = 1. / (1. - retv * qsat)
605	qsat = qsat * zcor
606	dqsat_dt = foede(t2, zdelta, zcvm5, qsat, zcor)
607	! correction lineaire pour conserver Te_th
608	! << Te = T2 + DT + RLvCp(qsatbef + dq/dTDT >>
609	denom = 1. + rlvcp * dqsat_dt
610	rnum = te2 - t2 - rlvcp * qsat
611	dt = rnum / denom
612	t2 = t2 + dt
613	END DO
614	q2 = min(qt_th(i), qsat)
615	ql2 = max(0., qt_th(i) - q2)
616	! jusqu'a PBLH y compris
617	IF (zm(i)<pblh(i)) THEN
618
619	! mais a PBLH, interpolation et complements
620	IF (zp(i)>=pblh(i)) THEN
621	reduc = (pblh(i) - zm(i)) / (zp(i) - zm(i))
622	spblh = s(i, k - 1) + reduc * (s(i, k) - s(i, k - 1))
623	! CAPE et EauLiq a pblH
624	cape(i) = kape(i) + reduc * (zp(i) - zm(i)) * rg * .5 / (tv2 + tv1) * max(0., (tv2 - tv1))
625	eauliq(i) = eauliq(i) + reduc * (paprs(i, k - 1) - paprs(i, k)) * ql2 / rg
626	! CTEI
627	the2 = (t2 + rlvcp * q2) / spblh
628	! T1 est en realite la Tl env (on a donc strict The1)
629	t1 = (t(i, k - 1) + reduc * (t(i, k) - t(i, k - 1)))
630	the1 = (t1 + rlvcp * q(i, k)) / spblh
631	! Calcul de la Cloud Top Entrainement Instability
632	! cf Mathieu Lahellec QJRMS (2005) Comments to DYCOMS-II
633	! saut a l'inversion :
634	delt_the = the1 - the2 ! negatif
635	delt_qt = q(i, k) - qt_th(i) ! negatif
636	d_qt(i) = -delt_qt
637	dlt_2(i) = .63 * delt_the - the2 * delt_qt
638	! init ctei(i)
639	ctei(i) = dlt_2(i)
640	IF (dlt_2(i)<-0.1) THEN
641	! integrale de Peter :
642	aa = delt_the - delt_qt * (rlvcp - retv * the2)
643	bb = (rlvcp - (1. + retv) * the2) * ql2
644	d_thv(i) = aa - bb
645	! approx de Xhi_s et de l'integrale Xint=ctei(i)
646	xhis(i) = bb / (aa - dlt_2(i))
647	! trmb1(i) = xhis
648	! trmb3(i) = dlt_2
649	xnull = bb / aa
650	IF (xhis(i)>0.1) THEN
651	ctei(i) = dlt_2(i) * xhis(i) + aa * (1. - xhis(i)) + bb * alog(xhis(i))
652	ELSE
653	ctei(i) = .5 * (dlt_2(i) + aa - bb)
654	END IF
655	IF (xnull>0.) THEN
656	posint(i) = aa - bb + bb * alog(xnull)
657	ELSE
658	posint(i) = 0.
659	END IF
660	ELSE
661	ctei(i) = 1.
662	posint(i) = 1.
663	END IF
664	check(i) = .FALSE.
665	omegafl(i) = .TRUE.
666	END IF ! end a pblh
667	IF (check(i)) eauliq(i) = eauliq(i) + (paprs(i, k) - paprs(i, k + 1)) * ql2 / rg
668	END IF
669
670	END IF ! Zsat
671
672	! KAPE : thermique / environnement
673	tv2 = t2 * (1. + retv * q2 - ql2)
674	! diag
675	! dTv21(i,k) = Tv2-Tv1
676	! Kape courante
677	kape(i) = kape(i) + (zp(i) - zm(i)) * rg * .5 / (tv2 + tv1) * max(0., (tv2 - tv1))
678	! Cin
679	IF (zcin(i) .AND. tv2 - tv1>0.) THEN
680	zcin(i) = .FALSE.
681	cin(i) = kin(i)
682	END IF
683	IF (.NOT. zcin(i) .AND. tv2 - tv1<0.) THEN
684	zcin(i) = .TRUE.
685	kin(i) = kin(i) + (zp(i) - zm(i)) * rg * .5 / (tv2 + tv1) * min(0., (tv2 - tv1))
686	END IF
687	IF (kape(i) + kin(i)<0.) THEN
688	omega(i) = zm(i)
689	! trmb3(i) = paprs(i,k)
690	omegafl(i) = .FALSE.
691	! diag
692	! PRINT*,'Tv2-Tv1 (k): ',i,(dTv21(i,j),j=1,k)
693	END IF
694	! CC EndIf !plcl
695	END IF ! check(i)
696	END DO
697	END DO ! end of level loop
698	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ end k loop
699
700	END SUBROUTINE hbtm2l

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format