source: LMDZ5/branches/testing/libf/phylmd/hbtm2l.F90 @ 2219

Last change on this file since 2219 was 2187, checked in by Laurent Fairhead, 10 years ago

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