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