CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
hedmf.f
1
4
7 module hedmf
8
9 use tridi_mod
10 use mfpbl_mod
11
12 contains
13
17 subroutine hedmf_init (hybedmf,moninq_fac,errmsg,errflg)
18 use machine, only : kind_phys
19 implicit none
20
21 logical, intent(in) :: hybedmf
22
23 real(kind=kind_phys), intent(in) :: moninq_fac
24 character(len=*), intent(out) :: errmsg
25 integer, intent(out) :: errflg
26 ! Initialize CCPP error handling variables
27 errmsg = ''
28 errflg = 0
29
30 ! Consistency checks
31 if (.not. hybedmf) then
32 errflg = 1
33 write(errmsg,'(*(a))') 'Logic error: hybedmf = .false.'
34 return
35 end if
36
37 if (moninq_fac == 0) then
38 errflg = 1
39 write(errmsg,'(*(a))') 'Logic error: moninq_fac == 0',
40 & ' is incompatible with hedmf'
41 end if
42 end subroutine hedmf_init
43
69 subroutine hedmf_run (im,km,ntrac,ntcw,dv,du,tau,rtg, &
70 & u1,v1,t1,q1,swh,hlw,xmu, &
71 & psk,rbsoil,zorl,u10m,v10m,fm,fh, &
72 & tsea,heat,evap,stress,spd1,kpbl, &
73 & prsi,del,prsl,prslk,phii,phil,delt,dspheat, &
74 & dusfc,dvsfc,dtsfc,dqsfc,hpbl,hgamt,hgamq,dkt,dku, &
75 & kinver,xkzm_m,xkzm_h,xkzm_s,lprnt,ipr, &
76 & xkzminv,moninq_fac,hurr_pbl,islimsk,var_ric, &
77 & coef_ric_l,coef_ric_s,ldiag3d,ntqv,rtg_ozone_index,ntoz, &
78 & dtend,dtidx,index_of_process_pbl,index_of_x_wind, &
79 & index_of_y_wind,index_of_temperature, &
80 & flag_for_pbl_generic_tend,errmsg,errflg)
81!
82 use machine , only : kind_phys
83 use funcphys , only : fpvs
84 !GJF: Note that sending these constants through the argument list
85 !results in regression test failures with "PROD" mode compilation
86 !flags (specifically, grav and cp)
87 use physcons, grav => con_g, cp => con_cp,
88 & hvap => con_hvap, fv => con_fvirt
89
90 implicit none
91!
92! arguments
93!
94 logical, intent(in) :: lprnt, hurr_pbl, ldiag3d
95 logical, intent(in) :: flag_for_pbl_generic_tend
96 integer, intent(in) :: ipr, islimsk(:)
97 integer, intent(in) :: im, km, ntrac, ntcw, kinver(:), ntoz
98 integer, intent(out) :: kpbl(:)
99
100!
101 real(kind=kind_phys), intent(in) :: delt, xkzm_m, xkzm_h, xkzm_s
102 real(kind=kind_phys), intent(in) :: xkzminv, moninq_fac, var_ric, &
103 & coef_ric_l, coef_ric_s
104 real(kind=kind_phys), intent(inout) :: dv(:,:), du(:,:), &
105 & tau(:,:), rtg(:,:,:)
106 ! dtend is only allocated if ldiag3d or qdiag3d are true
107 real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:)
108 integer, intent(in) :: dtidx(:,:)
109 integer, intent(in) :: index_of_x_wind, index_of_y_wind, &
110 & index_of_process_pbl, index_of_temperature, ntqv, rtg_ozone_index
111 real(kind=kind_phys), intent(in) :: &
112 & u1(:,:), v1(:,:), &
113 & t1(:,:), q1(:,:,:), &
114 & swh(:,:), hlw(:,:), &
115 & xmu(:), psk(:), &
116 & rbsoil(:), zorl(:), &
117 & u10m(:), v10m(:), &
118 & fm(:), fh(:), &
119 & tsea(:), &
120 & heat(:), evap(:), &
121 & stress(:), spd1(:)
122 real(kind=kind_phys), intent(in) :: &
123 & prsi(:,:), del(:,:), &
124 & prsl(:,:), prslk(:,:), &
125 & phii(:,:), phil(:,:)
126 real(kind=kind_phys), intent(out) :: &
127 & dusfc(:), dvsfc(:), &
128 & dtsfc(:), dqsfc(:), &
129 & hpbl(:)
130 real(kind=kind_phys), intent(out) :: &
131 & dkt(:,:), dku(:,:)
132 real(kind=kind_phys), intent(inout) :: &
133 & hgamt(:), hgamq(:)
134!
135 logical, intent(in) :: dspheat
136! flag for tke dissipative heating
137 character(len=*), intent(out) :: errmsg
138 integer, intent(out) :: errflg
139
140!
141! locals
142
143 integer i,iprt,is,iun,k,kk,km1,kmpbl,latd,lond
144 integer lcld(im),icld(im),kcld(im),krad(im)
145 integer kx1(im), kpblx(im)
146!
147! real(kind=kind_phys) betaq(im), betat(im), betaw(im),
148 real(kind=kind_phys) phih(im), phim(im), hpblx(im), &
149 & rbdn(im), rbup(im), &
150 & beta(im), sflux(im), &
151 & z0(im), crb(im), wstar(im), &
152 & zol(im), ustmin(im), ustar(im), &
153 & thermal(im),wscale(im), wscaleu(im)
154!
155 real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km), &
156 & qlx(im,km), thetae(im,km), &
157 & qtx(im,km), bf(im,km-1), diss(im,km), &
158 & radx(im,km-1), &
159 & govrth(im), hrad(im), &
160! & hradm(im), radmin(im), vrad(im), &
161 & radmin(im), vrad(im), &
162 & zd(im), zdd(im), thlvx1(im)
163!
164 real(kind=kind_phys) rdzt(im,km-1),dktx(im,km-1), &
165 & zi(im,km+1), zl(im,km), &
166 & xkzo(im,km-1), xkzmo(im,km-1), &
167 & cku(im,km-1), ckt(im,km-1), &
168 & ti(im,km-1), shr2(im,km-1), &
169 & al(im,km-1), ad(im,km), &
170 & au(im,km-1), a1(im,km), &
171 & a2(im,km*ntrac)
172!
173 real(kind=kind_phys) tcko(im,km), qcko(im,km,ntrac), &
174 & ucko(im,km), vcko(im,km), xmf(im,km)
175!
176 real(kind=kind_phys) prinv(im), rent(im)
177!
178 logical pblflg(im), sfcflg(im), scuflg(im), flg(im)
179 logical ublflg(im), pcnvflg(im)
180!
181! pcnvflg: true for convective(strongly unstable) pbl
182! ublflg: true for unstable but not convective(strongly unstable) pbl
183!
184 real(kind=kind_phys) aphi16, aphi5, bvf2, wfac,
185 & cfac, conq, cont, conw,
186 & dk, dkmax, dkmin,
187 & dq1, dsdz2, dsdzq, dsdzt,
188 & dsdzu, dsdzv,
189 & dsig, dt2, dthe1, dtodsd,
190 & dtodsu, dw2, dw2min,
191 & gamcrq, gamcrt, gocp,
192 & gravi, f0,
193 & prnum, prmax, prmin, pfac, crbcon,
194 & qmin, tdzmin, qtend, crbmin,crbmax,
195 & rbint, rdt, rdz, qlmin,
196 & ri, rimin, rl2, rlam, rlamun,
197 & rone, rzero, sfcfrac,
198 & spdk2, sri, zol1, zolcr, zolcru,
199 & robn, ttend,
200 & utend, vk, vk2,
201 & ust3, wst3,
202 & vtend, zfac, vpert, cteit,
203 & rentf1, rentf2, radfac,
204 & zfmin, zk, tem, tem1, tem2,
205 & xkzm, xkzmu,
206 & ptem, ptem1, ptem2, tx1(im), tx2(im)
207!
208 real(kind=kind_phys) zstblmax,h1, h2, qlcr, actei,
209 & cldtime
210 real :: ttend_fac
211
212 integer :: idtend1, idtend2
213
214 !! for hurricane application
215 real(kind=kind_phys) wspm(im,km-1)
216 integer kLOC ! RGF
217 real :: xDKU ! RGF
218
219 integer, parameter :: useshape=2!0-- no change, original ALPHA adjustment,1-- shape1, 2-- shape2(adjust above sfc)
220 real :: smax,ashape,sz2h, sksfc,skmax,ashape1,skminusk0, hmax
221cc
222 parameter(gravi=1.0/grav)
223 parameter(gocp=grav/cp)
224 parameter(cont=cp/grav,conq=hvap/grav,conw=1.0/grav) ! for del in pa
225! parameter(cont=1000.*cp/grav,conq=1000.*hvap/grav,conw=1000./grav) ! for del in kpa
226 parameter(rlam=30.0,vk=0.4,vk2=vk*vk)
227 parameter(prmin=0.25,prmax=4.,zolcr=0.2,zolcru=-0.5)
228 parameter(dw2min=0.0001,dkmin=0.0,dkmax=1000.,rimin=-100.)
229 parameter(crbcon=0.25,crbmin=0.15,crbmax=0.35)
230 parameter(wfac=7.0,cfac=6.5,pfac=2.0,sfcfrac=0.1)
231! parameter(qmin=1.e-8,xkzm=1.0,zfmin=1.e-8,aphi5=5.,aphi16=16.)
232 parameter(qmin=1.e-8, zfmin=1.e-8,aphi5=5.,aphi16=16.)
233 parameter(tdzmin=1.e-3,qlmin=1.e-12,f0=1.e-4)
234 parameter(h1=0.33333333,h2=0.66666667)
235! parameter(cldtime=500.,xkzminv=0.3)
236 parameter(cldtime=500.)
237! parameter(cldtime=500.,xkzmu=3.0,xkzminv=0.3)
238! parameter(gamcrt=3.,gamcrq=2.e-3,rlamun=150.0)
239 parameter(gamcrt=3.,gamcrq=0.,rlamun=150.0)
240 parameter(rentf1=0.2,rentf2=1.0,radfac=0.85)
241 parameter(iun=84)
242!
243! parameter (zstblmax = 2500., qlcr=1.0e-5)
244! parameter (zstblmax = 2500., qlcr=3.0e-5)
245! parameter (zstblmax = 2500., qlcr=3.5e-5)
246! parameter (zstblmax = 2500., qlcr=1.0e-4)
247 parameter(zstblmax = 2500., qlcr=3.5e-5)
248! parameter (actei = 0.23)
249 parameter(actei = 0.7)
250c
251c-----------------------------------------------------------------------
252c
253 601 format(1x,' moninp lat lon step hour ',3i6,f6.1)
254 602 format(1x,' k',' z',' t',' th',
255 1 ' tvh',' q',' u',' v',
256 2 ' sp')
257 603 format(1x,i5,8f9.1)
258 604 format(1x,' sfc',9x,f9.1,18x,f9.1)
259 605 format(1x,' k zl spd2 thekv the1v'
260 1 ,' thermal rbup')
261 606 format(1x,i5,6f8.2)
262 607 format(1x,' kpbl hpbl fm fh hgamt',
263 1 ' hgamq ws ustar cd ch')
264 608 format(1x,i5,9f8.2)
265 609 format(1x,' k pr dkt dku ',i5,3f8.2)
266 610 format(1x,' k pr dkt dku ',i5,3f8.2,' l2 ri t2',
267 1 ' sr2 ',2f8.2,2e10.2)
268! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269! Initialize CCPP error handling variables
270 errmsg = ''
271 errflg = 0
272
273! compute preliminary variables
274!
275! iprt = 0
276! if(iprt.eq.1) then
277!cc latd = 0
278! lond = 0
279! else
280!cc latd = 0
281! lond = 0
282! endif
283!
284 dt2 = delt
285 rdt = 1. / dt2
286 km1 = km - 1
287 kmpbl = km / 2
288
289 idtend1 = 0
290 idtend2 = 0
291
293 do k=1,km
294 do i=1,im
295 zi(i,k) = phii(i,k) * gravi
296 zl(i,k) = phil(i,k) * gravi
297 enddo
298 enddo
299 do i=1,im
300 zi(i,km+1) = phii(i,km+1) * gravi
301 enddo
303 do k = 1,km1
304 do i=1,im
305 rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k))
306 enddo
307 enddo
309 do i=1,im
310 kx1(i) = 1
311 tx1(i) = 1.0 / prsi(i,1)
312 tx2(i) = tx1(i)
313 enddo
315 do k = 1,km1
316 do i=1,im
317 xkzo(i,k) = 0.0
318 xkzmo(i,k) = 0.0
319 if (k < kinver(i)) then
320! vertical background diffusivity
321 ptem = prsi(i,k+1) * tx1(i)
322 tem1 = 1.0 - ptem
323 tem1 = tem1 * tem1 * 10.0
324 xkzo(i,k) = xkzm_h * min(1.0, exp(-tem1))
325
326! vertical background diffusivity for momentum
327 if (ptem >= xkzm_s) then
328 xkzmo(i,k) = xkzm_m
329 kx1(i) = k + 1
330 else
331 if (k == kx1(i) .and. k > 1) tx2(i) = 1.0 / prsi(i,k)
332 tem1 = 1.0 - prsi(i,k+1) * tx2(i)
333 tem1 = tem1 * tem1 * 5.0
334 xkzmo(i,k) = xkzm_m * min(1.0, exp(-tem1))
335 endif
336 endif
337 enddo
338 enddo
339! if (lprnt) then
340! print *,' xkzo=',(xkzo(ipr,k),k=1,km1)
341! print *,' xkzmo=',(xkzmo(ipr,k),k=1,km1)
342! endif
343!
344! diffusivity in the inversion layer is set to be xkzminv (m^2/s)
346 do k = 1,kmpbl
347 do i=1,im
348! if(zi(i,k+1) > 200..and.zi(i,k+1) < zstblmax) then
349 if(zi(i,k+1) > 250.) then
350 tem1 = (t1(i,k+1)-t1(i,k)) * rdzt(i,k)
351 if(tem1 > 1.e-5) then
352 xkzo(i,k) = min(xkzo(i,k),xkzminv)
353 endif
354 endif
355 enddo
356 enddo
358 do i = 1,im
359 z0(i) = 0.01 * zorl(i)
360 dusfc(i) = 0.
361 dvsfc(i) = 0.
362 dtsfc(i) = 0.
363 dqsfc(i) = 0.
364 wscale(i)= 0.
365 wscaleu(i)= 0.
366 kpbl(i) = 1
367 hpbl(i) = zi(i,1)
368 hpblx(i) = zi(i,1)
369 pblflg(i)= .true.
370 sfcflg(i)= .true.
371 if(rbsoil(i) > 0.) sfcflg(i) = .false.
372 ublflg(i)= .false.
373 pcnvflg(i)= .false.
374 scuflg(i)= .true.
375 if(scuflg(i)) then
376 radmin(i)= 0.
377 rent(i) = rentf1
378 hrad(i) = zi(i,1)
379! hradm(i) = zi(i,1)
380 krad(i) = 1
381 icld(i) = 0
382 lcld(i) = km1
383 kcld(i) = km1
384 zd(i) = 0.
385 endif
386 enddo
388 do k = 1,km
389 do i = 1,im
390 theta(i,k) = t1(i,k) * psk(i) / prslk(i,k)
391 qlx(i,k) = max(q1(i,k,ntcw),qlmin)
392 qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k)
393 ptem = qlx(i,k)
394 ptem1 = hvap*max(q1(i,k,1),qmin)/(cp*t1(i,k))
395 thetae(i,k)= theta(i,k)*(1.+ptem1)
396 thvx(i,k) = theta(i,k)*(1.+fv*max(q1(i,k,1),qmin)-ptem)
397 ptem2 = theta(i,k)-(hvap/cp)*ptem
398 thlvx(i,k) = ptem2*(1.+fv*qtx(i,k))
399 enddo
400 enddo
402 do k = 1,km
403 do i = 1,im
404 dku(i,k) = 0.
405 dkt(i,k) = 0.
406 enddo
407 enddo
408 do k = 1,km1
409 do i = 1,im
410 dktx(i,k) = 0.
411 cku(i,k) = 0.
412 ckt(i,k) = 0.
413 tem = zi(i,k+1)-zi(i,k)
414 radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k))
415 enddo
416 enddo
418 do i=1,im
419 flg(i) = scuflg(i)
420 enddo
421 do k = 1, km1
422 do i=1,im
423 if(flg(i).and.zl(i,k) >= zstblmax) then
424 lcld(i)=k
425 flg(i)=.false.
426 endif
427 enddo
428 enddo
429!
430! compute virtual potential temp gradient (bf) and winshear square
432 do k = 1, km1
433 do i = 1, im
434 rdz = rdzt(i,k)
435 bf(i,k) = (thvx(i,k+1)-thvx(i,k))*rdz
436 ti(i,k) = 2./(t1(i,k)+t1(i,k+1))
437 dw2 = (u1(i,k)-u1(i,k+1))**2
438 & + (v1(i,k)-v1(i,k+1))**2
439 shr2(i,k) = max(dw2,dw2min)*rdz*rdz
440 enddo
441 enddo
443 do i = 1,im
444 govrth(i) = grav/theta(i,1)
445 enddo
446!
447 do i=1,im
448 beta(i) = dt2 / (zi(i,2)-zi(i,1))
449 enddo
450!
451 do i=1,im
452 ustar(i) = sqrt(stress(i))
453 enddo
454!
455 do i = 1,im
456 sflux(i) = heat(i) + evap(i)*fv*theta(i,1)
457 if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false.
458 enddo
463! compute the pbl height
464!
465 if (.not. (hurr_pbl .and. moninq_fac < 0.0)) then
466 do i=1,im
467 flg(i) = .false.
468 rbup(i) = rbsoil(i)
469 !
470 if(pblflg(i)) then
471 thermal(i) = thvx(i,1)
472 crb(i) = crbcon
473 else
474 thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
475 tem = sqrt(u10m(i)**2+v10m(i)**2)
476 tem = max(tem, 1.)
477 robn = tem / (f0 * z0(i))
478 tem1 = 1.e-7 * robn
479 crb(i) = 0.16 * (tem1 ** (-0.18))
480 crb(i) = max(min(crb(i), crbmax), crbmin)
481 endif
482 enddo
483 else
484 do i=1,im
485 flg(i) = .false.
486 rbup(i) = rbsoil(i)
487
488 ! use variable Ri for all conditions
489 if(pblflg(i)) then
490 thermal(i) = thvx(i,1)
491 else
492 thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
493 endif
494 tem = sqrt(u10m(i)**2+v10m(i)**2)
495 tem = max(tem, 1.)
496 robn = tem / (f0 * z0(i))
497 tem1 = 1.e-7 * robn
498 crb(i) = crbcon
499 if (var_ric .eq. 1.) then
500 if (islimsk(i) .eq. 1) crb(i) = coef_ric_l*(tem1)**(-0.18)
501 if (islimsk(i) .eq. 0) crb(i) = coef_ric_s*(tem1)**(-0.18)
502 endif
503 crb(i) = max(min(crb(i), crbmax), crbmin)
504 enddo
505 endif
506
515 do k = 1, kmpbl
516 do i = 1, im
517 if(.not.flg(i)) then
518 rbdn(i) = rbup(i)
519 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
520 rbup(i) = (thvx(i,k)-thermal(i))*
521 & (grav*zl(i,k)/thvx(i,1))/spdk2
522 kpbl(i) = k
523 flg(i) = rbup(i) > crb(i)
524 endif
525 enddo
526 enddo
528 do i = 1,im
529 if(kpbl(i) > 1) then
530 k = kpbl(i)
531 if(rbdn(i) >= crb(i)) then
532 rbint = 0.
533 elseif(rbup(i) <= crb(i)) then
534 rbint = 1.
535 else
536 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
537 endif
538 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
539 if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1
540 else
541 hpbl(i) = zl(i,1)
542 kpbl(i) = 1
543 endif
544 kpblx(i) = kpbl(i)
545 hpblx(i) = hpbl(i)
546 enddo
547!
548! compute similarity parameters
564 do i=1,im
565 zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
566 if(sfcflg(i)) then
567 zol(i) = min(zol(i),-zfmin)
568 else
569 zol(i) = max(zol(i),zfmin)
570 endif
571 zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1)
572 if(sfcflg(i)) then
573! phim(i) = (1.-aphi16*zol1)**(-1./4.)
574! phih(i) = (1.-aphi16*zol1)**(-1./2.)
575 tem = 1.0 / (1. - aphi16*zol1)
576 phih(i) = sqrt(tem)
577 phim(i) = sqrt(phih(i))
578 else
579 phim(i) = 1. + aphi5*zol1
580 phih(i) = phim(i)
581 endif
582 wscale(i) = ustar(i)/phim(i)
583 ustmin(i) = ustar(i)/aphi5
584 wscale(i) = max(wscale(i),ustmin(i))
585 enddo
586 do i=1,im
587 if(pblflg(i)) then
588 if(zol(i) < zolcru .and. kpbl(i) > 1) then
589 pcnvflg(i) = .true.
590 else
591 ublflg(i) = .true.
592 endif
593 wst3 = govrth(i)*sflux(i)*hpbl(i)
594 wstar(i)= wst3**h1
595 ust3 = ustar(i)**3.
596 wscaleu(i) = (ust3+wfac*vk*wst3*sfcfrac)**h1
597 wscaleu(i) = max(wscaleu(i),ustmin(i))
598 endif
599 enddo
600!
601! compute counter-gradient mixing term for heat and moisture
604 do i = 1,im
605 if(ublflg(i)) then
606 hgamt(i) = min(cfac*heat(i)/wscaleu(i),gamcrt)
607 hgamq(i) = min(cfac*evap(i)/wscaleu(i),gamcrq)
608 vpert = hgamt(i) + hgamq(i)*fv*theta(i,1)
609 vpert = min(vpert,gamcrt)
610 thermal(i)= thermal(i)+max(vpert,0.)
611 hgamt(i) = max(hgamt(i),0.0)
612 hgamq(i) = max(hgamq(i),0.0)
613 endif
614 enddo
615!
616! enhance the pbl height by considering the thermal excess
618 do i=1,im
619 flg(i) = .true.
620 if(ublflg(i)) then
621 flg(i) = .false.
622 rbup(i) = rbsoil(i)
623 endif
624 enddo
625 do k = 2, kmpbl
626 do i = 1, im
627 if(.not.flg(i)) then
628 rbdn(i) = rbup(i)
629 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
630 rbup(i) = (thvx(i,k)-thermal(i))*
631 & (grav*zl(i,k)/thvx(i,1))/spdk2
632 kpbl(i) = k
633 flg(i) = rbup(i) > crb(i)
634 endif
635 enddo
636 enddo
637 do i = 1,im
638 if(ublflg(i)) then
639 k = kpbl(i)
640 if(rbdn(i) >= crb(i)) then
641 rbint = 0.
642 elseif(rbup(i) <= crb(i)) then
643 rbint = 1.
644 else
645 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
646 endif
647 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
648 if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1
649 if(kpbl(i) <= 1) then
650 ublflg(i) = .false.
651 pblflg(i) = .false.
652 endif
653 endif
654 enddo
655!
656! look for stratocumulus
659 do i = 1, im
660 flg(i)=scuflg(i)
661 enddo
662 do k = kmpbl,1,-1
663 do i = 1, im
664 if(flg(i) .and. k <= lcld(i)) then
665 if(qlx(i,k).ge.qlcr) then
666 kcld(i)=k
667 flg(i)=.false.
668 endif
669 endif
670 enddo
671 enddo
672 do i = 1, im
673 if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false.
674 enddo
676 do i = 1, im
677 flg(i)=scuflg(i)
678 enddo
679 do k = kmpbl,1,-1
680 do i = 1, im
681 if(flg(i) .and. k <= kcld(i)) then
682 if(qlx(i,k) >= qlcr) then
683 if(radx(i,k) < radmin(i)) then
684 radmin(i)=radx(i,k)
685 krad(i)=k
686 endif
687 else
688 flg(i)=.false.
689 endif
690 endif
691 enddo
692 enddo
693 do i = 1, im
694 if(scuflg(i) .and. krad(i) <= 1) scuflg(i)=.false.
695 if(scuflg(i) .and. radmin(i)>=0.) scuflg(i)=.false.
696 enddo
698 do i = 1, im
699 flg(i)=scuflg(i)
700 enddo
701 do k = kmpbl,2,-1
702 do i = 1, im
703 if(flg(i) .and. k <= krad(i)) then
704 if(qlx(i,k) >= qlcr) then
705 icld(i)=icld(i)+1
706 else
707 flg(i)=.false.
708 endif
709 endif
710 enddo
711 enddo
712 do i = 1, im
713 if(scuflg(i) .and. icld(i) < 1) scuflg(i)=.false.
714 enddo
716 do i = 1, im
717 if(scuflg(i)) then
718 hrad(i) = zi(i,krad(i)+1)
719! hradm(i)= zl(i,krad(i))
720 endif
721 enddo
722!
723 do i = 1, im
724 if(scuflg(i) .and. hrad(i)<zi(i,2)) scuflg(i)=.false.
725 enddo
727 do i = 1, im
728 if(scuflg(i)) then
729 k = krad(i)
730 tem = zi(i,k+1)-zi(i,k)
731 tem1 = cldtime*radmin(i)/tem
732 thlvx1(i) = thlvx(i,k)+tem1
733! if(thlvx1(i) > thlvx(i,k-1)) scuflg(i)=.false.
734 endif
735 enddo
737 do i = 1, im
738 flg(i)=scuflg(i)
739 enddo
740 do k = kmpbl,1,-1
741 do i = 1, im
742 if(flg(i) .and. k <= krad(i))then
743 if(thlvx1(i) <= thlvx(i,k))then
744 tem=zi(i,k+1)-zi(i,k)
745 zd(i)=zd(i)+tem
746 else
747 flg(i)=.false.
748 endif
749 endif
750 enddo
751 enddo
753 do i = 1, im
754 if(scuflg(i))then
755 kk = max(1, krad(i)+1-icld(i))
756 zdd(i) = hrad(i)-zi(i,kk)
757 endif
758 enddo
760 do i = 1, im
761 if(scuflg(i))then
762 zd(i) = max(zd(i),zdd(i))
763 zd(i) = min(zd(i),hrad(i))
764 tem = govrth(i)*zd(i)*(-radmin(i))
765 vrad(i)= tem**h1
766 endif
767 enddo
768!
769! compute inverse prandtl number
772 do i = 1, im
773 if(ublflg(i)) then
774 tem = phih(i)/phim(i)+cfac*vk*sfcfrac
775 else
776 tem = phih(i)/phim(i)
777 endif
778 prinv(i) = 1.0 / tem
779 prinv(i) = min(prinv(i),prmax)
780 prinv(i) = max(prinv(i),prmin)
781 enddo
782 do i = 1, im
783 if(zol(i) > zolcr) then
784 kpbl(i) = 1
785 endif
786 enddo
787
788
789!!! 20150915 WeiguoWang added alpha (moninq_fac) and wind-dependent modification of K by RGF
790! -------------------------------------------------------------------------------------
791! begin RGF modifications
792! this is version MOD05
793
794! RGF determine wspd at roughly 500 m above surface, or as close as possible,
795! reuse SPDK2
796! zi(i,k) is AGL, right? May not matter if applied only to water grid points
797 if(hurr_pbl .and. moninq_fac < 0.0) then
798 do i=1,im
799 spdk2 = 0.
800 wspm(i,1) = 0.
801 do k = 1, kmpbl ! kmpbl is like a max possible pbl height
802 if (zi(i,k) .le. 500. .and. zi(i,k+1) .gt. 500.) then ! find level bracketing 500 m
803 spdk2 = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)) ! wspd near 500 m
804 wspm(i,1) = spdk2/0.6 ! now the Km limit for 500 m. just store in K=1
805 wspm(i,2) = float(k) ! height of level at gridpoint i. store in K=2
806 endif
807 enddo !k
808 enddo ! i
809 endif ! hurr_pbl and moninq_fac < 0
810
811
812! compute diffusion coefficients below pbl
815 if (.not. (hurr_pbl .and. moninq_fac < 0.0)) then
816 do k = 1, kmpbl
817 do i=1,im
818 if(k < kpbl(i)) then
819! zfac = max((1.-(zi(i,k+1)-zl(i,1))/
820! 1 (hpbl(i)-zl(i,1))), zfmin)
821 zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin)
822 tem = zi(i,k+1) * (zfac**pfac) * moninq_fac ! lmh suggested by kg
823 if(pblflg(i)) then
824 tem1 = vk * wscaleu(i) * tem
825! dku(i,k) = xkzmo(i,k) + tem1
826! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i)
827 dku(i,k) = tem1
828 dkt(i,k) = tem1 * prinv(i)
829 else
830 tem1 = vk * wscale(i) * tem
831! dku(i,k) = xkzmo(i,k) + tem1
832! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i)
833 dku(i,k) = tem1
834 dkt(i,k) = tem1 * prinv(i)
835 endif
836 dku(i,k) = min(dku(i,k),dkmax)
837 dku(i,k) = max(dku(i,k),xkzmo(i,k))
838 dkt(i,k) = min(dkt(i,k),dkmax)
839 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
840 dktx(i,k)= dkt(i,k)
841 endif
842 enddo !i
843 enddo !k
844 else
845 !hurricane PBL case and moninq_fac < 0 (note that the i and k loop order has been switched)
846 do i=1, im
847 do k=1, kmpbl
848 if (k < kpbl(i)) then
849! zfac = max((1.-(zi(i,k+1)-zl(i,1))/
850! 1 (hpbl(i)-zl(i,1))), zfmin)
851 zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin)
852 tem = zi(i,k+1) * (zfac**pfac) * abs(moninq_fac)
853
854!!!! CHANGES FOR HEIGHT-DEPENDENT K ADJUSTMENT, WANG W
855 if (useshape .ge. 1) then
856 sz2h=(zi(i,k+1)-zl(i,1))/(hpbl(i)-zl(i,1))
857 sz2h=max(sz2h,zfmin)
858 sz2h=min(sz2h,1.0)
859 zfac=(1.0-sz2h)**pfac
860! smax=0.148 !! max value of this shape function
861 smax=0.148 !! max value of this shape function
862 hmax=0.333 !! roughly height if max K
863 skmax=hmax*(1.0-hmax)**pfac
864 sksfc=min(zi(i,2)/hpbl(i),0.05) ! surface layer top, 0.05H or ZI(2) (Zi(1)=0)
865 sksfc=sksfc*(1-sksfc)**pfac
866
867 zfac=max(zfac,zfmin)
868 ashape=max(abs(moninq_fac),0.2) ! should not be smaller than 0.2, otherwise too much adjustment(?)
869 if (useshape == 1) then
870 ashape=(1.0 - ((sz2h*zfac/smax)**0.25) *(1.0 - ashape))
871 tem = zi(i,k+1) * (zfac) * ashape
872 elseif (useshape == 2) then !only adjus K that is > K_surface_top
873 ashape1=1.0
874 if (skmax > sksfc) then
875 ashape1=(skmax*ashape-sksfc)/(skmax-sksfc)
876 endif
877 skminusk0 = zi(i,k+1)*zfac - hpbl(i)*sksfc
878 tem = zi(i,k+1) * (zfac) ! no adjustment
879 if (skminusk0 > 0) then ! only adjust K which is > surface top K
880 tem = skminusk0*ashape1 + hpbl(i)*sksfc
881 endif
882 endif ! useshape == 1 or 2
883 endif ! endif useshape>1
884!!!! END OF CHANGES , WANG W
885
886!!If alpha >= 0, this is the only modification of K
887! if alpha = -1, the above provides the first guess for DKU, based on assumption
888! alpha = +1
889! (other values of alpha < 0 can also be applied)
890! if alpha > 0, the above applies the alpha suppression factor and we are
891! finished
892
893 if(pblflg(i)) then
894 tem1 = vk * wscaleu(i) * tem
895! dku(i,k) = xkzmo(i,k) + tem1
896! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i)
897 dku(i,k) = tem1
898 dkt(i,k) = tem1 * prinv(i)
899 else
900 tem1 = vk * wscale(i) * tem
901! dku(i,k) = xkzmo(i,k) + tem1
902! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i)
903 dku(i,k) = tem1
904 dkt(i,k) = tem1 * prinv(i)
905 endif
906 dku(i,k) = min(dku(i,k),dkmax)
907 dku(i,k) = max(dku(i,k),xkzmo(i,k))
908 dkt(i,k) = min(dkt(i,k),dkmax)
909 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
910 dktx(i,k)= dkt(i,k)
911 endif !k < kpbl(i)
912 enddo !K loop
913
914! possible modification of first guess DKU, under certain conditions
915! (1) this applies only to columns over water
916 if (islimsk(i) .eq. 0) then ! sea only
917! (2) alpha test
918! if alpha < 0, find alpha for each column and do the loop again
919! if alpha > 0, we are finished
920
921!GJF: redundant check for moninq_fac < 0?
922 if (moninq_fac .lt. 0.) then ! variable alpha test
923! k-level of layer around 500 m
924 kloc = int(wspm(i,2))
925! print *,' kLOC ',kLOC,' KPBL ',KPBL(I)
926
927! (3) only do this IF KPBL(I) >= kLOC. Otherwise, we are finished, with DKU as
928! if alpha = +1
929 if(kpbl(i) .gt. kloc) then
930 xdku = dku(i,kloc) ! Km at k-level
931! (4) DKU check.
932! WSPM(i,1) is the KM cap for the 500-m level.
933! if DKU at 500-m level < WSPM(i,1), do not limit Km ANYWHERE. Alpha =
934! abs(alpha). No need to recalc.
935! if DKU at 500-m level > WSPM(i,1), then alpha = WSPM(i,1)/xDKU for entire
936! column
937 if(xdku .ge. wspm(i,1)) then ! ONLY if DKU at 500-m exceeds cap, otherwise already done
938 wspm(i,3) = wspm(i,1)/xdku ! ratio of cap to Km at k-level, store in WSPM(i,3)
939 !WSPM(i,4) = amin1(WSPM(I,3),1.0) ! this is new column alpha. cap at 1. ! should never be needed
940 wspm(i,4) = min(wspm(i,3),1.0) ! this is new column alpha. cap at 1. ! should never be needed
941 !! recalculate K capped by WSPM(i,1)
942 do k = 1, kmpbl
943 if(k < kpbl(i)) then
944! zfac = max((1.-(zi(i,k+1)-zl(i,1))/
945! 1 (hpbl(i)-zl(i,1))), zfmin)
946 zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin)
947 tem = zi(i,k+1) * (zfac**pfac) * wspm(i,4)
948!!! wang use different K shape, options!!!!!!!!!!!!!!!!!!!!!!!!!
949!!!! HANGES FOR HEIGHT-DEPENDENT K ADJUSTMENT, WANG W
950 if(useshape .ge. 1) then
951 sz2h=(zi(i,k+1)-zl(i,1))/(hpbl(i)-zl(i,1))
952 sz2h=max(sz2h,zfmin)
953 sz2h=min(sz2h,1.0)
954 zfac=(1.0-sz2h)**pfac
955 smax=0.148 !! max value of this shape function
956 hmax=0.333 !! roughly height if max K
957 skmax=hmax*(1.0-hmax)**pfac
958 sksfc=min(zi(i,2)/hpbl(i),0.05) ! surface layer top, 0.05H or ZI(2) (Zi(1)=0)
959 sksfc=sksfc*(1-sksfc)**pfac
960
961 zfac=max(zfac,zfmin)
962 ashape=max(wspm(i,4),0.2) !! adjustment coef should not smaller than 0.2
963 if(useshape ==1) then
964 ashape=(1.0 - ((sz2h*zfac/smax)**0.25)*
965 & (1.0 - ashape))
966 tem = zi(i,k+1) * (zfac) * ashape
967 elseif (useshape == 2) then !only adjus K that is > K_surface_top
968 ashape1=1.0
969 if (skmax > sksfc) then
970 ashape1=(skmax*ashape-sksfc)/(skmax-sksfc)
971 endif
972 skminusk0=zi(i,k+1)*zfac - hpbl(i)*sksfc
973 tem = zi(i,k+1) * (zfac) ! no adjustment
974 if (skminusk0 > 0) then ! only adjust K which is > surface top K
975 tem = skminusk0*ashape1 + hpbl(i)*sksfc
976 endif
977 endif ! endif useshape=1 or 2
978 endif ! endif useshape>1
979!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
980 if(pblflg(i)) then
981 tem1 = vk * wscaleu(i) * tem
982! dku(i,k) = xkzmo(i,k) + tem1
983! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i)
984 dku(i,k) = tem1
985 dkt(i,k) = tem1 * prinv(i)
986 else
987 tem1 = vk * wscale(i) * tem
988! dku(i,k) = xkzmo(i,k) + tem1
989! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i)
990 dku(i,k) = tem1
991 dkt(i,k) = tem1 * prinv(i)
992 endif !pblflg
993 dku(i,k) = min(dku(i,k),dkmax)
994 dku(i,k) = max(dku(i,k),xkzmo(i,k))
995 dkt(i,k) = min(dkt(i,k),dkmax)
996 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
997 dktx(i,k)= dkt(i,k)
998 endif ! k < kpbl(i)
999 enddo ! K loop
1000 endif ! xDKU .ge. wspm(i,1)
1001 endif ! kpbl(i) .ge. kLOC
1002 endif ! moninq_fac < 0 (GJF: redundant?)
1003 endif ! islimsk == 0
1004 enddo ! I loop
1005 endif ! not (hurr_pbl and moninq_fac < 0)
1006!
1007! compute diffusion coefficients based on local scheme above pbl
1042 do k = 1, km1
1043 do i=1,im
1044 if(k >= kpbl(i)) then
1045 bvf2 = grav*bf(i,k)*ti(i,k)
1046 ri = max(bvf2/shr2(i,k),rimin)
1047 zk = vk*zi(i,k+1)
1048 if(ri < 0.) then ! unstable regime
1049 rl2 = zk*rlamun/(rlamun+zk)
1050 dk = rl2*rl2*sqrt(shr2(i,k))
1051 sri = sqrt(-ri)
1052! dku(i,k) = xkzmo(i,k) + dk*(1+8.*(-ri)/(1+1.746*sri))
1053! dkt(i,k) = xkzo(i,k) + dk*(1+8.*(-ri)/(1+1.286*sri))
1054 dku(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
1055 dkt(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
1056 else ! stable regime
1057 rl2 = zk*rlam/(rlam+zk)
1058!! tem = rlam * sqrt(0.01*prsi(i,k))
1059!! rl2 = zk*tem/(tem+zk)
1060 dk = rl2*rl2*sqrt(shr2(i,k))
1061 tem1 = dk/(1+5.*ri)**2
1062!
1063 if(k >= kpblx(i)) then
1064 prnum = 1.0 + 2.1*ri
1065 prnum = min(prnum,prmax)
1066 else
1067 prnum = 1.0
1068 endif
1069! dku(i,k) = xkzmo(i,k) + tem1 * prnum
1070! dkt(i,k) = xkzo(i,k) + tem1
1071 dku(i,k) = tem1 * prnum
1072 dkt(i,k) = tem1
1073 endif
1074!
1075 dku(i,k) = min(dku(i,k),dkmax)
1076 dku(i,k) = max(dku(i,k),xkzmo(i,k))
1077 dkt(i,k) = min(dkt(i,k),dkmax)
1078 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
1079!
1080 endif
1081!
1082 enddo
1083 enddo
1084!
1085!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1086! compute components for mass flux mixing by large thermals
1089 do k = 1, km
1090 do i = 1, im
1091 if(pcnvflg(i)) then
1092 tcko(i,k) = t1(i,k)
1093 ucko(i,k) = u1(i,k)
1094 vcko(i,k) = v1(i,k)
1095 xmf(i,k) = 0.
1096 endif
1097 enddo
1098 enddo
1099 do kk = 1, ntrac
1100 do k = 1, km
1101 do i = 1, im
1102 if(pcnvflg(i)) then
1103 qcko(i,k,kk) = q1(i,k,kk)
1104 endif
1105 enddo
1106 enddo
1107 enddo
1109 call mfpbl(im,im,km,ntrac,dt2,pcnvflg,
1110 & zl,zi,thvx,q1,t1,u1,v1,hpbl,kpbl,
1111 & sflux,ustar,wstar,xmf,tcko,qcko,ucko,vcko)
1112!
1113!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1114! compute diffusion coefficients for cloud-top driven diffusion
1115! if the condition for cloud-top instability is met,
1116! increase entrainment flux at cloud top
1117!
1128 do i = 1, im
1129 if(scuflg(i)) then
1130 k = krad(i)
1131 tem = thetae(i,k) - thetae(i,k+1)
1132 tem1 = qtx(i,k) - qtx(i,k+1)
1133 if (tem > 0. .and. tem1 > 0.) then
1134 cteit= cp*tem/(hvap*tem1)
1135 if(cteit > actei) rent(i) = rentf2
1136 endif
1137 endif
1138 enddo
1139 do i = 1, im
1140 if(scuflg(i)) then
1141 k = krad(i)
1142 tem1 = max(bf(i,k),tdzmin)
1143 ckt(i,k) = -rent(i)*radmin(i)/tem1
1144 cku(i,k) = ckt(i,k)
1145 endif
1146 enddo
1147!
1148 do k = 1, kmpbl
1149 do i=1,im
1150 if(scuflg(i) .and. k < krad(i)) then
1151 tem1=hrad(i)-zd(i)
1152 tem2=zi(i,k+1)-tem1
1153 if(tem2 > 0.) then
1154 ptem= tem2/zd(i)
1155 if(ptem.ge.1.) ptem= 1.
1156 ptem= tem2*ptem*sqrt(1.-ptem)
1157 ckt(i,k) = radfac*vk*vrad(i)*ptem
1158 cku(i,k) = 0.75*ckt(i,k)
1159 ckt(i,k) = max(ckt(i,k),dkmin)
1160 ckt(i,k) = min(ckt(i,k),dkmax)
1161 cku(i,k) = max(cku(i,k),dkmin)
1162 cku(i,k) = min(cku(i,k),dkmax)
1163 endif
1164 endif
1165 enddo
1166 enddo
1167!
1168!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1169!
1171 if (.not. hurr_pbl) then
1172 do k = 1, kmpbl
1173 do i=1,im
1174 if(scuflg(i)) then
1175 dkt(i,k) = dkt(i,k)+ckt(i,k)
1176 dku(i,k) = dku(i,k)+cku(i,k)
1177 dkt(i,k) = min(dkt(i,k),dkmax)
1178 dku(i,k) = min(dku(i,k),dkmax)
1179 endif
1180 enddo
1181 enddo
1182 else
1183 do k = 1, kmpbl
1184 do i=1,im
1185 if(scuflg(i)) then
1186 !! if K needs to be adjusted by alpha, then no need to add this term
1187 if (.not. (hurr_pbl .and. moninq_fac < 0.0)) then
1188 dkt(i,k) = dkt(i,k)+ckt(i,k)
1189 dku(i,k) = dku(i,k)+cku(i,k)
1190 end if
1191 dkt(i,k) = min(dkt(i,k),dkmax)
1192 dku(i,k) = min(dku(i,k),dkmax)
1193 endif
1194 enddo
1195 enddo
1196 endif
1197!
1198! compute tridiagonal matrix elements for heat and moisture
1199!
1202 do i=1,im
1203 ad(i,1) = 1.
1204 a1(i,1) = t1(i,1) + beta(i) * heat(i)
1205 a2(i,1) = q1(i,1,1) + beta(i) * evap(i)
1206 enddo
1207
1208 if(ntrac >= 2) then
1209 do k = 2, ntrac
1210 is = (k-1) * km
1211 do i = 1, im
1212 a2(i,1+is) = q1(i,1,k)
1213 enddo
1214 enddo
1215 endif
1216!
1217 do k = 1,km1
1218 do i = 1,im
1219 dtodsd = dt2/del(i,k)
1220 dtodsu = dt2/del(i,k+1)
1221 dsig = prsl(i,k)-prsl(i,k+1)
1222 rdz = rdzt(i,k)
1223 tem1 = dsig * dkt(i,k) * rdz
1224 dsdz2 = tem1 * rdz
1225 au(i,k) = -dtodsd*dsdz2
1226 al(i,k) = -dtodsu*dsdz2
1227!
1228 if(pcnvflg(i) .and. k < kpbl(i)) then
1229 tem2 = dsig * rdz
1230 ptem = 0.5 * tem2 * xmf(i,k)
1231 ptem1 = dtodsd * ptem
1232 ptem2 = dtodsu * ptem
1233 ad(i,k) = ad(i,k)-au(i,k)-ptem1
1234 ad(i,k+1) = 1.-al(i,k)+ptem2
1235 au(i,k) = au(i,k)-ptem1
1236 al(i,k) = al(i,k)+ptem2
1237 ptem = tcko(i,k) + tcko(i,k+1)
1238 dsdzt = tem1 * gocp
1239 a1(i,k) = a1(i,k)+dtodsd*dsdzt-ptem1*ptem
1240 a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+ptem2*ptem
1241 ptem = qcko(i,k,1) + qcko(i,k+1,1)
1242 a2(i,k) = a2(i,k) - ptem1 * ptem
1243 a2(i,k+1) = q1(i,k+1,1) + ptem2 * ptem
1244 elseif(ublflg(i) .and. k < kpbl(i)) then
1245 ptem1 = dsig * dktx(i,k) * rdz
1246 tem = 1.0 / hpbl(i)
1247 dsdzt = tem1 * gocp - ptem1 * hgamt(i) * tem
1248 dsdzq = - ptem1 * hgamq(i) * tem
1249 ad(i,k) = ad(i,k)-au(i,k)
1250 ad(i,k+1) = 1.-al(i,k)
1251 a1(i,k) = a1(i,k)+dtodsd*dsdzt
1252 a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
1253 a2(i,k) = a2(i,k)+dtodsd*dsdzq
1254 a2(i,k+1) = q1(i,k+1,1)-dtodsu*dsdzq
1255 else
1256 ad(i,k) = ad(i,k)-au(i,k)
1257 ad(i,k+1) = 1.-al(i,k)
1258 dsdzt = tem1 * gocp
1259 a1(i,k) = a1(i,k)+dtodsd*dsdzt
1260 a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
1261 a2(i,k+1) = q1(i,k+1,1)
1262 endif
1263!
1264 enddo
1265 enddo
1266!
1267 if(ntrac >= 2) then
1268 do kk = 2, ntrac
1269 is = (kk-1) * km
1270 do k = 1, km1
1271 do i = 1, im
1272 if(pcnvflg(i) .and. k < kpbl(i)) then
1273 dtodsd = dt2/del(i,k)
1274 dtodsu = dt2/del(i,k+1)
1275 dsig = prsl(i,k)-prsl(i,k+1)
1276 tem = dsig * rdzt(i,k)
1277 ptem = 0.5 * tem * xmf(i,k)
1278 ptem1 = dtodsd * ptem
1279 ptem2 = dtodsu * ptem
1280 tem1 = qcko(i,k,kk) + qcko(i,k+1,kk)
1281 a2(i,k+is) = a2(i,k+is) - ptem1*tem1
1282 a2(i,k+1+is)= q1(i,k+1,kk) + ptem2*tem1
1283 else
1284 a2(i,k+1+is) = q1(i,k+1,kk)
1285 endif
1286 enddo
1287 enddo
1288 enddo
1289 endif
1290!
1291! solve tridiagonal problem for heat and moisture
1292!
1294 call tridin(im,km,ntrac,al,ad,au,a1,a2,au,a1,a2)
1295
1296!
1297! recover tendencies of heat and moisture
1298!
1300 do k = 1,km
1301 do i = 1,im
1302 ttend = (a1(i,k)-t1(i,k)) * rdt
1303 qtend = (a2(i,k)-q1(i,k,1))*rdt
1304 tau(i,k) = tau(i,k)+ttend
1305 rtg(i,k,1) = rtg(i,k,1)+qtend
1306 dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend
1307 dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend
1308 enddo
1309 enddo
1310 if(.not.flag_for_pbl_generic_tend) then
1311 idtend1=dtidx(index_of_temperature,index_of_process_pbl)
1312 idtend2=dtidx(ntqv+100,index_of_process_pbl)
1313 if(idtend1>=1) then
1314 do k = 1,km
1315 do i = 1,im
1316 ttend = (a1(i,k)-t1(i,k)) * rdt
1317 dtend(i,k,idtend1) = dtend(i,k,idtend1) + ttend*delt
1318 enddo
1319 enddo
1320 endif
1321 if(idtend2>=1) then
1322 do k = 1,km
1323 do i = 1,im
1324 qtend = (a2(i,k)-q1(i,k,1))*rdt
1325 dtend(i,k,idtend2) = dtend(i,k,idtend2) + qtend*delt
1326 enddo
1327 enddo
1328 endif
1329 endif
1330 if(ntrac >= 2) then
1331 do kk = 2, ntrac
1332 is = (kk-1) * km
1333 do k = 1, km
1334 do i = 1, im
1335 qtend = (a2(i,k+is)-q1(i,k,kk))*rdt
1336 rtg(i,k,kk) = rtg(i,k,kk)+qtend
1337 enddo
1338 enddo
1339 enddo
1340 if(.not.flag_for_pbl_generic_tend .and. ldiag3d .and. &
1341 & rtg_ozone_index>0) then
1342 idtend1 = dtidx(100+ntoz,index_of_process_pbl)
1343 if(idtend1>=1) then
1344 kk = rtg_ozone_index
1345 is = (kk-1) * km
1346 do k = 1, km
1347 do i = 1, im
1348 qtend = (a2(i,k+is)-q1(i,k,kk))
1349 dtend(i,k,idtend1) = dtend(i,k,idtend1)+qtend
1350 enddo
1351 enddo
1352 endif
1353 endif
1354 endif
1355!
1356! compute tke dissipation rate
1357!
1360 if(dspheat) then
1361!
1362 do k = 1,km1
1363 do i = 1,im
1364 diss(i,k) = dku(i,k)*shr2(i,k)-grav*ti(i,k)*dkt(i,k)*bf(i,k)
1365! diss(i,k) = dku(i,k)*shr2(i,k)
1366 enddo
1367 enddo
1368!
1369! add dissipative heating at the first model layer
1370!
1372 if (hurr_pbl .and. moninq_fac < 0.0) then
1373 ttend_fac = 0.7
1374 else
1375 ttend_fac = 0.5
1376 endif
1377
1378 do i = 1,im
1379 tem = govrth(i)*sflux(i)
1380 tem1 = tem + stress(i)*spd1(i)/zl(i,1)
1381 tem2 = 0.5 * (tem1+diss(i,1))
1382 tem2 = max(tem2, 0.)
1383 ttend = tem2 / cp
1384 tau(i,1) = tau(i,1)+ttend_fac*ttend
1385 enddo
1386!
1387! add dissipative heating above the first model layer
1388!
1389 do k = 2,km1
1390 do i = 1,im
1391 tem = 0.5 * (diss(i,k-1)+diss(i,k))
1392 tem = max(tem, 0.)
1393 ttend = tem / cp
1394 tau(i,k) = tau(i,k) + ttend_fac*ttend
1395 enddo
1396 enddo
1397!
1398 endif
1399!
1400! compute tridiagonal matrix elements for momentum
1401!
1404 do i=1,im
1405 ad(i,1) = 1.0 + beta(i) * stress(i) / spd1(i)
1406 a1(i,1) = u1(i,1)
1407 a2(i,1) = v1(i,1)
1408 enddo
1409!
1410 do k = 1,km1
1411 do i=1,im
1412 dtodsd = dt2/del(i,k)
1413 dtodsu = dt2/del(i,k+1)
1414 dsig = prsl(i,k)-prsl(i,k+1)
1415 rdz = rdzt(i,k)
1416 tem1 = dsig*dku(i,k)*rdz
1417 dsdz2 = tem1 * rdz
1418 au(i,k) = -dtodsd*dsdz2
1419 al(i,k) = -dtodsu*dsdz2
1420!
1421 if(pcnvflg(i) .and. k < kpbl(i)) then
1422 tem2 = dsig * rdz
1423 ptem = 0.5 * tem2 * xmf(i,k)
1424 ptem1 = dtodsd * ptem
1425 ptem2 = dtodsu * ptem
1426 ad(i,k) = ad(i,k)-au(i,k)-ptem1
1427 ad(i,k+1) = 1.-al(i,k)+ptem2
1428 au(i,k) = au(i,k)-ptem1
1429 al(i,k) = al(i,k)+ptem2
1430 ptem = ucko(i,k) + ucko(i,k+1)
1431 a1(i,k) = a1(i,k) - ptem1 * ptem
1432 a1(i,k+1) = u1(i,k+1) + ptem2 * ptem
1433 ptem = vcko(i,k) + vcko(i,k+1)
1434 a2(i,k) = a2(i,k) - ptem1 * ptem
1435 a2(i,k+1) = v1(i,k+1) + ptem2 * ptem
1436 else
1437 ad(i,k) = ad(i,k)-au(i,k)
1438 ad(i,k+1) = 1.-al(i,k)
1439 a1(i,k+1) = u1(i,k+1)
1440 a2(i,k+1) = v1(i,k+1)
1441 endif
1442!
1443 enddo
1444 enddo
1445
1446!
1447! solve tridiagonal problem for momentum
1448!
1449 call tridi2(im,km,al,ad,au,a1,a2,au,a1,a2)
1450!
1451! recover tendencies of momentum
1452!
1454 do k = 1,km
1455 do i = 1,im
1456 utend = (a1(i,k)-u1(i,k))*rdt
1457 vtend = (a2(i,k)-v1(i,k))*rdt
1458 du(i,k) = du(i,k) + utend
1459 dv(i,k) = dv(i,k) + vtend
1460 dusfc(i) = dusfc(i) + conw*del(i,k)*utend
1461 dvsfc(i) = dvsfc(i) + conw*del(i,k)*vtend
1462!
1463! for dissipative heating for ecmwf model
1464!
1465! tem1 = 0.5*(a1(i,k)+u1(i,k))
1466! tem2 = 0.5*(a2(i,k)+v1(i,k))
1467! diss(i,k) = -(tem1*utend+tem2*vtend)
1468! diss(i,k) = max(diss(i,k),0.)
1469! ttend = diss(i,k) / cp
1470! tau(i,k) = tau(i,k) + ttend
1471!
1472 enddo
1473 enddo
1474 if(.not.flag_for_pbl_generic_tend) then
1475 idtend1 = dtidx(index_of_x_wind,index_of_process_pbl)
1476 if(idtend1>=1) then
1477 do k = 1,km
1478 do i = 1,im
1479 utend = (a1(i,k)-u1(i,k))*rdt
1480 dtend(i,k,idtend1) = dtend(i,k,idtend1) + utend*delt
1481 enddo
1482 enddo
1483 endif
1484
1485 idtend2 = dtidx(index_of_y_wind,index_of_process_pbl)
1486 if(idtend2>=1) then
1487 do k = 1,km
1488 do i = 1,im
1489 vtend = (a2(i,k)-v1(i,k))*rdt
1490 dtend(i,k,idtend2) = dtend(i,k,idtend2) + vtend*delt
1491 enddo
1492 enddo
1493 endif
1494 endif
1495!
1496!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1497!
1498 do i = 1, im
1499 hpbl(i) = hpblx(i)
1500 kpbl(i) = kpblx(i)
1501 enddo
1502!
1503!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1504 return
1505 end subroutine hedmf_run
1508
1509 end module hedmf
subroutine hedmf_run(im, km, ntrac, ntcw, dv, du, tau, rtg, u1, v1, t1, q1, swh, hlw, xmu, psk, rbsoil, zorl, u10m, v10m, fm, fh, tsea, heat, evap, stress, spd1, kpbl, prsi, del, prsl, prslk, phii, phil, delt, dspheat, dusfc, dvsfc, dtsfc, dqsfc, hpbl, hgamt, hgamq, dkt, dku, kinver, xkzm_m, xkzm_h, xkzm_s, lprnt, ipr, xkzminv, moninq_fac, hurr_pbl, islimsk, var_ric, coef_ric_l, coef_ric_s, ldiag3d, ntqv, rtg_ozone_index, ntoz, dtend, dtidx, index_of_process_pbl, index_of_x_wind, index_of_y_wind, index_of_temperature, flag_for_pbl_generic_tend, errmsg, errflg)
Definition hedmf.f:81
subroutine mfpbl(im, ix, km, ntrac, delt, cnvflg, zl, zm, thvx, q1, t1, u1, v1, hpbl, kpbl, sflx, ustar, wstar, xmf, tcko, qcko, ucko, vcko)
This subroutine is used for calculating the mass flux and updraft properties.
Definition mfpbl.f:48
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 .
Definition tridi.f:95
subroutine tridi2(l, n, cl, cm, cu, r1, r2, au, a1, a2)
This subroutine ..
Definition tridi.f:50
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module contains the CCPP-compliant hybrid eddy-diffusivity mass-flux scheme.
Definition hedmf.f:7
This module contains some of the most frequently used math and physics constants for GCM models.
Definition physcons.F90:36