/[lmdze]/trunk/phylmd/Interface_surf/hbtm.f90
ViewVC logotype

Contents of /trunk/phylmd/Interface_surf/hbtm.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (show annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 6 months ago) by guez
Original Path: trunk/Sources/phylmd/hbtm.f
File size: 13477 byte(s)
Rename phisinit to phis in restart.nc: clearer, same name as Fortran variable.

In aaam_bud, use rlat and rlon from phyetat0_m instead of having these
module variables associated to actual arguments in physiq.

In clmain, too many wind variables make the procedure hard to
understand. Use yu(:knon, 1) and yv(:knon, 1) instead of u1lay(:knon)
and v1lay(:knon). Note that when yu(:knon, 1) and yv(:knon, 1) are
used as actual arguments, they are probably copied to new arrays since
the elements are not contiguous. Rename yu10m to wind10m because this
is the norm of wind vector, not its zonal component. Rename yustar to
ustar. Rename uzon and vmer to u1 and v1 since these are wind
components at first layer and u1 and v1 are the names of corresponding
dummy arguments in stdlevvar.

In clmain, rename yzlev to zlev.

In clmain, screenc, stdlevvar and coefcdrag, remove the code
corresponding to zxli true (not used in LMDZ either).

Subroutine ustarhb becomes a function. Simplifications using the fact
that zx_alf2 = 0 and zx_alf1 = 1 (discarding the possibility to change
this).

In procedure vdif_kcay, remove unused dummy argument plev. Remove
useless computations of sss and sssq.

In clouds_gno, exp(100.) would overflow in single precision. Set
maximum to exp(80.) instead.

In physiq, use u(:, 1) and v(:, 1) as arguments to phytrac instead of
creating ad hoc variables yu1 and yv1.

In stdlevvar, rename dummy argument u_10m to wind10m, following the
corresponding modification in clmain. Simplifications using the fact
that ok_pred = 0 and ok_corr = 1 (discarding the possibility to change
this).

