CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
satmedmfvdifq.F
1
2
6 module satmedmfvdifq
7 use mfpbltq_mod
8 use tridi_mod
9 use mfscuq_mod
10 contains
11
32 subroutine satmedmfvdifq_init (satmedmf, &
33 & isatmedmf,isatmedmf_vdifq, &
34 & errmsg,errflg)
35
36 logical, intent(in ) :: satmedmf
37 integer, intent(in) :: isatmedmf,isatmedmf_vdifq
38
39 character(len=*), intent(out) :: errmsg
40 integer, intent(out) :: errflg
41
42! Initialize CCPP error handling variables
43 errmsg = ''
44 errflg = 0
45
46! Consistency checks
47 if (.not. satmedmf) then
48 write(errmsg,fmt='(*(a))') 'Logic error: satmedmf = .false.'
49 errflg = 1
50 return
51 end if
52
53 if (.not. isatmedmf==isatmedmf_vdifq) then
54 write(errmsg,fmt='(*(a))') 'Logic error: satmedmfvdif is ',
55 & 'called, but isatmedmf/=isatmedmf_vdifq.'
56 errflg = 1
57 return
58 end if
59
60 end subroutine satmedmfvdifq_init
61
76 subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
77 & ntiw,ntke,grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, &
78 & dv,du,tdt,rtg,u1,v1,t1,q1,usfco,vsfco,icplocn2atm, &
79 & swh,hlw,xmu,garea,zvfun,sigmaf, &
80 & psk,rbsoil,zorl,u10m,v10m,fm,fh, &
81 & tsea,heat,evap,stress,spd1,kpbl, &
82 & prsi,del,prsl,prslk,phii,phil,delt, &
83 & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku, &
84 & kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr, &
85 & rlmx,elmx,sfc_rlm,tc_pbl, &
86 & ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind, &
87 & index_of_y_wind,index_of_process_pbl,gen_tend,ldiag3d, &
88 & errmsg,errflg)
89!
90 use machine , only : kind_phys
91 use funcphys , only : fpvs
92!
93 implicit none
94!
95!----------------------------------------------------------------------
96 integer, intent(in) :: im, km, ntrac, ntcw, ntrw, ntiw, &
97 & ntke, ntqv
98 integer, intent(in) :: sfc_rlm
99 integer, intent(in) :: tc_pbl
100 integer, intent(in) :: kinver(:)
101 integer, intent(out) :: kpbl(:)
102 logical, intent(in) :: gen_tend,ldiag3d
103!
104 real(kind=kind_phys), intent(in) :: grav,rd,cp,rv,hvap,hfus,fv, &
105 & eps,epsm1
106 real(kind=kind_phys), intent(in) :: delt, xkzm_m, xkzm_h, xkzm_s
107 real(kind=kind_phys), intent(in) :: dspfac, bl_upfr, bl_dnfr
108 real(kind=kind_phys), intent(in) :: rlmx, elmx
109 real(kind=kind_phys), intent(inout) :: dv(:,:), du(:,:), &
110 & tdt(:,:), rtg(:,:,:)
111 real(kind=kind_phys), intent(in) :: &
112 & u1(:,:), v1(:,:), &
113 & usfco(:), vsfco(:), &
114 & t1(:,:), q1(:,:,:), &
115 & swh(:,:), hlw(:,:), &
116 & xmu(:), garea(:), &
117 & zvfun(:), sigmaf(:), &
118 & psk(:), rbsoil(:), &
119 & zorl(:), tsea(:), &
120 & u10m(:), v10m(:), &
121 & fm(:), fh(:), &
122 & evap(:), heat(:), &
123 & stress(:), spd1(:), &
124 & prsi(:,:), del(:,:), &
125 & prsl(:,:), prslk(:,:), &
126 & phii(:,:), phil(:,:)
127 real(kind=kind_phys), intent(inout), dimension(:,:,:), optional ::&
128 & dtend
129 integer, intent(in) :: dtidx(:,:), index_of_temperature, &
130 & index_of_x_wind, index_of_y_wind, index_of_process_pbl
131 integer, intent(in) :: icplocn2atm
132 real(kind=kind_phys), intent(out) :: &
133 & dusfc(:), dvsfc(:), &
134 & dtsfc(:), dqsfc(:), &
135 & hpbl(:)
136 real(kind=kind_phys), intent(out) :: &
137 & dkt(:,:), dku(:,:)
138!
139 logical, intent(in) :: dspheat
140 character(len=*), intent(out) :: errmsg
141 integer, intent(out) :: errflg
142
143! flag for tke dissipative heating
144!
145!----------------------------------------------------------------------
146!***
147!*** local variables
148 real(kind=kind_phys) spd1_m
149!***
150 integer i,is,k,n,ndt,km1,kmpbl,kmscu,ntrac1,idtend
151 integer kps,kbx,kmx
152 integer lcld(im),kcld(im),krad(im),mrad(im)
153 integer kx1(im), kb1(im), kpblx(im)
154!
155 real(kind=kind_phys) tke(im,km), tkeh(im,km-1), e2(im,0:km)
156!
157 real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km),
158 & qlx(im,km), thetae(im,km),thlx(im,km),
159 & slx(im,km), svx(im,km), qtx(im,km),
160 & tvx(im,km), pix(im,km), radx(im,km-1),
161 & dkq(im,km-1),cku(im,km-1), ckt(im,km-1)
162!
163 real(kind=kind_phys) plyr(im,km), rhly(im,km), cfly(im,km),
164 & qstl(im,km)
165!
166 real(kind=kind_phys) dtdz1(im), gdx(im),
167 & phih(im), phim(im),
168 & phims(im), prn(im,km-1),
169 & rbdn(im), rbup(im), thermal(im),
170 & ustar(im), wstar(im), hpblx(im),
171 & ust3(im), wst3(im), rho_a(im),
172 & z0(im), crb(im), tkemean(im),
173 & hgamt(im), hgamq(im),
174 & wscale(im),vpert(im),
175 & zol(im), sflux(im),
176 & sumx(im), tx1(im), tx2(im)
177!
178 real(kind=kind_phys) radmin(im)
179!
180 real(kind=kind_phys) zi(im,km+1), zl(im,km), zm(im,km),
181 & xkzo(im,km-1),xkzmo(im,km-1),
182 & xkzm_hx(im), xkzm_mx(im), tkmnz(im,km-1),
183 & rdzt(im,km-1),rlmnz(im,km),
184 & al(im,km-1), ad(im,km), au(im,km-1),
185 & f1(im,km), f2(im,km*(ntrac-1))
186!
187 real(kind=kind_phys) elm(im,km), ele(im,km),
188 & ckz(im,km), chz(im,km),
189 & diss(im,km-1),prod(im,km-1),
190 & bf(im,km-1), shr2(im,km-1), wush(im,km),
191 & xlamue(im,km-1), xlamde(im,km-1),
192 & gotvx(im,km), rlam(im,km-1)
193!
194! variables for updrafts (thermals)
195!
196 real(kind=kind_phys) tcko(im,km), qcko(im,km,ntrac),
197 & ucko(im,km), vcko(im,km),
198 & buou(im,km), xmf(im,km)
199!
200! variables for stratocumulus-top induced downdrafts
201!
202 real(kind=kind_phys) tcdo(im,km), qcdo(im,km,ntrac),
203 & ucdo(im,km), vcdo(im,km),
204 & buod(im,km), xmfd(im,km)
205!
206! variables for Total Variation Diminishing (TVD) flux-limiter scheme
207!
208 real(kind=kind_phys) e_half(im,km-1), e_diff(im,0:km-1),
209 & q_half(im,km-1,ntrac-1),
210 & qh(im,km-1,ntrac-1),
211 & q_diff(im,0:km-1,ntrac-1)
212 real(kind=kind_phys) rrkp, phkp
213 real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im)
214 real(kind=kind_phys) sfcpbl(im), vez0fun(im)
215!
216 logical pblflg(im), sfcflg(im), flg(im)
217 logical scuflg(im), pcnvflg(im)
218 logical mlenflg
219!
220! pcnvflg: true for unstable pbl
221!
222 real(kind=kind_phys) aphi16, aphi5,
223 & wfac, cfac,
224 & gamcrt, gamcrq, sfcfrac,
225! & conq, cont, conw,
226 & dsdz2, dsdzt, dkmax,
227 & dsig, dt2, dtodsd,
228 & dtodsu, g, factor, dz,
229 & gocp, gravi, zol1, zolcru,
230 & buop, shrp, dtn,
231 & prnum, prmax, prmin, prtke,
232 & prscu, pr0, ri,
233 & dw2, dw2min, zk,
234 & elmfac, elefac, dspmax,
235 & alp, clwt, cql,
236 & f0, robn, crbmin, crbmax,
237 & es, qs, value, onemrh,
238 & cfh, gamma, elocp, el2orc,
239 & epsi, beta, chx, cqx,
240 & rdt, rdz, qmin, qlmin,
241 & rimin, rbcr, rbint, tdzmin,
242 & rlmn, rlmn0, rlmn1, rlmn2,
243 & ttend, utend, vtend, qtend,
244 & zfac, zfmin, vk, spdk2,
245 & tkmin, tkbmx, xkgdx,
246 & xkinv1, xkinv2,
247 & zlup, zldn, cs0, csmf,
248 & tem, tem1, tem2, tem3,
249 & ptem, ptem0, ptem1, ptem2
250!
251 real(kind=kind_phys) slfac
252!
253 real(kind=kind_phys) vegflo, vegfup, z0lo, z0up, vc0, zc0
254!
255 real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck
256!
257 real(kind=kind_phys) qlcr, zstblmax, hcrinv
258!
259 real(kind=kind_phys) h1
260
261 real(kind=kind_phys) bfac, mffac
262!!
263 parameter(bfac=100.)
264 parameter(wfac=7.0,cfac=4.5)
265 parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1)
266 parameter(vk=0.4,rimin=-100.,slfac=0.1)
267 parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3)
268 parameter(rlmn=30.,rlmn0=5.,rlmn1=5.,rlmn2=10.)
269 parameter(prmin=0.25,prmax=4.0)
270 parameter(pr0=1.0,prtke=1.0,prscu=0.67)
271 parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35)
272 parameter(tkmin=1.e-9,tkbmx=0.2,dspmax=10.0)
273 parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8)
274 parameter(aphi5=5.,aphi16=16.)
275 parameter(elmfac=1.0,elefac=1.0,cql=100.)
276 parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=1000.)
277 parameter(qlcr=3.5e-5,zstblmax=2500.)
278 parameter(xkinv1=0.4,xkinv2=0.3)
279 parameter(h1=0.33333333,hcrinv=250.)
280 parameter(vegflo=0.1,vegfup=1.0,z0lo=0.1,z0up=1.0)
281 parameter(vc0=1.0,zc0=1.0)
282 parameter(ck1=0.15,ch1=0.15)
283 parameter(cs0=0.4,csmf=0.5)
284 parameter(rchck=1.5,ndt=20)
285
286 if (tc_pbl == 0) then
287 ck0 = 0.4
288 ch0 = 0.4
289 ce0 = 0.4
290 else if (tc_pbl == 1) then
291 ck0 = 0.55
292 ch0 = 0.55
293 ce0 = 0.12
294 endif
295 gravi = 1.0 / grav
296 g = grav
297 gocp = g / cp
298! cont=cp/g
299! conq=hvap/g
300! conw=1.0/g ! for del in pa
301!! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) !kpa
302 elocp = hvap / cp
303 el2orc = hvap * hvap / (rv * cp)
304!
305!************************************************************************
306! Initialize CCPP error handling variables
307 errmsg = ''
308 errflg = 0
309
311 dt2 = delt
312 rdt = 1. / dt2
313!
314! the code is written assuming ntke=ntrac
315! if ntrac > ntke, the code needs to be modified
316!
317 ntrac1 = ntrac - 1
318 km1 = km - 1
319 kmpbl = km / 2
320 kmscu = km / 2
323 do k=1,km
324 do i=1,im
325 zi(i,k) = phii(i,k) * gravi
326 zl(i,k) = phil(i,k) * gravi
327 xmf(i,k) = 0.
328 xmfd(i,k) = 0.
329 buou(i,k) = 0.
330 buod(i,k) = 0.
331 wush(i,k) = 0.
332 ckz(i,k) = ck1
333 chz(i,k) = ch1
334 rlmnz(i,k) = rlmn0
335 enddo
336 enddo
337 do i=1,im
338 zi(i,km+1) = phii(i,km+1) * gravi
339 enddo
340 do k=1,km
341 do i=1,im
342 zm(i,k) = zi(i,k+1)
343 enddo
344 enddo
346 do i=1,im
347 gdx(i) = sqrt(garea(i))
348 enddo
351 do k=1,km
352 do i=1,im
353 tke(i,k) = max(q1(i,k,ntke), tkmin)
354 enddo
355 enddo
356 do k=1,km1
357 do i=1,im
358 tkeh(i,k) = 0.5 * (tke(i,k) + tke(i,k+1))
359 enddo
360 enddo
362 do k = 1,km1
363 do i=1,im
364 rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k))
365 prn(i,k) = pr0
366 enddo
367 enddo
368!
370
372
374
379!
380 do i=1,im
381 kx1(i) = 1
382 tx1(i) = 1.0 / prsi(i,1)
383 tx2(i) = tx1(i)
384 if(gdx(i) >= xkgdx) then
385 xkzm_hx(i) = xkzm_h
386 xkzm_mx(i) = xkzm_m
387 else
388 tem = gdx(i) / xkgdx
389 xkzm_hx(i) = xkzm_h * tem
390 xkzm_mx(i) = xkzm_m * tem
391 endif
392 enddo
393 do k = 1,km1
394 do i=1,im
395 xkzo(i,k) = 0.0
396 xkzmo(i,k) = 0.0
397 if (k < kinver(i)) then
398! minimum turbulent mixing length
399 ptem = prsi(i,k+1) * tx1(i)
400 tem1 = 1.0 - ptem
401 tem2 = tem1 * tem1 * 2.5
402 tem2 = min(1.0, exp(-tem2))
403 rlmnz(i,k)= rlmn * tem2
404 rlmnz(i,k)= max(rlmnz(i,k), rlmn0)
405! vertical background diffusivity
406 tem2 = tem1 * tem1 * 10.0
407 tem2 = min(1.0, exp(-tem2))
408 xkzo(i,k) = xkzm_hx(i) * tem2
409! vertical background diffusivity for
410! momentum
411 if (ptem >= xkzm_s) then
412 xkzmo(i,k) = xkzm_mx(i)
413 kx1(i) = k + 1
414 else
415 if (k == kx1(i) .and. k > 1) tx2(i) = 1.0 / prsi(i,k)
416 tem1 = 1.0 - prsi(i,k+1) * tx2(i)
417 tem1 = tem1 * tem1 * 5.0
418 xkzmo(i,k) = xkzm_mx(i) * min(1.0, exp(-tem1))
419 endif
420 endif
421 enddo
422 enddo
423!
425 do i = 1,im
426 z0(i) = 0.01 * zorl(i)
427 rho_a(i) = prsl(i,1)/(rd*t1(i,1)*(1.+fv*max(q1(i,1,1),qmin)))
428 dusfc(i) = 0.
429 dvsfc(i) = 0.
430 dtsfc(i) = 0.
431 dqsfc(i) = 0.
432 kpbl(i) = 1
433 hpbl(i) = 0.
434 kpblx(i) = 1
435 hpblx(i) = 0.
436 pblflg(i)= .true.
437 sfcflg(i)= .true.
438 if(rbsoil(i) > 0.) sfcflg(i) = .false.
439 pcnvflg(i)= .false.
440 scuflg(i)= .true.
441 if(scuflg(i)) then
442 radmin(i)= 0.
443 mrad(i) = km1
444 krad(i) = 1
445 lcld(i) = km1
446 kcld(i) = km1
447 endif
448 enddo
449!
453!
454 do i = 1,im
455 tem = (sigmaf(i) - vegflo) / (vegfup - vegflo)
456 tem = min(max(tem, 0.), 1.)
457 tem1 = sqrt(tem)
458 ptem = (z0(i) - z0lo) / (z0up - z0lo)
459 ptem = min(max(ptem, 0.), 1.)
460 vez0fun(i) = (1. + vc0 * tem1) * (1. + zc0 * ptem)
461 enddo
462!
465 do k=1,km
466 do i=1,im
467 pix(i,k) = psk(i) / prslk(i,k)
468 theta(i,k) = t1(i,k) * pix(i,k)
469 if(ntiw > 0) then
470 tem = max(q1(i,k,ntcw),qlmin)
471 tem1 = max(q1(i,k,ntiw),qlmin)
472 qlx(i,k) = tem + tem1
473 ptem = hvap*tem + (hvap+hfus)*tem1
474 slx(i,k) = cp * t1(i,k) + phil(i,k) - ptem
475 else
476 qlx(i,k) = max(q1(i,k,ntcw),qlmin)
477 slx(i,k) = cp * t1(i,k) + phil(i,k) - hvap*qlx(i,k)
478 endif
479 tem2 = 1.+fv*max(q1(i,k,1),qmin)-qlx(i,k)
480 thvx(i,k) = theta(i,k) * tem2
481 tvx(i,k) = t1(i,k) * tem2
482 qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k)
483 thlx(i,k) = theta(i,k) - pix(i,k)*elocp*qlx(i,k)
484 thlvx(i,k) = thlx(i,k) * (1. + fv * qtx(i,k))
485 svx(i,k) = cp * tvx(i,k)
486 ptem1 = elocp * pix(i,k) * max(q1(i,k,1),qmin)
487 thetae(i,k)= theta(i,k) + ptem1
488! gotvx(i,k) = g / tvx(i,k)
489 gotvx(i,k) = g / thvx(i,k)
490 enddo
491 enddo
492
495 do k = 1, km
496 do i = 1, im
497 plyr(i,k) = 0.01 * prsl(i,k) ! pa to mb (hpa)
498! --- ... compute relative humidity
499 es = 0.01 * fpvs(t1(i,k)) ! fpvs in pa
500 qs = max(qmin, eps * es / (plyr(i,k) + epsm1*es))
501 rhly(i,k) = max(0.0, min(1.0, max(qmin, q1(i,k,1))/qs))
502 qstl(i,k) = qs
503 enddo
504 enddo
505!
506 do k = 1, km
507 do i = 1, im
508 cfly(i,k) = 0.
509 clwt = 1.0e-6 * (plyr(i,k)*0.001)
510 if (qlx(i,k) > clwt) then
511 onemrh = max(1.e-10, 1.0-rhly(i,k))
512 tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0)
513 tem1 = cql / tem1
514 value = max(min( tem1*qlx(i,k), 50.0), 0.0)
515 tem2 = sqrt(sqrt(rhly(i,k)))
516 cfly(i,k) = min(max(tem2*(1.0-exp(-value)), 0.0), 1.0)
517 endif
518 enddo
519 enddo
520!
522!
523 do k = 1, km1
524 do i = 1, im
525 tem = 0.5 * (svx(i,k) + svx(i,k+1))
526 tem1 = 0.5 * (t1(i,k) + t1(i,k+1))
527 tem2 = 0.5 * (qstl(i,k) + qstl(i,k+1))
528 cfh = min(cfly(i,k+1),0.5*(cfly(i,k)+cfly(i,k+1)))
529 alp = g / tem
530 gamma = el2orc * tem2 / (tem1**2)
531 epsi = tem1 / elocp
532 beta = (1. + gamma*epsi*(1.+fv)) / (1. + gamma)
533 chx = cfh * alp * beta + (1. - cfh) * alp
534 cqx = cfh * alp * hvap * (beta - epsi)
535 cqx = cqx + (1. - cfh) * fv * g
536 ptem1 = (slx(i,k+1)-slx(i,k))*rdzt(i,k)
537 ptem2 = (qtx(i,k+1)-qtx(i,k))*rdzt(i,k)
538 bf(i,k) = chx * ptem1 + cqx * ptem2
539 enddo
540 enddo
541!
542!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
543!
546 do k=1,km
547 do i=1,im
548 dku(i,k) = 0.
549 dkt(i,k) = 0.
550 enddo
551 enddo
552 do k=1,km1
553 do i=1,im
554 dkq(i,k) = 0.
555 cku(i,k) = 0.
556 ckt(i,k) = 0.
557 tem = zi(i,k+1)-zi(i,k)
558 radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k))
559 enddo
560 enddo
564 do i = 1,im
565 sflux(i) = heat(i) + evap(i)*fv*theta(i,1)
566 if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false.
567 enddo
568!
586 do i = 1,im
587 if(pblflg(i)) then
588! thermal(i) = thvx(i,1)
589 thermal(i) = thlvx(i,1)
590 crb(i) = rbcr
591 else
592 thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
593 tem = sqrt(u10m(i)**2+v10m(i)**2)
594 tem = max(tem, 1.)
595 robn = tem / (f0 * z0(i))
596 tem1 = 1.e-7 * robn
597 crb(i) = 0.16 * (tem1 ** (-0.18))
598 crb(i) = max(min(crb(i), crbmax), crbmin)
599 endif
600 enddo
602 do i=1,im
603 dtdz1(i) = dt2 / (zi(i,2)-zi(i,1))
604 enddo
605!
606 do i=1,im
607 ustar(i) = sqrt(stress(i))
608 enddo
609!
612!
613 do k = 1, km1
614 do i = 1, im
615 rdz = rdzt(i,k)
616! bf(i,k) = gotvx(i,k)*(thvx(i,k+1)-thvx(i,k))*rdz
617 dw2 = (u1(i,k)-u1(i,k+1))**2
618 & + (v1(i,k)-v1(i,k+1))**2
619 shr2(i,k) = max(dw2,dw2min)*rdz*rdz
620 enddo
621 enddo
622!
623! Find first quess pbl height based on bulk richardson number (mrf pbl scheme)
624! and also for diagnostic purpose
625!
626 do i=1,im
627 flg(i) = .false.
628 rbup(i) = rbsoil(i)
629 enddo
635 do k = 1, kmpbl
636 do i = 1, im
637 if(.not.flg(i)) then
638 rbdn(i) = rbup(i)
639 if (tc_pbl == 0) then
640 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
641! rbup(i) = (thvx(i,k)-thermal(i))*
642! & (g*zl(i,k)/thvx(i,1))/spdk2
643 rbup(i) = (thlvx(i,k)-thermal(i))*
644 & (g*zl(i,k)/thlvx(i,1))/spdk2
645 else if (tc_pbl == 1) then
646 spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.)
647 & + bfac*ustar(i)**2
648 rbup(i) = (thlvx(i,k)-thermal(i))*
649 & (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2
650 endif
651 kpblx(i) = k
652 flg(i) = rbup(i) > crb(i)
653 endif
654 enddo
655 enddo
659 do i = 1,im
660 if(kpblx(i) > 1) then
661 k = kpblx(i)
662 if(rbdn(i) >= crb(i)) then
663 rbint = 0.
664 elseif(rbup(i) <= crb(i)) then
665 rbint = 1.
666 else
667 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
668 endif
669 hpblx(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
670 if(hpblx(i) < zi(i,kpblx(i))) kpblx(i)=kpblx(i)-1
671 else
672 hpblx(i) = zl(i,1)
673 kpblx(i) = 1
674 endif
675 hpbl(i) = hpblx(i)
676 kpbl(i) = kpblx(i)
677 if(kpbl(i) <= 1) pblflg(i)=.false.
678 enddo
679!
680! update thermal at a level of slfac*hpbl for unstable pbl
681!
682 do i=1,im
683 sfcpbl(i) = slfac * hpbl(i)
684 kb1(i) = 1
685 flg(i) = .false.
686 if(pblflg(i)) then
687 flg(i) = .true.
688 endif
689 enddo
690 do k = 2, kmpbl
691 do i=1,im
692 if (flg(i) .and. zl(i,k) <= sfcpbl(i)) then
693 kb1(i) = k
694 else
695 flg(i) = .false.
696 endif
697 enddo
698 enddo
699 do i=1,im
700 if(pblflg(i)) kb1(i)=min(kb1(i),kpbl(i))
701 enddo
702!
703! re-compute pbl height with the updated thermal
704!
705 do i=1,im
706 flg(i) = .true.
707 if(pblflg(i) .and. kb1(i) > 1) then
708 flg(i) = .false.
709 rbup(i) = rbsoil(i)
710! thermal(i) = thvx(i,kb1(i))
711 thermal(i) = thlvx(i,kb1(i))
712 kpblx(i) = kb1(i)
713 hpblx(i) = zl(i,kb1(i))
714 endif
715 enddo
716 do k = 2, kmpbl
717 do i = 1, im
718 if(.not.flg(i) .and. k > kb1(i)) then
719 rbdn(i) = rbup(i)
720 if (tc_pbl == 0) then
721 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
722! rbup(i) = (thvx(i,k)-thermal(i))*
723! & (g*zl(i,k)/thvx(i,1))/spdk2
724 rbup(i) = (thlvx(i,k)-thermal(i))*
725 & (g*zl(i,k)/thlvx(i,1))/spdk2
726 else if (tc_pbl == 1) then
727 spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.)
728 & + bfac*ustar(i)**2
729 rbup(i) = (thlvx(i,k)-thermal(i))*
730 & (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2
731 endif
732 kpblx(i) = k
733 flg(i) = rbup(i) > crb(i)
734 endif
735 enddo
736 enddo
737 do i = 1,im
738 if(pblflg(i) .and. kb1(i) > 1) then
739 k = kpblx(i)
740 if(rbdn(i) >= crb(i)) then
741 rbint = 0.
742 elseif(rbup(i) <= crb(i)) then
743 rbint = 1.
744 else
745 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
746 endif
747 hpblx(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
748 if(hpblx(i) < zi(i,kpblx(i))) kpblx(i)=kpblx(i)-1
749 hpbl(i) = hpblx(i)
750 kpbl(i) = kpblx(i)
751 endif
752 enddo
753!
755!
756 do i = 1, im
757 sumx(i) = 0.
758 tkemean(i) = 0.
759 enddo
760 do k = 1, kmpbl
761 do i = 1, im
762 if(k < kpbl(i)) then
763 dz = zi(i,k+1) - zi(i,k)
764 tkemean(i) = tkemean(i) + tke(i,k) * dz
765 sumx(i) = sumx(i) + dz
766 endif
767 enddo
768 enddo
769 do i = 1, im
770 if(tkemean(i) > 0. .and. sumx(i) > 0.) then
771 tkemean(i) = tkemean(i) / sumx(i)
772 endif
773 enddo
774!
777!
778 kps = max(kmpbl, kmscu)
779 do k = 2, kps
780 do i = 1, im
781 dz = zi(i,k+1) - zi(i,k)
782 tem = (0.5*(u1(i,k-1)-u1(i,k+1))/dz)**2
783 tem1 = tem+(0.5*(v1(i,k-1)-v1(i,k+1))/dz)**2
784 wush(i,k) = csmf * sqrt(tem1)
785 enddo
786 enddo
787!
797 do i=1,im
798 zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
799 if(sfcflg(i)) then
800 zol(i) = min(zol(i),-zfmin)
801 else
802 zol(i) = max(zol(i),zfmin)
803 endif
815 zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1)
816 if(sfcflg(i)) then
817 tem = 1.0 / (1. - aphi16*zol1)
818 phih(i) = sqrt(tem)
819 phim(i) = sqrt(phih(i))
820 tem1 = 1.0 / (1. - aphi16*zol(i))
821 phims(i) = sqrt(sqrt(tem1))
822 else
823 phim(i) = 1. + aphi5*zol1
824 phih(i) = phim(i)
825 phims(i) = 1. + aphi5*zol(i)
826 endif
827 enddo
828!
844 do i=1,im
845 if(pblflg(i)) then
846 if(zol(i) < zolcru) then
847 pcnvflg(i) = .true.
848 endif
849 wst3(i) = gotvx(i,1)*sflux(i)*hpbl(i)
850 wstar(i)= wst3(i)**h1
851 ust3(i) = ustar(i)**3.
852 wscale(i)=(ust3(i)+wfac*vk*wst3(i)*sfcfrac)**h1
853 ptem = ustar(i)/aphi5
854 wscale(i) = max(wscale(i),ptem)
855 endif
856 enddo
857!
860!
861 do i = 1,im
862 if(pcnvflg(i)) then
863 hgamt(i) = heat(i)/wscale(i)
864 hgamq(i) = evap(i)/wscale(i)
865 vpert(i) = hgamt(i) + hgamq(i)*fv*theta(i,1)
866 vpert(i) = max(vpert(i),0.)
867 tem = min(cfac*vpert(i),gamcrt)
868 thermal(i)= thermal(i) + tem
869 endif
870 enddo
871!
872! enhance the pbl height by considering the thermal excess
873! (overshoot pbl top)
874!
875 do i=1,im
876 flg(i) = .true.
877 if(pcnvflg(i)) then
878 flg(i) = .false.
879 rbup(i) = rbsoil(i)
880 endif
881 enddo
882 do k = 2, kmpbl
883 do i = 1, im
884 if(.not.flg(i) .and. k > kb1(i)) then
885 rbdn(i) = rbup(i)
886 if (tc_pbl == 0) then
887 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
888 rbup(i) = (thlvx(i,k)-thermal(i))*
889 & (g*zl(i,k)/thlvx(i,1))/spdk2
890 else if (tc_pbl == 1) then
891 spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.)
892 & + bfac*ustar(i)**2
893 rbup(i) = (thlvx(i,k)-thermal(i))*
894 & (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2
895 endif
896 kpbl(i) = k
897 flg(i) = rbup(i) > crb(i)
898 endif
899 enddo
900 enddo
901 do i = 1,im
902 if(pcnvflg(i)) then
903 k = kpbl(i)
904 if(rbdn(i) >= crb(i)) then
905 rbint = 0.
906 elseif(rbup(i) <= crb(i)) then
907 rbint = 1.
908 else
909 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
910 endif
911 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
912 if(hpbl(i) < zi(i,kpbl(i))) then
913 kpbl(i) = kpbl(i) - 1
914 endif
915 if(kpbl(i) <= 1) then
916 pcnvflg(i) = .false.
917 pblflg(i) = .false.
918 endif
919 endif
920 enddo
921!
922!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
923! look for stratocumulus
928 do i=1,im
929 flg(i) = scuflg(i)
930 enddo
931 do k = 1, km1
932 do i=1,im
933 if(flg(i).and.zl(i,k) >= zstblmax) then
934 lcld(i)=k
935 flg(i)=.false.
936 endif
937 enddo
938 enddo
939 do i = 1, im
940 flg(i)=scuflg(i)
941 enddo
942 do k = kmscu,1,-1
943 do i = 1, im
944 if(flg(i) .and. k <= lcld(i)) then
945 if(qlx(i,k) >= qlcr) then
946 kcld(i)=k
947 flg(i)=.false.
948 endif
949 endif
950 enddo
951 enddo
952 do i = 1, im
953 if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false.
954 enddo
960 do i = 1, im
961 flg(i)=scuflg(i)
962 enddo
963 do k = kmscu,1,-1
964 do i = 1, im
965 if(flg(i) .and. k <= kcld(i)) then
966 if(qlx(i,k) >= qlcr) then
967 if(radx(i,k) < radmin(i)) then
968 radmin(i)=radx(i,k)
969 krad(i)=k
970 endif
971 else
972 flg(i)=.false.
973 endif
974 endif
975 enddo
976 enddo
977 do i = 1, im
978 if(scuflg(i) .and. krad(i) <= 1) scuflg(i)=.false.
979 if(scuflg(i) .and. radmin(i)>=0.) scuflg(i)=.false.
980 enddo
981!
982!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
987 do k = 1, km
988 do i = 1, im
989 if(pcnvflg(i)) then
990 tcko(i,k) = t1(i,k)
991 ucko(i,k) = u1(i,k)
992 vcko(i,k) = v1(i,k)
993 endif
994 if(scuflg(i)) then
995 tcdo(i,k) = t1(i,k)
996 ucdo(i,k) = u1(i,k)
997 vcdo(i,k) = v1(i,k)
998 endif
999 enddo
1000 enddo
1001 do n = 1, ntrac1
1002 do k = 1, km
1003 do i = 1, im
1004 if(pcnvflg(i)) then
1005 qcko(i,k,n) = q1(i,k,n)
1006 endif
1007 if(scuflg(i)) then
1008 qcdo(i,k,n) = q1(i,k,n)
1009 endif
1010 enddo
1011 enddo
1012 enddo
1015 call mfpbltq(im,im,km,kmpbl,ntcw,ntrac1,dt2,
1016 & pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx,
1017 & gdx,hpbl,kpbl,vpert,buou,wush,tkemean,vez0fun,xmf,
1018 & tcko,qcko,ucko,vcko,xlamue,bl_upfr)
1021 call mfscuq(im,im,km,kmscu,ntcw,ntrac1,dt2,
1022 & scuflg,zl,zm,q1,t1,u1,v1,plyr,pix,
1023 & thlx,thvx,thlvx,gdx,thetae,
1024 & krad,mrad,radmin,buod,wush,tkemean,vez0fun,xmfd,
1025 & tcdo,qcdo,ucdo,vcdo,xlamde,bl_dnfr)
1026
1027 if (tc_pbl == 1) then
1029 do i=1,im
1030 if(zol(i) > -0.5) then
1031 do k = 1, km
1032 xmf(i,k) = 0.0
1033 end do
1034 end if
1035 end do
1037 do i = 1,im
1038 tem = sqrt(u10m(i)**2+v10m(i)**2)
1039 mffac = (1. - min(max(tem - 20.0, 0.0), 10.0)/10.)
1040 do k = 1, km
1041 xmf(i,k) = xmf(i,k)*mffac
1042 xmfd(i,k) = xmfd(i,k)*mffac
1043 enddo
1044 enddo
1045 endif
1046!
1047!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1048
1050 do k = 1, kmpbl
1051 do i = 1, im
1052 if(k < kpbl(i)) then
1053 tem = phih(i)/phim(i)
1054 ptem = sfcfrac*hpbl(i)
1055 tem1 = max(zi(i,k+1)-ptem, 0.)
1056 tem2 = tem1 / (hpbl(i) - ptem)
1057 if(pcnvflg(i)) then
1058 tem = min(tem, pr0)
1059 prn(i,k) = tem + (pr0 - tem) * tem2
1060 else
1061 tem = max(tem, pr0)
1062 prn(i,k) = tem
1063 endif
1064 prn(i,k) = min(prn(i,k),prmax)
1065 prn(i,k) = max(prn(i,k),prmin)
1066!
1067 ckz(i,k) = ck0 + (ck1 - ck0) * tem2
1068 ckz(i,k) = max(min(ckz(i,k), ck0), ck1)
1069 chz(i,k) = ch0 + (ch1 - ch0) * tem2
1070 chz(i,k) = max(min(chz(i,k), ch0), ch1)
1071!
1072 endif
1073 enddo
1074 enddo
1075!
1076! Above a threshold height (hcrinv), the background vertical
1077! diffusivities & mixing length
1078! in the inversion layers are set to much smaller values (xkinv1 &
1079! rlmn1)
1080!
1081! Below the threshold height (hcrinv), the background vertical
1082! diffusivities & mixing length
1083! in the inversion layers are increased with increasing roughness
1084! length & vegetation fraction
1085!
1086 do k = 1,km1
1087 do i=1,im
1088 if(zi(i,k+1) > hcrinv) then
1089 tem1 = tvx(i,k+1)-tvx(i,k)
1090 if(tem1 >= 0.) then
1091 xkzo(i,k) = min(xkzo(i,k), xkinv1)
1092 xkzmo(i,k) = min(xkzmo(i,k), xkinv1)
1093 rlmnz(i,k) = min(rlmnz(i,k), rlmn1)
1094 endif
1095 else
1096 tem1 = tvx(i,k+1)-tvx(i,k)
1097 if(tem1 > 0.) then
1098 ptem = xkzo(i,k) * zvfun(i)
1099 xkzo(i,k) = min(max(ptem, xkinv2), xkzo(i,k))
1100 ptem = xkzmo(i,k) * zvfun(i)
1101 xkzmo(i,k) = min(max(ptem, xkinv2), xkzmo(i,k))
1102 ptem = rlmnz(i,k) * zvfun(i)
1103 rlmnz(i,k) = min(max(ptem, rlmn2), rlmnz(i,k))
1104 endif
1105 endif
1106 enddo
1107 enddo
1108 do k = 2,km1
1109 do i=1,im
1110 rlmnz(i,k) = 0.5 * (rlmnz(i,k-1) + rlmnz(i,k))
1111 enddo
1112 enddo
1113!
1114!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1116!
1117 do k = 1, km1
1118 do i = 1, im
1119 zlup = 0.0
1120 mlenflg = .true.
1121 e2(i,k) = max(2.*tke(i,k), 0.001)
1122 do n = k, km1
1123 if(mlenflg) then
1124 dz = zl(i,n+1) - zl(i,n)
1125 tem1 = 2.*gotvx(i,n+1)*(thvx(i,k)-thvx(i,n+1))
1126 tem2 = cs0*sqrt(e2(i,n))*sqrt(shr2(i,n))
1127 e2(i,n+1) = e2(i,n) + (tem1 - tem2) * dz
1128 zlup = zlup + dz
1129 if(e2(i,n+1) < 0.) then
1130 ptem = e2(i,n+1) / (e2(i,n+1) - e2(i,n))
1131 zlup = zlup - ptem * dz
1132 zlup = max(zlup, 0.)
1133 mlenflg = .false.
1134 endif
1135 endif
1136 enddo
1137 zldn = 0.0
1138 mlenflg = .true.
1139 do n = k, 1, -1
1140 if(mlenflg) then
1141 if(n == 1) then
1142 dz = zl(i,1)
1143 tem = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
1144 tem1 = 2.*gotvx(i,n)*(tem-thvx(i,k))
1145 tem2 = ustar(i)*phims(i)/(vk*dz)
1146 tem2 = cs0*sqrt(e2(i,n))*tem2
1147 e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz
1148 else
1149 dz = zl(i,n) - zl(i,n-1)
1150 tem1 = 2.*gotvx(i,n-1)*(thvx(i,n-1)-thvx(i,k))
1151 tem2 = cs0*sqrt(e2(i,n))*sqrt(shr2(i,n-1))
1152 e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz
1153 endif
1154 zldn = zldn + dz
1155 if(e2(i,n-1) < 0.) then
1156 ptem = e2(i,n-1) / (e2(i,n-1) - e2(i,n))
1157 zldn = zldn - ptem * dz
1158 zldn = max(zldn, 0.)
1159 mlenflg = .false.
1160 endif
1161 endif
1162 enddo
1163!
1164 tem = 0.5 * (zi(i,k+1)-zi(i,k))
1165 tem1 = min(tem, rlmnz(i,k))
1178!
1179! Following Rodier et. al (2017), environmental wind shear effect on
1180! mixing length was included.
1181!
1182 ptem2 = min(zlup,zldn)
1183 rlam(i,k) = elmfac * ptem2
1184 rlam(i,k) = max(rlam(i,k), tem1)
1185 rlam(i,k) = min(rlam(i,k), rlmx)
1186!
1187 ptem2 = sqrt(zlup*zldn)
1188 ele(i,k) = elefac * ptem2
1189 ele(i,k) = max(ele(i,k), tem1)
1190 ele(i,k) = min(ele(i,k), elmx)
1191!
1192 enddo
1193 enddo
1196 do k = 1, km1
1197 do i = 1, im
1198 tem = vk * zl(i,k)
1199 if (zol(i) < 0.) then
1200 ptem = 1. - 100. * zol(i)
1201 ptem1 = ptem**0.2
1202 zk = tem * ptem1
1203 elseif (zol(i) >= 1.) then
1204 zk = tem / 3.7
1205 else
1206 ptem = 1. + 2.7 * zol(i)
1207 zk = tem / ptem
1208 endif
1209
1210 if (tc_pbl == 0) then
1211 elm(i,k) = zk*rlam(i,k)/(rlam(i,k)+zk)
1213 if ( sfc_rlm == 1 ) then
1214 if ( sfcflg(i) .and.
1215 & zl(i,k) < min(100.0,hpbl(i)*0.05) ) elm(i,k)=zk
1216 endif
1217 else if (tc_pbl == 1) then
1218 ! new blending method for mixing length
1219 elm(i,k) = sqrt( 1.0/( 1.0/(zk**2)+1.0/(rlam(i,k)**2) ) )
1220 endif
1221
1222!
1223 if(k == 1) elm(i,k)=zk
1224!
1225 dz = zi(i,k+1) - zi(i,k)
1226 tem = max(gdx(i),dz)
1227 elm(i,k) = min(elm(i,k), tem)
1228
1229 if (tc_pbl == 0) then
1230 ele(i,k) = min(ele(i,k), tem)
1231 else if (tc_pbl == 1) then
1232 ele(i,k) = elm(i,k)
1233 endif
1234!
1235 enddo
1236 enddo
1237 do i = 1, im
1238 elm(i,km) = elm(i,km1)
1239 ele(i,km) = ele(i,km1)
1240 enddo
1241!
1242!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1245!
1246 do k = 1, km1
1247 do i = 1, im
1248 tem = 0.5 * (elm(i,k) + elm(i,k+1))
1249 tem = tem * sqrt(tkeh(i,k))
1250 ri = max(bf(i,k)/shr2(i,k),rimin)
1251 if(k < kpbl(i)) then
1252 if(pcnvflg(i)) then
1253 dku(i,k) = ckz(i,k) * tem
1254 dkt(i,k) = dku(i,k) / prn(i,k)
1255 else
1256 if(ri < 0.) then ! unstable regime
1257 dku(i,k) = ckz(i,k) * tem
1258 dkt(i,k) = dku(i,k) / prn(i,k)
1259 else ! stable regime
1260 dkt(i,k) = chz(i,k) * tem
1261 dku(i,k) = dkt(i,k) * prn(i,k)
1262 endif
1263 endif
1264 else
1265 if(ri < 0.) then ! unstable regime
1266 dku(i,k) = ck1 * tem
1267 dkt(i,k) = rchck * dku(i,k)
1268 else ! stable regime
1269 dkt(i,k) = ch1 * tem
1270 prnum = 1.0 + 2.1 * ri
1271 prnum = min(prnum,prmax)
1272 dku(i,k) = dkt(i,k) * prnum
1273 endif
1274 endif
1275!
1276 if(scuflg(i)) then
1277 if(k >= mrad(i) .and. k < krad(i)) then
1278 tem1 = ckz(i,k) * tem
1279 ptem1 = tem1 / prscu
1280 dku(i,k) = max(dku(i,k), tem1)
1281 dkt(i,k) = max(dkt(i,k), ptem1)
1282 endif
1283 endif
1284!
1285 dkq(i,k) = prtke * dkt(i,k)
1286!
1287 dkt(i,k) = min(dkt(i,k),dkmax)
1288 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
1289 dkq(i,k) = min(dkq(i,k),dkmax)
1290 dkq(i,k) = max(dkq(i,k),xkzo(i,k))
1291 dku(i,k) = min(dku(i,k),dkmax)
1292 dku(i,k) = max(dku(i,k),xkzmo(i,k))
1293!
1294 enddo
1295 enddo
1298!
1299 do k = 1, km1
1300 do i = 1, im
1301 if(k == 1) then
1302 tem = ckz(i,1)
1303 tem1 = 0.5 * xkzmo(i,1)
1304 else
1305 tem = 0.5 * (ckz(i,k-1) + ckz(i,k))
1306 tem1 = 0.5 * (xkzmo(i,k-1) + xkzmo(i,k))
1307 endif
1308 ptem = tem1 / (tem * elm(i,k))
1309 tkmnz(i,k) = ptem * ptem
1310 tkmnz(i,k) = min(tkmnz(i,k), tkbmx)
1311 tkmnz(i,k) = max(tkmnz(i,k), tkmin)
1312 enddo
1313 enddo
1314!
1315!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1318!
1319 do k = 1, km1
1320 do i = 1, im
1321 if (k == 1) then
1322 tem = -dkt(i,1) * bf(i,1)
1323! if(pcnvflg(i)) then
1324! ptem1 = xmf(i,1) * buou(i,1)
1325! else
1326 ptem1 = 0.
1327! endif
1328 if(scuflg(i) .and. mrad(i) == 1) then
1329 ptem2 = xmfd(i,1) * buod(i,1)
1330 else
1331 ptem2 = 0.
1332 endif
1333 tem = tem + ptem1 + ptem2
1334 buop = 0.5 * (gotvx(i,1) * sflux(i) + tem)
1335!
1336 tem1 = dku(i,1) * shr2(i,1)
1337!
1338 tem = (u1(i,2)-u1(i,1))*rdzt(i,1)
1339! if(pcnvflg(i)) then
1340! ptem = xmf(i,1) * tem
1341! ptem1 = 0.5 * ptem * (u1(i,2)-ucko(i,2))
1342! else
1343 ptem1 = 0.
1344! endif
1345 if(scuflg(i) .and. mrad(i) == 1) then
1346 ptem = ucdo(i,1)+ucdo(i,2)-u1(i,1)-u1(i,2)
1347 ptem = 0.5 * tem * xmfd(i,1) * ptem
1348 else
1349 ptem = 0.
1350 endif
1351 ptem1 = ptem1 + ptem
1352!
1353 tem = (v1(i,2)-v1(i,1))*rdzt(i,1)
1354! if(pcnvflg(i)) then
1355! ptem = xmf(i,1) * tem
1356! ptem2 = 0.5 * ptem * (v1(i,2)-vcko(i,2))
1357! else
1358 ptem2 = 0.
1359! endif
1360 if(scuflg(i) .and. mrad(i) == 1) then
1361 ptem = vcdo(i,1)+vcdo(i,2)-v1(i,1)-v1(i,2)
1362 ptem = 0.5 * tem * xmfd(i,1) * ptem
1363 else
1364 ptem = 0.
1365 endif
1366 ptem2 = ptem2 + ptem
1367!
1368 tem2 = stress(i)*ustar(i)*phims(i)/(vk*zl(i,1))
1369 shrp = 0.5 * (tem1 + ptem1 + ptem2 + tem2)
1370 else
1371 tem1 = -dkt(i,k-1) * bf(i,k-1)
1372 tem2 = -dkt(i,k) * bf(i,k)
1373 tem = 0.5 * (tem1 + tem2)
1374 if(pcnvflg(i) .and. k <= kpbl(i)) then
1375 ptem = 0.5 * (xmf(i,k-1) + xmf(i,k))
1376 ptem1 = ptem * buou(i,k)
1377 else
1378 ptem1 = 0.
1379 endif
1380 if(scuflg(i)) then
1381 if(k >= mrad(i) .and. k < krad(i)) then
1382 ptem0 = 0.5 * (xmfd(i,k-1) + xmfd(i,k))
1383 ptem2 = ptem0 * buod(i,k)
1384 else
1385 ptem2 = 0.
1386 endif
1387 else
1388 ptem2 = 0.
1389 endif
1390 buop = tem + ptem1 + ptem2
1391!
1392 tem1 = dku(i,k-1) * shr2(i,k-1)
1393 tem2 = dku(i,k) * shr2(i,k)
1394 tem = 0.5 * (tem1 + tem2)
1395 tem1 = (u1(i,k+1)-u1(i,k))*rdzt(i,k)
1396 tem2 = (u1(i,k)-u1(i,k-1))*rdzt(i,k-1)
1397 if(pcnvflg(i) .and. k <= kpbl(i)) then
1398 ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2
1399 ptem1 = 0.5 * ptem * (u1(i,k)-ucko(i,k))
1400 else
1401 ptem1 = 0.
1402 endif
1403 if(scuflg(i)) then
1404 if(k >= mrad(i) .and. k < krad(i)) then
1405 ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2
1406 ptem2 = 0.5 * ptem0 * (ucdo(i,k)-u1(i,k))
1407 else
1408 ptem2 = 0.
1409 endif
1410 else
1411 ptem2 = 0.
1412 endif
1413 shrp = tem + ptem1 + ptem2
1414 tem1 = (v1(i,k+1)-v1(i,k))*rdzt(i,k)
1415 tem2 = (v1(i,k)-v1(i,k-1))*rdzt(i,k-1)
1416 if(pcnvflg(i) .and. k <= kpbl(i)) then
1417 ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2
1418 ptem1 = 0.5 * ptem * (v1(i,k)-vcko(i,k))
1419 else
1420 ptem1 = 0.
1421 endif
1422 if(scuflg(i)) then
1423 if(k >= mrad(i) .and. k < krad(i)) then
1424 ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2
1425 ptem2 = 0.5 * ptem0 * (vcdo(i,k)-v1(i,k))
1426 else
1427 ptem2 = 0.
1428 endif
1429 else
1430 ptem2 = 0.
1431 endif
1432 shrp = shrp + ptem1 + ptem2
1433 endif
1434 prod(i,k) = buop + shrp
1435 enddo
1436 enddo
1437!
1438!----------------------------------------------------------------------
1440!
1441 dtn = dt2 / float(ndt)
1442 do n = 1, ndt
1443 do k = 1,km1
1444 do i=1,im
1445 tem = sqrt(tke(i,k))
1446 ptem = ce0 / ele(i,k)
1447 diss(i,k) = ptem * tke(i,k) * tem
1448 tem1 = prod(i,k) + tke(i,k) / dtn
1449 diss(i,k)=max(min(diss(i,k), tem1), 0.)
1450 tke(i,k) = tke(i,k) + dtn * (prod(i,k)-diss(i,k))
1451! tke(i,k) = max(tke(i,k), tkmin)
1452 tke(i,k) = max(tke(i,k), tkmnz(i,k))
1453 enddo
1454 enddo
1455 enddo
1456!
1458!
1459 do k = 1, km
1460 do i = 1, im
1461 if(pcnvflg(i)) then
1462 qcko(i,k,ntke) = tke(i,k)
1463 endif
1464 if(scuflg(i)) then
1465 qcdo(i,k,ntke) = tke(i,k)
1466 endif
1467 enddo
1468 enddo
1469 do k = 2, kmpbl
1470 do i = 1, im
1471 if (pcnvflg(i) .and. k <= kpbl(i)) then
1472 dz = zl(i,k) - zl(i,k-1)
1473 tem = 0.5 * xlamue(i,k-1) * dz
1474 factor = 1. + tem
1475 qcko(i,k,ntke)=((1.-tem)*qcko(i,k-1,ntke)+tem*
1476 & (tke(i,k)+tke(i,k-1)))/factor
1477 endif
1478 enddo
1479 enddo
1480 do k = kmscu, 1, -1
1481 do i = 1, im
1482 if (scuflg(i) .and. k < krad(i)) then
1483 if(k >= mrad(i)) then
1484 dz = zl(i,k+1) - zl(i,k)
1485 tem = 0.5 * xlamde(i,k) * dz
1486 factor = 1. + tem
1487 qcdo(i,k,ntke)=((1.-tem)*qcdo(i,k+1,ntke)+tem*
1488 & (tke(i,k)+tke(i,k+1)))/factor
1489 endif
1490 endif
1491 enddo
1492 enddo
1493!
1494!--------------------------------------------------------
1495! compute variables for TVD flux-limiter scheme
1496! on environmental subsidence and uplifting
1497!
1498 kps = max(kmpbl, kmscu)
1499!
1500! for moisture and tracers including hydrometeors
1501!
1502 do n=1,ntrac1
1503 do k=1,kps
1504 do i=1,im
1505 qh(i,k,n) = 0.5 * (q1(i,k,n)+q1(i,k+1,n))
1506 enddo
1507 enddo
1508 do k=1,kps
1509 do i=1,im
1510 q_diff(i,k,n) = q1(i,k,n) - q1(i,k+1,n)
1511 enddo
1512 enddo
1513 do i=1,im
1514 if(q1(i,1,n) >= 0.) then
1515 q_diff(i,0,n) = max(0.,2.*q1(i,1,n)-q1(i,2,n))-
1516 & q1(i,1,n)
1517 else
1518 q_diff(i,0,n) = min(0.,2.*q1(i,1,n)-q1(i,2,n))-
1519 & q1(i,1,n)
1520 endif
1521 enddo
1522 enddo
1523!
1524 do n = 1, ntrac1
1525!
1526 do k = 1, kps
1527 do i = 1, im
1528 kmx = max(kpbl(i), krad(i))
1529 q_half(i,k,n) = qh(i,k,n)
1530 if((pcnvflg(i) .or. scuflg(i)) .and. k < kmx) then
1531 tem = 0.
1532 if(pcnvflg(i) .and. k < kpbl(i)) then
1533 tem = xmf(i,k)
1534 endif
1535 if(scuflg(i) .and.
1536 & (k >= mrad(i) .and. k < krad(i))) then
1537 tem = tem - xmfd(i,k)
1538 endif
1539 if(tem > 0.) then
1540 rrkp = 0.
1541 if(abs(q_diff(i,k,n)) > 1.e-22)
1542 & rrkp = q_diff(i,k+1,n) / q_diff(i,k,n)
1543 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
1544 q_half(i,k,n) = q1(i,k+1,n) +
1545 & phkp*(qh(i,k,n)-q1(i,k+1,n))
1546 elseif (tem < 0.) then
1547 rrkp = 0.
1548 if(abs(q_diff(i,k,n)) > 1.e-22)
1549 & rrkp = q_diff(i,k-1,n) / q_diff(i,k,n)
1550 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
1551 q_half(i,k,n) = q1(i,k,n) +
1552 & phkp*(qh(i,k,n)-q1(i,k,n))
1553 endif
1554 endif
1555 enddo
1556 enddo
1557!
1558 enddo
1559!
1560! for tke
1561!
1562 do k=1,kps
1563 do i=1,im
1564 tkeh(i,k) = 0.5 * (tke(i,k)+tke(i,k+1))
1565 enddo
1566 enddo
1567 do k=1,kps
1568 do i=1,im
1569 e_diff(i,k) = tke(i,k) - tke(i,k+1)
1570 enddo
1571 enddo
1572 do i=1,im
1573 if(tke(i,1) >= 0.) then
1574 e_diff(i,0) = max(0.,2.*tke(i,1)-tke(i,2))-
1575 & tke(i,1)
1576 else
1577 e_diff(i,0) = min(0.,2.*tke(i,1)-tke(i,2))-
1578 & tke(i,1)
1579 endif
1580 enddo
1581!
1582 do k = 1, kps
1583 do i = 1, im
1584 kmx = max(kpbl(i), krad(i))
1585 e_half(i,k) = tkeh(i,k)
1586 if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
1587 tem = 0.
1588 if(pcnvflg(i) .and. k < kpbl(i)) then
1589 tem = xmf(i,k)
1590 endif
1591 if(scuflg(i) .and.
1592 & (k >= mrad(i) .and. k < krad(i))) then
1593 tem = tem - xmfd(i,k)
1594 endif
1595 if(tem > 0.) then
1596 rrkp = 0.
1597 if(abs(e_diff(i,k)) > 1.e-22)
1598 & rrkp = e_diff(i,k+1) / e_diff(i,k)
1599 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
1600 e_half(i,k) = tke(i,k+1) +
1601 & phkp*(tkeh(i,k)-tke(i,k+1))
1602 elseif (tem < 0.) then
1603 rrkp = 0.
1604 if(abs(e_diff(i,k)) > 1.e-22)
1605 & rrkp = e_diff(i,k-1) / e_diff(i,k)
1606 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
1607 e_half(i,k) = tke(i,k) +
1608 & phkp*(tkeh(i,k)-tke(i,k))
1609 endif
1610 endif
1611 enddo
1612 enddo
1613!
1614!----------------------------------------------------------------------
1616!
1617 do i=1,im
1618 ad(i,1) = 1.0
1619 f1(i,1) = tke(i,1)
1620 enddo
1621!
1622 do k = 1,km1
1623 do i=1,im
1624 dtodsd = dt2/del(i,k)
1625 dtodsu = dt2/del(i,k+1)
1626 dsig = prsl(i,k)-prsl(i,k+1)
1627 rdz = rdzt(i,k)
1628 tem1 = dsig * dkq(i,k) * rdz
1629 dsdz2 = tem1 * rdz
1630 au(i,k) = -dtodsd*dsdz2
1631 al(i,k) = -dtodsu*dsdz2
1632 ad(i,k) = ad(i,k)-au(i,k)
1633 ad(i,k+1)= 1.-al(i,k)
1634 tem2 = dsig * rdz
1635!
1636 if(pcnvflg(i) .and. k < kpbl(i)) then
1637 ptem = 0.5 * tem2 * xmf(i,k)
1638 ptem1 = dtodsd * ptem
1639 ptem2 = dtodsu * ptem
1640 ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke)
1641 f1(i,k) = f1(i,k) - ptem * ptem1
1642 f1(i,k+1) = tke(i,k+1) + ptem * ptem2
1643 else
1644 f1(i,k+1) = tke(i,k+1)
1645 endif
1646!
1647 if(scuflg(i)) then
1648 if(k >= mrad(i) .and. k < krad(i)) then
1649 ptem = 0.5 * tem2 * xmfd(i,k)
1650 ptem1 = dtodsd * ptem
1651 ptem2 = dtodsu * ptem
1652 ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke)
1653 f1(i,k) = f1(i,k) + ptem * ptem1
1654 f1(i,k+1) = f1(i,k+1) - ptem * ptem2
1655 endif
1656 endif
1657!
1658 kmx = max(kpbl(i), krad(i))
1659 if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
1660 ptem = tem2 * (xmf(i,k) - xmfd(i,k))
1661 ptem1 = dtodsd * ptem
1662 ptem2 = dtodsu * ptem
1663 f1(i,k) = f1(i,k) + e_half(i,k) * ptem1
1664 f1(i,k+1) = f1(i,k+1) - e_half(i,k) * ptem2
1665 endif
1666!
1667 enddo
1668 enddo
1669c
1671c
1672 call tridit(im,km,1,al,ad,au,f1,au,f1)
1673!
1674! Negative TKE is set to zero after borrowing it from positive
1675! values within the mass-flux transport layers
1676!
1677 do i = 1,im
1678 tsumn(i) = 0.
1679 tsump(i) = 0.
1680 rtnp(i) = 1.
1681 enddo
1682 do k = 1,kps
1683 do i = 1,im
1684 if(pcnvflg(i) .and. scuflg(i)) then
1685 kbx = 1
1686 kmx = max(kpbl(i), krad(i))
1687 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
1688 kbx = 1
1689 kmx = kpbl(i)
1690 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
1691 kbx = mrad(i)
1692 kmx = krad(i)
1693 endif
1694 if((pcnvflg(i) .or. scuflg(i)) .and.
1695 & (k >= kbx .and. k <= kmx)) then
1696 tem = f1(i,k) * del(i,k) * gravi
1697 if(f1(i,k) < 0.) tsumn(i) = tsumn(i) + tem
1698 if(f1(i,k) > 0.) tsump(i) = tsump(i) + tem
1699 endif
1700 enddo
1701 enddo
1702 do i = 1,im
1703 if(pcnvflg(i) .or. scuflg(i)) then
1704 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
1705 if(tsump(i) > abs(tsumn(i))) then
1706 rtnp(i) = tsumn(i) / tsump(i)
1707 else
1708 rtnp(i) = tsump(i) / tsumn(i)
1709 endif
1710 endif
1711 endif
1712 enddo
1713 do k = 1,kps
1714 do i = 1,im
1715 if(pcnvflg(i) .and. scuflg(i)) then
1716 kbx = 1
1717 kmx = max(kpbl(i), krad(i))
1718 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
1719 kbx = 1
1720 kmx = kpbl(i)
1721 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
1722 kbx = mrad(i)
1723 kmx = krad(i)
1724 endif
1725 if((pcnvflg(i) .or. scuflg(i)) .and.
1726 & (k >= kbx .and. k <= kmx)) then
1727 if(rtnp(i) < 0.) then
1728 if(tsump(i) > abs(tsumn(i))) then
1729 if(f1(i,k) < 0.) f1(i,k) = 0.
1730 if(f1(i,k) > 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k)
1731 else
1732 if(f1(i,k) < 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k)
1733 if(f1(i,k) > 0.) f1(i,k) = 0.
1734 endif
1735 endif
1736 endif
1737 enddo
1738 enddo
1739!
1740! To remove negative TKEs which were leaked out of the mass-flux transport layers
1741! by eddy diffusion or potential negative TKEs from the diffusion scheme,
1742! positive TKEs are borrowed again now from the entire layers
1743!
1744 do i = 1,im
1745 tsumn(i) = 0.
1746 tsump(i) = 0.
1747 rtnp(i) = 1.
1748 enddo
1749 do k = 1,km
1750 do i = 1,im
1751 tem = f1(i,k) * del(i,k) * gravi
1752 if(f1(i,k) < 0.) tsumn(i) = tsumn(i) + tem
1753 if(f1(i,k) > 0.) tsump(i) = tsump(i) + tem
1754 enddo
1755 enddo
1756 do i = 1,im
1757 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
1758 if(tsump(i) > abs(tsumn(i))) then
1759 rtnp(i) = tsumn(i) / tsump(i)
1760 else
1761 rtnp(i) = tsump(i) / tsumn(i)
1762 endif
1763 endif
1764 enddo
1765 do k = 1,km
1766 do i = 1,im
1767 if(rtnp(i) < 0.) then
1768 if(tsump(i) > abs(tsumn(i))) then
1769 if(f1(i,k) < 0.) f1(i,k) = 0.
1770 if(f1(i,k) > 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k)
1771 else
1772 if(f1(i,k) < 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k)
1773 if(f1(i,k) > 0.) f1(i,k) = 0.
1774 endif
1775 endif
1776 enddo
1777 enddo
1778c
1780c
1781 do k = 1,km
1782 do i = 1,im
1783! f1(i,k) = max(f1(i,k), tkmin)
1784 qtend = (f1(i,k)-q1(i,k,ntke))*rdt
1785 rtg(i,k,ntke) = rtg(i,k,ntke)+qtend
1786 enddo
1787 enddo
1788 if(ldiag3d) then
1789 idtend = dtidx(ntke+100,index_of_process_pbl)
1790 if(idtend>0) then
1791 dtend(1:im,1:km,idtend) = dtend(1:im,1:km,idtend) + &
1792 & (f1(1:im,1:km)-q1(1:im,1:km,ntke))*rdt
1793 endif
1794 endif
1795c
1797c
1798 do i=1,im
1799 ad(i,1) = 1.
1800 f1(i,1) = t1(i,1) + dtdz1(i) * heat(i)
1801 f2(i,1) = q1(i,1,1) + dtdz1(i) * evap(i)
1802 enddo
1803 if(ntrac1 >= 2) then
1804 do n = 2, ntrac1
1805 is = (n-1) * km
1806 do i = 1, im
1807 f2(i,1+is) = q1(i,1,n)
1808 enddo
1809 enddo
1810 endif
1811c
1812 do k = 1,km1
1813 do i = 1,im
1814 dtodsd = dt2/del(i,k)
1815 dtodsu = dt2/del(i,k+1)
1816 dsig = prsl(i,k)-prsl(i,k+1)
1817 rdz = rdzt(i,k)
1818 tem1 = dsig * dkt(i,k) * rdz
1819 dsdzt = tem1 * gocp
1820 dsdz2 = tem1 * rdz
1821 au(i,k) = -dtodsd*dsdz2
1822 al(i,k) = -dtodsu*dsdz2
1823 ad(i,k) = ad(i,k)-au(i,k)
1824 ad(i,k+1)= 1.-al(i,k)
1825 tem2 = dsig * rdz
1826!
1827 if(pcnvflg(i) .and. k < kpbl(i)) then
1828 ptem = 0.5 * tem2 * xmf(i,k)
1829 ptem1 = dtodsd * ptem
1830 ptem2 = dtodsu * ptem
1831 tem = t1(i,k) + t1(i,k+1)
1832 ptem = tcko(i,k) + tcko(i,k+1)
1833 f1(i,k) = f1(i,k)+dtodsd*dsdzt-(ptem-tem)*ptem1
1834 f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+(ptem-tem)*ptem2
1835 ptem = qcko(i,k,1) + qcko(i,k+1,1)
1836 f2(i,k) = f2(i,k) - ptem * ptem1
1837 f2(i,k+1) = q1(i,k+1,1) + ptem * ptem2
1838 else
1839 f1(i,k) = f1(i,k)+dtodsd*dsdzt
1840 f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
1841 f2(i,k+1) = q1(i,k+1,1)
1842 endif
1843!
1844 if(scuflg(i)) then
1845 if(k >= mrad(i) .and. k < krad(i)) then
1846 ptem = 0.5 * tem2 * xmfd(i,k)
1847 ptem1 = dtodsd * ptem
1848 ptem2 = dtodsu * ptem
1849 ptem = tcdo(i,k) + tcdo(i,k+1)
1850 tem = t1(i,k) + t1(i,k+1)
1851 f1(i,k) = f1(i,k) + (ptem - tem) * ptem1
1852 f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2
1853 ptem = qcdo(i,k,1) + qcdo(i,k+1,1)
1854 f2(i,k) = f2(i,k) + ptem * ptem1
1855 f2(i,k+1) = f2(i,k+1) - ptem * ptem2
1856 endif
1857 endif
1858!
1859 kmx = max(kpbl(i), krad(i))
1860 if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
1861 ptem = tem2 * (xmf(i,k) - xmfd(i,k))
1862 ptem1 = dtodsd * ptem
1863 ptem2 = dtodsu * ptem
1864 f2(i,k) = f2(i,k) + q_half(i,k,1) * ptem1
1865 f2(i,k+1) = f2(i,k+1) - q_half(i,k,1) * ptem2
1866 endif
1867!
1868 enddo
1869 enddo
1870!
1871 if(ntrac1 >= 2) then
1872 do n = 2, ntrac1
1873 is = (n-1) * km
1874 do k = 1, km1
1875 do i = 1, im
1876 dtodsd = dt2/del(i,k)
1877 dtodsu = dt2/del(i,k+1)
1878 dsig = prsl(i,k)-prsl(i,k+1)
1879 tem2 = dsig * rdzt(i,k)
1880!
1881 if(pcnvflg(i) .and. k < kpbl(i)) then
1882 ptem = 0.5 * tem2 * xmf(i,k)
1883 ptem1 = dtodsd * ptem
1884 ptem2 = dtodsu * ptem
1885 ptem = qcko(i,k,n) + qcko(i,k+1,n)
1886 f2(i,k+is) = f2(i,k+is) - ptem * ptem1
1887 f2(i,k+1+is)= q1(i,k+1,n) + ptem * ptem2
1888 else
1889 f2(i,k+1+is) = q1(i,k+1,n)
1890 endif
1891!
1892 if(scuflg(i)) then
1893 if(k >= mrad(i) .and. k < krad(i)) then
1894 ptem = 0.5 * tem2 * xmfd(i,k)
1895 ptem1 = dtodsd * ptem
1896 ptem2 = dtodsu * ptem
1897 ptem = qcdo(i,k,n) + qcdo(i,k+1,n)
1898 f2(i,k+is) = f2(i,k+is) + ptem * ptem1
1899 f2(i,k+1+is)= f2(i,k+1+is) - ptem * ptem2
1900 endif
1901 endif
1902!
1903 kmx = max(kpbl(i), krad(i))
1904 if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
1905 ptem = tem2 * (xmf(i,k) - xmfd(i,k))
1906 ptem1 = dtodsd * ptem
1907 ptem2 = dtodsu * ptem
1908 f2(i,k+is) = f2(i,k+is) + q_half(i,k,n) * ptem1
1909 f2(i,k+1+is) = f2(i,k+1+is) - q_half(i,k,n) * ptem2
1910 endif
1911!
1912 enddo
1913 enddo
1914 enddo
1915 endif
1916c
1918c
1919 call tridin(im,km,ntrac1,al,ad,au,f1,f2,au,f1,f2)
1920!
1921! Negative moisture is set to zero after borrowing it from
1922! positive values within the mass-flux transport layers
1923!
1924 do i = 1,im
1925 tsumn(i) = 0.
1926 tsump(i) = 0.
1927 rtnp(i) = 1.
1928 enddo
1929 do k = 1,kps
1930 do i = 1,im
1931 if(pcnvflg(i) .and. scuflg(i)) then
1932 kbx = 1
1933 kmx = max(kpbl(i), krad(i))
1934 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
1935 kbx = 1
1936 kmx = kpbl(i)
1937 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
1938 kbx = mrad(i)
1939 kmx = krad(i)
1940 endif
1941 if((pcnvflg(i) .or. scuflg(i)) .and.
1942 & (k >= kbx .and. k <= kmx)) then
1943 tem = f2(i,k) * del(i,k) * gravi
1944 if(f2(i,k) < 0.) tsumn(i) = tsumn(i) + tem
1945 if(f2(i,k) > 0.) tsump(i) = tsump(i) + tem
1946 endif
1947 enddo
1948 enddo
1949 do i = 1,im
1950 if(pcnvflg(i) .or. scuflg(i)) then
1951 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
1952 if(tsump(i) > abs(tsumn(i))) then
1953 rtnp(i) = tsumn(i) / tsump(i)
1954 else
1955 rtnp(i) = tsump(i) / tsumn(i)
1956 endif
1957 endif
1958 endif
1959 enddo
1960 do k = 1,kps
1961 do i = 1,im
1962 if(pcnvflg(i) .and. scuflg(i)) then
1963 kbx = 1
1964 kmx = max(kpbl(i), krad(i))
1965 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
1966 kbx = 1
1967 kmx = kpbl(i)
1968 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
1969 kbx = mrad(i)
1970 kmx = krad(i)
1971 endif
1972 if((pcnvflg(i) .or. scuflg(i)) .and.
1973 & (k >= kbx .and. k <= kmx)) then
1974 if(rtnp(i) < 0.) then
1975 if(tsump(i) > abs(tsumn(i))) then
1976 if(f2(i,k) < 0.) f2(i,k) = 0.
1977 if(f2(i,k) > 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k)
1978 else
1979 if(f2(i,k) < 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k)
1980 if(f2(i,k) > 0.) f2(i,k) = 0.
1981 endif
1982 endif
1983 endif
1984 enddo
1985 enddo
1986!
1987! To remove negative moistures which were leaked out of the mass-flux transport layers
1988! by eddy diffusion or potential negative moistures from the diffusion scheme
1989! especially due to downward surface latent heat flux during nighttime,
1990! positive moistures are borrowed again now from the entire layers
1991!
1992 do i = 1,im
1993 tsumn(i) = 0.
1994 tsump(i) = 0.
1995 rtnp(i) = 1.
1996 enddo
1997 do k = 1,km
1998 do i = 1,im
1999 tem = f2(i,k) * del(i,k) * gravi
2000 if(f2(i,k) < 0.) tsumn(i) = tsumn(i) + tem
2001 if(f2(i,k) > 0.) tsump(i) = tsump(i) + tem
2002 enddo
2003 enddo
2004 do i = 1,im
2005 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2006 if(tsump(i) > abs(tsumn(i))) then
2007 rtnp(i) = tsumn(i) / tsump(i)
2008 else
2009 rtnp(i) = tsump(i) / tsumn(i)
2010 endif
2011 endif
2012 enddo
2013 do k = 1,km
2014 do i = 1,im
2015 if(rtnp(i) < 0.) then
2016 if(tsump(i) > abs(tsumn(i))) then
2017 if(f2(i,k) < 0.) f2(i,k) = 0.
2018 if(f2(i,k) > 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k)
2019 else
2020 if(f2(i,k) < 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k)
2021 if(f2(i,k) > 0.) f2(i,k) = 0.
2022 endif
2023 endif
2024 enddo
2025 enddo
2026!
2027! Negative hydrometeors & tracers are set to zero after
2028! borrowing them from positive values within the mass-flux
2029! transport layers
2030!
2031! For the negative liquid water, first borrow water from vapor
2032! and then borrow it from the other layers if there is still
2033! negative water
2034!
2035 if(ntrac1 >= 2) then
2036 is = (ntcw-1) * km
2037 do k = 1,kps
2038 do i = 1,im
2039 if(pcnvflg(i) .and. scuflg(i)) then
2040 kbx = 1
2041 kmx = max(kpbl(i), krad(i))
2042 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2043 kbx = 1
2044 kmx = kpbl(i)
2045 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2046 kbx = mrad(i)
2047 kmx = krad(i)
2048 endif
2049 if((pcnvflg(i) .or. scuflg(i)) .and.
2050 & (k >= kbx .and. k <= kmx)) then
2051 if(f2(i,k+is) < 0.) then
2052 tem = f2(i,k) + f2(i,k+is)
2053 if(tem >= 0.0) then
2054 f2(i,k) = tem
2055 f1(i,k) = f1(i,k) - elocp * f2(i,k+is)
2056 f2(i,k+is) = 0.
2057 elseif (f2(i,k) > 0.0) then
2058 f2(i,k+is) = tem
2059 f1(i,k) = f1(i,k) + elocp * f2(i,k)
2060 f2(i,k) = 0.
2061 endif
2062 endif
2063 endif
2064 enddo
2065 enddo
2066 endif
2067!
2068! For the negative rain water, first borrow water from vapor
2069! and then borrow it from the other layers if there is still
2070! negative water
2071!
2072 if(ntrac1 >= 2 .and. ntrw > 0) then
2073 is = (ntrw-1) * km
2074 do k = 1,kps
2075 do i = 1,im
2076 if(pcnvflg(i) .and. scuflg(i)) then
2077 kbx = 1
2078 kmx = max(kpbl(i), krad(i))
2079 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2080 kbx = 1
2081 kmx = kpbl(i)
2082 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2083 kbx = mrad(i)
2084 kmx = krad(i)
2085 endif
2086 if((pcnvflg(i) .or. scuflg(i)) .and.
2087 & (k >= kbx .and. k <= kmx)) then
2088 if(f2(i,k+is) < 0.) then
2089 tem = f2(i,k) + f2(i,k+is)
2090 if(tem >= 0.0) then
2091 f2(i,k) = tem
2092 f1(i,k) = f1(i,k) - elocp * f2(i,k+is)
2093 f2(i,k+is) = 0.
2094 elseif (f2(i,k) > 0.0) then
2095 f2(i,k+is) = tem
2096 f1(i,k) = f1(i,k) + elocp * f2(i,k)
2097 f2(i,k) = 0.
2098 endif
2099 endif
2100 endif
2101 enddo
2102 enddo
2103 endif
2104!
2105 if(ntrac1 >= 2) then
2106 do n = 2, ntrac1
2107 is = (n-1) * km
2108!
2109 do i = 1,im
2110 tsumn(i) = 0.
2111 tsump(i) = 0.
2112 rtnp(i) = 1.
2113 enddo
2114 do k = 1,kps
2115 do i = 1,im
2116 if(pcnvflg(i) .and. scuflg(i)) then
2117 kbx = 1
2118 kmx = max(kpbl(i), krad(i))
2119 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2120 kbx = 1
2121 kmx = kpbl(i)
2122 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2123 kbx = mrad(i)
2124 kmx = krad(i)
2125 endif
2126 if((pcnvflg(i) .or. scuflg(i)) .and.
2127 & (k >= kbx .and. k <= kmx)) then
2128 tem = f2(i,k+is) * del(i,k) * gravi
2129 if(f2(i,k+is) < 0.) tsumn(i) = tsumn(i) + tem
2130 if(f2(i,k+is) > 0.) tsump(i) = tsump(i) + tem
2131 endif
2132 enddo
2133 enddo
2134 do i = 1,im
2135 if(pcnvflg(i) .or. scuflg(i)) then
2136 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2137 if(tsump(i) > abs(tsumn(i))) then
2138 rtnp(i) = tsumn(i) / tsump(i)
2139 else
2140 rtnp(i) = tsump(i) / tsumn(i)
2141 endif
2142 endif
2143 endif
2144 enddo
2145 do k = 1,kps
2146 do i = 1,im
2147 if(pcnvflg(i) .and. scuflg(i)) then
2148 kbx = 1
2149 kmx = max(kpbl(i), krad(i))
2150 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2151 kbx = 1
2152 kmx = kpbl(i)
2153 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2154 kbx = mrad(i)
2155 kmx = krad(i)
2156 endif
2157 if((pcnvflg(i) .or. scuflg(i)) .and.
2158 & (k >= kbx .and. k <= kmx)) then
2159 if(rtnp(i) < 0.) then
2160 if(tsump(i) > abs(tsumn(i))) then
2161 if(f2(i,k+is)<0.) f2(i,k+is)=0.
2162 if(f2(i,k+is)>0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is)
2163 else
2164 if(f2(i,k+is)<0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is)
2165 if(f2(i,k+is)>0.) f2(i,k+is)=0.
2166 endif
2167 endif
2168 endif
2169 enddo
2170 enddo
2171!
2172! To remove negative hydrometeors & tracers which were leaked out of the mass-flux transport layers
2173! by eddy diffusion or potential negative hydrometeors & tracers from the diffusion scheme
2174! especially due to downward surface fluxes during nighttime,
2175! positive hydrometeors & tracers are borrowed again now from the entire layers
2176!
2177 do i = 1,im
2178 tsumn(i) = 0.
2179 tsump(i) = 0.
2180 rtnp(i) = 1.
2181 enddo
2182 do k = 1,km
2183 do i = 1,im
2184 tem = f2(i,k+is) * del(i,k) * gravi
2185 if(f2(i,k+is) < 0.) tsumn(i) = tsumn(i) + tem
2186 if(f2(i,k+is) > 0.) tsump(i) = tsump(i) + tem
2187 enddo
2188 enddo
2189 do i = 1,im
2190 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2191 if(tsump(i) > abs(tsumn(i))) then
2192 rtnp(i) = tsumn(i) / tsump(i)
2193 else
2194 rtnp(i) = tsump(i) / tsumn(i)
2195 endif
2196 endif
2197 enddo
2198 do k = 1,km
2199 do i = 1,im
2200 if(rtnp(i) < 0.) then
2201 if(tsump(i) > abs(tsumn(i))) then
2202 if(f2(i,k+is)<0.) f2(i,k+is)=0.
2203 if(f2(i,k+is)>0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is)
2204 else
2205 if(f2(i,k+is)<0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is)
2206 if(f2(i,k+is)>0.) f2(i,k+is)=0.
2207 endif
2208 endif
2209 enddo
2210 enddo
2211!
2212 enddo
2213 endif
2214c
2216c
2217 do k = 1,km
2218 do i = 1,im
2219 ttend = (f1(i,k)-t1(i,k))*rdt
2220 qtend = (f2(i,k)-q1(i,k,1))*rdt
2221 tdt(i,k) = tdt(i,k)+ttend
2222 rtg(i,k,1) = rtg(i,k,1)+qtend
2223! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend
2224! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend
2225 enddo
2226 enddo
2227!
2228 do i = 1,im
2229 dtsfc(i) = rho_a(i) * cp * heat(i)
2230 dqsfc(i) = rho_a(i) * hvap * evap(i)
2231 enddo
2232!
2233 if(ldiag3d .and. .not. gen_tend) then
2234 idtend = dtidx(index_of_temperature,index_of_process_pbl)
2235 if(idtend>=1) then
2236 do k = 1,km
2237 do i = 1,im
2238 ttend = (f1(i,k)-t1(i,k))*rdt
2239 dtend(i,k,idtend) = dtend(i,k,idtend)+ttend*delt
2240 enddo
2241 enddo
2242 endif
2243 ! Send tendencies just for QV; other tracers are below.
2244 idtend = dtidx(100+ntqv,index_of_process_pbl)
2245 if(idtend>=1) then
2246 do k = 1,km
2247 do i = 1,im
2248 qtend = (f2(i,k)-q1(i,k,1))*rdt
2249 dtend(i,k,idtend) = dtend(i,k,idtend)+qtend*delt
2250 enddo
2251 enddo
2252 endif
2253 endif
2254!
2255 if(ntrac1 >= 2) then
2256 do n = 2, ntrac1
2257 is = (n-1) * km
2258 do k = 1, km
2259 do i = 1, im
2260 qtend = (f2(i,k+is)-q1(i,k,n))*rdt
2261 rtg(i,k,n) = rtg(i,k,n)+qtend
2262 enddo
2263 enddo
2264 enddo
2265 if(ldiag3d .and. .not. gen_tend) then
2266 ! Send tendencies for all tracers that were selected.
2267 do n = 2, ntrac1
2268 is = (n-1) * km
2269 idtend = dtidx(n+100,index_of_process_pbl)
2270 if(idtend>=1) then
2271 if(n/=ntke) then
2272 do k = 1, km
2273 do i = 1, im
2274 qtend = (f2(i,k+is)-q1(i,k,n))*rdt
2275 dtend(i,k,idtend) = dtend(i,k,idtend)+qtend*delt
2276 enddo
2277 enddo
2278 endif
2279 endif
2280 enddo
2281 endif
2282 endif
2283!
2285!
2286 if(dspheat) then
2287 do k = 1,km1
2288 do i = 1,im
2289! tem = min(diss(i,k), dspmax)
2290! ttend = tem / cp
2291 ttend = diss(i,k) / cp
2292 tdt(i,k) = tdt(i,k) + dspfac * ttend
2293 enddo
2294 enddo
2295 if(ldiag3d .and. .not. gen_tend) then
2296 idtend = dtidx(index_of_temperature,index_of_process_pbl)
2297 if(idtend>=1) then
2298 do k = 1,km1
2299 do i = 1,im
2300 ttend = diss(i,k) / cp
2301 dtend(i,k,idtend) = dtend(i,k,idtend)+dspfac*ttend*delt
2302 enddo
2303 enddo
2304 endif
2305 endif
2306 endif
2307c
2309c
2310 do i=1,im
2311 ad(i,1) = 1.0 + dtdz1(i) * stress(i) / spd1(i)
2312 f1(i,1) = u1(i,1)
2313 f2(i,1) = v1(i,1)
2314 enddo
2315c
2316 do k = 1,km1
2317 do i=1,im
2318 dtodsd = dt2/del(i,k)
2319 dtodsu = dt2/del(i,k+1)
2320 dsig = prsl(i,k)-prsl(i,k+1)
2321 rdz = rdzt(i,k)
2322 tem1 = dsig * dku(i,k) * rdz
2323 dsdz2 = tem1*rdz
2324 au(i,k) = -dtodsd*dsdz2
2325 al(i,k) = -dtodsu*dsdz2
2326 ad(i,k) = ad(i,k)-au(i,k)
2327 ad(i,k+1)= 1.-al(i,k)
2328 tem2 = dsig * rdz
2329!
2330 if(pcnvflg(i) .and. k < kpbl(i)) then
2331 ptem = 0.5 * tem2 * xmf(i,k)
2332 ptem1 = dtodsd * ptem
2333 ptem2 = dtodsu * ptem
2334 tem = u1(i,k) + u1(i,k+1)
2335 ptem = ucko(i,k) + ucko(i,k+1)
2336 f1(i,k) = f1(i,k) - (ptem - tem) * ptem1
2337 f1(i,k+1) = u1(i,k+1) + (ptem - tem) * ptem2
2338 tem = v1(i,k) + v1(i,k+1)
2339 ptem = vcko(i,k) + vcko(i,k+1)
2340 f2(i,k) = f2(i,k) - (ptem - tem) * ptem1
2341 f2(i,k+1) = v1(i,k+1) + (ptem - tem) * ptem2
2342 else
2343 f1(i,k+1) = u1(i,k+1)
2344 f2(i,k+1) = v1(i,k+1)
2345 endif
2346!
2347 if(scuflg(i)) then
2348 if(k >= mrad(i) .and. k < krad(i)) then
2349 ptem = 0.5 * tem2 * xmfd(i,k)
2350 ptem1 = dtodsd * ptem
2351 ptem2 = dtodsu * ptem
2352 tem = u1(i,k) + u1(i,k+1)
2353 ptem = ucdo(i,k) + ucdo(i,k+1)
2354 f1(i,k) = f1(i,k) + (ptem - tem) *ptem1
2355 f1(i,k+1) = f1(i,k+1) - (ptem - tem) *ptem2
2356 tem = v1(i,k) + v1(i,k+1)
2357 ptem = vcdo(i,k) + vcdo(i,k+1)
2358 f2(i,k) = f2(i,k) + (ptem - tem) * ptem1
2359 f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2
2360 endif
2361 endif
2362!
2363 enddo
2364 enddo
2365c
2367c
2368 call tridi2(im,km,al,ad,au,f1,f2,au,f1,f2)
2369c
2371c
2372 do k = 1,km
2373 do i = 1,im
2374 utend = (f1(i,k)-u1(i,k))*rdt
2375 vtend = (f2(i,k)-v1(i,k))*rdt
2376 du(i,k) = du(i,k)+utend
2377 dv(i,k) = dv(i,k)+vtend
2378! dusfc(i) = dusfc(i)+conw*del(i,k)*utend
2379! dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend
2380 enddo
2381 enddo
2382 do i = 1,im
2383 if(icplocn2atm == 0) then
2384 dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i)
2385 dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i)
2386 else if (icplocn2atm ==1) then
2387 spd1_m=sqrt( (u1(i,1)-usfco(i))**2+(v1(i,1)-vsfco(i))**2 )
2388 dusfc(i) = -1.*rho_a(i)*stress(i)*(u1(i,1)-usfco(i))/spd1_m
2389 dvsfc(i) = -1.*rho_a(i)*stress(i)*(v1(i,1)-vsfco(i))/spd1_m
2390 endif
2391 enddo
2392!
2393 if(ldiag3d .and. .not. gen_tend) then
2394 idtend = dtidx(index_of_x_wind,index_of_process_pbl)
2395 if(idtend>=1) then
2396 do k = 1,km
2397 do i = 1,im
2398 utend = (f1(i,k)-u1(i,k))*rdt
2399 dtend(i,k,idtend) = dtend(i,k,idtend) + utend*delt
2400 enddo
2401 enddo
2402 endif
2403
2404 idtend = dtidx(index_of_y_wind,index_of_process_pbl)
2405 if(idtend>=1) then
2406 do k = 1,km
2407 do i = 1,im
2408 vtend = (f2(i,k)-v1(i,k))*rdt
2409 dtend(i,k,idtend) = dtend(i,k,idtend) + vtend*delt
2410 enddo
2411 enddo
2412 endif
2413 endif
2414!
2415!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2417!
2418 do i = 1, im
2419 hpbl(i) = hpblx(i)
2420 kpbl(i) = kpblx(i)
2421 enddo
2422!
2423!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2424 return
2425 end subroutine satmedmfvdifq_run
2427 end module satmedmfvdifq
mfpbltq_mod::mfpbltq
subroutine mfpbltq(im, ix, km, kmpbl, ntcw, ntrac1, delt,
This subroutine computes mass flux and updraft parcel properties for thermals driven by surface heati...
Definition mfpbltq.f:17
satmedmfvdifq::satmedmfvdifq_run
subroutine satmedmfvdifq_run(im, km, ntrac, ntcw, ntrw, ntiw, ntke, grav, rd, cp, rv, hvap, hfus, fv, eps, epsm1, dv, du, tdt, rtg, u1, v1, t1, q1, usfco, vsfco, icplocn2atm, swh, hlw, xmu, garea, zvfun, sigmaf, psk, rbsoil, zorl, u10m, v10m, fm, fh, tsea, heat, evap, stress, spd1, kpbl, prsi, del, prsl, prslk, phii, phil, delt, dspheat, dusfc, dvsfc, dtsfc, dqsfc, hpbl, dkt, dku, kinver, xkzm_m, xkzm_h, xkzm_s, dspfac, bl_upfr, bl_dnfr, rlmx, elmx, sfc_rlm, tc_pbl, ntqv, dtend, dtidx, index_of_temperature, index_of_x_wind, index_of_y_wind, index_of_process_pbl, gen_tend, ldiag3d, errmsg, errflg)
Definition satmedmfvdifq.F:89
tridi_mod::tridit
subroutine tridit(l, n, nt, cl, cm, cu, rt, au, at)
This subroutine solves tridiagonal problem for TKE.
Definition tridi.f:167
mfscuq_mod::mfscuq
subroutine mfscuq(im, ix, km, kmscu, ntcw, ntrac1, delt,
This subroutine computes mass flux and downdraft parcel properties for stratocumulus-top-driven turbu...
Definition mfscuq.f:15
tridi_mod::tridin
subroutine tridin(l, n, nt, cl, cm, cu, r1, r2, au, a1, a2)
Routine to solve the tridiagonal system to calculate u- and v-momentum at ; part of two-part process ...
Definition tridi.f:98
satmedmfvdifq::satmedmfvdifq_init
subroutine satmedmfvdifq_init(satmedmf, isatmedmf, isatmedmf_vdifq, errmsg, errflg)
Definition satmedmfvdifq.F:35
tridi_mod::tridi2
subroutine tridi2(l, n, cl, cm, cu, r1, r2, au, a1, a2)
This subroutine ..
Definition tridi.f:53
mfpbltq_mod
This module contains the subroutine that calculates mass flux and updraft parcel properties for therm...
Definition mfpbltq.f:9
mfscuq_mod
This module contains the mass flux and downdraft parcel properties parameterization for stratocumulus...
Definition mfscuq.f:6
satmedmfvdifq
This file contains the CCPP-compliant SATMEDMF scheme (updated version) which computes subgrid vertic...
Definition satmedmfvdifq.F:6
tridi_mod
This module contains routine to compute tridiagonal matrix elements for TKE, heat,...
Definition tridi.f:6