source: LMDZ6/trunk/libf/phylmd/hbtm2l.f90 @ 5319

Last change on this file since 5319 was 5285, checked in by abarral, 10 days ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 25.5 KB
Line 
1
2! $Header$
3
4SUBROUTINE hbtm2l(knon, paprs, pplay, t2m, t10m, q2m, q10m, ustar, flux_t, flux_q, u, v, t, q, pblh, therm, plcl, cape, &
5    cin, eauliq, ctei, d_qt, d_thv, dlt_2, xhis, posint, omega, diagok)
6  USE yomcst_mod_h
7  USE dimphy
8  USE yoethf_mod_h
9  IMPLICIT NONE
10
11  ! ***************************************************************
12  ! *                                                             *
13  ! * HBTM2L   D'apres Holstag&Boville et Troen&Mahrt             *
14  ! *                 JAS 47              BLM                     *
15  ! * Algorithmes These Anne Mathieu                              *
16  ! * Critere d'Entrainement Peter Duynkerke (JAS 50)             *
17  ! * written by  : Anne MATHIEU & Alain LAHELLEC, 22/11/99       *
18  ! * features : implem. exces Mathieu                            *
19  ! ***************************************************************
20  ! * mods : decembre 99 passage th a niveau plus bas. voir fixer *
21  ! ***************************************************************
22  ! * fin therm a la HBTM passage a forme Mathieu 12/09/2001      *
23  ! ***************************************************************
24  ! AM Fev 2003 Adaptation a LMDZ version couplee                 *
25  ! *
26  ! Pour le moment on fait passer en argument les grdeurs de surface :
27  ! flux, t,q2m, t,q10m, on va utiliser systematiquement les grdeurs a 2m ms
28  ! on garde la possibilite de changer si besoin est (jusqu'a present la
29  ! forme de HB avec le 1er niveau modele etait conservee)       *
30  ! ***************************************************************
31  ! * re-ecriture complete Alain Mars 2012 dans LMDZ5V5           *
32  ! ***************************************************************
33  REAL rlvcp, reps
34  ! Arguments:
35
36  INTEGER knon ! nombre de points a calculer
37  ! AM
38  REAL t2m(klon), t10m(klon) ! temperature a 2 et 10m
39  REAL q2m(klon), q10m(klon) ! q a 2 et 10m
40  REAL ustar(klon)
41  REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
42  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
43  REAL flux_t(klon, klev), flux_q(klon, klev) ! Flux
44  REAL u(klon, klev) ! vitesse U (m/s)
45  REAL v(klon, klev) ! vitesse V (m/s)
46  REAL t(klon, klev) ! temperature (K)
47  REAL q(klon, klev) ! vapeur d'eau (kg/kg)
48
49  INTEGER isommet
50  REAL vk
51  PARAMETER (vk=0.35) ! Von Karman => passer a .41 ! cf U.Olgstrom
52  REAL ricr
53  PARAMETER (ricr=0.4)
54  REAL fak
55  PARAMETER (fak=8.5) ! b calcul du Prandtl et de dTetas
56  REAL fakn
57  PARAMETER (fakn=7.2) ! a
58  REAL onet
59  PARAMETER (onet=1.0/3.0)
60  REAL betam
61  PARAMETER (betam=15.0) ! pour Phim / h dans la S.L stable
62  REAL betah
63  PARAMETER (betah=15.0)
64  REAL betas
65  PARAMETER (betas=5.0) ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
66  REAL sffrac
67  PARAMETER (sffrac=0.1) ! S.L. = z/h < .1
68  REAL binm
69  PARAMETER (binm=betam*sffrac)
70  REAL binh
71  PARAMETER (binh=betah*sffrac)
72
73  REAL q_star, t_star
74  REAL b1, b2, b212, b2sr ! Lambert correlations T' q' avec T* q*
75  PARAMETER (b1=70., b2=20.) ! b1 entre 70 et 100
76
77  REAL z(klon, klev)
78  ! AM
79  REAL zref, dt0
80  PARAMETER (zref=2.) ! Niveau de ref a 2m
81  PARAMETER (dt0=0.1) ! convergence do while
82
83  INTEGER i, k, j
84  REAL khfs(klon) ! surface kinematic heat flux [mK/s]
85  REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
86  REAL heatv(klon) ! surface virtual heat flux
87  REAL rhino(klon, klev) ! bulk Richardon no. mais en Theta_v
88  LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
89  LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
90  LOGICAL omegafl(klon) ! flag de prolongement cape pour pt Omega
91  REAL obklen(klon) ! Monin-Obukhov lengh
92
93  REAL pblh(klon) ! PBL H               (m)
94  REAL therm(klon) ! exces du thermique  (K)
95  REAL plcl(klon) ! Lifted Cnd Level (Pa)
96  REAL cape(klon) ! Cape
97  REAL cin(klon) ! Inhibition
98  REAL eauliq(klon) ! Eau Liqu integree
99  REAL ctei(klon) ! Cld Top Entr. Instab.
100  REAL d_qt(klon) ! Saut de qT a l'inversion
101  REAL d_thv(klon) !         Theta_e
102  REAL dlt_2(klon) ! Ordonnee a gauche de courbe de melange
103  REAL xhis(klon) ! fraction de melange pour flottab nulle
104  REAL posint(klon) ! partie positive de l'int. de Peter
105  REAL omega(klon) ! point ultime de l'ascention du thermique
106  REAL diagok(klon) ! pour traiter les sous-mailles sans info
107  ! Algorithme thermique
108  REAL s(klon, klev) ! [P/Po]^Kappa milieux couches
109  REAL th_th(klon) ! potential temperature of thermal
110  REAL the_th(klon) ! equivalent potential temperature of thermal
111  REAL qt_th(klon) ! total water of thermal
112  REAL tbef(klon) ! T thermique niveau ou calcul precedent
113  LOGICAL zsat(klon) ! le thermique est sature
114  LOGICAL zcin(klon) ! calcul d'inhibition
115  REAL kape(klon) ! Cape locale
116  REAL kin(klon) ! Cin locale
117  ! calcul de CTEI etc
118  REAL the1, the2, aa, bb, zthvd, zthvu, qsat, chi, rh, zxt, zdu2
119  REAL rnum, denom, th1, th2, tv1, tv2, thv1, thv2, ql1, ql2, dt
120  REAL dqsat_dt, qsat2, qt1, q1, q2, t1, t2, tl1, te2, xnull, delt_the
121  REAL delt_qt, quadsat, spblh, reduc
122  ! diag      REAL dTv21(klon,klev)
123
124  REAL phiminv(klon) ! inverse phi function for momentum
125  REAL phihinv(klon) ! inverse phi function for heat
126  REAL wm(klon) ! turbulent velocity scale for momentum
127  REAL zm(klon) ! current level height
128  REAL zp(klon) ! current level height + one level up
129  REAL zcor, zdelta, zcvm5
130  REAL fac, pblmin
131  REAL missing_val
132
133  include "FCTTRE.h"
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+.608*qcls)
201
202    ! ;Pour avoir w'theta',
203    ! ; il faut diviser par ro.Cp
204    ! Cp=Cpd*(1+0.84*qcls)
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+.608*xle*tcls
210    ! ***************************************************
211    ! dif khfs est deja w't'_v / heatv(i) = khfs(i) + RETV*zxt*kqfs(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)*fcsv*z(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      ! endif else if (z(ind) gt 20) then begin
316      ! wm=(u_star^3+c1*w_star^3)^(1/3.)
317      ! endif
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(20*100).T*.q* + (.608*Tv)^2*20.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.*.608*q_10(i,j))*(fcs(i,j)/wm(i,j))^2+$
330      ! (.608*Thetav_10(i,j))^2*b2*(xle(i,j)/wm(i,j))^2+$
331      ! 2.*.608*thetav_10(i,j)*sqrt(b1*b2)*(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      ! endif 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+fac*ustar(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.*RETV*qT_th(i))*t_star**2
486      ! trmb2(i) = (RETV*Th_th(i))**2*b2*q_star*q_star
487      ! trmb3(i) = 2.*RETV*Th_th(i)*b212*q_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+RETV*qT_th(i))-t(i,1)*(1+RETV*q(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/dT*DT >>
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/dT*DT >>
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  RETURN
700END SUBROUTINE hbtm2l
Note: See TracBrowser for help on using the repository browser.