1 module HBTM_m
2
3 IMPLICIT none
4
5 contains
6
7 SUBROUTINE HBTM(paprs, pplay, t2m, q2m, ustar, flux_t, flux_q, u, v, t, q, &
8 pblh, cape, EauLiq, ctei, pblT, therm, trmb1, trmb2, trmb3, plcl)
9
10 ! D'apr\'es Holstag et Boville et Troen et Mahrt
11 ! JAS 47 BLM
12
13 ! Algorithme th\'ese Anne Mathieu. Crit\'ere d'entra\^inement
14 ! Peter Duynkerke (JAS 50). Written by: Anne MATHIEU and Alain
15 ! LAHELLEC, 22nd November 1999.
16
17 ! Modifications : d\'ecembre 99 passage th \`a niveau plus bas. Voir fixer
18 ! la prise du th \`a z/Lambda = -.2 (max Ray)
19 ! Autre algorithme : entra\^inement ~ Theta + v =constante
20 ! mais comment ? The ?
21 ! On peut fixer q \`a 0.7 qsat (cf. non adiabatique) d'où T2 et The2.
22 ! Voir aussi KE pblh = niveau The_e ou l = env.
23
24 ! Adaptation \`a LMDZ version coupl\'ee. Pour le moment on fait
25 ! passer en argument les grandeurs de surface : flux, t, q2m. On
26 ! va utiliser syst\'ematiquement les grandeurs \`a 2 m mais on
27 ! garde la possibilit\'e de changer si besoin (jusqu'\`a pr\'esent
28 ! la forme de HB avec le premier niveau mod\`ele \'etait
29 ! conserv\'ee).
30
31 USE dimphy, ONLY: klev, klon
32 USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rtt
33 USE yoethf_m, ONLY: r2es, rvtmp2
34 USE fcttre, ONLY: foeew
35
36 ! Arguments:
37
38 ! pression a inter-couche (Pa)
39 REAL, intent(in):: paprs(klon, klev+1)
40 ! pression au milieu de couche (Pa)
41 REAL, intent(in):: pplay(klon, klev)
42 REAL, intent(in):: t2m(klon) ! temperature a 2 m
43 ! q a 2 et 10m
44 REAL, intent(in):: q2m(klon)
45 REAL, intent(in):: ustar(:) ! (knon)
46 REAL, intent(in):: flux_t(:), flux_q(:) ! (knon) flux à la surface
47
48 REAL, intent(in):: u(klon, klev) ! vitesse U (m/s)
49 REAL, intent(in):: v(klon, klev) ! vitesse V (m/s)
50
51 ! temperature (K)
52 REAL, intent(in):: t(klon, klev)
53 ! vapeur d'eau (kg/kg)
54 REAL, intent(in):: q(klon, klev)
55
56 REAL, intent(out):: pblh(:) ! (knon)
57 ! Cape du thermique
58 REAL Cape(klon)
59 ! Eau liqu integr du thermique
60 REAL EauLiq(klon)
61 ! Critere d'instab d'entrainmt des nuages de
62 REAL ctei(klon)
63 REAL pblT(klon)
64 ! thermal virtual temperature excess
65 REAL therm(klon)
66 REAL trmb1(klon), trmb2(klon), trmb3(klon)
67 REAL plcl(klon)
68
69 ! Local:
70
71 INTEGER knon ! nombre de points a calculer
72 INTEGER isommet
73 ! limite max sommet pbl
74 PARAMETER (isommet=klev)
75 REAL vk
76 ! Von Karman => passer a .41 ! cf U.Olgstrom
77 PARAMETER (vk=0.35)
78 REAL ricr
79 PARAMETER (ricr=0.4)
80 ! a
81 REAL onet
82 PARAMETER (onet=1.0/3.0)
83 REAL zkmin
84 PARAMETER (zkmin=0.01)
85 REAL betam
86 ! pour Phim / h dans la S.L stable
87 PARAMETER (betam=15.0)
88 ! z/OBL<>1
89 REAL sffrac
90 ! S.L. = z/h < .1
91 PARAMETER (sffrac=0.1)
92 REAL binm
93 PARAMETER (binm=betam*sffrac)
94
95 REAL q_star, t_star
96 ! Lambert correlations T' q' avec T* q*
97 REAL b1, b2, b212, b2sr
98 PARAMETER (b1=70., b2=20.)
99
100 REAL z(klon, klev)
101
102 REAL zref
103 ! Niveau de ref a 2m peut eventuellement
104 PARAMETER (zref=2.)
105 ! etre choisi a 10m
106
107 INTEGER i, k
108 REAL zxt
109 ! surface kinematic heat flux [mK/s]
110 REAL khfs(klon)
111 ! sfc kinematic constituent flux [m/s]
112 REAL kqfs(klon)
113 ! surface virtual heat flux
114 REAL heatv(klon)
115 ! bulk Richardon no. mais en Theta_v
116 REAL rhino(klon, klev)
117 ! pts w/unstbl pbl (positive virtual ht flx)
118 LOGICAL unstbl(klon)
119 LOGICAL check(klon) ! Richardson number > critical
120 ! flag de prolongerment cape pour pt Omega
121 LOGICAL omegafl(klon)
122
123 ! Monin-Obukhov lengh
124 REAL obklen(klon)
125
126 REAL zdu2
127 ! Algorithme thermique
128 REAL s(klon, klev) ! [P/Po]^Kappa milieux couches
129 ! total water of thermal
130 REAL qT_th(klon)
131 ! T thermique niveau precedent
132 REAL qsatbef(klon)
133 ! le thermique est sature
134 LOGICAL Zsat(klon)
135 REAL zthvd, zthvu, qqsat
136 REAL t2
137
138 ! inverse phi function for momentum
139 REAL phiminv(klon)
140 ! turbulent velocity scale for momentum
141 REAL wm(klon)
142 ! current level height + one level up
143 REAL zp(klon)
144 REAL zcor
145
146 REAL pblmin
147
148 !-----------------------------------------------------------------
149
150 knon = size(pblh)
151
152 ! initialisations
153 q_star = 0
154 t_star = 0
155
156 b212=sqrt(b1*b2)
157 b2sr=sqrt(b2)
158
159 ! Calculer les hauteurs de chaque couche
160 ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
161 ! pourquoi ne pas utiliser Phi/RG ?
162 DO i = 1, knon
163 z(i, 1) = RD * t(i, 1) / (0.5*(paprs(i, 1)+pplay(i, 1))) &
164 * (paprs(i, 1)-pplay(i, 1)) / RG
165 s(i, 1) = (pplay(i, 1)/paprs(i, 1))**RKappa
166 ENDDO
167 ! s(k) = [pplay(k)/ps]^kappa
168 ! + + + + + + + + + pplay <-> s(k) t dp=pplay(k-1)-pplay(k)
169 ! ----------------- paprs <-> sig(k)
170 ! + + + + + + + + + pplay <-> s(k-1)
171 ! + + + + + + + + + pplay <-> s(1) t dp=paprs-pplay z(1)
172 ! ----------------- paprs <-> sig(1)
173
174 DO k = 2, klev
175 DO i = 1, knon
176 z(i, k) = z(i, k-1) &
177 + RD * 0.5*(t(i, k-1)+t(i, k)) / paprs(i, k) &
178 * (pplay(i, k-1)-pplay(i, k)) / RG
179 s(i, k) = (pplay(i, k) / paprs(i, 1))**RKappa
180 ENDDO
181 ENDDO
182
183 ! Determination des grandeurs de surface
184 DO i = 1, knon
185 ! Niveau de ref choisi a 2m
186 zxt = t2m(i)
187
188 ! convention >0 vers le bas ds lmdz
189 khfs(i) = - flux_t(i)*zxt*Rd / (RCPD*paprs(i, 1))
190 kqfs(i) = - flux_q(i)*zxt*Rd / paprs(i, 1)
191 ! verifier que khfs et kqfs sont bien de la forme w'l'
192 heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
193 ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
194 ! Theta et qT du thermique sans exces (interpolin vers surf)
195 ! chgt de niveau du thermique (jeudi 30/12/1999)
196 ! (interpolation lineaire avant integration phi_h)
197 qT_th(i) = q2m(i)
198 ENDDO
199
200 DO i = 1, knon
201 ! Global Richardson
202 rhino(i, 1) = 0.0
203 check(i) = .TRUE.
204 ! on initialise pblh a l'altitude du 1er niv
205 pblh(i) = z(i, 1)
206 plcl(i) = 6000.
207 ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
208 obklen(i) = -t(i, 1)*ustar(i)**3/(RG*vk*heatv(i))
209 trmb1(i) = 0.
210 trmb2(i) = 0.
211 trmb3(i) = 0.
212 ENDDO
213
214 ! PBL height calculation: Search for level of pbl. Scan upward
215 ! until the Richardson number between the first level and the
216 ! current level exceeds the "critical" value. (bonne idee Nu de
217 ! separer le Ric et l'exces de temp du thermique)
218 DO k = 2, isommet
219 DO i = 1, knon
220 IF (check(i)) THEN
221 ! pourquoi / niveau 1 (au lieu du sol) et le terme en u*^2 ?
222 zdu2 = u(i, k)**2+v(i, k)**2
223 zdu2 = max(zdu2, 1.0e-20)
224 ! Theta_v environnement
225 zthvd=t(i, k)/s(i, k)*(1.+RETV*q(i, k))
226
227 ! therm Theta_v sans exces (avec hypothese fausse de H&B, sinon,
228 ! passer par Theta_e et virpot)
229 zthvu = T2m(i)*(1.+RETV*qT_th(i))
230 ! Le Ri par Theta_v
231 ! On a nveau de ref a 2m ???
232 rhino(i, k) = (z(i, k)-zref)*RG*(zthvd-zthvu) &
233 /(zdu2*0.5*(zthvd+zthvu))
234
235 IF (rhino(i, k).GE.ricr) THEN
236 pblh(i) = z(i, k-1) + (z(i, k-1)-z(i, k)) * &
237 (ricr-rhino(i, k-1))/(rhino(i, k-1)-rhino(i, k))
238 ! test04
239 pblh(i) = pblh(i) + 100.
240 pblT(i) = t(i, k-1) + (t(i, k)-t(i, k-1)) * &
241 (pblh(i)-z(i, k-1))/(z(i, k)-z(i, k-1))
242 check(i) = .FALSE.
243 ENDIF
244 ENDIF
245 ENDDO
246 ENDDO
247
248 ! Set pbl height to maximum value where computation exceeds number of
249 ! layers allowed
250 DO i = 1, knon
251 if (check(i)) pblh(i) = z(i, isommet)
252 ENDDO
253
254 ! Improve estimate of pbl height for the unstable points.
255 ! Find unstable points (sensible heat flux is upward):
256 DO i = 1, knon
257 IF (heatv(i) > 0.) THEN
258 unstbl(i) = .TRUE.
259 check(i) = .TRUE.
260 ELSE
261 unstbl(i) = .FALSE.
262 check(i) = .FALSE.
263 ENDIF
264 ENDDO
265
266 ! For the unstable case, compute velocity scale and the
267 ! convective temperature excess:
268 DO i = 1, knon
269 IF (check(i)) THEN
270 phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
271
272 ! CALCUL DE wm
273 ! Ici on considerera que l'on est dans la couche de surf jusqu'a 100
274 ! On prend svt couche de surface=0.1*h mais on ne connait pas h
275 ! Dans la couche de surface
276 wm(i)= ustar(i)*phiminv(i)
277
278 ! forme Mathieu :
279 q_star = kqfs(i)/wm(i)
280 t_star = khfs(i)/wm(i)
281
282 therm(i) = sqrt( b1*(1.+2.*RETV*qT_th(i))*t_star**2 &
283 + (RETV*T2m(i))**2*b2*q_star**2 &
284 + max(0., 2.*RETV*T2m(i)*b212*q_star*t_star))
285
286 ! Theta et qT du thermique (forme H&B) avec exces
287 ! (attention, on ajoute therm(i) qui est virtuelle ...)
288 ! pourquoi pas sqrt(b1)*t_star ?
289 qT_th(i) = qT_th(i) + b2sr*q_star
290 ! new on differre le calcul de Theta_e
291 rhino(i, 1) = 0.
292 ENDIF
293 ENDDO
294
295 ! Improve pblh estimate for unstable conditions using the
296 ! convective temperature excess :
297 DO k = 2, isommet
298 DO i = 1, knon
299 IF (check(i)) THEN
300 zdu2 = u(i, k)**2 + v(i, k)**2
301 zdu2 = max(zdu2, 1e-20)
302 ! Theta_v environnement
303 zthvd=t(i, k)/s(i, k)*(1.+RETV*q(i, k))
304
305 ! et therm Theta_v (avec hypothese de constance de H&B,
306 zthvu = T2m(i)*(1.+RETV*qT_th(i)) + therm(i)
307
308 ! Le Ri par Theta_v
309 ! Niveau de ref 2m
310 rhino(i, k) = (z(i, k)-zref)*RG*(zthvd-zthvu) &
311 /(zdu2*0.5*(zthvd+zthvu))
312
313 IF (rhino(i, k).GE.ricr) THEN
314 pblh(i) = z(i, k-1) + (z(i, k-1)-z(i, k)) * &
315 (ricr-rhino(i, k-1))/(rhino(i, k-1)-rhino(i, k))
316 ! test04
317 pblh(i) = pblh(i) + 100.
318 pblT(i) = t(i, k-1) + (t(i, k)-t(i, k-1)) * &
319 (pblh(i)-z(i, k-1))/(z(i, k)-z(i, k-1))
320 check(i) = .FALSE.
321 ENDIF
322 ENDIF
323 ENDDO
324 ENDDO
325
326 ! Set pbl height to maximum value where computation exceeds number of
327 ! layers allowed
328 DO i = 1, knon
329 if (check(i)) pblh(i) = z(i, isommet)
330 ENDDO
331
332 ! PBL height must be greater than some minimum mechanical mixing depth
333 ! Several investigators have proposed minimum mechanical mixing depth
334 ! relationships as a function of the local friction velocity, u*. We
335 ! make use of a linear relationship of the form h = c u* where c=700.
336 ! The scaling arguments that give rise to this relationship most often
337 ! represent the coefficient c as some constant over the local coriolis
338 ! parameter. Here we make use of the experimental results of Koracin
339 ! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
340 ! where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
341 ! latitude value for f so that c = 0.07/f = 700.
342 DO i = 1, knon
343 pblmin = 700. * ustar(i)
344 pblh(i) = MAX(pblh(i), pblmin)
345 ! par exemple :
346 pblT(i) = t(i, 2) + (t(i, 3)-t(i, 2)) * &
347 (pblh(i)-z(i, 2))/(z(i, 3)-z(i, 2))
348 ENDDO
349
350 ! pblh is now available; do preparation for diffusivity calculation:
351 DO i = 1, knon
352 check(i) = .TRUE.
353 Zsat(i) = .FALSE.
354 ! omegafl utilise pour prolongement CAPE
355 omegafl(i) = .FALSE.
356 Cape(i) = 0.
357 EauLiq(i) = 0.
358 CTEI(i) = 0.
359
360 ! Do additional preparation for unstable cases only, set temperature
361 ! and moisture perturbations depending on stability.
362 ! Remarque : les formule sont prises dans leur forme CS
363 IF (unstbl(i)) THEN
364 ! Niveau de ref du thermique
365 zxt=(T2m(i)-zref*0.5*RG/RCPD/(1.+RVTMP2*qT_th(i))) &
366 *(1.+RETV*qT_th(i))
367 phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
368 wm(i) = ustar(i)*phiminv(i)
369 ENDIF
370 ENDDO
371
372 ! Main level loop to compute the diffusivities and
373 ! counter-gradient terms:
374 loop_level: DO k = 2, isommet
375 ! Find levels within boundary layer:
376 DO i = 1, knon
377 zp(i) = z(i, k)
378 IF (zkmin == 0. .AND. zp(i) > pblh(i)) zp(i) = pblh(i)
379 ENDDO
380
381 ! For all layers, compute integral info and CTEI
382 DO i = 1, knon
383 if (check(i) .or. omegafl(i)) then
384 if (.not. Zsat(i)) then
385 T2 = T2m(i) * s(i, k)
386 ! thermodyn functions
387 qqsat= r2es * FOEEW(T2, RTT >= T2) / pplay(i, k)
388 qqsat=MIN(0.5, qqsat)
389 zcor=1./(1.-retv*qqsat)
390 qqsat=qqsat*zcor
391
392 if (qqsat < qT_th(i)) then
393 ! on calcule lcl
394 if (k == 2) then
395 plcl(i) = z(i, k)
396 else
397 plcl(i) = z(i, k-1) + (z(i, k-1)-z(i, k)) &
398 * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)
399 endif
400 Zsat(i) = .true.
401 endif
402 endif
403 qsatbef(i) = qqsat
404 ! cette ligne a deja ete faite normalement ?
405 endif
406 ENDDO
407 end DO loop_level
408
409 END SUBROUTINE HBTM
410
411 end module HBTM_m

Properties

Name Value
svn:executable

  ViewVC Help
Powered by ViewVC 1.1.21