GFS Operational Physics Documentation  Revision: 81451
mfpbl.f
Go to the documentation of this file.
1 
3 
38  subroutine mfpbl(im,ix,km,ntrac,delt,cnvflg,
39  & zl,zm,thvx,q1,t1,u1,v1,hpbl,kpbl,
40  & sflx,ustar,wstar,xmf,tcko,qcko,ucko,vcko)
41 !
42  use machine , only : kind_phys
43  use physcons, grav => con_g, cp => con_cp
44 !
45  implicit none
46 !
47  integer im, ix, km, ntrac
48 ! &, me
49  integer kpbl(im)
50  logical cnvflg(im)
51  real(kind=kind_phys) delt
52  real(kind=kind_phys) q1(ix,km,ntrac), t1(ix,km),
53  & u1(ix,km), v1(ix,km),
54  & thvx(im,km),
55  & zl(im,km), zm(im,km+1),
56  & hpbl(im), sflx(im), ustar(im),
57  & wstar(im), xmf(im,km),
58  & tcko(im,km),qcko(im,km,ntrac),
59  & ucko(im,km),vcko(im,km)
60 !
61 c local variables and arrays
62 !
63  integer i, j, k, n, kmpbl
64 !
65  real(kind=kind_phys) dt2, dz, ce0,
66  & h1, factor, gocp,
67  & g, c1, d1,
68  & b1, f1, bb1, bb2,
69  & alp, a1, qmin, zfmin,
70  & xmmx, rbint, tau,
71 ! & rbint, tau,
72  & tem, tem1, tem2,
73  & ptem, ptem1, ptem2,
74  & pgcon
75 !
76  real(kind=kind_phys) sigw1(im), usws3(im), xlamax(im),
77  & rbdn(im), rbup(im), delz(im)
78 !
79  real(kind=kind_phys) wu2(im,km), xlamue(im,km),
80  & thvu(im,km), zi(im,km),
81  & buo(im,km)
82 !
83  logical totflg, flg(im)
84 !
85 c physical parameters
86  parameter(g=grav)
87  parameter(gocp=g/cp)
88 ! parameter(ce0=0.37,qmin=1.e-8,alp=1.0,pgcon=0.55)
89  parameter(ce0=0.38,qmin=1.e-8,alp=1.0,pgcon=0.55)
90  parameter(a1=0.08,b1=0.5,f1=0.15,c1=0.3,d1=2.58,tau=500.)
91  parameter(zfmin=1.e-8,h1=0.33333333)
92 !
93 c-----------------------------------------------------------------------
94 !
95 !************************************************************************
96 !
97  kmpbl = km/2 + 1
98  dt2 = delt
100  totflg = .true.
101  do i=1,im
102  totflg = totflg .and. (.not. cnvflg(i))
103  enddo
104  if(totflg) return
105 !!
106  do k = 1, km
107  do i=1,im
108  if (cnvflg(i)) then
109  zi(i,k) = zm(i,k+1)
110  endif
111  enddo
112  enddo
115  do i=1,im
116  if(cnvflg(i)) then
117  k = kpbl(i) / 2
118  k = max(k, 1)
119  delz(i) = zl(i,k+1) - zl(i,k)
120  xlamax(i) = ce0 / delz(i)
121  endif
122  enddo
123  do k = 1, kmpbl
124  do i=1,im
125  if(cnvflg(i)) then
126  if(k < kpbl(i)) then
127  ptem = 1./(zi(i,k)+delz(i))
128  tem = max((hpbl(i)-zi(i,k)+delz(i)) ,delz(i))
129  ptem1 = 1./tem
130  xlamue(i,k) = ce0 * (ptem+ptem1)
131  else
132  xlamue(i,k) = xlamax(i)
133  endif
134  endif
135  enddo
136  enddo
137 c
138 c compute thermal excess
139 c
141  do i=1,im
142  if(cnvflg(i)) then
143  tem = zl(i,1)/hpbl(i)
144  usws3(i) = (ustar(i)/wstar(i))**3.
145  tem1 = usws3(i) + 0.6*tem
146  tem2 = max((1.-tem), zfmin)
147  ptem = (tem1**h1) * sqrt(tem2)
148  sigw1(i) = 1.3 * ptem * wstar(i)
149  ptem1 = alp * sflx(i) / sigw1(i)
150  thvu(i,1) = thvx(i,1) + ptem1
151  buo(i,1) = g * (thvu(i,1)/thvx(i,1)-1.)
152  endif
153  enddo
154 c
155 c compute potential temperature and buoyancy for updraft air parcel
156 c
162  do k = 2, kmpbl
163  do i=1,im
164  if(cnvflg(i)) then
165  dz = zl(i,k) - zl(i,k-1)
166  tem = xlamue(i,k-1) * dz
167  ptem = 2. + tem
168  ptem1 = (2. - tem) / ptem
169  tem1 = tem * (thvx(i,k)+thvx(i,k-1)) / ptem
170  thvu(i,k) = ptem1 * thvu(i,k-1) + tem1
171  buo(i,k) = g * (thvu(i,k)/thvx(i,k)-1.)
172  endif
173  enddo
174  enddo
175 c
176 c compute updraft velocity square(wu2)
177 c
183 ! tem = 1.-2.*f1
184 ! bb1 = 2. * b1 / tem
185 ! bb2 = 2. / tem
186 ! from soares et al. (2004,qjrms)
187 ! bb1 = 2.
188 ! bb2 = 4.
189 !
190 ! from bretherton et al. (2004, mwr)
191 ! bb1 = 4.
192 ! bb2 = 2.
193 !
194 ! from our tuning
195  bb1 = 1.8
196  bb2 = 3.5
197 !
198  do i = 1, im
199  if(cnvflg(i)) then
200 !
201 ! tem = zi(i,1)/hpbl(i)
202 ! tem1 = usws3(i) + 0.6*tem
203 ! tem2 = max((1.-tem), zfmin)
204 ! ptem = (tem1**h1) * sqrt(tem2)
205 ! ptem1 = 1.3 * ptem * wstar(i)
206 ! wu2(i,1) = d1*d1*ptem1*ptem1
207 !
208  dz = zi(i,1)
209  tem = 0.5*bb1*xlamue(i,1)*dz
210  tem1 = bb2 * buo(i,1) * dz
211  ptem1 = 1. + tem
212  wu2(i,1) = tem1 / ptem1
213 !
214  endif
215  enddo
216  do k = 2, kmpbl
217  do i = 1, im
218  if(cnvflg(i)) then
219  dz = zi(i,k) - zi(i,k-1)
220  tem = 0.25*bb1*(xlamue(i,k)+xlamue(i,k-1))*dz
221  tem1 = bb2 * buo(i,k) * dz
222  ptem = (1. - tem) * wu2(i,k-1)
223  ptem1 = 1. + tem
224  wu2(i,k) = (ptem + tem1) / ptem1
225  endif
226  enddo
227  enddo
228 c
229 c update pbl height as the height where updraft velocity vanishes
230 c
233  do i=1,im
234  flg(i) = .true.
235  if(cnvflg(i)) then
236  flg(i) = .false.
237  rbup(i) = wu2(i,1)
238  endif
239  enddo
240  do k = 2, kmpbl
241  do i = 1, im
242  if(.not.flg(i)) then
243  rbdn(i) = rbup(i)
244  rbup(i) = wu2(i,k)
245  kpbl(i) = k
246  flg(i) = rbup(i).le.0.
247  endif
248  enddo
249  enddo
250  do i = 1,im
251  if(cnvflg(i)) then
252  k = kpbl(i)
253  if(rbdn(i) <= 0.) then
254  rbint = 0.
255  elseif(rbup(i) >= 0.) then
256  rbint = 1.
257  else
258  rbint = rbdn(i)/(rbdn(i)-rbup(i))
259  endif
260  hpbl(i) = zi(i,k-1) + rbint*(zi(i,k)-zi(i,k-1))
261  endif
262  enddo
263 c
265  do i=1,im
266  if(cnvflg(i)) then
267  k = kpbl(i) / 2
268  k = max(k, 1)
269  delz(i) = zl(i,k+1) - zl(i,k)
270  xlamax(i) = ce0 / delz(i)
271  endif
272  enddo
273 c
274 c update entrainment rate
275 c
276 ! do k = 1, kmpbl
277 ! do i=1,im
278 ! if(cnvflg(i)) then
279 ! if(k < kpbl(i)) then
280 ! tem = tau * sqrt(wu2(i,k))
281 ! tem1 = 1. / tem
282 ! ptem = ce0 / zi(i,k)
283 ! xlamue(i,k) = max(tem1, ptem)
284 ! else
285 ! xlamue(i,k) = xlamax(i)
286 ! endif
287 ! endif
288 ! enddo
289 ! enddo
290 !
291  do k = 1, kmpbl
292  do i=1,im
293  if(cnvflg(i)) then
294  if(k < kpbl(i)) then
295  ptem = 1./(zi(i,k)+delz(i))
296  tem = max((hpbl(i)-zi(i,k)+delz(i)) ,delz(i))
297  ptem1 = 1./tem
298  xlamue(i,k) = ce0 * (ptem+ptem1)
299  else
300  xlamue(i,k) = xlamax(i)
301  endif
302  endif
303  enddo
304  enddo
305 c
306 c updraft mass flux as a function of sigmaw
307 c (0.3*sigmaw[square root of vertical turbulence variance])
308 c
310 ! do k = 1, kmpbl
311 ! do i=1,im
312 ! if(cnvflg(i) .and. k < kpbl(i)) then
313 ! tem = zi(i,k)/hpbl(i)
314 ! tem1 = usws3(i) + 0.6*tem
315 ! tem2 = max((1.-tem), zfmin)
316 ! ptem = (tem1**h1) * sqrt(tem2)
317 ! ptem1 = 1.3 * ptem * wstar(i)
318 ! xmf(i,k) = c1 * ptem1
319 ! endif
320 ! enddo
321 ! enddo
322 c
323 c updraft mass flux as a function of updraft velocity profile
324 c
330  do k = 1, kmpbl
331  do i = 1, im
332  if (cnvflg(i) .and. k < kpbl(i)) then
333  xmf(i,k) = a1 * sqrt(wu2(i,k))
334  dz = zl(i,k+1) - zl(i,k)
335  xmmx = dz / dt2
336  xmf(i,k) = min(xmf(i,k),xmmx)
337  endif
338  enddo
339  enddo
340 c
341 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
342 c compute updraft property
343 c
357  do k = 2, kmpbl
358  do i = 1, im
359  if (cnvflg(i) .and. k <= kpbl(i)) then
360  dz = zl(i,k) - zl(i,k-1)
361  tem = 0.5 * xlamue(i,k-1) * dz
362  factor = 1. + tem
363  ptem = tem + pgcon
364  ptem1= tem - pgcon
365 !
366  tcko(i,k) = ((1.-tem)*tcko(i,k-1)+tem*
367  & (t1(i,k)+t1(i,k-1))-gocp*dz)/factor
368  ucko(i,k) = ((1.-tem)*ucko(i,k-1)+ptem*u1(i,k)
369  & +ptem1*u1(i,k-1))/factor
370  vcko(i,k) = ((1.-tem)*vcko(i,k-1)+ptem*v1(i,k)
371  & +ptem1*v1(i,k-1))/factor
372  endif
373  enddo
374  enddo
375  do n = 1, ntrac
376  do k = 2, kmpbl
377  do i = 1, im
378  if (cnvflg(i) .and. k <= kpbl(i)) then
379  dz = zl(i,k) - zl(i,k-1)
380  tem = 0.5 * xlamue(i,k-1) * dz
381  factor = 1. + tem
382 
383  qcko(i,k,n) = ((1.-tem)*qcko(i,k-1,n)+tem*
384  & (q1(i,k,n)+q1(i,k-1,n)))/factor
385  endif
386  enddo
387  enddo
388  enddo
389 !
390  return
391  end
392 
physcons::con_g
real(kind=kind_phys), parameter con_g
gravity ( )
Definition: physcons.f:59
mfpbl
subroutine mfpbl(im, ix, km, ntrac, delt, cnvflg, zl, zm, thvx, q1, t1, u1, v1, hpbl, kpbl, sflx, ustar, wstar, xmf, tcko, qcko, ucko, vcko)
This subroutine is used for calculating the mass flux and updraft properties.
Definition: mfpbl.f:41
physcons::con_cp
real(kind=kind_phys), parameter con_cp
spec heat air at p ( )
Definition: physcons.f:80