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