/[lmdze]/trunk/Sources/phylmd/hbtm.f
ViewVC logotype

Annotation of /trunk/Sources/phylmd/hbtm.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (hide annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 7 months ago) by guez
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 guez 37 module HBTM_m
2 guez 3
3 guez 37 IMPLICIT none
4 guez 3
5 guez 37 contains
6 guez 3
7 guez 206 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 guez 3
10 guez 149 ! D'apr\'es Holstag et Boville et Troen et Mahrt
11 guez 37 ! JAS 47 BLM
12 guez 3
13 guez 186 ! 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 guez 3
17 guez 186 ! 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 guez 3
24 guez 186 ! 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 guez 3
31 guez 71 USE dimphy, ONLY: klev, klon
32 guez 178 USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rtt
33 guez 49 USE yoethf_m, ONLY: r2es, rvtmp2
34     USE fcttre, ONLY: foeew
35    
36 guez 37 ! Arguments:
37 guez 3
38 guez 186 ! 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 guez 37 REAL, intent(in):: t2m(klon) ! temperature a 2 m
43     ! q a 2 et 10m
44 guez 186 REAL, intent(in):: q2m(klon)
45 guez 227 REAL, intent(in):: ustar(:) ! (knon)
46 guez 206 REAL, intent(in):: flux_t(:), flux_q(:) ! (knon) flux à la surface
47    
48 guez 227 REAL, intent(in):: u(klon, klev) ! vitesse U (m/s)
49     REAL, intent(in):: v(klon, klev) ! vitesse V (m/s)
50    
51 guez 37 ! temperature (K)
52 guez 186 REAL, intent(in):: t(klon, klev)
53 guez 37 ! vapeur d'eau (kg/kg)
54 guez 186 REAL, intent(in):: q(klon, klev)
55 guez 3
56 guez 186 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 guez 206
71     INTEGER knon ! nombre de points a calculer
72 guez 37 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 guez 178 REAL zthvd, zthvu, qqsat
136 guez 149 REAL t2
137 guez 37
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 guez 149 REAL zcor
145 guez 37
146 guez 178 REAL pblmin
147 guez 37
148     !-----------------------------------------------------------------
149    
150 guez 206 knon = size(pblh)
151    
152 guez 37 ! 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 guez 206 khfs(i) = - flux_t(i)*zxt*Rd / (RCPD*paprs(i, 1))
190     kqfs(i) = - flux_q(i)*zxt*Rd / paprs(i, 1)
191 guez 37 ! 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 guez 3 phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
368 guez 37 wm(i) = ustar(i)*phiminv(i)
369     ENDIF
370     ENDDO
371 guez 3
372 guez 37 ! Main level loop to compute the diffusivities and
373     ! counter-gradient terms:
374 guez 49 loop_level: DO k = 2, isommet
375 guez 37 ! 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 guez 49 if (check(i) .or. omegafl(i)) then
384     if (.not. Zsat(i)) then
385 guez 37 T2 = T2m(i) * s(i, k)
386     ! thermodyn functions
387 guez 103 qqsat= r2es * FOEEW(T2, RTT >= T2) / pplay(i, k)
388 guez 37 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 guez 3 endif
406 guez 37 ENDDO
407 guez 49 end DO loop_level
408 guez 37
409     END SUBROUTINE HBTM
410    
411     end module HBTM_m

Properties

Name Value
svn:executable

  ViewVC Help
Powered by ViewVC 1.1.21