source: LMDZ5/trunk/libf/phylmd/hbtm.F90 @ 2103

Last change on this file since 2103 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

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