CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
mfscu.f
1
4 module mfscu_mod
5 contains
11 subroutine mfscu(im,ix,km,kmscu,ntcw,ntrac1,delt, &
12 & cnvflg,zl,zm,q1,t1,u1,v1,plyr,pix, &
13 & thlx,thvx,thlvx,gdx,thetae,radj, &
14 & krad,mrad,radmin,buo,xmfd, &
15 & tcdo,qcdo,ucdo,vcdo,xlamde)
16!
17 use machine , only : kind_phys
18 use funcphys , only : fpvs
19 use physcons, grav => con_g, cp => con_cp &
20 &, rv => con_rv, hvap => con_hvap &
21 &, fv => con_fvirt &
22 &, eps => con_eps, epsm1 => con_epsm1
23!
24 implicit none
25!
26 integer im, ix, km, kmscu, ntcw, ntrac1
27! &, me
28 integer krad(im), mrad(im)
29!
30 logical cnvflg(im)
31 real(kind=kind_phys) delt
32 real(kind=kind_phys) q1(ix,km,ntrac1),t1(ix,km), &
33 & u1(ix,km), v1(ix,km), &
34 & plyr(im,km), pix(im,km), &
35 & thlx(im,km), &
36 & thvx(im,km), thlvx(im,km), &
37 & gdx(im), radj(im), &
38 & zl(im,km), zm(im,km), &
39 & thetae(im,km), radmin(im), &
40 & buo(im,km), xmfd(im,km), &
41 & tcdo(im,km), qcdo(im,km,ntrac1), &
42 & ucdo(im,km), vcdo(im,km), &
43 & xlamde(im,km-1)
44!
45! local variables and arrays
46!
47!
48 integer i,j,indx, k, n, kk, ndc
49 integer krad1(im), mradx(im), mrady(im)
50!
51 real(kind=kind_phys) dt2, dz, ce0, cm,
52 & gocp, factor, g, tau,
53 & b1, f1, bb1, bb2,
54 & a1, a2, a11, a22,
55 & cteit, pgcon,
56 & qmin, qlmin,
57 & xmmx, tem, tem1, tem2,
58 & ptem, ptem1, ptem2
59!
60 real(kind=kind_phys) elocp, el2orc, qs, es,
61 & tld, gamma, qld, thdn,
62 & thvd, dq
63!
64 real(kind=kind_phys) wd2(im,km), thld(im,km),
65 & qtx(im,km), qtd(im,km),
66 & thlvd(im), hrad(im),
67 & xlamdem(im,km-1), ra1(im), ra2(im)
68!
69 real(kind=kind_phys) xlamavg(im), sigma(im),
70 & scaldfunc(im), sumx(im)
71!
72 logical totflg, flg(im)
73!
74 real(kind=kind_phys) actei, cldtime
75!
76c physical parameters
77 parameter(g=grav)
78 parameter(gocp=g/cp)
79 parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
80 parameter(ce0=0.4,cm=1.0,pgcon=0.55)
81 parameter(qmin=1.e-8,qlmin=1.e-12)
82 parameter(b1=0.45,f1=0.15)
83 parameter(a1=0.12,a2=0.5)
84 parameter(a11=0.2,a22=1.0)
85 parameter(cldtime=500.)
86 parameter(actei = 0.7)
87! parameter(actei = 0.23)
88!
89!************************************************************************
90!!
91 totflg = .true.
92 do i=1,im
93 totflg = totflg .and. (.not. cnvflg(i))
94 enddo
95 if(totflg) return
96!
97 dt2 = delt
98!!
99 do k = 1, km
100 do i=1,im
101 if(cnvflg(i)) then
102 buo(i,k) = 0.
103 wd2(i,k) = 0.
104 qtx(i,k) = q1(i,k,1) + q1(i,k,ntcw)
105 endif
106 enddo
107 enddo
108!
109 do i = 1, im
110 if(cnvflg(i)) then
111 hrad(i) = zm(i,krad(i))
112 krad1(i) = krad(i)-1
113 endif
114 enddo
115!
116 do i = 1, im
117 if(cnvflg(i)) then
118 k = krad(i)
119 tem = zm(i,k+1)-zm(i,k)
120 tem1 = cldtime*radmin(i)/tem
121 tem1 = max(tem1, -3.0)
122 thld(i,k)= thlx(i,k) + tem1
123 qtd(i,k) = qtx(i,k)
124 thlvd(i) = thlvx(i,k) + tem1
125 buo(i,k) = - g * tem1 / thvx(i,k)
126 endif
127 enddo
128!
130!
131 do i=1,im
132 if(cnvflg(i)) then
133 ra1(i) = a1
134 ra2(i) = a11
135 endif
136 enddo
137!
140!
141 do i = 1, im
142 if(cnvflg(i)) then
143 k = krad(i)
144 tem = thetae(i,k) - thetae(i,k+1)
145 tem1 = qtx(i,k) - qtx(i,k+1)
146 if (tem > 0. .and. tem1 > 0.) then
147 cteit= cp*tem/(hvap*tem1)
148 if(cteit > actei) then
149 ra1(i) = a2
150 ra2(i) = a22
151 endif
152 endif
153 endif
154 enddo
155!
157!
158 do i = 1, im
159 if(cnvflg(i)) then
160 radj(i) = -ra2(i) * radmin(i)
161 endif
162 enddo
163!
165!
166 do i = 1, im
167 flg(i) = cnvflg(i)
168 mrad(i) = krad(i)
169 enddo
170 do k = kmscu,1,-1
171 do i = 1, im
172 if(flg(i) .and. k < krad(i)) then
173 if(thlvd(i) <= thlvx(i,k)) then
174 mrad(i) = k
175 else
176 flg(i)=.false.
177 endif
178 endif
179 enddo
180 enddo
181 do i=1,im
182 if (cnvflg(i)) then
183 kk = krad(i)-mrad(i)
184 if(kk < 1) cnvflg(i)=.false.
185 endif
186 enddo
187!!
188 totflg = .true.
189 do i=1,im
190 totflg = totflg .and. (.not. cnvflg(i))
191 enddo
192 if(totflg) return
193!!
194!
196!
197 do k = 1, kmscu
198 do i=1,im
199 if(cnvflg(i)) then
200 dz = zl(i,k+1) - zl(i,k)
201 if(k >= mrad(i) .and. k < krad(i)) then
202 if(mrad(i) == 1) then
203 ptem = 1./(zm(i,k)+dz)
204 else
205 ptem = 1./(zm(i,k)-zm(i,mrad(i)-1)+dz)
206 endif
207 tem = max((hrad(i)-zm(i,k)+dz) ,dz)
208 ptem1 = 1./tem
209 xlamde(i,k) = ce0 * (ptem+ptem1)
210 else
211 xlamde(i,k) = ce0 / dz
212 endif
213 xlamdem(i,k) = cm * xlamde(i,k)
214 endif
215 enddo
216 enddo
217!
219!
220 do k = kmscu,1,-1
221 do i=1,im
222 if(cnvflg(i) .and. k < krad(i)) then
223 dz = zl(i,k+1) - zl(i,k)
224 tem = 0.5 * xlamde(i,k) * dz
225 factor = 1. + tem
226!
227 thld(i,k) = ((1.-tem)*thld(i,k+1)+tem*
228 & (thlx(i,k)+thlx(i,k+1)))/factor
229 qtd(i,k) = ((1.-tem)*qtd(i,k+1)+tem*
230 & (qtx(i,k)+qtx(i,k+1)))/factor
231!
232 tld = thld(i,k) / pix(i,k)
233 es = 0.01 * fpvs(tld) ! fpvs in pa
234 qs = max(qmin, eps * es / (plyr(i,k)+epsm1*es))
235 dq = qtd(i,k) - qs
236!
237 if (dq > 0.) then
238 gamma = el2orc * qs / (tld**2)
239 qld = dq / (1. + gamma)
240 qtd(i,k) = qs + qld
241 tem1 = 1. + fv * qs - qld
242 thdn = thld(i,k) + pix(i,k) * elocp * qld
243 thvd = thdn * tem1
244 else
245 tem1 = 1. + fv * qtd(i,k)
246 thvd = thld(i,k) * tem1
247 endif
248 buo(i,k) = g * (1. - thvd / thvx(i,k))
249!
250 endif
251 enddo
252 enddo
253!
255!
256! tem = 1.-2.*f1
257! bb1 = 2. * b1 / tem
258! bb2 = 2. / tem
259! from Soares et al. (2004,QJRMS)
260! bb1 = 2.
261! bb2 = 4.
262!
263! from Bretherton et al. (2004, MWR)
264! bb1 = 4.
265! bb2 = 2.
266!
267! from our tuning
268 bb1 = 2.0
269 bb2 = 4.0
270!
271 do i = 1, im
272 if(cnvflg(i)) then
273 k = krad1(i)
274 dz = zm(i,k+1) - zm(i,k)
275! tem = 0.25*bb1*(xlamde(i,k)+xlamde(i,k+1))*dz
276 tem = 0.5*bb1*xlamde(i,k)*dz
277 tem1 = bb2 * buo(i,k+1) * dz
278 ptem1 = 1. + tem
279 wd2(i,k) = tem1 / ptem1
280 endif
281 enddo
282 do k = kmscu,1,-1
283 do i = 1, im
284 if(cnvflg(i) .and. k < krad1(i)) then
285 dz = zm(i,k+1) - zm(i,k)
286 tem = 0.25*bb1*(xlamde(i,k)+xlamde(i,k+1))*dz
287 tem1 = bb2 * buo(i,k+1) * dz
288 ptem = (1. - tem) * wd2(i,k+1)
289 ptem1 = 1. + tem
290 wd2(i,k) = (ptem + tem1) / ptem1
291 endif
292 enddo
293 enddo
294c
295 do i = 1, im
296 flg(i) = cnvflg(i)
297 mrady(i) = mrad(i)
298 if(flg(i)) mradx(i) = krad(i)
299 enddo
300 do k = kmscu,1,-1
301 do i = 1, im
302 if(flg(i) .and. k < krad(i)) then
303 if(wd2(i,k) > 0.) then
304 mradx(i) = k
305 else
306 flg(i)=.false.
307 endif
308 endif
309 enddo
310 enddo
311!
312 do i = 1,im
313 if(cnvflg(i)) then
314 if(mrad(i) < mradx(i)) then
315 mrad(i) = mradx(i)
316 endif
317 endif
318 enddo
319!
320 do i=1,im
321 if (cnvflg(i)) then
322 kk = krad(i)-mrad(i)
323 if(kk < 1) cnvflg(i)=.false.
324 endif
325 enddo
326!!
327 totflg = .true.
328 do i=1,im
329 totflg = totflg .and. (.not. cnvflg(i))
330 enddo
331 if(totflg) return
332!!
333!
335!
336 do k = 1, kmscu
337 do i=1,im
338 if(cnvflg(i) .and. mrady(i) < mradx(i)) then
339 dz = zl(i,k+1) - zl(i,k)
340 if(k >= mrad(i) .and. k < krad(i)) then
341 if(mrad(i) == 1) then
342 ptem = 1./(zm(i,k)+dz)
343 else
344 ptem = 1./(zm(i,k)-zm(i,mrad(i)-1)+dz)
345 endif
346 tem = max((hrad(i)-zm(i,k)+dz) ,dz)
347 ptem1 = 1./tem
348 xlamde(i,k) = ce0 * (ptem+ptem1)
349 else
350 xlamde(i,k) = ce0 / dz
351 endif
352 xlamdem(i,k) = cm * xlamde(i,k)
353 endif
354 enddo
355 enddo
356!
358!
359 do i = 1, im
360 xlamavg(i) = 0.
361 sumx(i) = 0.
362 enddo
363 do k = kmscu, 1, -1
364 do i = 1, im
365 if(cnvflg(i) .and.
366 & (k >= mrad(i) .and. k < krad(i))) then
367 dz = zl(i,k+1) - zl(i,k)
368 xlamavg(i) = xlamavg(i) + xlamde(i,k) * dz
369 sumx(i) = sumx(i) + dz
370 endif
371 enddo
372 enddo
373 do i = 1, im
374 if(cnvflg(i)) then
375 xlamavg(i) = xlamavg(i) / sumx(i)
376 endif
377 enddo
378!
380!
381 do k = kmscu, 1, -1
382 do i = 1, im
383 if(cnvflg(i) .and.
384 & (k >= mrad(i) .and. k < krad(i))) then
385 if(wd2(i,k) > 0.) then
386 tem = sqrt(wd2(i,k))
387 else
388 tem = 0.
389 endif
390 xmfd(i,k) = ra1(i) * tem
391 endif
392 enddo
393 enddo
394!
397!
398 do i = 1, im
399 if(cnvflg(i)) then
400 tem = 0.2 / xlamavg(i)
401 tem1 = 3.14 * tem * tem
402 sigma(i) = tem1 / (gdx(i) * gdx(i))
403 sigma(i) = max(sigma(i), 0.001)
404 sigma(i) = min(sigma(i), 0.999)
405 endif
406 enddo
407!
410!
411 do i = 1, im
412 if(cnvflg(i)) then
413 if (sigma(i) > ra1(i)) then
414 scaldfunc(i) = (1.-sigma(i)) * (1.-sigma(i))
415 scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.)
416 else
417 scaldfunc(i) = 1.0
418 endif
419 endif
420 enddo
421!
423!
424 do k = kmscu, 1, -1
425 do i = 1, im
426 if(cnvflg(i) .and.
427 & (k >= mrad(i) .and. k < krad(i))) then
428 xmfd(i,k) = scaldfunc(i) * xmfd(i,k)
429 dz = zl(i,k+1) - zl(i,k)
430 xmmx = dz / dt2
431 xmfd(i,k) = min(xmfd(i,k),xmmx)
432 endif
433 enddo
434 enddo
435!
436!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
438!
439 do i = 1, im
440 if(cnvflg(i)) then
441 k = krad(i)
442 thld(i,k)= thlx(i,k)
443 endif
444 enddo
445!
446! do i = 1, im
447! if(cnvflg(i)) then
448! k = krad(i)
449! ptem1 = max(qcdo(i,k,ntcw), 0.)
450! tld = thld(i,k) / pix(i,k)
451! tcdo(i,k) = tld + elocp * ptem1
452! qcdo(i,k,1) = qcdo(i,k,1)+0.2*qcdo(i,k,1)
453! qcdo(i,k,ntcw) = qcdo(i,k,ntcw)+0.2*qcdo(i,k,ntcw)
454! endif
455! enddo
456!
457 do k = kmscu,1,-1
458 do i=1,im
459 if(cnvflg(i) .and.
460 & (k >= mrad(i) .and. k < krad(i))) then
461 dz = zl(i,k+1) - zl(i,k)
462 tem = 0.5 * xlamde(i,k) * dz
463 factor = 1. + tem
464!
465 thld(i,k) = ((1.-tem)*thld(i,k+1)+tem*
466 & (thlx(i,k)+thlx(i,k+1)))/factor
467 qtd(i,k) = ((1.-tem)*qtd(i,k+1)+tem*
468 & (qtx(i,k)+qtx(i,k+1)))/factor
469!
470 tld = thld(i,k) / pix(i,k)
471 es = 0.01 * fpvs(tld) ! fpvs in pa
472 qs = max(qmin, eps * es / (plyr(i,k)+epsm1*es))
473 dq = qtd(i,k) - qs
474!
475 if (dq > 0.) then
476 gamma = el2orc * qs / (tld**2)
477 qld = dq / (1. + gamma)
478 qtd(i,k) = qs + qld
479 qcdo(i,k,1) = qs
480 qcdo(i,k,ntcw) = qld
481 tcdo(i,k) = tld + elocp * qld
482 else
483 qcdo(i,k,1) = qtd(i,k)
484 qcdo(i,k,ntcw) = 0.
485 tcdo(i,k) = tld
486 endif
487!
488 endif
489 enddo
490 enddo
491!
492 do k = kmscu, 1, -1
493 do i = 1, im
494 if (cnvflg(i) .and. k < krad(i)) then
495 if(k >= mrad(i)) then
496 dz = zl(i,k+1) - zl(i,k)
497 tem = 0.5 * xlamdem(i,k) * dz
498 factor = 1. + tem
499 ptem = tem - pgcon
500 ptem1= tem + pgcon
501!
502 ucdo(i,k) = ((1.-tem)*ucdo(i,k+1)+ptem*u1(i,k+1)
503 & +ptem1*u1(i,k))/factor
504 vcdo(i,k) = ((1.-tem)*vcdo(i,k+1)+ptem*v1(i,k+1)
505 & +ptem1*v1(i,k))/factor
506 endif
507 endif
508 enddo
509 enddo
510!
511 if(ntcw > 2) then
512!
513 do n = 2, ntcw-1
514 do k = kmscu, 1, -1
515 do i = 1, im
516 if (cnvflg(i) .and. k < krad(i)) then
517 if(k >= mrad(i)) then
518 dz = zl(i,k+1) - zl(i,k)
519 tem = 0.5 * xlamde(i,k) * dz
520 factor = 1. + tem
521!
522 qcdo(i,k,n) = ((1.-tem)*qcdo(i,k+1,n)+tem*
523 & (q1(i,k,n)+q1(i,k+1,n)))/factor
524 endif
525 endif
526 enddo
527 enddo
528 enddo
529!
530 endif
531!
532 ndc = ntrac1 - ntcw
533!
534 if(ndc > 0) then
535!
536 do n = ntcw+1, ntrac1
537 do k = kmscu, 1, -1
538 do i = 1, im
539 if (cnvflg(i) .and. k < krad(i)) then
540 if(k >= mrad(i)) then
541 dz = zl(i,k+1) - zl(i,k)
542 tem = 0.5 * xlamde(i,k) * dz
543 factor = 1. + tem
544!
545 qcdo(i,k,n) = ((1.-tem)*qcdo(i,k+1,n)+tem*
546 & (q1(i,k,n)+q1(i,k+1,n)))/factor
547 endif
548 endif
549 enddo
550 enddo
551 enddo
552!
553 endif
554!
555 return
556 end
557
558 end module mfscu_mod
subroutine mfscu(im, ix, km, kmscu, ntcw, ntrac1, delt, cnvflg, zl, zm, q1, t1, u1, v1, plyr, pix, thlx, thvx, thlvx, gdx, thetae, radj, krad, mrad, radmin, buo, xmfd, tcdo, qcdo, ucdo, vcdo, xlamde)
This subroutine computes mass flux and downdraft parcel properties for stratocumulus-top-driven turbu...
Definition mfscu.f:16