← Back to Documentation Home
CCPP SciDoc v8.0.0  v8.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
satmedmfvdifq_can.F
1
2
6
7!! if(tte_edmf=.true.), the TKE-EDMF parameterization becomes
8!! TTE(total turbulent energy)-based moist (TTE-EDMF) parameterization
9!!
10 module satmedmfvdifq_can_mod
11
12 contains
13
50 subroutine satmedmfvdifq_can(im,km,ntrac,ntcw,ntrw,
51 & ntiw,ntke,grav,pi,rd,cp,rv,hvap,hfus,fv,eps,epsm1,
52!The following three variables are for SA-3D-TKE
53 & def_1,def_2,def_3,sa3dtke,dku3d_h,dku3d_e,
54 & dv,du,tdt,rtg,u1,v1,t1,q1,usfco,vsfco,use_oceanuv,
55 & swh,hlw,xmu,garea,zvfun,sigmaf,
56 & psk,rbsoil,zorl,u10m,v10m,fm,fh,
57 & tsea,heat,evap,stress,spd1,kpbl,
58 & prsi,del,prsl,prslk,phii,phil,delt,tte_edmf,
59 & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku,tkeh,
60 & dkt_can, dku_can,
61 & kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr,
62 & rlmx,elmx,sfc_rlm,tc_pbl,use_lpt,
63 & do_canopy, cplaqm,
64 & ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind,
65 & index_of_y_wind,index_of_process_pbl,gen_tend,ldiag3d,
66 & errmsg,errflg)
67!IVAI: aux arrays
68! & naux2d,naux3d,aux2d,aux3d)
69
70!
71 use machine , only : kind_phys
72 use funcphys , only : fpvs
73 use mfpbltq_mod
74 use tridi_mod
75 use mfscuq_mod
76!
77 !PCC_CANOPY_utilities
78 use canopy_utils_mod
79!
80 implicit none
81!
82!----------------------------------------------------------------------
83 integer, intent(in) :: im, km, ntrac, ntcw, ntrw, ntiw,
84 & ntke, ntqv
85 integer, intent(in) :: sfc_rlm
86 integer, intent(in) :: tc_pbl
87 integer, intent(in) :: use_lpt
88 integer, intent(in) :: kinver(:)
89 integer, intent(out) :: kpbl(:)
90 logical, intent(in) :: gen_tend,ldiag3d
91!
92 real(kind=kind_phys), intent(in) :: grav,pi,rd,cp,rv,hvap,hfus,fv,
93 & eps,epsm1
94 real(kind=kind_phys), intent(in) :: delt, xkzm_m, xkzm_h, xkzm_s
95 real(kind=kind_phys), intent(in) :: dspfac, bl_upfr, bl_dnfr
96 real(kind=kind_phys), intent(in) :: rlmx, elmx
97!PCC CANOPY------------------------------------
98 logical, intent(in) :: do_canopy, cplaqm
99 real(kind=kind_phys), optional, intent(in) ::
100 & dkt_can(:,:), dku_can(:,:)
101 !----------------------------------------------
102 real(kind=kind_phys), intent(inout) :: dv(:,:), du(:,:),
103 & tdt(:,:), rtg(:,:,:), tkeh(:,:)
104 real(kind=kind_phys), intent(in) ::
105 & u1(:,:), v1(:,:),
106 & usfco(:), vsfco(:),
107 & t1(:,:), q1(:,:,:),
108!The following two variables are for SA-3D-TKE
109 & def_1(:,:), def_2(:,:), def_3(:,:),
110 & swh(:,:), hlw(:,:),
111 & xmu(:), garea(:),
112 & zvfun(:), sigmaf(:),
113 & psk(:), rbsoil(:),
114 & zorl(:), tsea(:),
115 & u10m(:), v10m(:),
116 & fm(:), fh(:),
117 & evap(:), heat(:),
118 & stress(:), spd1(:),
119 & prsi(:,:), del(:,:),
120 & prsl(:,:), prslk(:,:),
121 & phii(:,:), phil(:,:)
122 real(kind=kind_phys), intent(inout), dimension(:,:,:), optional ::
123 & dtend
124 integer, intent(in) :: dtidx(:,:), index_of_temperature,
125 & index_of_x_wind, index_of_y_wind, index_of_process_pbl
126 logical, intent(in) :: use_oceanuv
127 real(kind=kind_phys), intent(out) ::
128 & dusfc(:), dvsfc(:),
129 & dtsfc(:), dqsfc(:),
130 & hpbl(:)
131 real(kind=kind_phys), intent(out) ::
132 & dkt(:,:), dku(:,:)
133
134!
135 logical, intent(in) :: sa3dtke !flag for SA-3D-TKE scheme
136!
137! flag for tke dissipative heating
138 logical, intent(in) :: dspheat
139! flag for TTE-EDMF scheme
140 logical, intent(in) :: tte_edmf
141!
142 character(len=*), intent(out) :: errmsg
143 integer, intent(out) :: errflg
144
145!For passing dku to the dyn_core (SA-3D-TKE scheme)
146 real(kind=kind_phys), intent(out) ::
147 & dku3d_h(:,:),dku3d_e(:,:)
148
149!
150!----------------------------------------------------------------------
151!***
152!*** local variables
153 real(kind=kind_phys) spd1_m
154!***
155 integer i,is,k,n,ndt,km1,kmpbl,kmscu,ntrac1,idtend
156 integer kps,kbx,kmx
157 integer lcld(im),kcld(im),krad(im),mrad(im)
158 integer kx1(im), kb1(im), kpblx(im)
159!
160 real(kind=kind_phys) te(im,km), tei(im,km-1), tke(im,km),
161 & tteh(im,km), tesq(im,km-1),e2(im,0:km)
162!
163 real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km),
164 & qlx(im,km), thetae(im,km),thlx(im,km),
165 & slx(im,km), svx(im,km), qtx(im,km),
166 & tvx(im,km), pix(im,km), radx(im,km-1),
167 & dkq(im,km-1)
168!
169 real(kind=kind_phys) plyr(im,km), rhly(im,km), cfly(im,km),
170 & qstl(im,km)
171!
172 real(kind=kind_phys) dtdz1(im), gdx(im),
173 & phih(im), phim(im), phihs(im),
174 & phims(im), prn(im,km-1),
175 & rbdn(im), rbup(im), thermal(im),
176 & ustar(im), wstar(im), hpblx(im),
177 & ust3(im), wst3(im), rho_a(im),
178 & z0(im), crb(im), tkemean(im),
179 & hgamt(im), hgamq(im),
180 & wscale(im),vpert(im), thvs(im),
181 & zol(im), sflux(im), ris(im),
182 & sumx(im), tx1(im), tx2(im)
183!
184 real(kind=kind_phys) radmin(im)
185!
186 real(kind=kind_phys) zi(im,km+1), zl(im,km), zm(im,km),
187 & xkzo(im,km-1),xkzmo(im,km-1),
188 & xkzm_hx(im), xkzm_mx(im), tkmnz(im,km-1),
189 & rdzt(im,km-1),rlmnz(im,km),
190 & al(im,km-1), ad(im,km), au(im,km-1),
191 & f1(im,km), f2(im,km*(ntrac-1))
192!
193 real(kind=kind_phys) elm(im,km), ele(im,km),
194 & ckz(im,km), chz(im,km),
195 & diss(im,km-1),prod(im,km-1),
196 & bf(im,km-1), shr2(im,km-1), wush(im,km),
197 & xlamue(im,km-1), xlamde(im,km-1),
198 & gotvx(im,km), rlam(im,km-1)
199!
200! variables for updrafts (thermals)
201!
202 real(kind=kind_phys) tcko(im,km), qcko(im,km,ntrac),
203 & ucko(im,km), vcko(im,km),
204 & buou(im,km), xmf(im,km)
205!
206! variables for stratocumulus-top induced downdrafts
207!
208 real(kind=kind_phys) tcdo(im,km), qcdo(im,km,ntrac),
209 & ucdo(im,km), vcdo(im,km),
210 & buod(im,km), xmfd(im,km)
211!
212! variables for Total Variation Diminishing (TVD) flux-limiter scheme
213!
214 real(kind=kind_phys) e_half(im,km-1), e_diff(im,0:km-1),
215 & q_half(im,km-1,ntrac-1),
216 & qh(im,km-1,ntrac-1),
217 & q_diff(im,0:km-1,ntrac-1)
218 real(kind=kind_phys) rrkp, phkp
219 real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im)
220 real(kind=kind_phys) sfcpbl(im), vez0fun(im)
221!
222 logical pblflg(im), sfcflg(im), flg(im)
223 logical scuflg(im), pcnvflg(im)
224 logical mlenflg
225!
226! pcnvflg: true for unstable pbl
227!
228 real(kind=kind_phys) aphi16, aphi5,
229 & wfac, cfac,
230 & gamcrt, gamcrq, sfcfrac,
231! & conq, cont, conw,
232 & dsdz2, dsdzt, dkmax,
233 & dsig, dt2, dtodsd,
234 & dtodsu, g, factor, dz,
235 & gocp, gravi, zol1, zolcru,
236 & buop, shrp, dtn,
237 & prnum, prmax, prmin, prtke,
238 & prscu, pr0, ri,
239 & dw2, dw2min, zk,
240 & elmfac, elefac, dspmax,
241 & alp, clwt, cql,
242 & f0, robn, crbmin, crbmax,
243 & es, qs, value, onemrh,
244 & cfh, gamma, elocp, el2orc,
245 & epsi, beta, chx, cqx,
246 & rdt, rdz, qmin, qlmin,
247 & rimin, rbcr, rbint, tdzmin,
248 & rlmn, rlmn0, rlmn1, rlmn2,
249 & ttend, utend, vtend, qtend,
250 & zfac, zfmin, vk, spdk2,
251 & tkmin, tkbmx, disste, xkgdx,
252 & xkinv1, xkinv2,
253 & zlup, zldn, cs0, csmf,
254 & tem, tem1, tem2, tem3,
255 & ptem, ptem0, ptem1, ptem2
256!
257!The following variables are for SA-3D-TKE
258 integer kk
259 real(kind=kind_phys) thetal(im,km),dku_les(im,km),dkt_les(im,km),
260 & elmh(im,km),ele_les(im,km),pftke(im),
261 & dkq_les(im,km),pfl(im),pfdx(im),
262 & dku_h(im,km),dkq_h(im,km),
263 & elmhfac,elmhmx,ckh,elm_les,
264 & cpl1,cpl2,cpl3,cpl4,cpl5,cpl6,
265 & cptke1,cptke2,cptke3
266 integer ktkemax(im)
267 real(kind=kind_phys) tkemax(im),scl(im)
268 real(kind=kind_phys) sclmax,sclmin,dkmaxles
269! end of SA-3D-TKE variables
270!
271 real(kind=kind_phys) slfac
272!
273 real(kind=kind_phys) vegflo, vegfup, z0lo, z0up, vc0, zc0
274!
275 real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck
276!
277!
278 real(kind=kind_phys) epotte
279!
280 real(kind=kind_phys) qlcr, zstblmax, hcrinv
281!
282 real(kind=kind_phys) h1
283
284 real(kind=kind_phys) bfac, mffac
285
286 real(kind=kind_phys) qice(im,km),qliq(im,km)
287
288!PCC_CANOPY------------------------------------
289 integer kount !IVAI
290!----------------------------------------------
291
292!IVAI
293! integer, intent(in) :: naux2d,naux3d
294! real(kind_phys), intent(inout) :: aux2d(:,:)
295! real(kind_phys), intent(inout) :: aux3d(:,:,:)
296!IVAI
297
298!!
299 parameter(bfac=100.)
300 parameter(wfac=7.0)
301 parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1)
302 parameter(vk=0.4,rimin=-100.,slfac=0.1)
303 parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3)
304 parameter(rlmn=30.,rlmn0=5.,rlmn1=5.,rlmn2=10.)
305 parameter(prmin=0.25)
306 parameter(pr0=1.0,prtke=1.0)
307 parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35)
308 parameter(tkmin=1.e-9,tkbmx=0.2,dspmax=10.0)
309 parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8)
310 parameter(aphi5=5.,aphi16=16.)
311 parameter(elmfac=1.0,elefac=1.0,cql=100.)
312 parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=1000.)
313 parameter(qlcr=3.5e-5,zstblmax=2500.)
314 parameter(xkinv1=0.15,xkinv2=0.3)
315 parameter(h1=0.33333333,hcrinv=250.)
316 parameter(vegflo=0.1,vegfup=1.0,z0lo=0.1,z0up=1.0)
317 parameter(vc0=1.0,zc0=1.0)
318 parameter(cs0=0.4,csmf=0.5)
319 parameter(rchck=1.5,ndt=20)
320 !The following variables are for SA-3D-TKE
321 parameter(cpl1=0.280,cpl2=0.870,cpl3=0.913)
322 parameter(cpl4=0.153,cpl5=0.278,cpl6=0.720)
323 parameter(cptke1=0.07,cptke2=0.142,cptke3=0.071)
324 parameter(dkmaxles=300.0,sclmin=500.,sclmax=2500.)
325 parameter(elmhfac=1.5,elmhmx=1000.,ckh=0.4)
326!
327!PCC_CANOPY------------------------------------
328 if (do_canopy) then
329! Initialize layer thickness in non-canopy columns
330 rdzt(:,:) = 0.0 ! IVAI
331 endif
332!----------------------------------------------
333 if (tc_pbl == 0) then
334 ck0 = 0.4
335 ch0 = 0.4
336 ce0 = 0.4
337 else if (tc_pbl == 1) then
338 ck0 = 0.55
339 ch0 = 0.55
340 ce0 = 0.12
341 endif
342!
343 if(tte_edmf) then
344 cfac = 3.0
345 prmax = 6.0
346 prscu = 0.4
347 ck1 = 0.16
348 ch1 = 0.16
349 else
350 cfac = 4.5
351 prmax = 4.0
352 prscu = 0.67
353 ck1 = 0.15
354 ch1 = 0.15
355 endif
356!
357 gravi = 1.0 / grav
358 g = grav
359 gocp = g / cp
360! cont=cp/g
361! conq=hvap/g
362! conw=1.0/g ! for del in pa
363!! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) !kpa
364 elocp = hvap / cp
365 el2orc = hvap * hvap / (rv * cp)
366!
367!************************************************************************
368! Initialize CCPP error handling variables
369 errmsg = ''
370 errflg = 0
371
373 dt2 = delt
374 rdt = 1. / dt2
375!
376! the code is written assuming ntke=ntrac
377! if ntrac > ntke, the code needs to be modified
378!
379 ntrac1 = ntrac - 1
380 km1 = km - 1
381 kmpbl = km / 2
382 kmscu = km / 2
385 do k=1,km
386 do i=1,im
387 zi(i,k) = phii(i,k) * gravi
388 zl(i,k) = phil(i,k) * gravi
389 xmf(i,k) = 0.
390 xmfd(i,k) = 0.
391 buou(i,k) = 0.
392 buod(i,k) = 0.
393 wush(i,k) = 0.
394 ckz(i,k) = ck1
395 chz(i,k) = ch1
396 rlmnz(i,k) = rlmn0
397 enddo
398 enddo
399 do i=1,im
400 zi(i,km+1) = phii(i,km+1) * gravi
401 enddo
402 do k=1,km
403 do i=1,im
404 zm(i,k) = zi(i,k+1)
405 enddo
406 enddo
408 do i=1,im
409 gdx(i) = sqrt(garea(i))
410 enddo
413 do k=1,km
414 do i=1,im
415 te(i,k) = max(q1(i,k,ntke), tkmin)
416 tkeh(i,k) = 0
417 tteh(i,k) = 0
418 enddo
419 enddo
420 if(tte_edmf) then
421 do k=1,km1
422 do i=1,im
423 tteh(i,k) = 0.5 * (te(i,k) + te(i,k+1))
424 enddo
425 enddo
426 else
427 do k = 1, km
428 do i = 1, im
429 tke(i,k) = te(i,k)
430 enddo
431 enddo
432 do k=1,km1
433 do i=1,im
434 tkeh(i,k) = 0.5 * (tke(i,k) + tke(i,k+1))
435 enddo
436 enddo
437 endif
439 do k = 1,km1
440 do i=1,im
441 rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k))
442 prn(i,k) = pr0
443 enddo
444 enddo
445!
447
449
451
456!
457 do i=1,im
458 kx1(i) = 1
459 tx1(i) = 1.0 / prsi(i,1)
460 tx2(i) = tx1(i)
461 if(gdx(i) >= xkgdx) then
462 xkzm_hx(i) = xkzm_h
463 xkzm_mx(i) = xkzm_m
464 else
465 tem = gdx(i) / xkgdx
466 xkzm_hx(i) = xkzm_h * tem
467 xkzm_mx(i) = xkzm_m * tem
468 endif
469 enddo
470 do k = 1,km1
471 do i=1,im
472 xkzo(i,k) = 0.0
473 xkzmo(i,k) = 0.0
474 if (k < kinver(i)) then
475! minimum turbulent mixing length
476 ptem = prsi(i,k+1) * tx1(i)
477 tem1 = 1.0 - ptem
478 tem2 = tem1 * tem1 * 2.5
479 tem2 = min(1.0, exp(-tem2))
480 rlmnz(i,k)= rlmn * tem2
481 rlmnz(i,k)= max(rlmnz(i,k), rlmn0)
482! vertical background diffusivity
483 tem2 = tem1 * tem1 * 10.0
484 tem2 = min(1.0, exp(-tem2))
485 xkzo(i,k) = xkzm_hx(i) * tem2
486! vertical background diffusivity for
487! momentum
488 if (ptem >= xkzm_s) then
489 xkzmo(i,k) = xkzm_mx(i)
490 kx1(i) = k + 1
491 else
492 if (k == kx1(i) .and. k > 1) tx2(i) = 1.0 / prsi(i,k)
493 tem1 = 1.0 - prsi(i,k+1) * tx2(i)
494 tem1 = tem1 * tem1 * 5.0
495 xkzmo(i,k) = xkzm_mx(i) * min(1.0, exp(-tem1))
496 endif
497 endif
498 enddo
499 enddo
500!
502 do i = 1,im
503 z0(i) = 0.01 * zorl(i)
504 rho_a(i) = prsl(i,1)/(rd*t1(i,1)*(1.+fv*max(q1(i,1,1),qmin)))
505 dusfc(i) = 0.
506 dvsfc(i) = 0.
507 dtsfc(i) = 0.
508 dqsfc(i) = 0.
509 kpbl(i) = 1
510 hpbl(i) = 0.
511 kpblx(i) = 1
512 hpblx(i) = 0.
513 pfl(i)=1.0
514 pftke(i)=1.0
515 pblflg(i)= .true.
516 sfcflg(i)= .true.
517 if(rbsoil(i) > 0.) sfcflg(i) = .false.
518 pcnvflg(i)= .false.
519 scuflg(i)= .true.
520 if(scuflg(i)) then
521 radmin(i)= 0.
522 mrad(i) = km1
523 krad(i) = 1
524 lcld(i) = km1
525 kcld(i) = km1
526 endif
527 enddo
528!
532!
533 do i = 1,im
534 tem = (sigmaf(i) - vegflo) / (vegfup - vegflo)
535 tem = min(max(tem, 0.), 1.)
536 tem1 = sqrt(tem)
537 ptem = (z0(i) - z0lo) / (z0up - z0lo)
538 ptem = min(max(ptem, 0.), 1.)
539 vez0fun(i) = (1. + vc0 * tem1) * (1. + zc0 * ptem)
540 enddo
541!
544 do k=1,km
545 do i=1,im
546 pix(i,k) = psk(i) / prslk(i,k)
547 theta(i,k) = t1(i,k) * pix(i,k)
548 qice(i,k) = 0.0
549 qliq(i,k) = 0.0
550 if(ntiw > 0) then
551 tem = max(q1(i,k,ntcw),qlmin)
552 qliq(i,k) = tem
553 if(sa3dtke) then
554 tem1=max(q1(i,k,ntiw)+q1(i,k,5)+q1(i,k,6),qlmin) !for SA-3D-TKE
555 qice(i,k) = tem1
556 else
557 tem1=max(q1(i,k,ntiw),qlmin)
558 qice(i,k) = tem1
559 endif
560 qlx(i,k) = tem + tem1
561 ptem = hvap*tem + (hvap+hfus)*tem1
562 slx(i,k) = cp * t1(i,k) + phil(i,k) - ptem
563 else
564 qlx(i,k) = max(q1(i,k,ntcw),qlmin)
565 slx(i,k) = cp * t1(i,k) + phil(i,k) - hvap*qlx(i,k)
566 qliq(i,k) = qlx(i,k)
567 endif
568 tem2 = 1.+fv*max(q1(i,k,1),qmin)-qlx(i,k)
569 thvx(i,k) = theta(i,k) * tem2
570 tvx(i,k) = t1(i,k) * tem2
571 qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k)
572 thlx(i,k) = theta(i,k) - pix(i,k)*elocp*qlx(i,k)
573 thlvx(i,k) = thlx(i,k) * (1. + fv * qtx(i,k))
574 svx(i,k) = cp * tvx(i,k)
575 ptem1 = elocp * pix(i,k) * max(q1(i,k,1),qmin)
576 thetae(i,k)= theta(i,k) + ptem1
577! gotvx(i,k) = g / tvx(i,k)
578 gotvx(i,k) = g / thvx(i,k)
579 enddo
580 enddo
581!
584 do k = 1, km
585 do i = 1, im
586 plyr(i,k) = 0.01 * prsl(i,k) ! pa to mb (hpa)
587! --- ... compute relative humidity
588 es = 0.01 * fpvs(t1(i,k)) ! fpvs in pa
589 qs = max(qmin, eps * es / (plyr(i,k) + epsm1*es))
590 rhly(i,k) = max(0.0, min(1.0, max(qmin, q1(i,k,1))/qs))
591 qstl(i,k) = qs
592 enddo
593 enddo
594!
595 do k = 1, km
596 do i = 1, im
597 cfly(i,k) = 0.
598 clwt = 1.0e-6 * (plyr(i,k)*0.001)
599 if (qlx(i,k) > clwt) then
600 onemrh = max(1.e-10, 1.0-rhly(i,k))
601 tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0)
602 tem1 = cql / tem1
603 value = max(min( tem1*qlx(i,k), 50.0), 0.0)
604 tem2 = sqrt(sqrt(rhly(i,k)))
605 cfly(i,k) = min(max(tem2*(1.0-exp(-value)), 0.0), 1.0)
606 endif
607 enddo
608 enddo
609!
611!
612 do k = 1, km1
613 do i = 1, im
614 tem = 0.5 * (svx(i,k) + svx(i,k+1))
615 tem1 = 0.5 * (t1(i,k) + t1(i,k+1))
616 tem2 = 0.5 * (qstl(i,k) + qstl(i,k+1))
617 cfh = min(cfly(i,k+1),0.5*(cfly(i,k)+cfly(i,k+1)))
618 alp = g / tem
619 gamma = el2orc * tem2 / (tem1**2)
620 epsi = tem1 / elocp
621 beta = (1. + gamma*epsi*(1.+fv)) / (1. + gamma)
622 chx = cfh * alp * beta + (1. - cfh) * alp
623 cqx = cfh * alp * hvap * (beta - epsi)
624 cqx = cqx + (1. - cfh) * fv * g
625 ptem1 = (slx(i,k+1)-slx(i,k))*rdzt(i,k)
626 ptem2 = (qtx(i,k+1)-qtx(i,k))*rdzt(i,k)
627 bf(i,k) = chx * ptem1 + cqx * ptem2
628 enddo
629 enddo
630!
631!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
632!
635 do k=1,km
636 do i=1,im
637 dku(i,k) = 0.
638 dkt(i,k) = 0.
639 enddo
640 enddo
641 do k=1,km1
642 do i=1,im
643 dkq(i,k) = 0.
644 tem = zi(i,k+1)-zi(i,k)
645 radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k))
646 enddo
647 enddo
651 do i = 1,im
652 sflux(i) = heat(i) + evap(i)*fv*theta(i,1)
653 if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false.
654 enddo
655!
673 do i = 1,im
674 thvs(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
675 if(pblflg(i)) then
676! thermal(i) = thvx(i,1)
677 thermal(i) = thlvx(i,1)
678 crb(i) = rbcr
679 else
680 thermal(i) = thvs(i)
681 tem = sqrt(u10m(i)**2+v10m(i)**2)
682 tem = max(tem, 1.)
683 robn = tem / (f0 * z0(i))
684 tem1 = 1.e-7 * robn
685 crb(i) = 0.16 * (tem1 ** (-0.18))
686 crb(i) = max(min(crb(i), crbmax), crbmin)
687 endif
688 enddo
690 do i=1,im
691 dtdz1(i) = dt2 / (zi(i,2)-zi(i,1))
692 enddo
693!
694 do i=1,im
695 ustar(i) = sqrt(stress(i))
696 enddo
697!
700!
701 do k = 1, km1
702 do i = 1, im
703 rdz = rdzt(i,k)
704! bf(i,k) = gotvx(i,k)*(thvx(i,k+1)-thvx(i,k))*rdz
705 dw2 = (u1(i,k)-u1(i,k+1))**2
706 & + (v1(i,k)-v1(i,k+1))**2
707 shr2(i,k) = max(dw2,dw2min)*rdz*rdz
708 enddo
709 enddo
710!
711! Find first quess pbl height based on bulk richardson number (mrf pbl scheme)
712! and also for diagnostic purpose
713!
714 do i=1,im
715 flg(i) = .false.
716 rbup(i) = rbsoil(i)
717 enddo
723 do k = 1, kmpbl
724 do i = 1, im
725 if(.not.flg(i)) then
726 rbdn(i) = rbup(i)
727 if (tc_pbl == 0) then
728 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
729! rbup(i) = (thvx(i,k)-thermal(i))*
730! & (g*zl(i,k)/thvx(i,1))/spdk2
731 rbup(i) = (thlvx(i,k)-thermal(i))*
732 & (g*zl(i,k)/thlvx(i,1))/spdk2
733 else if (tc_pbl == 1) then
734 spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.)
735 & + bfac*ustar(i)**2
736 rbup(i) = (thlvx(i,k)-thermal(i))*
737 & (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2
738 endif
739 kpblx(i) = k
740 flg(i) = rbup(i) > crb(i)
741 endif
742 enddo
743 enddo
747 do i = 1,im
748 if(kpblx(i) > 1) then
749 k = kpblx(i)
750 if(rbdn(i) >= crb(i)) then
751 rbint = 0.
752 elseif(rbup(i) <= crb(i)) then
753 rbint = 1.
754 else
755 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
756 endif
757 hpblx(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
758 if(hpblx(i) < zi(i,kpblx(i))) kpblx(i)=kpblx(i)-1
759 else
760 hpblx(i) = zl(i,1)
761 kpblx(i) = 1
762 endif
763 hpbl(i) = hpblx(i)
764 kpbl(i) = kpblx(i)
765 if(kpbl(i) <= 1) pblflg(i)=.false.
766 enddo
767!
768! update thermal at a level of slfac*hpbl for unstable pbl
769!
770 do i=1,im
771 sfcpbl(i) = slfac * hpbl(i)
772 kb1(i) = 1
773 flg(i) = .false.
774 if(pblflg(i)) then
775 flg(i) = .true.
776 endif
777 enddo
778 do k = 2, kmpbl
779 do i=1,im
780 if (flg(i) .and. zl(i,k) <= sfcpbl(i)) then
781 kb1(i) = k
782 else
783 flg(i) = .false.
784 endif
785 enddo
786 enddo
787 do i=1,im
788 if(pblflg(i)) kb1(i)=min(kb1(i),kpbl(i))
789 enddo
790!
791! re-compute pbl height with the updated thermal
792!
793 do i=1,im
794 flg(i) = .true.
795 if(pblflg(i) .and. kb1(i) > 1) then
796 flg(i) = .false.
797 rbup(i) = rbsoil(i)
798! thermal(i) = thvx(i,kb1(i))
799 thermal(i) = thlvx(i,kb1(i))
800 kpblx(i) = kb1(i)
801 hpblx(i) = zl(i,kb1(i))
802 endif
803 enddo
804 do k = 2, kmpbl
805 do i = 1, im
806 if(.not.flg(i) .and. k > kb1(i)) then
807 rbdn(i) = rbup(i)
808 if (tc_pbl == 0) then
809 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
810! rbup(i) = (thvx(i,k)-thermal(i))*
811! & (g*zl(i,k)/thvx(i,1))/spdk2
812 rbup(i) = (thlvx(i,k)-thermal(i))*
813 & (g*zl(i,k)/thlvx(i,1))/spdk2
814 else if (tc_pbl == 1) then
815 spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.)
816 & + bfac*ustar(i)**2
817 rbup(i) = (thlvx(i,k)-thermal(i))*
818 & (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2
819 endif
820 kpblx(i) = k
821 flg(i) = rbup(i) > crb(i)
822 endif
823 enddo
824 enddo
825 do i = 1,im
826 if(pblflg(i) .and. kb1(i) > 1) then
827 k = kpblx(i)
828 if(rbdn(i) >= crb(i)) then
829 rbint = 0.
830 elseif(rbup(i) <= crb(i)) then
831 rbint = 1.
832 else
833 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
834 endif
835 hpblx(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
836 if(hpblx(i) < zi(i,kpblx(i))) kpblx(i)=kpblx(i)-1
837 hpbl(i) = hpblx(i)
838 kpbl(i) = kpblx(i)
839 endif
840 enddo
841!
842 if(.not.tte_edmf) then
843!
845!
846 do i = 1, im
847 sumx(i) = 0.
848 tkemean(i) = 0.
849 enddo
850 do k = 1, kmpbl
851 do i = 1, im
852 if(k < kpbl(i)) then
853 dz = zi(i,k+1) - zi(i,k)
854 tkemean(i) = tkemean(i) + tke(i,k) * dz
855 sumx(i) = sumx(i) + dz
856 endif
857 enddo
858 enddo
859 do i = 1, im
860 if(tkemean(i) > 0. .and. sumx(i) > 0.) then
861 tkemean(i) = tkemean(i) / sumx(i)
862 endif
863 enddo
864!
865 endif
866!
869!
870 kps = max(kmpbl, kmscu)
871 do k = 2, kps
872 do i = 1, im
873 dz = zi(i,k+1) - zi(i,k)
874 tem = (0.5*(u1(i,k-1)-u1(i,k+1))/dz)**2
875 tem1 = tem+(0.5*(v1(i,k-1)-v1(i,k+1))/dz)**2
876 wush(i,k) = csmf * sqrt(tem1)
877 enddo
878 enddo
879!
889 do i=1,im
890 zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
891 if(sfcflg(i)) then
892 zol(i) = min(zol(i),-zfmin)
893 else
894 zol(i) = max(zol(i),zfmin)
895 endif
907 zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1)
908 if(sfcflg(i)) then
909 tem = 1.0 / (1. - aphi16*zol1)
910 phih(i) = sqrt(tem)
911 phim(i) = sqrt(phih(i))
912 tem1 = 1.0 / (1. - aphi16*zol(i))
913 phihs(i) = sqrt(tem1)
914 phims(i) = sqrt(phihs(i))
915 else
916 phim(i) = 1. + aphi5*zol1
917 phih(i) = phim(i)
918 phims(i) = 1. + aphi5*zol(i)
919 phihs(i) = phims(i)
920 endif
921 enddo
922!
938 do i=1,im
939 if(pblflg(i)) then
940 if(zol(i) < zolcru) then
941 pcnvflg(i) = .true.
942 endif
943 wst3(i) = gotvx(i,1)*sflux(i)*hpbl(i)
944 wstar(i)= wst3(i)**h1
945 ust3(i) = ustar(i)**3.
946 wscale(i)=(ust3(i)+wfac*vk*wst3(i)*sfcfrac)**h1
947 ptem = ustar(i)/aphi5
948 wscale(i) = max(wscale(i),ptem)
949 endif
950 enddo
951!
954!
955 do i = 1,im
956 if(pcnvflg(i)) then
957 hgamt(i) = heat(i)/wscale(i)
958 hgamq(i) = evap(i)/wscale(i)
959 vpert(i) = hgamt(i) + hgamq(i)*fv*theta(i,1)
960 vpert(i) = max(vpert(i),0.)
961 tem = min(cfac*vpert(i),gamcrt)
962 thermal(i)= thermal(i) + tem
963 endif
964 enddo
965!
966! enhance the pbl height by considering the thermal excess
967! (overshoot pbl top)
968!
969 do i=1,im
970 flg(i) = .true.
971 if(pcnvflg(i)) then
972 flg(i) = .false.
973 rbup(i) = rbsoil(i)
974 endif
975 enddo
976 do k = 2, kmpbl
977 do i = 1, im
978 if(.not.flg(i) .and. k > kb1(i)) then
979 rbdn(i) = rbup(i)
980 if (tc_pbl == 0) then
981 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
982 rbup(i) = (thlvx(i,k)-thermal(i))*
983 & (g*zl(i,k)/thlvx(i,1))/spdk2
984 else if (tc_pbl == 1) then
985 spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.)
986 & + bfac*ustar(i)**2
987 rbup(i) = (thlvx(i,k)-thermal(i))*
988 & (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2
989 endif
990 kpbl(i) = k
991 flg(i) = rbup(i) > crb(i)
992 endif
993 enddo
994 enddo
995 do i = 1,im
996 if(pcnvflg(i)) then
997 k = kpbl(i)
998 if(rbdn(i) >= crb(i)) then
999 rbint = 0.
1000 elseif(rbup(i) <= crb(i)) then
1001 rbint = 1.
1002 else
1003 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
1004 endif
1005 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
1006 if(hpbl(i) < zi(i,kpbl(i))) then
1007 kpbl(i) = kpbl(i) - 1
1008 endif
1009 if(kpbl(i) <= 1) then
1010 pcnvflg(i) = .false.
1011 pblflg(i) = .false.
1012 endif
1013 endif
1014 enddo
1015!
1016!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1017! compute tke using tte & ri for TTE-EDMF
1018!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1019!
1020 if(tte_edmf) then
1021!
1022 do i = 1, im
1023 tem = phims(i) * phims(i)
1024 ris(i) = zol(i) * phihs(i) / tem
1025 ris(i) = max(ris(i), rimin)
1026 enddo
1027 do k = 1, km1
1028 do i = 1, im
1029 ptem = sfcfrac*hpbl(i)
1030 if (zl(i,k) <= ptem) then
1031 ri = ris(i)
1032 else
1033 if(k == 1) then
1034 tem = gotvx(i,1) * (thlvx(i,1)-thvs(i))
1035 tem1 = tem / zl(i,1)
1036 tem1 = 0.5 * (tem1 + bf(i,1))
1037 ptem = max((u1(i,1)**2+v1(i,1)**2), 1.)
1038 ptem1 = ptem / (zl(i,1) * zl(i,1))
1039 ptem1 = 0.5 * (ptem1 + shr2(i,1))
1040 ri = max(tem1/ptem1, rimin)
1041 else
1042 tem1 = 0.5 * (bf(i,k-1) + bf(i,k))
1043 ptem1 = 0.5 * (shr2(i,k-1) + shr2(i,k))
1044 ri = max(tem1/ptem1, rimin)
1045 endif
1046 endif
1047 if(ri < 0) then
1048 tem = 2. * ri - pr0
1049 epotte = ri / tem
1050 else
1051 tem = pr0 + 3. * ri
1052 epotte = ri / tem
1053 endif
1054 tke(i,k) = te(i,k) * (1. - epotte)
1055 enddo
1056 enddo
1057 do i=1,im
1058 tke(i,km) = tke(i,km1)
1059 enddo
1060!
1062!
1063 do i = 1, im
1064 sumx(i) = 0.
1065 tkemean(i) = 0.
1066 enddo
1067 do k = 1, kmpbl
1068 do i = 1, im
1069 if(k < kpbl(i)) then
1070 dz = zi(i,k+1) - zi(i,k)
1071 tkemean(i) = tkemean(i) + tke(i,k) * dz
1072 sumx(i) = sumx(i) + dz
1073 endif
1074 enddo
1075 enddo
1076 do i = 1, im
1077 if(tkemean(i) > 0. .and. sumx(i) > 0.) then
1078 tkemean(i) = tkemean(i) / sumx(i)
1079 endif
1080 enddo
1081!
1082 endif ! end of if(tte_edmf)
1083!
1084!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1085! look for stratocumulus
1090 do i=1,im
1091 flg(i) = scuflg(i)
1092 enddo
1093 do k = 1, km1
1094 do i=1,im
1095 if(flg(i).and.zl(i,k) >= zstblmax) then
1096 lcld(i)=k
1097 flg(i)=.false.
1098 endif
1099 enddo
1100 enddo
1101 do i = 1, im
1102 flg(i)=scuflg(i)
1103 enddo
1104 do k = kmscu,1,-1
1105 do i = 1, im
1106 if(flg(i) .and. k <= lcld(i)) then
1107 if(qlx(i,k) >= qlcr) then
1108 kcld(i)=k
1109 flg(i)=.false.
1110 endif
1111 endif
1112 enddo
1113 enddo
1114 do i = 1, im
1115 if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false.
1116 enddo
1122 do i = 1, im
1123 flg(i)=scuflg(i)
1124 enddo
1125 do k = kmscu,1,-1
1126 do i = 1, im
1127 if(flg(i) .and. k <= kcld(i)) then
1128 if(qlx(i,k) >= qlcr) then
1129 if(radx(i,k) < radmin(i)) then
1130 radmin(i)=radx(i,k)
1131 krad(i)=k
1132 endif
1133 else
1134 flg(i)=.false.
1135 endif
1136 endif
1137 enddo
1138 enddo
1139 do i = 1, im
1140 if(scuflg(i) .and. krad(i) <= 1) scuflg(i)=.false.
1141 if(scuflg(i) .and. radmin(i)>=0.) scuflg(i)=.false.
1142 enddo
1143!
1144!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1149 do k = 1, km
1150 do i = 1, im
1151 if(pcnvflg(i)) then
1152 tcko(i,k) = t1(i,k)
1153 ucko(i,k) = u1(i,k)
1154 vcko(i,k) = v1(i,k)
1155 endif
1156 if(scuflg(i)) then
1157 tcdo(i,k) = t1(i,k)
1158 ucdo(i,k) = u1(i,k)
1159 vcdo(i,k) = v1(i,k)
1160 endif
1161 enddo
1162 enddo
1163 do n = 1, ntrac1
1164 do k = 1, km
1165 do i = 1, im
1166 if(pcnvflg(i)) then
1167 qcko(i,k,n) = q1(i,k,n)
1168 endif
1169 if(scuflg(i)) then
1170 qcdo(i,k,n) = q1(i,k,n)
1171 endif
1172 enddo
1173 enddo
1174 enddo
1177 call mfpbltq(im,im,km,kmpbl,ntcw,ntrac1,dt2,
1178 & pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx,
1179 & gdx,hpbl,kpbl,vpert,buou,wush,tkemean,vez0fun,xmf,
1180 & tcko,qcko,ucko,vcko,xlamue,bl_upfr)
1183 call mfscuq(im,im,km,kmscu,ntcw,ntrac1,dt2,
1184 & scuflg,zl,zm,q1,t1,u1,v1,plyr,pix,
1185 & thlx,thvx,thlvx,gdx,thetae,
1186 & krad,mrad,radmin,buod,wush,tkemean,vez0fun,xmfd,
1187 & tcdo,qcdo,ucdo,vcdo,xlamde,bl_dnfr)
1188
1189 if (tc_pbl == 1) then
1191 do i=1,im
1192 if(zol(i) > -0.5) then
1193 do k = 1, km
1194 xmf(i,k) = 0.0
1195 end do
1196 end if
1197 end do
1199 do i = 1,im
1200 tem = sqrt(u10m(i)**2+v10m(i)**2)
1201 mffac = (1. - min(max(tem - 20.0, 0.0), 10.0)/10.)
1202 do k = 1, km
1203 xmf(i,k) = xmf(i,k)*mffac
1204 xmfd(i,k) = xmfd(i,k)*mffac
1205 enddo
1206 enddo
1207 endif
1208!
1209!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1210
1212 do k = 1, kmpbl
1213 do i = 1, im
1214 if(k < kpbl(i)) then
1215 tem = phih(i)/phim(i)
1216 ptem = sfcfrac*hpbl(i)
1217 tem1 = max(zi(i,k+1)-ptem, 0.)
1218 tem2 = tem1 / (hpbl(i) - ptem)
1219 if(pcnvflg(i)) then
1220 tem = min(tem, pr0)
1221 prn(i,k) = tem + (pr0 - tem) * tem2
1222 else
1223 tem = max(tem, pr0)
1224 prn(i,k) = tem
1225 endif
1226 prn(i,k) = min(prn(i,k),prmax)
1227 prn(i,k) = max(prn(i,k),prmin)
1228!
1229 ckz(i,k) = ck0 + (ck1 - ck0) * tem2
1230 ckz(i,k) = max(min(ckz(i,k), ck0), ck1)
1231 chz(i,k) = ch0 + (ch1 - ch0) * tem2
1232 chz(i,k) = max(min(chz(i,k), ch0), ch1)
1233!
1234 endif
1235 enddo
1236 enddo
1237!
1238! Above a threshold height (hcrinv), the background vertical
1239! diffusivities & mixing length
1240! in the inversion layers are set to much smaller values (xkinv1 &
1241! rlmn1)
1242!
1243! Below the threshold height (hcrinv), the background vertical
1244! diffusivities & mixing length
1245! in the inversion layers are increased with increasing roughness
1246! length & vegetation fraction
1247!
1248 do k = 1,km1
1249 do i=1,im
1250 if(zi(i,k+1) > hcrinv) then
1251 tem1 = tvx(i,k+1)-tvx(i,k)
1252 if(tem1 >= 0.) then
1253 xkzo(i,k) = min(xkzo(i,k), xkinv1)
1254 xkzmo(i,k) = min(xkzmo(i,k), xkinv1)
1255 rlmnz(i,k) = min(rlmnz(i,k), rlmn1)
1256 endif
1257 else
1258 tem1 = tvx(i,k+1)-tvx(i,k)
1259 if(tem1 > 0.) then
1260 ptem = xkzo(i,k) * zvfun(i)
1261 xkzo(i,k) = min(max(ptem, xkinv2), xkzo(i,k))
1262 ptem = xkzmo(i,k) * zvfun(i)
1263 xkzmo(i,k) = min(max(ptem, xkinv2), xkzmo(i,k))
1264 ptem = rlmnz(i,k) * zvfun(i)
1265 rlmnz(i,k) = min(max(ptem, rlmn2), rlmnz(i,k))
1266 endif
1267 endif
1268 enddo
1269 enddo
1270 do k = 2,km1
1271 do i=1,im
1272 rlmnz(i,k) = 0.5 * (rlmnz(i,k-1) + rlmnz(i,k))
1273 enddo
1274 enddo
1275!
1276!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1278!
1279 do k = 1, km1
1280 do i = 1, im
1281 zlup = 0.0
1282 mlenflg = .true.
1283 e2(i,k) = max(2.*tke(i,k), 0.001)
1284 do n = k, km1
1285 if(mlenflg) then
1286 dz = zl(i,n+1) - zl(i,n)
1287 tem1 = 2.*gotvx(i,n+1)*(thvx(i,k)-thvx(i,n+1))
1288 tem2 = cs0*sqrt(e2(i,n))*sqrt(shr2(i,n))
1289 e2(i,n+1) = e2(i,n) + (tem1 - tem2) * dz
1290 zlup = zlup + dz
1291 if(e2(i,n+1) < 0.) then
1292 ptem = e2(i,n+1) / (e2(i,n+1) - e2(i,n))
1293 zlup = zlup - ptem * dz
1294 zlup = max(zlup, 0.)
1295 mlenflg = .false.
1296 endif
1297 endif
1298 enddo
1299 zldn = 0.0
1300 mlenflg = .true.
1301 do n = k, 1, -1
1302 if(mlenflg) then
1303 if(n == 1) then
1304 dz = zl(i,1)
1305 tem = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
1306 tem1 = 2.*gotvx(i,n)*(tem-thvx(i,k))
1307 tem2 = ustar(i)*phims(i)/(vk*dz)
1308 tem2 = cs0*sqrt(e2(i,n))*tem2
1309 e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz
1310 else
1311 dz = zl(i,n) - zl(i,n-1)
1312 tem1 = 2.*gotvx(i,n-1)*(thvx(i,n-1)-thvx(i,k))
1313 tem2 = cs0*sqrt(e2(i,n))*sqrt(shr2(i,n-1))
1314 e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz
1315 endif
1316 zldn = zldn + dz
1317 if(e2(i,n-1) < 0.) then
1318 ptem = e2(i,n-1) / (e2(i,n-1) - e2(i,n))
1319 zldn = zldn - ptem * dz
1320 zldn = max(zldn, 0.)
1321 mlenflg = .false.
1322 endif
1323 endif
1324 enddo
1325!
1326 tem = 0.5 * (zi(i,k+1)-zi(i,k))
1327 tem1 = min(tem, rlmnz(i,k))
1340!
1341! Following Rodier et. al (2017), environmental wind shear effect on
1342! mixing length was included.
1343!
1344 ptem2 = min(zlup,zldn)
1345 rlam(i,k) = elmfac * ptem2
1346 rlam(i,k) = max(rlam(i,k), tem1)
1347 rlam(i,k) = min(rlam(i,k), rlmx)
1348!
1349 ptem2 = sqrt(zlup*zldn)
1350 ele(i,k) = elefac * ptem2
1351 ele(i,k) = max(ele(i,k), tem1)
1352 elmh(i,k)= elmhfac * ele(i,k)
1353 ele(i,k) = min(ele(i,k), elmx)
1354 elmh(i,k)= min(elmh(i,k), elmhmx)
1355!
1356 enddo
1357 enddo
1360 do k = 1, km1
1361 do i = 1, im
1362 tem = vk * zl(i,k)
1363 if (zol(i) < 0.) then
1364 ptem = 1. - 100. * zol(i)
1365 ptem1 = ptem**0.2
1366 zk = tem * ptem1
1367 elseif (zol(i) >= 1.) then
1368 zk = tem / 3.7
1369 else
1370 ptem = 1. + 2.7 * zol(i)
1371 zk = tem / ptem
1372 endif
1373
1374 if (tc_pbl == 0) then
1375 elm(i,k) = zk*rlam(i,k)/(rlam(i,k)+zk)
1377 if ( sfc_rlm == 1 ) then
1378 if ( sfcflg(i) .and.
1379 & zl(i,k) < min(100.0,hpbl(i)*0.05) ) elm(i,k)=zk
1380 endif
1381 else if (tc_pbl == 1) then
1382 ! new blending method for mixing length
1383 elm(i,k) = sqrt( 1.0/( 1.0/(zk**2)+1.0/(rlam(i,k)**2) ) )
1384 endif
1385
1386!
1387 if(k == 1) elm(i,k)=zk
1388!
1389 dz = zi(i,k+1) - zi(i,k)
1390 tem = max(gdx(i),dz)
1391 elm(i,k) = min(elm(i,k), tem)
1392
1393 if (tc_pbl == 0) then
1394 ele(i,k) = min(ele(i,k), tem)
1395 else if (tc_pbl == 1) then
1396 ele(i,k) = elm(i,k)
1397 endif
1398!
1399 enddo
1400 enddo
1401 do i = 1, im
1402 elm(i,km) = elm(i,km1)
1403 ele(i,km) = ele(i,km1)
1404 elmh(i,km)= elmh(i,km1)
1405 enddo
1406!
1407!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1410!
1411 if(tte_edmf) then
1412!
1413 do k = 1, km1
1414 do i = 1, im
1415 ptem = sfcfrac*hpbl(i)
1416 if (zi(i,k+1) <= ptem) then
1417 ri = ris(i)
1418 else
1419 ri = max(bf(i,k)/shr2(i,k),rimin)
1420 endif
1421 if(ri < 0) then
1422 tem = 2. * ri - pr0
1423 epotte = ri / tem
1424 else
1425 tem = pr0 + 3. * ri
1426 epotte = ri / tem
1427 endif
1428 tkeh(i,k) = tteh(i,k) * (1. - epotte)
1429 tesq(i,k) = tkeh(i,k) / sqrt(tteh(i,k))
1430 enddo
1431 enddo
1432!
1433 else
1434!
1435 do k = 1, km1
1436 do i = 1, im
1437 tesq(i,k) = sqrt(tkeh(i,k))
1438 enddo
1439 enddo
1440!
1441 endif
1442!
1443 do k = 1, km1
1444 do i = 1, im
1445 tem = 0.5 * (elm(i,k) + elm(i,k+1))
1446 tem = tem * tesq(i,k)
1447 ri = max(bf(i,k)/shr2(i,k),rimin)
1448 if(k < kpbl(i)) then
1449 if(pcnvflg(i)) then
1450 dku(i,k) = ckz(i,k) * tem
1451 dkt(i,k) = dku(i,k) / prn(i,k)
1452 else
1453 if(ri < 0.) then ! unstable regime
1454 dku(i,k) = ckz(i,k) * tem
1455 dkt(i,k) = dku(i,k) / prn(i,k)
1456 else ! stable regime
1457 dkt(i,k) = chz(i,k) * tem
1458 dku(i,k) = dkt(i,k) * prn(i,k)
1459 endif
1460 endif
1461 else
1462 if(ri < 0.) then ! unstable regime
1463 dku(i,k) = ck1 * tem
1464 dkt(i,k) = rchck * dku(i,k)
1465 else ! stable regime
1466 dkt(i,k) = ch1 * tem
1467 prnum = pr0 + 2.1 * ri
1468 prnum = min(prnum,prmax)
1469 dku(i,k) = dkt(i,k) * prnum
1470 endif
1471 endif
1472!
1473 if(scuflg(i)) then
1474 if(k >= mrad(i) .and. k < krad(i)) then
1475 if(tte_edmf) then
1476 tem1 = ck0 * tem
1477 else
1478 tem1 = ckz(i,k) * tem
1479 endif
1480 ptem1 = tem1 / prscu
1481 dku(i,k) = max(dku(i,k), tem1)
1482 dkt(i,k) = max(dkt(i,k), ptem1)
1483 endif
1484 endif
1485!
1486 dkq(i,k) = prtke * dkt(i,k)
1487!
1488 dkt(i,k) = min(dkt(i,k),dkmax)
1489 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
1490 dkq(i,k) = min(dkq(i,k),dkmax)
1491 dkq(i,k) = max(dkq(i,k),xkzo(i,k))
1492 dku(i,k) = min(dku(i,k),dkmax)
1493 dku(i,k) = max(dku(i,k),xkzmo(i,k))
1494!
1495 enddo
1496 enddo
1497!
1498!The following is for SA-3D-TKE
1499 if(sa3dtke) then
1500! 1. compute LES component of km, kh, and kq (Deardorff 1980)
1501! calculate thetal
1502 do k=1,km
1503 do i=1,im
1504 pix(i,k) = psk(i) / prslk(i,k)
1505 theta(i,k) = t1(i,k) * pix(i,k)
1506 tem=theta(i,k)/t1(i,k)
1507 if(ntiw > 0) then
1508 tem1=max(q1(i,k,ntcw),qlmin)+
1509 & max(q1(i,k,ntiw)+q1(i,k,5)+q1(i,k,6),qlmin)
1510 thetal(i,k)=theta(i,k)-(hvap+hfus)/cp*tem*tem1
1511 else
1512 tem1=max(q1(i,k,ntcw),qlmin)
1513 thetal(i,k)=theta(i,k)-hvap/cp*tem*tem1
1514 endif
1515 enddo
1516 enddo
1517
1518 do k=1,km
1519 do i=1,im
1520 dku_les(i,k) = 0.
1521 dkt_les(i,k) = 0.
1522 dkq_les(i,k) = 0.
1523 enddo
1524 enddo
1525!
1526! eddy diffusivities at model interface (zm level) in LES scale
1527!
1528 do k = 1, km1
1529 do i = 1, im
1530 tem=gotvx(i,k)*(thetal(i,k+1)-thetal(i,k))*rdzt(i,k)
1531 dz = zl(i,k+1) - zl(i,k)
1532 tem1=(garea(i)*dz)**h1
1533! calculate LES mixing length
1534 if(tem > 0.0) then
1535 elm_les=0.76*sqrt(tke(i,k))/sqrt(tem)
1536 elm_les=min(elm_les,tem1)
1537 else
1538 elm_les=tem1
1539 endif
1540! calculate km, kh, and kq for LES
1541 dku_les(i,k)=0.1*elm_les*sqrt(tkeh(i,k))
1542 dkt_les(i,k)=(1.0+2.0*elm_les/tem1)*dku_les(i,k)
1543 dkq_les(i,k)=dkt_les(i,k)
1544 dku_les(i,k) = min(dku_les(i,k),dkmaxles)
1545 dkt_les(i,k) = min(dkt_les(i,k),dkmaxles)
1546 dkq_les(i,k) = min(dkq_les(i,k),dkmaxles)
1547 enddo
1548 enddo
1549!
1550! calculate blending coefficients for km, kt, kq, and nonlocal mixing
1551! finding scale of large eddies from TKE
1552 do i=1,im
1553 tkemax(i) = tke(i,1)
1554 ktkemax(i) = 1
1555 enddo
1556 do k = 2, kmpbl
1557 do i = 1, im
1558 if(tke(i,k) > tkemax(i)) then
1559 tkemax(i) = tke(i,k)
1560 ktkemax(i) = k
1561 endif
1562 enddo
1563 enddo
1564 do i=1,im
1565 flg(i) = .true.
1566 scl(i) = 0.
1567 if(zl(i,ktkemax(i)) > sclmax) then
1568 flg(i) = .false.
1569 scl(i) = sclmin
1570 endif
1571 enddo
1572 do k = 1, kmpbl
1573 do i = 1, im
1574 if(flg(i) .and. k > ktkemax(i)) then
1575 scl(i) = zl(i,k)
1576 tem = 0.5*tkemax(i)
1577 if(tke(i,k) < tem) flg(i) = .false.
1578 endif
1579 enddo
1580 enddo
1581 do i=1,im
1582 scl(i)=max(scl(i), sclmin)
1583 scl(i)=min(scl(i), sclmax)
1584 scl(i)=max(scl(i), hpbl(i))
1585 pfdx(i)=gdx(i)/scl(i)
1586 enddo
1587!
1588 do i = 1, im
1589! partition function for local fluxes
1590 pfl(i)=cpl1*(pfdx(i)**2+cpl2*pfdx(i)**0.5-cpl3)/
1591 & (pfdx(i)**2+cpl4*pfdx(i)**0.5+cpl5)+cpl6
1592 pfl(i)=min(max(pfl(i),0.0),1.0)
1593! partition function for TKE
1594 pftke(i)=(pfdx(i)**2+cptke1*pfdx(i)**(2./3.))/
1595 & (pfdx(i)**2+cptke2*pfdx(i)**(2./3.)+cptke3)
1596 pftke(i)=min(max(pftke(i),0.0),1.0)
1597 enddo
1598!
1599! blending LES and MS components of vertical km,kt, and kq
1600!
1601 do k = 1,km1
1602 do i=1,im
1603 dkq(i,k)=(1.0-pfl(i))*dkq_les(i,k)+pfl(i)*dkq(i,k)
1604 dkt(i,k)=(1.0-pfl(i))*dkt_les(i,k)+pfl(i)*dkt(i,k)
1605 dku(i,k)=(1.0-pfl(i))*dku_les(i,k)+pfl(i)*dku(i,k)
1606 enddo
1607 enddo
1608!
1609! 2. compute MS horizontal km
1610!
1611 do k = 1, km
1612 do i = 1, im
1613 dku_h(i,k)=ckh*elmh(i,k)*sqrt(tke(i,k))
1614 dkq_h(i,k)=dku_h(i,k)
1615 enddo
1616 enddo
1617!
1618! eddy diffusivities at model layer (zl level) in LES scale
1619!
1620 do k = 1, km1
1621 do i = 1, im
1622 if(k > 1) then
1623 dz = zl(i,k+1) - zl(i,k-1)
1624 tem=gotvx(i,k)*(thetal(i,k+1)-thetal(i,k-1))/dz
1625 else
1626 dz = zl(i,k+1)
1627 tem=gotvx(i,k)*(thetal(i,k+1)-thvs(i))/dz
1628 endif
1629 dz = zi(i,k+1) - zi(i,k)
1630 tem1=(garea(i)*dz)**h1
1631! calculate LES mixing length
1632 if(tem > 0.0) then
1633 elm_les=0.76*sqrt(tke(i,k))/sqrt(tem)
1634 elm_les=min(elm_les,tem1)
1635 else
1636 elm_les=tem1
1637 endif
1638 ele_les(i,k)=elm_les
1639! calculate km, kh, and kq for LES
1640 dku_les(i,k)=0.1*elm_les*sqrt(tke(i,k))
1641 dkq_les(i,k)=(1.0+2.0*elm_les/tem1)*dku_les(i,k)
1642 dku_les(i,k) = min(dku_les(i,k),dkmaxles)
1643 dkq_les(i,k) = min(dkq_les(i,k),dkmaxles)
1644 enddo
1645 enddo
1646!
1647 do k = 1,km1
1648 do i=1,im
1649 dku_h(i,k)=(1.0-pfl(i))*dku_les(i,k)+pfl(i)*dku_h(i,k)
1650 dkq_h(i,k)=(1.0-pfl(i))*dkq_les(i,k)+pfl(i)*dkq_h(i,k)
1651 enddo
1652 enddo
1653 do i = 1, im
1654 dku_h(i,km)=dku_h(i,km1)
1655 dkq_h(i,km)=dkq_h(i,km1)
1656 enddo
1657!
1658 endif !sa3dtke
1659
1660!PCC_CANOPY------------------------------------
1661 kount=0 !IVAI
1662 if (do_canopy .and. cplaqm) then
1663
1664!IVAI
1665! Output 3D pbl diags
1666! aux3d(:,:,5) = dku(:,:) ! Out
1667! aux3d(:,:,3) = dkt(:,:) ! Out
1668
1669! 3-Layer Sub-Canopy effect
1670 if (present(dku_can)) dku(:,1:km) = dku_can(:,1:km)
1671 if (present(dkt_can)) then
1672 dkt(:,1:km) = dkt_can(:,1:km)
1673 dkq(:,1:km-1) = prtke * dkt_can(:,1:km-1)
1674 endif
1675
1676 endif !do_canopy .and. cplaqm
1677
1680!
1681 do k = 1, km1
1682 do i = 1, im
1683 if(k == 1) then
1684 tem = ckz(i,1)
1685 tem1 = 0.5 * xkzmo(i,1)
1686 else
1687 tem = 0.5 * (ckz(i,k-1) + ckz(i,k))
1688 tem1 = 0.5 * (xkzmo(i,k-1) + xkzmo(i,k))
1689 endif
1690 ptem = tem1 / (tem * elm(i,k))
1691 tkmnz(i,k) = ptem * ptem
1692 tkmnz(i,k) = min(tkmnz(i,k), tkbmx)
1693 tkmnz(i,k) = max(tkmnz(i,k), tkmin)
1694 enddo
1695 enddo
1696!
1697!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1700!
1701 do k = 1, km1
1702 do i = 1, im
1703 if (k == 1) then
1704 tem = -dkt(i,1) * bf(i,1)
1705! if(pcnvflg(i)) then
1706! ptem1 = xmf(i,1) * buou(i,1)
1707! else
1708 ptem1 = 0.
1709! endif
1710 if(scuflg(i) .and. mrad(i) == 1) then
1711 ptem2 = xmfd(i,1) * buod(i,1)
1712 else
1713 ptem2 = 0.
1714 endif
1715 tem = tem + ptem1 + ptem2
1716 buop = 0.5 * (gotvx(i,1) * sflux(i) + tem)
1717!
1718 if(sa3dtke) then
1719 tem = 2. * dku_h(i,1)
1720 tem1 = dku(i,1)*def_1(i,1)+tem*def_2(i,1)
1721 else
1722 tem1 = dku(i,1) * shr2(i,1)
1723 endif
1724!
1725 tem = (u1(i,2)-u1(i,1))*rdzt(i,1)
1726! if(pcnvflg(i)) then
1727! ptem = xmf(i,1) * tem
1728! ptem1 = 0.5 * ptem * (u1(i,2)-ucko(i,2))
1729! else
1730 ptem1 = 0.
1731! endif
1732 if(scuflg(i) .and. mrad(i) == 1) then
1733 ptem = ucdo(i,1)+ucdo(i,2)-u1(i,1)-u1(i,2)
1734 ptem = 0.5 * tem * xmfd(i,1) * ptem
1735 else
1736 ptem = 0.
1737 endif
1738 ptem1 = ptem1 + ptem
1739!
1740 tem = (v1(i,2)-v1(i,1))*rdzt(i,1)
1741! if(pcnvflg(i)) then
1742! ptem = xmf(i,1) * tem
1743! ptem2 = 0.5 * ptem * (v1(i,2)-vcko(i,2))
1744! else
1745 ptem2 = 0.
1746! endif
1747 if(scuflg(i) .and. mrad(i) == 1) then
1748 ptem = vcdo(i,1)+vcdo(i,2)-v1(i,1)-v1(i,2)
1749 ptem = 0.5 * tem * xmfd(i,1) * ptem
1750 else
1751 ptem = 0.
1752 endif
1753 ptem2 = ptem2 + ptem
1754!
1755 tem2 = stress(i)*ustar(i)*phims(i)/(vk*zl(i,1))
1756 shrp = 0.5 * (tem1 + ptem1 + ptem2 + tem2)
1757 else
1758 tem1 = -dkt(i,k-1) * bf(i,k-1)
1759 tem2 = -dkt(i,k) * bf(i,k)
1760 tem = 0.5 * (tem1 + tem2)
1761 if(pcnvflg(i) .and. k <= kpbl(i)) then
1762 ptem = 0.5 * (xmf(i,k-1) + xmf(i,k))
1763 ptem1 = ptem * buou(i,k)
1764 else
1765 ptem1 = 0.
1766 endif
1767 if(scuflg(i)) then
1768 if(k >= mrad(i) .and. k < krad(i)) then
1769 ptem0 = 0.5 * (xmfd(i,k-1) + xmfd(i,k))
1770 ptem2 = ptem0 * buod(i,k)
1771 else
1772 ptem2 = 0.
1773 endif
1774 else
1775 ptem2 = 0.
1776 endif
1777 buop = tem + ptem1 + ptem2
1778!
1779 if(sa3dtke) then
1780! obtaining 3d shear production from dycore
1781 tem2 = 2.*dku_h(i,k)
1782 tem1 = dku(i,k-1)*def_1(i,k-1)
1783 tem2 = dku(i,k)*def_1(i,k)+tem2*def_2(i,k)
1784 else
1785 tem1 = dku(i,k-1) * shr2(i,k-1)
1786 tem2 = dku(i,k) * shr2(i,k)
1787 endif
1788 tem = 0.5 * (tem1 + tem2)
1789 tem1 = (u1(i,k+1)-u1(i,k))*rdzt(i,k)
1790 tem2 = (u1(i,k)-u1(i,k-1))*rdzt(i,k-1)
1791 if(pcnvflg(i) .and. k <= kpbl(i)) then
1792 ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2
1793 ptem1 = 0.5 * ptem * (u1(i,k)-ucko(i,k))
1794 else
1795 ptem1 = 0.
1796 endif
1797 if(scuflg(i)) then
1798 if(k >= mrad(i) .and. k < krad(i)) then
1799 ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2
1800 ptem2 = 0.5 * ptem0 * (ucdo(i,k)-u1(i,k))
1801 else
1802 ptem2 = 0.
1803 endif
1804 else
1805 ptem2 = 0.
1806 endif
1807 shrp = tem + ptem1 + ptem2
1808 tem1 = (v1(i,k+1)-v1(i,k))*rdzt(i,k)
1809 tem2 = (v1(i,k)-v1(i,k-1))*rdzt(i,k-1)
1810 if(pcnvflg(i) .and. k <= kpbl(i)) then
1811 ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2
1812 ptem1 = 0.5 * ptem * (v1(i,k)-vcko(i,k))
1813 else
1814 ptem1 = 0.
1815 endif
1816 if(scuflg(i)) then
1817 if(k >= mrad(i) .and. k < krad(i)) then
1818 ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2
1819 ptem2 = 0.5 * ptem0 * (vcdo(i,k)-v1(i,k))
1820 else
1821 ptem2 = 0.
1822 endif
1823 else
1824 ptem2 = 0.
1825 endif
1826 shrp = shrp + ptem1 + ptem2
1827 endif
1828 if(tte_edmf) then
1829 if(buop > 0.) then
1830 prod(i,k) = 2. * buop + shrp
1831 else
1832 prod(i,k) = shrp
1833 endif
1834 else
1835 prod(i,k) = buop + shrp
1836 endif
1837 enddo
1838 enddo
1839!
1840!----------------------------------------------------------------------
1842!
1843 if(sa3dtke) then
1844!The following is for SA-3D-TKE
1845 dtn = dt2 / float(ndt)
1846 do n = 1, ndt
1847 do k = 1,km1
1848 do i=1,im
1849 tem = sqrt(te(i,k))
1850! calculating 3D TKE transport and pressure correlation
1851 ptem1 = ce0 / ele(i,k)
1852 dz = zi(i,k+1) - zi(i,k)
1853 tem1=(garea(i)*dz)**h1
1854 tem2=0.19+0.51*ele_les(i,k)/tem1
1855 ptem2= tem2 / ele_les(i,k)
1856 ptem=(1.0-pftke(i))*ptem2+pftke(i)*ptem1
1857 disste = ptem * te(i,k) * tem
1858 tem1 = prod(i,k) + te(i,k) / dtn
1859 disste=max(min(disste, tem1), 0.)
1860 if(.not. tte_edmf) diss(i,k) = disste
1861! tem=2.0*def_3(i,k)
1862 tem=def_3(i,k)
1863! tem=min(tem,1.0)
1864 te(i,k) = te(i,k) + dtn * (prod(i,k)-disste+tem)
1865! te(i,k) = max(te(i,k), tkmin)
1866 te(i,k) = max(te(i,k), tkmnz(i,k))
1867 enddo
1868 enddo
1869 enddo
1870 else
1871 dtn = dt2 / float(ndt)
1872 do n = 1, ndt
1873 do k = 1,km1
1874 do i=1,im
1875 tem = sqrt(te(i,k))
1876 ptem = ce0 / ele(i,k)
1877 disste = ptem * te(i,k) * tem
1878 tem1 = prod(i,k) + te(i,k) / dtn
1879 disste = max(min(disste, tem1), 0.)
1880 if(.not. tte_edmf) diss(i,k) = disste
1881 te(i,k) = te(i,k) + dtn * (prod(i,k)-disste)
1882 te(i,k) = max(te(i,k), tkmnz(i,k))
1883! te(i,k) = max(te(i,k), tkmin)
1884 enddo
1885 enddo
1886 enddo
1887 endif !sa3dtke
1888!
1889! TKE dissipation for dissipative heating computation in TTE-EDMF
1890!
1891 if(tte_edmf) then
1892 do k = 1, km1
1893 do i = 1, im
1894 tem = sqrt(tke(i,k))
1895 if(sa3dtke) then
1896 ptem1 = ce0 / ele(i,k)
1897 dz = zi(i,k+1) - zi(i,k)
1898 tem1=(garea(i)*dz)**h1
1899 tem2=0.19+0.51*ele_les(i,k)/tem1
1900 ptem2= tem2 / ele_les(i,k)
1901 ptem=(1.0-pftke(i))*ptem2+pftke(i)*ptem1
1902 diss(i,k) = ptem * tke(i,k) * tem
1903 else
1904 ptem = ce0 / ele(i,k)
1905 diss(i,k) = ptem * tke(i,k) * tem
1906 endif
1907 enddo
1908 enddo
1909 endif
1910!
1912!
1913 do k = 1, km
1914 do i = 1, im
1915 if(pcnvflg(i)) then
1916 qcko(i,k,ntke) = te(i,k)
1917 endif
1918 if(scuflg(i)) then
1919 qcdo(i,k,ntke) = te(i,k)
1920 endif
1921 enddo
1922 enddo
1923 do k = 2, kmpbl
1924 do i = 1, im
1925 if (pcnvflg(i) .and. k <= kpbl(i)) then
1926 dz = zl(i,k) - zl(i,k-1)
1927 tem = 0.5 * xlamue(i,k-1) * dz
1928 factor = 1. + tem
1929 qcko(i,k,ntke)=((1.-tem)*qcko(i,k-1,ntke)+tem*
1930 & (te(i,k)+te(i,k-1)))/factor
1931 endif
1932 enddo
1933 enddo
1934 do k = kmscu, 1, -1
1935 do i = 1, im
1936 if (scuflg(i) .and. k < krad(i)) then
1937 if(k >= mrad(i)) then
1938 dz = zl(i,k+1) - zl(i,k)
1939 tem = 0.5 * xlamde(i,k) * dz
1940 factor = 1. + tem
1941 qcdo(i,k,ntke)=((1.-tem)*qcdo(i,k+1,ntke)+tem*
1942 & (te(i,k)+te(i,k+1)))/factor
1943 endif
1944 endif
1945 enddo
1946 enddo
1947!
1948!--------------------------------------------------------
1949! compute variables for TVD flux-limiter scheme
1950! on environmental subsidence and uplifting
1951!
1952 kps = max(kmpbl, kmscu)
1953!
1954! for moisture and tracers including hydrometeors
1955!
1956 do n=1,ntrac1
1957 do k=1,kps
1958 do i=1,im
1959 qh(i,k,n) = 0.5 * (q1(i,k,n)+q1(i,k+1,n))
1960 enddo
1961 enddo
1962 do k=1,kps
1963 do i=1,im
1964 q_diff(i,k,n) = q1(i,k,n) - q1(i,k+1,n)
1965 enddo
1966 enddo
1967 do i=1,im
1968 if(q1(i,1,n) >= 0.) then
1969 q_diff(i,0,n) = max(0.,2.*q1(i,1,n)-q1(i,2,n))-
1970 & q1(i,1,n)
1971 else
1972 q_diff(i,0,n) = min(0.,2.*q1(i,1,n)-q1(i,2,n))-
1973 & q1(i,1,n)
1974 endif
1975 enddo
1976 enddo
1977!
1978 do n = 1, ntrac1
1979!
1980 do k = 1, kps
1981 do i = 1, im
1982 kmx = max(kpbl(i), krad(i))
1983 q_half(i,k,n) = qh(i,k,n)
1984 if((pcnvflg(i) .or. scuflg(i)) .and. k < kmx) then
1985 tem = 0.
1986 if(pcnvflg(i) .and. k < kpbl(i)) then
1987 tem = xmf(i,k)
1988 endif
1989 if(scuflg(i) .and.
1990 & (k >= mrad(i) .and. k < krad(i))) then
1991 tem = tem - xmfd(i,k)
1992 endif
1993 if(tem > 0.) then
1994 rrkp = 0.
1995 if(abs(q_diff(i,k,n)) > 1.e-22)
1996 & rrkp = q_diff(i,k+1,n) / q_diff(i,k,n)
1997 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
1998 q_half(i,k,n) = q1(i,k+1,n) +
1999 & phkp*(qh(i,k,n)-q1(i,k+1,n))
2000 elseif (tem < 0.) then
2001 rrkp = 0.
2002 if(abs(q_diff(i,k,n)) > 1.e-22)
2003 & rrkp = q_diff(i,k-1,n) / q_diff(i,k,n)
2004 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
2005 q_half(i,k,n) = q1(i,k,n) +
2006 & phkp*(qh(i,k,n)-q1(i,k,n))
2007 endif
2008 endif
2009 enddo
2010 enddo
2011!
2012 enddo
2013!
2014! for TKE or TTE
2015!
2016 do k=1,kps
2017 do i=1,im
2018 tei(i,k) = 0.5 * (te(i,k)+te(i,k+1))
2019 enddo
2020 enddo
2021
2022 do k=1,kps
2023 do i=1,im
2024 e_diff(i,k) = te(i,k) - te(i,k+1)
2025 enddo
2026 enddo
2027 do i=1,im
2028 if(te(i,1) >= 0.) then
2029 e_diff(i,0) = max(0.,2.*te(i,1)-te(i,2))-
2030 & te(i,1)
2031 else
2032 e_diff(i,0) = min(0.,2.*te(i,1)-te(i,2))-
2033 & te(i,1)
2034 endif
2035 enddo
2036!
2037 do k = 1, kps
2038 do i = 1, im
2039 kmx = max(kpbl(i), krad(i))
2040 e_half(i,k) = tei(i,k)
2041 if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
2042 tem = 0.
2043 if(pcnvflg(i) .and. k < kpbl(i)) then
2044 tem = xmf(i,k)
2045 endif
2046 if(scuflg(i) .and.
2047 & (k >= mrad(i) .and. k < krad(i))) then
2048 tem = tem - xmfd(i,k)
2049 endif
2050 if(tem > 0.) then
2051 rrkp = 0.
2052 if(abs(e_diff(i,k)) > 1.e-22)
2053 & rrkp = e_diff(i,k+1) / e_diff(i,k)
2054 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
2055 e_half(i,k) = te(i,k+1) +
2056 & phkp*(tei(i,k)-te(i,k+1))
2057 elseif (tem < 0.) then
2058 rrkp = 0.
2059 if(abs(e_diff(i,k)) > 1.e-22)
2060 & rrkp = e_diff(i,k-1) / e_diff(i,k)
2061 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
2062 e_half(i,k) = te(i,k) +
2063 & phkp*(tei(i,k)-te(i,k))
2064 endif
2065 endif
2066 enddo
2067 enddo
2068!
2069!----------------------------------------------------------------------
2071!
2072 do i=1,im
2073 ad(i,1) = 1.0
2074 f1(i,1) = te(i,1)
2075 enddo
2076!
2077 do k = 1,km1
2078 do i=1,im
2079 dtodsd = dt2/del(i,k)
2080 dtodsu = dt2/del(i,k+1)
2081 dsig = prsl(i,k)-prsl(i,k+1)
2082 rdz = rdzt(i,k)
2083 tem1 = dsig * dkq(i,k) * rdz
2084 dsdz2 = tem1 * rdz
2085 au(i,k) = -dtodsd*dsdz2
2086 al(i,k) = -dtodsu*dsdz2
2087 ad(i,k) = ad(i,k)-au(i,k)
2088 ad(i,k+1)= 1.-al(i,k)
2089 tem2 = dsig * rdz
2090!
2091 if(pcnvflg(i) .and. k < kpbl(i)) then
2092 ptem = 0.5 * tem2 * xmf(i,k)
2093 ptem1 = dtodsd * ptem
2094 ptem2 = dtodsu * ptem
2095 ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke)
2096 f1(i,k) = f1(i,k) - ptem * ptem1
2097 f1(i,k+1) = te(i,k+1) + ptem * ptem2
2098 else
2099 f1(i,k+1) = te(i,k+1)
2100 endif
2101!
2102 if(scuflg(i)) then
2103 if(k >= mrad(i) .and. k < krad(i)) then
2104 ptem = 0.5 * tem2 * xmfd(i,k)
2105 ptem1 = dtodsd * ptem
2106 ptem2 = dtodsu * ptem
2107 ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke)
2108 f1(i,k) = f1(i,k) + ptem * ptem1
2109 f1(i,k+1) = f1(i,k+1) - ptem * ptem2
2110 endif
2111 endif
2112!
2113 kmx = max(kpbl(i), krad(i))
2114 if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
2115 ptem = tem2 * (xmf(i,k) - xmfd(i,k))
2116 ptem1 = dtodsd * ptem
2117 ptem2 = dtodsu * ptem
2118 f1(i,k) = f1(i,k) + e_half(i,k) * ptem1
2119 f1(i,k+1) = f1(i,k+1) - e_half(i,k) * ptem2
2120 endif
2121!
2122 enddo
2123 enddo
2124c
2126c
2127 call tridit(im,km,1,al,ad,au,f1,au,f1)
2128!
2129! Negative TKE or TTE are set to zero after borrowing it from positive
2130! values within the mass-flux transport layers
2131!
2132 do i = 1,im
2133 tsumn(i) = 0.
2134 tsump(i) = 0.
2135 rtnp(i) = 1.
2136 enddo
2137 do k = 1,kps
2138 do i = 1,im
2139 if(pcnvflg(i) .and. scuflg(i)) then
2140 kbx = 1
2141 kmx = max(kpbl(i), krad(i))
2142 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2143 kbx = 1
2144 kmx = kpbl(i)
2145 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2146 kbx = mrad(i)
2147 kmx = krad(i)
2148 endif
2149 if((pcnvflg(i) .or. scuflg(i)) .and.
2150 & (k >= kbx .and. k <= kmx)) then
2151 tem = f1(i,k) * del(i,k) * gravi
2152 if(f1(i,k) < 0.) tsumn(i) = tsumn(i) + tem
2153 if(f1(i,k) > 0.) tsump(i) = tsump(i) + tem
2154 endif
2155 enddo
2156 enddo
2157 do i = 1,im
2158 if(pcnvflg(i) .or. scuflg(i)) then
2159 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2160 if(tsump(i) > abs(tsumn(i))) then
2161 rtnp(i) = tsumn(i) / tsump(i)
2162 else
2163 rtnp(i) = tsump(i) / tsumn(i)
2164 endif
2165 endif
2166 endif
2167 enddo
2168 do k = 1,kps
2169 do i = 1,im
2170 if(pcnvflg(i) .and. scuflg(i)) then
2171 kbx = 1
2172 kmx = max(kpbl(i), krad(i))
2173 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2174 kbx = 1
2175 kmx = kpbl(i)
2176 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2177 kbx = mrad(i)
2178 kmx = krad(i)
2179 endif
2180 if((pcnvflg(i) .or. scuflg(i)) .and.
2181 & (k >= kbx .and. k <= kmx)) then
2182 if(rtnp(i) < 0.) then
2183 if(tsump(i) > abs(tsumn(i))) then
2184 if(f1(i,k) < 0.) f1(i,k) = 0.
2185 if(f1(i,k) > 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k)
2186 else
2187 if(f1(i,k) < 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k)
2188 if(f1(i,k) > 0.) f1(i,k) = 0.
2189 endif
2190 endif
2191 endif
2192 enddo
2193 enddo
2194!
2195! To remove negative TKEs or TTEs which were leaked out of the mass-flux transport layers
2196! by eddy diffusion or potential negative TKEs or TTEs from the diffusion scheme,
2197! positive TKEs or TTEs are borrowed again now from the entire layers
2198!
2199 do i = 1,im
2200 tsumn(i) = 0.
2201 tsump(i) = 0.
2202 rtnp(i) = 1.
2203 enddo
2204 do k = 1,km
2205 do i = 1,im
2206 tem = f1(i,k) * del(i,k) * gravi
2207 if(f1(i,k) < 0.) tsumn(i) = tsumn(i) + tem
2208 if(f1(i,k) > 0.) tsump(i) = tsump(i) + tem
2209 enddo
2210 enddo
2211 do i = 1,im
2212 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2213 if(tsump(i) > abs(tsumn(i))) then
2214 rtnp(i) = tsumn(i) / tsump(i)
2215 else
2216 rtnp(i) = tsump(i) / tsumn(i)
2217 endif
2218 endif
2219 enddo
2220 do k = 1,km
2221 do i = 1,im
2222 if(rtnp(i) < 0.) then
2223 if(tsump(i) > abs(tsumn(i))) then
2224 if(f1(i,k) < 0.) f1(i,k) = 0.
2225 if(f1(i,k) > 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k)
2226 else
2227 if(f1(i,k) < 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k)
2228 if(f1(i,k) > 0.) f1(i,k) = 0.
2229 endif
2230 endif
2231 enddo
2232 enddo
2233c
2235c
2236 do k = 1,km
2237 do i = 1,im
2238! f1(i,k) = max(f1(i,k), tkmin)
2239 qtend = (f1(i,k)-q1(i,k,ntke))*rdt
2240 rtg(i,k,ntke) = rtg(i,k,ntke)+qtend
2241 enddo
2242 enddo
2243 if(ldiag3d .and. present(dtend)) then
2244 idtend = dtidx(ntke+100,index_of_process_pbl)
2245 if(idtend>0) then
2246 dtend(1:im,1:km,idtend) = dtend(1:im,1:km,idtend) +
2247 & (f1(1:im,1:km)-q1(1:im,1:km,ntke))*rdt
2248 endif
2249 endif
2250c
2252c
2253 do i=1,im
2254 ad(i,1) = 1.
2255 f1(i,1) = t1(i,1) + dtdz1(i) * heat(i)
2256 f2(i,1) = q1(i,1,1) + dtdz1(i) * evap(i)
2257 enddo
2258 if(ntrac1 >= 2) then
2259 do n = 2, ntrac1
2260 is = (n-1) * km
2261 do i = 1, im
2262 f2(i,1+is) = q1(i,1,n)
2263 enddo
2264 enddo
2265 endif
2266c
2267 do k = 1,km1
2268 do i = 1,im
2269 dtodsd = dt2/del(i,k)
2270 dtodsu = dt2/del(i,k+1)
2271 dsig = prsl(i,k)-prsl(i,k+1)
2272 rdz = rdzt(i,k)
2273 tem1 = dsig * dkt(i,k) * rdz
2274 dsdzt = tem1 * gocp
2275 if (use_lpt > 0) then
2276 dsdzt = dsdzt-tem1*elocp*(qliq(i,k+1)-qliq(i,k))*rdz
2277 & -(1+0.33/2.5)*tem1*elocp*(qice(i,k+1)-qice(i,k))*rdz
2278 endif
2279 dsdz2 = tem1 * rdz
2280 au(i,k) = -dtodsd*dsdz2
2281 al(i,k) = -dtodsu*dsdz2
2282 ad(i,k) = ad(i,k)-au(i,k)
2283 ad(i,k+1)= 1.-al(i,k)
2284 tem2 = dsig * rdz
2285!
2286 if(pcnvflg(i) .and. k < kpbl(i)) then
2287 ptem = 0.5 * tem2 * xmf(i,k)
2288 ptem1 = dtodsd * ptem
2289 ptem2 = dtodsu * ptem
2290 tem = t1(i,k) + t1(i,k+1)
2291 ptem = tcko(i,k) + tcko(i,k+1)
2292 f1(i,k) = f1(i,k)+dtodsd*dsdzt-(ptem-tem)*ptem1
2293 f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+(ptem-tem)*ptem2
2294 ptem = qcko(i,k,1) + qcko(i,k+1,1)
2295 f2(i,k) = f2(i,k) - ptem * ptem1
2296 f2(i,k+1) = q1(i,k+1,1) + ptem * ptem2
2297 else
2298 f1(i,k) = f1(i,k)+dtodsd*dsdzt
2299 f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
2300 f2(i,k+1) = q1(i,k+1,1)
2301 endif
2302!
2303 if(scuflg(i)) then
2304 if(k >= mrad(i) .and. k < krad(i)) then
2305 ptem = 0.5 * tem2 * xmfd(i,k)
2306 ptem1 = dtodsd * ptem
2307 ptem2 = dtodsu * ptem
2308 ptem = tcdo(i,k) + tcdo(i,k+1)
2309 tem = t1(i,k) + t1(i,k+1)
2310 f1(i,k) = f1(i,k) + (ptem - tem) * ptem1
2311 f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2
2312 ptem = qcdo(i,k,1) + qcdo(i,k+1,1)
2313 f2(i,k) = f2(i,k) + ptem * ptem1
2314 f2(i,k+1) = f2(i,k+1) - ptem * ptem2
2315 endif
2316 endif
2317!
2318 kmx = max(kpbl(i), krad(i))
2319 if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
2320 ptem = tem2 * (xmf(i,k) - xmfd(i,k))
2321 ptem1 = dtodsd * ptem
2322 ptem2 = dtodsu * ptem
2323 f2(i,k) = f2(i,k) + q_half(i,k,1) * ptem1
2324 f2(i,k+1) = f2(i,k+1) - q_half(i,k,1) * ptem2
2325 endif
2326!
2327 enddo
2328 enddo
2329!
2330 if(ntrac1 >= 2) then
2331 do n = 2, ntrac1
2332 is = (n-1) * km
2333 do k = 1, km1
2334 do i = 1, im
2335 dtodsd = dt2/del(i,k)
2336 dtodsu = dt2/del(i,k+1)
2337 dsig = prsl(i,k)-prsl(i,k+1)
2338 tem2 = dsig * rdzt(i,k)
2339!
2340 if(pcnvflg(i) .and. k < kpbl(i)) then
2341 ptem = 0.5 * tem2 * xmf(i,k)
2342 ptem1 = dtodsd * ptem
2343 ptem2 = dtodsu * ptem
2344 ptem = qcko(i,k,n) + qcko(i,k+1,n)
2345 f2(i,k+is) = f2(i,k+is) - ptem * ptem1
2346 f2(i,k+1+is)= q1(i,k+1,n) + ptem * ptem2
2347 else
2348 f2(i,k+1+is) = q1(i,k+1,n)
2349 endif
2350!
2351 if(scuflg(i)) then
2352 if(k >= mrad(i) .and. k < krad(i)) then
2353 ptem = 0.5 * tem2 * xmfd(i,k)
2354 ptem1 = dtodsd * ptem
2355 ptem2 = dtodsu * ptem
2356 ptem = qcdo(i,k,n) + qcdo(i,k+1,n)
2357 f2(i,k+is) = f2(i,k+is) + ptem * ptem1
2358 f2(i,k+1+is)= f2(i,k+1+is) - ptem * ptem2
2359 endif
2360 endif
2361!
2362 kmx = max(kpbl(i), krad(i))
2363 if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
2364 ptem = tem2 * (xmf(i,k) - xmfd(i,k))
2365 ptem1 = dtodsd * ptem
2366 ptem2 = dtodsu * ptem
2367 f2(i,k+is) = f2(i,k+is) + q_half(i,k,n) * ptem1
2368 f2(i,k+1+is) = f2(i,k+1+is) - q_half(i,k,n) * ptem2
2369 endif
2370!
2371 enddo
2372 enddo
2373 enddo
2374 endif
2375c
2377c
2378 call tridin(im,km,ntrac1,al,ad,au,f1,f2,au,f1,f2)
2379!
2380! Negative moisture is set to zero after borrowing it from
2381! positive values within the mass-flux transport layers
2382!
2383 do i = 1,im
2384 tsumn(i) = 0.
2385 tsump(i) = 0.
2386 rtnp(i) = 1.
2387 enddo
2388 do k = 1,kps
2389 do i = 1,im
2390 if(pcnvflg(i) .and. scuflg(i)) then
2391 kbx = 1
2392 kmx = max(kpbl(i), krad(i))
2393 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2394 kbx = 1
2395 kmx = kpbl(i)
2396 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2397 kbx = mrad(i)
2398 kmx = krad(i)
2399 endif
2400 if((pcnvflg(i) .or. scuflg(i)) .and.
2401 & (k >= kbx .and. k <= kmx)) then
2402 tem = f2(i,k) * del(i,k) * gravi
2403 if(f2(i,k) < 0.) tsumn(i) = tsumn(i) + tem
2404 if(f2(i,k) > 0.) tsump(i) = tsump(i) + tem
2405 endif
2406 enddo
2407 enddo
2408 do i = 1,im
2409 if(pcnvflg(i) .or. scuflg(i)) then
2410 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2411 if(tsump(i) > abs(tsumn(i))) then
2412 rtnp(i) = tsumn(i) / tsump(i)
2413 else
2414 rtnp(i) = tsump(i) / tsumn(i)
2415 endif
2416 endif
2417 endif
2418 enddo
2419 do k = 1,kps
2420 do i = 1,im
2421 if(pcnvflg(i) .and. scuflg(i)) then
2422 kbx = 1
2423 kmx = max(kpbl(i), krad(i))
2424 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2425 kbx = 1
2426 kmx = kpbl(i)
2427 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2428 kbx = mrad(i)
2429 kmx = krad(i)
2430 endif
2431 if((pcnvflg(i) .or. scuflg(i)) .and.
2432 & (k >= kbx .and. k <= kmx)) then
2433 if(rtnp(i) < 0.) then
2434 if(tsump(i) > abs(tsumn(i))) then
2435 if(f2(i,k) < 0.) f2(i,k) = 0.
2436 if(f2(i,k) > 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k)
2437 else
2438 if(f2(i,k) < 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k)
2439 if(f2(i,k) > 0.) f2(i,k) = 0.
2440 endif
2441 endif
2442 endif
2443 enddo
2444 enddo
2445!
2446! To remove negative moistures which were leaked out of the mass-flux transport layers
2447! by eddy diffusion or potential negative moistures from the diffusion scheme
2448! especially due to downward surface latent heat flux during nighttime,
2449! positive moistures are borrowed again now from the entire layers
2450!
2451 do i = 1,im
2452 tsumn(i) = 0.
2453 tsump(i) = 0.
2454 rtnp(i) = 1.
2455 enddo
2456 do k = 1,km
2457 do i = 1,im
2458 tem = f2(i,k) * del(i,k) * gravi
2459 if(f2(i,k) < 0.) tsumn(i) = tsumn(i) + tem
2460 if(f2(i,k) > 0.) tsump(i) = tsump(i) + tem
2461 enddo
2462 enddo
2463 do i = 1,im
2464 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2465 if(tsump(i) > abs(tsumn(i))) then
2466 rtnp(i) = tsumn(i) / tsump(i)
2467 else
2468 rtnp(i) = tsump(i) / tsumn(i)
2469 endif
2470 endif
2471 enddo
2472 do k = 1,km
2473 do i = 1,im
2474 if(rtnp(i) < 0.) then
2475 if(tsump(i) > abs(tsumn(i))) then
2476 if(f2(i,k) < 0.) f2(i,k) = 0.
2477 if(f2(i,k) > 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k)
2478 else
2479 if(f2(i,k) < 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k)
2480 if(f2(i,k) > 0.) f2(i,k) = 0.
2481 endif
2482 endif
2483 enddo
2484 enddo
2485!
2486! Negative hydrometeors & tracers are set to zero after
2487! borrowing them from positive values within the mass-flux
2488! transport layers
2489!
2490! For the negative liquid water, first borrow water from vapor
2491! and then borrow it from the other layers if there is still
2492! negative water
2493!
2494 if(ntrac1 >= 2) then
2495 is = (ntcw-1) * km
2496 do k = 1,kps
2497 do i = 1,im
2498 if(pcnvflg(i) .and. scuflg(i)) then
2499 kbx = 1
2500 kmx = max(kpbl(i), krad(i))
2501 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2502 kbx = 1
2503 kmx = kpbl(i)
2504 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2505 kbx = mrad(i)
2506 kmx = krad(i)
2507 endif
2508 if((pcnvflg(i) .or. scuflg(i)) .and.
2509 & (k >= kbx .and. k <= kmx)) then
2510 if(f2(i,k+is) < 0.) then
2511 tem = f2(i,k) + f2(i,k+is)
2512 if(tem >= 0.0) then
2513 f2(i,k) = tem
2514 f1(i,k) = f1(i,k) - elocp * f2(i,k+is)
2515 f2(i,k+is) = 0.
2516 elseif (f2(i,k) > 0.0) then
2517 f2(i,k+is) = tem
2518 f1(i,k) = f1(i,k) + elocp * f2(i,k)
2519 f2(i,k) = 0.
2520 endif
2521 endif
2522 endif
2523 enddo
2524 enddo
2525 endif
2526!
2527! For the negative rain water, first borrow water from vapor
2528! and then borrow it from the other layers if there is still
2529! negative water
2530!
2531 if(ntrac1 >= 2 .and. ntrw > 0) then
2532 is = (ntrw-1) * km
2533 do k = 1,kps
2534 do i = 1,im
2535 if(pcnvflg(i) .and. scuflg(i)) then
2536 kbx = 1
2537 kmx = max(kpbl(i), krad(i))
2538 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2539 kbx = 1
2540 kmx = kpbl(i)
2541 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2542 kbx = mrad(i)
2543 kmx = krad(i)
2544 endif
2545 if((pcnvflg(i) .or. scuflg(i)) .and.
2546 & (k >= kbx .and. k <= kmx)) then
2547 if(f2(i,k+is) < 0.) then
2548 tem = f2(i,k) + f2(i,k+is)
2549 if(tem >= 0.0) then
2550 f2(i,k) = tem
2551 f1(i,k) = f1(i,k) - elocp * f2(i,k+is)
2552 f2(i,k+is) = 0.
2553 elseif (f2(i,k) > 0.0) then
2554 f2(i,k+is) = tem
2555 f1(i,k) = f1(i,k) + elocp * f2(i,k)
2556 f2(i,k) = 0.
2557 endif
2558 endif
2559 endif
2560 enddo
2561 enddo
2562 endif
2563!
2564 if(ntrac1 >= 2) then
2565 do n = 2, ntrac1
2566 is = (n-1) * km
2567!
2568 do i = 1,im
2569 tsumn(i) = 0.
2570 tsump(i) = 0.
2571 rtnp(i) = 1.
2572 enddo
2573 do k = 1,kps
2574 do i = 1,im
2575 if(pcnvflg(i) .and. scuflg(i)) then
2576 kbx = 1
2577 kmx = max(kpbl(i), krad(i))
2578 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2579 kbx = 1
2580 kmx = kpbl(i)
2581 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2582 kbx = mrad(i)
2583 kmx = krad(i)
2584 endif
2585 if((pcnvflg(i) .or. scuflg(i)) .and.
2586 & (k >= kbx .and. k <= kmx)) then
2587 tem = f2(i,k+is) * del(i,k) * gravi
2588 if(f2(i,k+is) < 0.) tsumn(i) = tsumn(i) + tem
2589 if(f2(i,k+is) > 0.) tsump(i) = tsump(i) + tem
2590 endif
2591 enddo
2592 enddo
2593 do i = 1,im
2594 if(pcnvflg(i) .or. scuflg(i)) then
2595 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2596 if(tsump(i) > abs(tsumn(i))) then
2597 rtnp(i) = tsumn(i) / tsump(i)
2598 else
2599 rtnp(i) = tsump(i) / tsumn(i)
2600 endif
2601 endif
2602 endif
2603 enddo
2604 do k = 1,kps
2605 do i = 1,im
2606 if(pcnvflg(i) .and. scuflg(i)) then
2607 kbx = 1
2608 kmx = max(kpbl(i), krad(i))
2609 elseif(pcnvflg(i) .and. .not. scuflg(i)) then
2610 kbx = 1
2611 kmx = kpbl(i)
2612 elseif(.not. pcnvflg(i) .and. scuflg(i)) then
2613 kbx = mrad(i)
2614 kmx = krad(i)
2615 endif
2616 if((pcnvflg(i) .or. scuflg(i)) .and.
2617 & (k >= kbx .and. k <= kmx)) then
2618 if(rtnp(i) < 0.) then
2619 if(tsump(i) > abs(tsumn(i))) then
2620 if(f2(i,k+is)<0.) f2(i,k+is)=0.
2621 if(f2(i,k+is)>0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is)
2622 else
2623 if(f2(i,k+is)<0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is)
2624 if(f2(i,k+is)>0.) f2(i,k+is)=0.
2625 endif
2626 endif
2627 endif
2628 enddo
2629 enddo
2630!
2631! To remove negative hydrometeors & tracers which were leaked out of the mass-flux transport layers
2632! by eddy diffusion or potential negative hydrometeors & tracers from the diffusion scheme
2633! especially due to downward surface fluxes during nighttime,
2634! positive hydrometeors & tracers are borrowed again now from the entire layers
2635!
2636 do i = 1,im
2637 tsumn(i) = 0.
2638 tsump(i) = 0.
2639 rtnp(i) = 1.
2640 enddo
2641 do k = 1,km
2642 do i = 1,im
2643 tem = f2(i,k+is) * del(i,k) * gravi
2644 if(f2(i,k+is) < 0.) tsumn(i) = tsumn(i) + tem
2645 if(f2(i,k+is) > 0.) tsump(i) = tsump(i) + tem
2646 enddo
2647 enddo
2648 do i = 1,im
2649 if(tsump(i) > 0. .and. tsumn(i) < 0.) then
2650 if(tsump(i) > abs(tsumn(i))) then
2651 rtnp(i) = tsumn(i) / tsump(i)
2652 else
2653 rtnp(i) = tsump(i) / tsumn(i)
2654 endif
2655 endif
2656 enddo
2657 do k = 1,km
2658 do i = 1,im
2659 if(rtnp(i) < 0.) then
2660 if(tsump(i) > abs(tsumn(i))) then
2661 if(f2(i,k+is)<0.) f2(i,k+is)=0.
2662 if(f2(i,k+is)>0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is)
2663 else
2664 if(f2(i,k+is)<0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is)
2665 if(f2(i,k+is)>0.) f2(i,k+is)=0.
2666 endif
2667 endif
2668 enddo
2669 enddo
2670!
2671 enddo
2672 endif
2673c
2675c
2676 do k = 1,km
2677 do i = 1,im
2678 ttend = (f1(i,k)-t1(i,k))*rdt
2679 qtend = (f2(i,k)-q1(i,k,1))*rdt
2680 tdt(i,k) = tdt(i,k)+ttend
2681 rtg(i,k,1) = rtg(i,k,1)+qtend
2682! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend
2683! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend
2684 enddo
2685 enddo
2686!
2687 do i = 1,im
2688 dtsfc(i) = rho_a(i) * cp * heat(i)
2689 dqsfc(i) = rho_a(i) * hvap * evap(i)
2690 enddo
2691!
2692 if(ldiag3d .and. .not. gen_tend .and. present(dtend)) then
2693 idtend = dtidx(index_of_temperature,index_of_process_pbl)
2694 if(idtend>=1) then
2695 do k = 1,km
2696 do i = 1,im
2697 ttend = (f1(i,k)-t1(i,k))*rdt
2698 dtend(i,k,idtend) = dtend(i,k,idtend)+ttend*delt
2699 enddo
2700 enddo
2701 endif
2702 ! Send tendencies just for QV; other tracers are below.
2703 idtend = dtidx(100+ntqv,index_of_process_pbl)
2704 if(idtend>=1) then
2705 do k = 1,km
2706 do i = 1,im
2707 qtend = (f2(i,k)-q1(i,k,1))*rdt
2708 dtend(i,k,idtend) = dtend(i,k,idtend)+qtend*delt
2709 enddo
2710 enddo
2711 endif
2712 endif
2713!
2714 if(ntrac1 >= 2) then
2715 do n = 2, ntrac1
2716 is = (n-1) * km
2717 do k = 1, km
2718 do i = 1, im
2719 qtend = (f2(i,k+is)-q1(i,k,n))*rdt
2720 rtg(i,k,n) = rtg(i,k,n)+qtend
2721 enddo
2722 enddo
2723 enddo
2724 if(ldiag3d .and. .not. gen_tend .and. present(dtend)) then
2725 ! Send tendencies for all tracers that were selected.
2726 do n = 2, ntrac1
2727 is = (n-1) * km
2728 idtend = dtidx(n+100,index_of_process_pbl)
2729 if(idtend>=1) then
2730 if(n/=ntke) then
2731 do k = 1, km
2732 do i = 1, im
2733 qtend = (f2(i,k+is)-q1(i,k,n))*rdt
2734 dtend(i,k,idtend) = dtend(i,k,idtend)+qtend*delt
2735 enddo
2736 enddo
2737 endif
2738 endif
2739 enddo
2740 endif
2741!IVAI
2742! aux3d(:,:, 7) = rtg(:,:, ntoz) ! ntoz=7 "o3mr" GFS
2743
2744! aux3d(:,:, 5) = rtg(:,:, 11) ! n=11 "no3"
2745! aux3d(:,:, 5) = rtg(:,:, 9 ) ! n=9 "no"
2746! aux3d(:,:, 3) = rtg(:,:, 10) ! n=10 "o3"
2747! aux3d(:,:, 1) = rtg(:,:, 8 ) ! n=8 "no2"
2748!IVAI
2749 endif
2750!
2752!
2753 if(dspheat) then
2754 do k = 1,km1
2755 do i = 1,im
2756! tem = min(diss(i,k), dspmax)
2757! ttend = tem / cp
2758 ttend = diss(i,k) / cp
2759 tdt(i,k) = tdt(i,k) + dspfac * ttend
2760 enddo
2761 enddo
2762 if(ldiag3d .and. .not. gen_tend .and. present(dtend)) then
2763 idtend = dtidx(index_of_temperature,index_of_process_pbl)
2764 if(idtend>=1) then
2765 do k = 1,km1
2766 do i = 1,im
2767 ttend = diss(i,k) / cp
2768 dtend(i,k,idtend) = dtend(i,k,idtend)+dspfac*ttend*delt
2769 enddo
2770 enddo
2771 endif
2772 endif
2773 endif
2774c
2776c
2777 do i=1,im
2778 ad(i,1) = 1.0 + dtdz1(i) * stress(i) / spd1(i)
2779 f1(i,1) = u1(i,1)
2780 f2(i,1) = v1(i,1)
2781 enddo
2782c
2783 do k = 1,km1
2784 do i=1,im
2785 dtodsd = dt2/del(i,k)
2786 dtodsu = dt2/del(i,k+1)
2787 dsig = prsl(i,k)-prsl(i,k+1)
2788 rdz = rdzt(i,k)
2789 tem1 = dsig * dku(i,k) * rdz
2790 dsdz2 = tem1*rdz
2791 au(i,k) = -dtodsd*dsdz2
2792 al(i,k) = -dtodsu*dsdz2
2793 ad(i,k) = ad(i,k)-au(i,k)
2794 ad(i,k+1)= 1.-al(i,k)
2795 tem2 = dsig * rdz
2796!
2797 if(pcnvflg(i) .and. k < kpbl(i)) then
2798 ptem = 0.5 * tem2 * xmf(i,k)
2799 ptem1 = dtodsd * ptem
2800 ptem2 = dtodsu * ptem
2801 tem = u1(i,k) + u1(i,k+1)
2802 ptem = ucko(i,k) + ucko(i,k+1)
2803 f1(i,k) = f1(i,k) - (ptem - tem) * ptem1
2804 f1(i,k+1) = u1(i,k+1) + (ptem - tem) * ptem2
2805 tem = v1(i,k) + v1(i,k+1)
2806 ptem = vcko(i,k) + vcko(i,k+1)
2807 f2(i,k) = f2(i,k) - (ptem - tem) * ptem1
2808 f2(i,k+1) = v1(i,k+1) + (ptem - tem) * ptem2
2809 else
2810 f1(i,k+1) = u1(i,k+1)
2811 f2(i,k+1) = v1(i,k+1)
2812 endif
2813!
2814 if(scuflg(i)) then
2815 if(k >= mrad(i) .and. k < krad(i)) then
2816 ptem = 0.5 * tem2 * xmfd(i,k)
2817 ptem1 = dtodsd * ptem
2818 ptem2 = dtodsu * ptem
2819 tem = u1(i,k) + u1(i,k+1)
2820 ptem = ucdo(i,k) + ucdo(i,k+1)
2821 f1(i,k) = f1(i,k) + (ptem - tem) *ptem1
2822 f1(i,k+1) = f1(i,k+1) - (ptem - tem) *ptem2
2823 tem = v1(i,k) + v1(i,k+1)
2824 ptem = vcdo(i,k) + vcdo(i,k+1)
2825 f2(i,k) = f2(i,k) + (ptem - tem) * ptem1
2826 f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2
2827 endif
2828 endif
2829!
2830 enddo
2831 enddo
2832c
2834c
2835 call tridi2(im,km,al,ad,au,f1,f2,au,f1,f2)
2836c
2838c
2839 do k = 1,km
2840 do i = 1,im
2841 utend = (f1(i,k)-u1(i,k))*rdt
2842 vtend = (f2(i,k)-v1(i,k))*rdt
2843 du(i,k) = du(i,k)+utend
2844 dv(i,k) = dv(i,k)+vtend
2845! dusfc(i) = dusfc(i)+conw*del(i,k)*utend
2846! dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend
2847 enddo
2848 enddo
2849 if (use_oceanuv) then
2850 do i = 1,im
2851 spd1_m=sqrt( (u1(i,1)-usfco(i))**2+(v1(i,1)-vsfco(i))**2 )
2852 dusfc(i) = -1.*rho_a(i)*stress(i)*(u1(i,1)-usfco(i))/spd1_m
2853 dvsfc(i) = -1.*rho_a(i)*stress(i)*(v1(i,1)-vsfco(i))/spd1_m
2854 enddo
2855 else
2856 do i = 1,im
2857 dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i)
2858 dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i)
2859 enddo
2860 endif
2861!
2862 if(ldiag3d .and. .not. gen_tend .and. present(dtend)) then
2863 idtend = dtidx(index_of_x_wind,index_of_process_pbl)
2864 if(idtend>=1) then
2865 do k = 1,km
2866 do i = 1,im
2867 utend = (f1(i,k)-u1(i,k))*rdt
2868 dtend(i,k,idtend) = dtend(i,k,idtend) + utend*delt
2869 enddo
2870 enddo
2871 endif
2872
2873 idtend = dtidx(index_of_y_wind,index_of_process_pbl)
2874 if(idtend>=1) then
2875 do k = 1,km
2876 do i = 1,im
2877 vtend = (f2(i,k)-v1(i,k))*rdt
2878 dtend(i,k,idtend) = dtend(i,k,idtend) + vtend*delt
2879 enddo
2880 enddo
2881 endif
2882 endif
2883!
2884!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2886!
2887 do i = 1, im
2888 hpbl(i) = hpblx(i)
2889 kpbl(i) = kpblx(i)
2890 enddo
2891 if(sa3dtke) then
2892 do k = 1, km
2893 do i = 1, im
2894 dku3d_h(i,k) = dku_h(i,k) ! pass dku3d_h to dyn_core
2895 dku3d_e(i,k) = dkq_h(i,k) ! pass dku3d_e to dyn_core
2896 enddo
2897 enddo
2898 endif !sa3dtke
2899!
2900!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2901 return
2902 end subroutine satmedmfvdifq_can
2904 end module satmedmfvdifq_can_mod
satmedmfvdifq_can_mod::satmedmfvdifq_can
subroutine satmedmfvdifq_can(im, km, ntrac, ntcw, ntrw,
Definition satmedmfvdifq_can.F:51
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
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
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_can_mod
This file contains the CCPP-compliant SATMEDMF scheme (updated version) which computes subgrid vertic...
Definition satmedmfvdifq_can.F:10
tridi_mod
This module contains routine to compute tridiagonal matrix elements for TKE, heat,...
Definition tridi.f:6