CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
satmedmfvdif.F
1
5
7 use tridi_mod
8 use mfscu_mod
9 use mfpblt_mod
10 contains
11
15 subroutine satmedmfvdif_init (satmedmf,
16 & isatmedmf,isatmedmf_vdif,
17 & errmsg,errflg)
18
19 logical, intent(in) :: satmedmf
20 integer, intent(in) :: isatmedmf,isatmedmf_vdif
21 character(len=*), intent(out) :: errmsg
22 integer, intent(out) :: errflg
23
24! Initialize CCPP error handling variables
25 errmsg = ''
26 errflg = 0
27
28! Consistency checks
29 if (.not. satmedmf) then
30 write(errmsg,fmt='(*(a))') 'Logic error: satmedmf = .false.'
31 errflg = 1
32 return
33 end if
34
35 if (.not. isatmedmf==isatmedmf_vdif) then
36 write(errmsg,fmt='(*(a))') 'Logic error: satmedmfvdif is ',
37 & 'called, but isatmedmf/=isatmedmf_vdif.'
38 errflg = 1
39 return
40 end if
41
42 end subroutine satmedmfvdif_init
43
62 subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, &
63 & grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, &
64 & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea, &
65 & psk,rbsoil,zorl,u10m,v10m,fm,fh, &
66 & tsea,heat,evap,stress,spd1,kpbl, &
67 & prsi,del,prsl,prslk,phii,phil,delt, &
68 & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl, &
69 & kinver,xkzm_m,xkzm_h,xkzm_s, &
70 & index_of_temperature,index_of_x_wind,index_of_y_wind, &
71 & index_of_process_pbl,ntqv,ntoz,dtend,dtidx, &
72 & gen_tend,ldiag3d,errmsg,errflg)
73!
74 use machine , only : kind_phys
75 use funcphys , only : fpvs
76!
77 implicit none
78!
79!----------------------------------------------------------------------
80 integer, intent(in) :: im, km, ntrac, ntcw, ntiw, ntke
81 integer, intent(in) :: kinver(:)
82 integer, intent(out) :: kpbl(:)
83!
84 logical, intent(in) :: gen_tend, ldiag3d
85 real(kind=kind_phys), intent(inout), dimension(:,:,:), optional ::&
86 & dtend
87 integer, intent(in) :: index_of_temperature,index_of_x_wind, &
88 & index_of_y_wind, ntqv, ntoz, dtidx(:,:), index_of_process_pbl
89!
90 real(kind=kind_phys), intent(in) :: grav,rd,cp,rv,hvap,hfus,fv, &
91 & eps,epsm1
92 real(kind=kind_phys), intent(in) :: delt, xkzm_m, xkzm_h, xkzm_s
93 real(kind=kind_phys), intent(inout) :: dv(:,:), du(:,:), &
94 & tdt(:,:), rtg(:,:,:)
95 real(kind=kind_phys), intent(in) :: &
96 & u1(:,:), v1(:,:), &
97 & t1(:,:), q1(:,:,:), &
98 & swh(:,:), hlw(:,:), &
99 & xmu(:), garea(:), &
100 & psk(:), rbsoil(:), &
101 & zorl(:), tsea(:), &
102 & u10m(:), v10m(:), &
103 & fm(:), fh(:), &
104 & evap(:), heat(:), &
105 & stress(:), spd1(:), &
106 & prsi(:,:), del(:,:), &
107 & prsl(:,:), prslk(:,:), &
108 & phii(:,:), phil(:,:)
109 real(kind=kind_phys), intent(out) :: &
110 & dusfc(:), dvsfc(:), &
111 & dtsfc(:), dqsfc(:), &
112 & hpbl(:)
113!
114 logical, intent(in) :: dspheat
115 character(len=*), intent(out) :: errmsg
116 integer, intent(out) :: errflg
117
118! flag for tke dissipative heating
119!
120!----------------------------------------------------------------------
121!***
122!*** local variables
123!***
124 integer i,is,k,kk,n,km1,kmpbl,kmscu,ntrac1
125 integer lcld(im),kcld(im),krad(im),mrad(im)
126 integer kx1(im), kpblx(im)
127!
128 real(kind=kind_phys) tke(im,km), tkeh(im,km-1)
129!
130 real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km),
131 & qlx(im,km), thetae(im,km),thlx(im,km),
132 & slx(im,km), svx(im,km), qtx(im,km),
133 & tvx(im,km), pix(im,km), radx(im,km-1),
134 & dku(im,km-1),dkt(im,km-1), dkq(im,km-1),
135 & cku(im,km-1),ckt(im,km-1)
136!
137 real(kind=kind_phys) plyr(im,km), rhly(im,km), cfly(im,km),
138 & qstl(im,km)
139!
140 real(kind=kind_phys) dtdz1(im), gdx(im),
141 & phih(im), phim(im), prn(im,km-1),
142 & rbdn(im), rbup(im), thermal(im),
143 & ustar(im), wstar(im), hpblx(im),
144 & ust3(im), wst3(im),
145 & z0(im), crb(im),
146 & hgamt(im), hgamq(im),
147 & wscale(im),vpert(im),
148 & zol(im), sflux(im), radj(im),
149 & tx1(im), tx2(im)
150!
151 real(kind=kind_phys) radmin(im)
152!
153 real(kind=kind_phys) zi(im,km+1), zl(im,km), zm(im,km),
154 & xkzo(im,km-1),xkzmo(im,km-1),
155 & xkzm_hx(im), xkzm_mx(im),
156 & rdzt(im,km-1),
157 & al(im,km-1), ad(im,km), au(im,km-1),
158 & f1(im,km), f2(im,km*(ntrac-1))
159!
160 real(kind=kind_phys) elm(im,km), ele(im,km), rle(im,km-1),
161 & ckz(im,km), chz(im,km),
162 & diss(im,km-1),prod(im,km-1),
163 & bf(im,km-1), shr2(im,km-1),
164 & xlamue(im,km-1), xlamde(im,km-1),
165 & gotvx(im,km), rlam(im,km-1)
166!
167! variables for updrafts (thermals)
168!
169 real(kind=kind_phys) tcko(im,km), qcko(im,km,ntrac),
170 & ucko(im,km), vcko(im,km),
171 & buou(im,km), xmf(im,km)
172!
173! variables for stratocumulus-top induced downdrafts
174!
175 real(kind=kind_phys) tcdo(im,km), qcdo(im,km,ntrac),
176 & ucdo(im,km), vcdo(im,km),
177 & buod(im,km), xmfd(im,km)
178!
179 logical pblflg(im), sfcflg(im), flg(im)
180 logical scuflg(im), pcnvflg(im)
181 logical mlenflg
182!
183! pcnvflg: true for unstable pbl
184!
185 real(kind=kind_phys) aphi16, aphi5,
186 & wfac, cfac,
187 & gamcrt, gamcrq, sfcfrac,
188 & conq, cont, conw,
189 & dsdz2, dsdzt, dkmax,
190 & dsig, dt2, dtodsd,
191 & dtodsu, g, factor, dz,
192 & gocp, gravi, zol1, zolcru,
193 & buop, shrp, dtn, cdtn,
194 & prnum, prmax, prmin, prtke,
195 & prscu, dw2, dw2min, zk,
196 & elmfac, elefac, dspmax,
197 & alp, clwt, cql,
198 & f0, robn, crbmin, crbmax,
199 & es, qs, value, onemrh,
200 & cfh, gamma, elocp, el2orc,
201 & epsi, beta, chx, cqx,
202 & rdt, rdz, qmin, qlmin,
203 & ri, rimin,
204 & rbcr, rbint, tdzmin,
205 & rlmn, rlmx, elmx,
206 & ttend, utend, vtend, qtend,
207 & zfac, zfmin, vk, spdk2,
208 & tkmin, xkzinv, dspfac, xkgdx,
209 & zlup, zldn, bsum,
210 & tem, tem1, tem2,
211 & ptem, ptem0, ptem1, ptem2
212!
213 real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck
214!
215 real(kind=kind_phys) qlcr, zstblmax
216!
217 real(kind=kind_phys) h1
218 integer :: idtend
219!!
220 parameter(wfac=7.0,cfac=4.5)
221 parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1)
222 parameter(vk=0.4,rimin=-100.)
223 parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3)
224 parameter(rlmn=30.,rlmx=500.,elmx=500.)
225 parameter(prmin=0.25,prmax=4.0,prtke=1.0,prscu=0.67)
226 parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35)
227 parameter(tkmin=1.e-9,dspfac=0.5,dspmax=10.0)
228 parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8)
229 parameter(aphi5=5.,aphi16=16.)
230 parameter(elmfac=1.0,elefac=1.0,cql=100.)
231 parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=25000.)
232 parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.15)
233 parameter(h1=0.33333333)
234 parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15,ce0=0.4)
235 parameter(rchck=1.5,cdtn=25.)
236
237 gravi=1.0/grav
238 g=grav
239 gocp=g/cp
240 cont=cp/g
241 conq=hvap/g
242 conw=1.0/g ! for del in pa
243! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) !kpa
244 elocp=hvap/cp
245 el2orc=hvap*hvap/(rv*cp)
246
247!
248!************************************************************************
249! Initialize CCPP error handling variables
250 errmsg = ''
251 errflg = 0
252
254 dt2 = delt
255 rdt = 1. / dt2
256!
257! the code is written assuming ntke=ntrac
258! if ntrac > ntke, the code needs to be modified
259!
260 ntrac1 = ntrac - 1
261 km1 = km - 1
262 kmpbl = km / 2
263 kmscu = km / 2
266 do k=1,km
267 do i=1,im
268 zi(i,k) = phii(i,k) * gravi
269 zl(i,k) = phil(i,k) * gravi
270 xmf(i,k) = 0.
271 xmfd(i,k) = 0.
272 buou(i,k) = 0.
273 buod(i,k) = 0.
274 ckz(i,k) = ck1
275 chz(i,k) = ch1
276 enddo
277 enddo
278 do i=1,im
279 zi(i,km+1) = phii(i,km+1) * gravi
280 enddo
281 do k=1,km
282 do i=1,im
283 zm(i,k) = zi(i,k+1)
284 enddo
285 enddo
287 do i=1,im
288 gdx(i) = sqrt(garea(i))
289 enddo
292 do k=1,km
293 do i=1,im
294 tke(i,k) = max(q1(i,k,ntke), tkmin)
295 enddo
296 enddo
297 do k=1,km1
298 do i=1,im
299 tkeh(i,k) = 0.5 * (tke(i,k) + tke(i,k+1))
300 enddo
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 prn(i,k) = 1.0
307 enddo
308 enddo
309!
317
318 do i=1,im
319 kx1(i) = 1
320 tx1(i) = 1.0 / prsi(i,1)
321 tx2(i) = tx1(i)
322 if(gdx(i) >= xkgdx) then
323 xkzm_hx(i) = xkzm_h
324 xkzm_mx(i) = xkzm_m
325 else
326 tem = 1. / (xkgdx - 5.)
327 tem1 = (xkzm_h - 0.01) * tem
328 tem2 = (xkzm_m - 0.01) * tem
329 ptem = gdx(i) - 5.
330 xkzm_hx(i) = 0.01 + tem1 * ptem
331 xkzm_mx(i) = 0.01 + tem2 * ptem
332 endif
333 enddo
334 do k = 1,km1
335 do i=1,im
336 xkzo(i,k) = 0.0
337 xkzmo(i,k) = 0.0
338 if (k < kinver(i)) then
339! vertical background diffusivity
340 ptem = prsi(i,k+1) * tx1(i)
341 tem1 = 1.0 - ptem
342 tem1 = tem1 * tem1 * 10.0
343 xkzo(i,k) = xkzm_hx(i) * min(1.0, exp(-tem1))
344! vertical background diffusivity for momentum
345 if (ptem >= xkzm_s) then
346 xkzmo(i,k) = xkzm_mx(i)
347 kx1(i) = k + 1
348 else
349 if (k == kx1(i) .and. k > 1) tx2(i) = 1.0 / prsi(i,k)
350 tem1 = 1.0 - prsi(i,k+1) * tx2(i)
351 tem1 = tem1 * tem1 * 5.0
352 xkzmo(i,k) = xkzm_mx(i) * min(1.0, exp(-tem1))
353 endif
354 endif
355 enddo
356 enddo
357!
359 do i = 1,im
360 z0(i) = 0.01 * zorl(i)
361 dusfc(i) = 0.
362 dvsfc(i) = 0.
363 dtsfc(i) = 0.
364 dqsfc(i) = 0.
365 kpbl(i) = 1
366 hpbl(i) = 0.
367 kpblx(i) = 1
368 hpblx(i) = 0.
369 pblflg(i)= .true.
370 sfcflg(i)= .true.
371 if(rbsoil(i) > 0.) sfcflg(i) = .false.
372 pcnvflg(i)= .false.
373 scuflg(i)= .true.
374 if(scuflg(i)) then
375 radmin(i)= 0.
376 mrad(i) = km1
377 krad(i) = 1
378 lcld(i) = km1
379 kcld(i) = km1
380 endif
381 enddo
382!
385 do k=1,km
386 do i=1,im
387 pix(i,k) = psk(i) / prslk(i,k)
388 theta(i,k) = t1(i,k) * pix(i,k)
389 if(ntiw > 0) then
390 tem = max(q1(i,k,ntcw),qlmin)
391 tem1 = max(q1(i,k,ntiw),qlmin)
392 qlx(i,k) = tem + tem1
393 ptem = hvap*tem + (hvap+hfus)*tem1
394 slx(i,k) = cp * t1(i,k) + phil(i,k) - ptem
395 else
396 qlx(i,k) = max(q1(i,k,ntcw),qlmin)
397 slx(i,k) = cp * t1(i,k) + phil(i,k) - hvap*qlx(i,k)
398 endif
399 tem2 = 1.+fv*max(q1(i,k,1),qmin)-qlx(i,k)
400 thvx(i,k) = theta(i,k) * tem2
401 tvx(i,k) = t1(i,k) * tem2
402 qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k)
403 thlx(i,k) = theta(i,k) - pix(i,k)*elocp*qlx(i,k)
404 thlvx(i,k) = thlx(i,k) * (1. + fv * qtx(i,k))
405 svx(i,k) = cp * tvx(i,k)
406 ptem1 = elocp * pix(i,k) * max(q1(i,k,1),qmin)
407 thetae(i,k)= theta(i,k) + ptem1
408 gotvx(i,k) = g / tvx(i,k)
409 enddo
410 enddo
411!
414!
415 do k = 1,km1
416 do i=1,im
417 tem1 = (tvx(i,k+1)-tvx(i,k)) * rdzt(i,k)
418 if(tem1 > 1.e-5) then
419 xkzo(i,k) = min(xkzo(i,k),xkzinv)
420 xkzmo(i,k) = min(xkzmo(i,k),xkzinv)
421 endif
422 enddo
423 enddo
424!
427 do k = 1, km
428 do i = 1, im
429 plyr(i,k) = 0.01 * prsl(i,k) ! pa to mb (hpa)
430! --- ... compute relative humidity
431 es = 0.01 * fpvs(t1(i,k)) ! fpvs in pa
432 qs = max(qmin, eps * es / (plyr(i,k) + epsm1*es))
433 rhly(i,k) = max(0.0, min(1.0, max(qmin, q1(i,k,1))/qs))
434 qstl(i,k) = qs
435 enddo
436 enddo
437!
438 do k = 1, km
439 do i = 1, im
440 cfly(i,k) = 0.
441 clwt = 1.0e-6 * (plyr(i,k)*0.001)
442 if (qlx(i,k) > clwt) then
443 onemrh= max(1.e-10, 1.0-rhly(i,k))
444 tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0)
445 tem1 = cql / tem1
446 value = max(min( tem1*qlx(i,k), 50.0), 0.0)
447 tem2 = sqrt(sqrt(rhly(i,k)))
448 cfly(i,k) = min(max(tem2*(1.0-exp(-value)), 0.0), 1.0)
449 endif
450 enddo
451 enddo
452!
454!
455 do k = 1, km1
456 do i = 1, im
457 tem = 0.5 * (svx(i,k) + svx(i,k+1))
458 tem1 = 0.5 * (t1(i,k) + t1(i,k+1))
459 tem2 = 0.5 * (qstl(i,k) + qstl(i,k+1))
460 cfh = min(cfly(i,k+1),0.5*(cfly(i,k)+cfly(i,k+1)))
461 alp = g / tem
462 gamma = el2orc * tem2 / (tem1**2)
463 epsi = tem1 / elocp
464 beta = (1. + gamma*epsi*(1.+fv)) / (1. + gamma)
465 chx = cfh * alp * beta + (1. - cfh) * alp
466 cqx = cfh * alp * hvap * (beta - epsi)
467 cqx = cqx + (1. - cfh) * fv * g
468 ptem1 = (slx(i,k+1)-slx(i,k))*rdzt(i,k)
469 ptem2 = (qtx(i,k+1)-qtx(i,k))*rdzt(i,k)
470 bf(i,k) = chx * ptem1 + cqx * ptem2
471 enddo
472 enddo
473!
474!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
475!
478 do k=1,km1
479 do i=1,im
480 dku(i,k) = 0.
481 dkt(i,k) = 0.
482 dkq(i,k) = 0.
483 cku(i,k) = 0.
484 ckt(i,k) = 0.
485 tem = zi(i,k+1)-zi(i,k)
486 radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k))
487 enddo
488 enddo
492 do i = 1,im
493 sflux(i) = heat(i) + evap(i)*fv*theta(i,1)
494 if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false.
495 enddo
496!
511!
512 do i = 1,im
513 if(pblflg(i)) then
514! thermal(i) = thvx(i,1)
515 thermal(i) = thlvx(i,1)
516 crb(i) = rbcr
517 else
518 thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
519 tem = sqrt(u10m(i)**2+v10m(i)**2)
520 tem = max(tem, 1.)
521 robn = tem / (f0 * z0(i))
522 tem1 = 1.e-7 * robn
523 crb(i) = 0.16 * (tem1 ** (-0.18))
524 crb(i) = max(min(crb(i), crbmax), crbmin)
525 endif
526 enddo
527!
529 do i=1,im
530 dtdz1(i) = dt2 / (zi(i,2)-zi(i,1))
531 enddo
532!
533 do i=1,im
534 ustar(i) = sqrt(stress(i))
535 enddo
536!
539!
540 do k = 1, km1
541 do i = 1, im
542 rdz = rdzt(i,k)
543! bf(i,k) = gotvx(i,k)*(thvx(i,k+1)-thvx(i,k))*rdz
544 dw2 = (u1(i,k)-u1(i,k+1))**2
545 & + (v1(i,k)-v1(i,k+1))**2
546 shr2(i,k) = max(dw2,dw2min)*rdz*rdz
547 enddo
548 enddo
549!
550! Find pbl height based on bulk Richardson number (mrf pbl scheme)
551! and also for diagnostic purpose
552!
553!
554 do i=1,im
555 flg(i) = .false.
556 rbup(i) = rbsoil(i)
557 enddo
563 do k = 1, kmpbl
564 do i = 1, im
565 if(.not.flg(i)) then
566 rbdn(i) = rbup(i)
567 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
568! rbup(i) = (thvx(i,k)-thermal(i))*
569! & (g*zl(i,k)/thvx(i,1))/spdk2
570 rbup(i) = (thlvx(i,k)-thermal(i))*
571 & (g*zl(i,k)/thlvx(i,1))/spdk2
572 kpblx(i) = k
573 flg(i) = rbup(i) > crb(i)
574 endif
575 enddo
576 enddo
580 do i = 1,im
581 if(kpblx(i) > 1) then
582 k = kpblx(i)
583 if(rbdn(i) >= crb(i)) then
584 rbint = 0.
585 elseif(rbup(i) <= crb(i)) then
586 rbint = 1.
587 else
588 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
589 endif
590 hpblx(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
591 if(hpblx(i) < zi(i,kpblx(i))) kpblx(i)=kpblx(i)-1
592 else
593 hpblx(i) = zl(i,1)
594 kpblx(i) = 1
595 endif
596 hpbl(i) = hpblx(i)
597 kpbl(i) = kpblx(i)
598 if(kpbl(i) <= 1) pblflg(i)=.false.
599 enddo
600!
610 do i=1,im
611 zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
612 if(sfcflg(i)) then
613 zol(i) = min(zol(i),-zfmin)
614 else
615 zol(i) = max(zol(i),zfmin)
616 endif
629 zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1)
630 if(sfcflg(i)) then
631 tem = 1.0 / (1. - aphi16*zol1)
632 phih(i) = sqrt(tem)
633 phim(i) = sqrt(phih(i))
634 else
635 phim(i) = 1. + aphi5*zol1
636 phih(i) = phim(i)
637 endif
638 enddo
639!
655 do i=1,im
656 if(pblflg(i)) then
657 if(zol(i) < zolcru) then
658 pcnvflg(i) = .true.
659 endif
660 wst3(i) = gotvx(i,1)*sflux(i)*hpbl(i)
661 wstar(i)= wst3(i)**h1
662 ust3(i) = ustar(i)**3.
663 wscale(i)=(ust3(i)+wfac*vk*wst3(i)*sfcfrac)**h1
664 ptem = ustar(i)/aphi5
665 wscale(i) = max(wscale(i),ptem)
666 endif
667 enddo
668!
670!
671 do i = 1,im
672 if(pcnvflg(i)) then
673 hgamt(i) = heat(i)/wscale(i)
674 hgamq(i) = evap(i)/wscale(i)
675 vpert(i) = hgamt(i) + hgamq(i)*fv*theta(i,1)
676 vpert(i) = max(vpert(i),0.)
677 tem = min(cfac*vpert(i),gamcrt)
678 thermal(i)= thermal(i) + tem
679 endif
680 enddo
681!
685! (overshoot pbl top)
686!
687 do i=1,im
688 flg(i) = .true.
689 if(pcnvflg(i)) then
690 flg(i) = .false.
691 rbup(i) = rbsoil(i)
692 endif
693 enddo
694 do k = 2, kmpbl
695 do i = 1, im
696 if(.not.flg(i)) then
697 rbdn(i) = rbup(i)
698 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
699 rbup(i) = (thlvx(i,k)-thermal(i))*
700 & (g*zl(i,k)/thlvx(i,1))/spdk2
701 kpbl(i) = k
702 flg(i) = rbup(i) > crb(i)
703 endif
704 enddo
705 enddo
706 do i = 1,im
707 if(pcnvflg(i)) then
708 k = kpbl(i)
709 if(rbdn(i) >= crb(i)) then
710 rbint = 0.
711 elseif(rbup(i) <= crb(i)) then
712 rbint = 1.
713 else
714 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
715 endif
716 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
717 if(hpbl(i) < zi(i,kpbl(i))) then
718 kpbl(i) = kpbl(i) - 1
719 endif
720 if(kpbl(i) <= 1) then
721 pcnvflg(i) = .false.
722 pblflg(i) = .false.
723 endif
724 endif
725 enddo
726!
727!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
728! look for stratocumulus
733 do i=1,im
734 flg(i) = scuflg(i)
735 enddo
736 do k = 1, km1
737 do i=1,im
738 if(flg(i).and.zl(i,k) >= zstblmax) then
739 lcld(i)=k
740 flg(i)=.false.
741 endif
742 enddo
743 enddo
744 do i = 1, im
745 flg(i)=scuflg(i)
746 enddo
747 do k = kmscu,1,-1
748 do i = 1, im
749 if(flg(i) .and. k <= lcld(i)) then
750 if(qlx(i,k) >= qlcr) then
751 kcld(i)=k
752 flg(i)=.false.
753 endif
754 endif
755 enddo
756 enddo
757 do i = 1, im
758 if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false.
759 enddo
765 do i = 1, im
766 flg(i)=scuflg(i)
767 enddo
768 do k = kmscu,1,-1
769 do i = 1, im
770 if(flg(i) .and. k <= kcld(i)) then
771 if(qlx(i,k) >= qlcr) then
772 if(radx(i,k) < radmin(i)) then
773 radmin(i)=radx(i,k)
774 krad(i)=k
775 endif
776 else
777 flg(i)=.false.
778 endif
779 endif
780 enddo
781 enddo
782 do i = 1, im
783 if(scuflg(i) .and. krad(i) <= 1) scuflg(i)=.false.
784 if(scuflg(i) .and. radmin(i)>=0.) scuflg(i)=.false.
785 enddo
786!
787!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
792 do k = 1, km
793 do i = 1, im
794 if(pcnvflg(i)) then
795 tcko(i,k) = t1(i,k)
796 ucko(i,k) = u1(i,k)
797 vcko(i,k) = v1(i,k)
798 endif
799 if(scuflg(i)) then
800 tcdo(i,k) = t1(i,k)
801 ucdo(i,k) = u1(i,k)
802 vcdo(i,k) = v1(i,k)
803 endif
804 enddo
805 enddo
806 do kk = 1, ntrac1
807 do k = 1, km
808 do i = 1, im
809 if(pcnvflg(i)) then
810 qcko(i,k,kk) = q1(i,k,kk)
811 endif
812 if(scuflg(i)) then
813 qcdo(i,k,kk) = q1(i,k,kk)
814 endif
815 enddo
816 enddo
817 enddo
820 call mfpblt(im,im,km,kmpbl,ntcw,ntrac1,dt2,
821 & pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx,
822 & gdx,hpbl,kpbl,vpert,buou,xmf,
823 & tcko,qcko,ucko,vcko,xlamue)
826 call mfscu(im,im,km,kmscu,ntcw,ntrac1,dt2,
827 & scuflg,zl,zm,q1,t1,u1,v1,plyr,pix,
828 & thlx,thvx,thlvx,gdx,thetae,radj,
829 & krad,mrad,radmin,buod,xmfd,
830 & tcdo,qcdo,ucdo,vcdo,xlamde)
831!
832!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
842 do k = 1, kmpbl
843 do i = 1, im
844 if(k < kpbl(i)) then
845 tem = phih(i)/phim(i)
846 ptem = -3.*(max(zi(i,k+1)-sfcfrac*hpbl(i),0.))**2.
847 & /hpbl(i)**2.
848 if(pcnvflg(i)) then
849 prn(i,k) = 1. + (tem-1.)*exp(ptem)
850 else
851 prn(i,k) = tem
852 endif
853 prn(i,k) = min(prn(i,k),prmax)
854 prn(i,k) = max(prn(i,k),prmin)
855!
856 ckz(i,k) = ck1 + (ck0-ck1)*exp(ptem)
857 ckz(i,k) = min(ckz(i,k),ck0)
858 ckz(i,k) = max(ckz(i,k),ck1)
859 chz(i,k) = ch1 + (ch0-ch1)*exp(ptem)
860 chz(i,k) = min(chz(i,k),ch0)
861 chz(i,k) = max(chz(i,k),ch1)
862 endif
863 enddo
864 enddo
865
866!
867!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
869!
870 do k = 1, km1
871 do i = 1, im
872 zlup = 0.0
873 bsum = 0.0
874 mlenflg = .true.
875 do n = k, km1
876 if(mlenflg) then
877 dz = zl(i,n+1) - zl(i,n)
878 ptem = gotvx(i,n)*(thvx(i,n+1)-thvx(i,k))*dz
879! ptem = gotvx(i,n)*(thlvx(i,n+1)-thlvx(i,k))*dz
880 bsum = bsum + ptem
881 zlup = zlup + dz
882 if(bsum >= tke(i,k)) then
883 if(ptem >= 0.) then
884 tem2 = max(ptem, zfmin)
885 else
886 tem2 = min(ptem, -zfmin)
887 endif
888 ptem1 = (bsum - tke(i,k)) / tem2
889 zlup = zlup - ptem1 * dz
890 zlup = max(zlup, 0.)
891 mlenflg = .false.
892 endif
893 endif
894 enddo
895 zldn = 0.0
896 bsum = 0.0
897 mlenflg = .true.
898 do n = k, 1, -1
899 if(mlenflg) then
900 if(n == 1) then
901 dz = zl(i,1)
902 tem1 = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
903 else
904 dz = zl(i,n) - zl(i,n-1)
905 tem1 = thvx(i,n-1)
906! tem1 = thlvx(i,n-1)
907 endif
908 ptem = gotvx(i,n)*(thvx(i,k)-tem1)*dz
909! ptem = gotvx(i,n)*(thlvx(i,k)-tem1)*dz
910 bsum = bsum + ptem
911 zldn = zldn + dz
912 if(bsum >= tke(i,k)) then
913 if(ptem >= 0.) then
914 tem2 = max(ptem, zfmin)
915 else
916 tem2 = min(ptem, -zfmin)
917 endif
918 ptem1 = (bsum - tke(i,k)) / tem2
919 zldn = zldn - ptem1 * dz
920 zldn = max(zldn, 0.)
921 mlenflg = .false.
922 endif
923 endif
924 enddo
925!
926 tem = 0.5 * (zi(i,k+1)-zi(i,k))
927 tem1 = min(tem, rlmn)
940 ptem2 = min(zlup,zldn)
941 rlam(i,k) = elmfac * ptem2
942 rlam(i,k) = max(rlam(i,k), tem1)
943 rlam(i,k) = min(rlam(i,k), rlmx)
944!
945 ptem2 = sqrt(zlup*zldn)
946 ele(i,k) = elefac * ptem2
947 ele(i,k) = max(ele(i,k), tem1)
948 ele(i,k) = min(ele(i,k), elmx)
949!
950 enddo
951 enddo
954 do k = 1, km1
955 do i = 1, im
956 tem = vk * zl(i,k)
957 if (zol(i) < 0.) then
958 ptem = 1. - 100. * zol(i)
959 ptem1 = ptem**0.2
960 zk = tem * ptem1
961 elseif (zol(i) >= 1.) then
962 zk = tem / 3.7
963 else
964 ptem = 1. + 2.7 * zol(i)
965 zk = tem / ptem
966 endif
967 elm(i,k) = zk*rlam(i,k)/(rlam(i,k)+zk)
968!
969 dz = zi(i,k+1) - zi(i,k)
970 tem = max(gdx(i),dz)
971 elm(i,k) = min(elm(i,k), tem)
972 ele(i,k) = min(ele(i,k), tem)
973!
974 enddo
975 enddo
976 do i = 1, im
977 elm(i,km) = elm(i,km1)
978 ele(i,km) = ele(i,km1)
979 enddo
980!
982!
983 do k = 1, km1
984 do i = 1, im
985 tem = 0.5 * (elm(i,k) + elm(i,k+1))
986 tem = tem * sqrt(tkeh(i,k))
987 if(k < kpbl(i)) then
988 if(pblflg(i)) then
989 dku(i,k) = ckz(i,k) * tem
990 dkt(i,k) = dku(i,k) / prn(i,k)
991 else
992 dkt(i,k) = chz(i,k) * tem
993 dku(i,k) = dkt(i,k) * prn(i,k)
994 endif
995 else
996 ri = max(bf(i,k)/shr2(i,k),rimin)
997 if(ri < 0.) then ! unstable regime
998 dku(i,k) = ck1 * tem
999 dkt(i,k) = rchck * dku(i,k)
1000 else ! stable regime
1001 dkt(i,k) = ch1 * tem
1002 prnum = 1.0 + 2.1*ri
1003 prnum = min(prnum,prmax)
1004 dku(i,k) = dkt(i,k) * prnum
1005 endif
1006 endif
1007!
1008 if(scuflg(i)) then
1009 if(k >= mrad(i) .and. k < krad(i)) then
1010 tem1 = ckz(i,k) * tem
1011 ptem1 = tem1 / prscu
1012 dku(i,k) = max(dku(i,k), tem1)
1013 dkt(i,k) = max(dkt(i,k), ptem1)
1014 endif
1015 endif
1016!
1017 dkq(i,k) = prtke * dkt(i,k)
1018!
1019 dkt(i,k) = min(dkt(i,k),dkmax)
1020 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
1021 dkq(i,k) = min(dkq(i,k),dkmax)
1022 dkq(i,k) = max(dkq(i,k),xkzo(i,k))
1023 dku(i,k) = min(dku(i,k),dkmax)
1024 dku(i,k) = max(dku(i,k),xkzmo(i,k))
1025!
1026 enddo
1027 enddo
1028!
1029 do i = 1, im
1030 if(scuflg(i)) then
1031 k = krad(i)
1032 tem = bf(i,k) / gotvx(i,k)
1033 tem1 = max(tem, tdzmin)
1034 ptem = radj(i) / tem1
1035 dkt(i,k) = dkt(i,k) + ptem
1036 dku(i,k) = dku(i,k) + ptem
1037 dkq(i,k) = dkq(i,k) + ptem
1038 endif
1039 enddo
1040!
1041!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1043!
1044 do k = 1, km1
1045 do i = 1, im
1046 if (k == 1) then
1047 tem = -dkt(i,1) * bf(i,1)
1048! if(pcnvflg(i)) then
1049! ptem1 = xmf(i,1) * buou(i,1)
1050! else
1051 ptem1 = 0.
1052! endif
1053 if(scuflg(i) .and. mrad(i) == 1) then
1054 ptem2 = xmfd(i,1) * buod(i,1)
1055 else
1056 ptem2 = 0.
1057 endif
1058 tem = tem + ptem1 + ptem2
1059 buop = 0.5 * (gotvx(i,1) * sflux(i) + tem)
1060!
1061 tem1 = dku(i,1) * shr2(i,1)
1062!
1063 tem = (u1(i,2)-u1(i,1))*rdzt(i,1)
1064! if(pcnvflg(i)) then
1065! ptem = xmf(i,1) * tem
1066! ptem1 = 0.5 * ptem * (u1(i,2)-ucko(i,2))
1067! else
1068 ptem1 = 0.
1069! endif
1070 if(scuflg(i) .and. mrad(i) == 1) then
1071 ptem = ucdo(i,1)+ucdo(i,2)-u1(i,1)-u1(i,2)
1072 ptem = 0.5 * tem * xmfd(i,1) * ptem
1073 else
1074 ptem = 0.
1075 endif
1076 ptem1 = ptem1 + ptem
1077!
1078 tem = (v1(i,2)-v1(i,1))*rdzt(i,1)
1079! if(pcnvflg(i)) then
1080! ptem = xmf(i,1) * tem
1081! ptem2 = 0.5 * ptem * (v1(i,2)-vcko(i,2))
1082! else
1083 ptem2 = 0.
1084! endif
1085 if(scuflg(i) .and. mrad(i) == 1) then
1086 ptem = vcdo(i,1)+vcdo(i,2)-v1(i,1)-v1(i,2)
1087 ptem = 0.5 * tem * xmfd(i,1) * ptem
1088 else
1089 ptem = 0.
1090 endif
1091 ptem2 = ptem2 + ptem
1092!
1093! tem2 = stress(i)*spd1(i)/zl(i,1)
1094 tem2 = stress(i)*ustar(i)*phim(i)/(vk*zl(i,1))
1095 shrp = 0.5 * (tem1 + ptem1 + ptem2 + tem2)
1096 else
1097 tem1 = -dkt(i,k-1) * bf(i,k-1)
1098 tem2 = -dkt(i,k) * bf(i,k)
1099 tem = 0.5 * (tem1 + tem2)
1100 if(pcnvflg(i) .and. k <= kpbl(i)) then
1101 ptem = 0.5 * (xmf(i,k-1) + xmf(i,k))
1102 ptem1 = ptem * buou(i,k)
1103 else
1104 ptem1 = 0.
1105 endif
1106 if(scuflg(i)) then
1107 if(k >= mrad(i) .and. k < krad(i)) then
1108 ptem0 = 0.5 * (xmfd(i,k-1) + xmfd(i,k))
1109 ptem2 = ptem0 * buod(i,k)
1110 else
1111 ptem2 = 0.
1112 endif
1113 else
1114 ptem2 = 0.
1115 endif
1116 buop = tem + ptem1 + ptem2
1117!
1118 tem1 = dku(i,k-1) * shr2(i,k-1)
1119 tem2 = dku(i,k) * shr2(i,k)
1120 tem = 0.5 * (tem1 + tem2)
1121 tem1 = (u1(i,k+1)-u1(i,k))*rdzt(i,k)
1122 tem2 = (u1(i,k)-u1(i,k-1))*rdzt(i,k-1)
1123 if(pcnvflg(i) .and. k <= kpbl(i)) then
1124 ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2
1125 ptem1 = 0.5 * ptem * (u1(i,k)-ucko(i,k))
1126 else
1127 ptem1 = 0.
1128 endif
1129 if(scuflg(i)) then
1130 if(k >= mrad(i) .and. k < krad(i)) then
1131 ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2
1132 ptem2 = 0.5 * ptem0 * (ucdo(i,k)-u1(i,k))
1133 else
1134 ptem2 = 0.
1135 endif
1136 else
1137 ptem2 = 0.
1138 endif
1139 shrp = tem + ptem1 + ptem2
1140 tem1 = (v1(i,k+1)-v1(i,k))*rdzt(i,k)
1141 tem2 = (v1(i,k)-v1(i,k-1))*rdzt(i,k-1)
1142 if(pcnvflg(i) .and. k <= kpbl(i)) then
1143 ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2
1144 ptem1 = 0.5 * ptem * (v1(i,k)-vcko(i,k))
1145 else
1146 ptem1 = 0.
1147 endif
1148 if(scuflg(i)) then
1149 if(k >= mrad(i) .and. k < krad(i)) then
1150 ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2
1151 ptem2 = 0.5 * ptem0 * (vcdo(i,k)-v1(i,k))
1152 else
1153 ptem2 = 0.
1154 endif
1155 else
1156 ptem2 = 0.
1157 endif
1158 shrp = shrp + ptem1 + ptem2
1159 endif
1160 prod(i,k) = buop + shrp
1161 enddo
1162 enddo
1163!
1164!----------------------------------------------------------------------
1166!
1167 do k = 1,km1
1168 do i=1,im
1169 rle(i,k) = ce0 / ele(i,k)
1170 enddo
1171 enddo
1172 kk = max(nint(dt2/cdtn), 1)
1173 dtn = dt2 / float(kk)
1174 do n = 1, kk
1175 do k = 1,km1
1176 do i=1,im
1177 tem = sqrt(tke(i,k))
1178 diss(i,k) = rle(i,k) * tke(i,k) * tem
1179 tem1 = prod(i,k) + tke(i,k) / dtn
1180 diss(i,k)=max(min(diss(i,k), tem1), 0.)
1181 tke(i,k) = tke(i,k) + dtn * (prod(i,k)-diss(i,k))
1182 tke(i,k) = max(tke(i,k), tkmin)
1183 enddo
1184 enddo
1185 enddo
1186!
1188!
1189 do k = 1, km
1190 do i = 1, im
1191 if(pcnvflg(i)) then
1192 qcko(i,k,ntke) = tke(i,k)
1193 endif
1194 if(scuflg(i)) then
1195 qcdo(i,k,ntke) = tke(i,k)
1196 endif
1197 enddo
1198 enddo
1199 do k = 2, kmpbl
1200 do i = 1, im
1201 if (pcnvflg(i) .and. k <= kpbl(i)) then
1202 dz = zl(i,k) - zl(i,k-1)
1203 tem = 0.5 * xlamue(i,k-1) * dz
1204 factor = 1. + tem
1205 qcko(i,k,ntke)=((1.-tem)*qcko(i,k-1,ntke)+tem*
1206 & (tke(i,k)+tke(i,k-1)))/factor
1207 endif
1208 enddo
1209 enddo
1210 do k = kmscu, 1, -1
1211 do i = 1, im
1212 if (scuflg(i) .and. k < krad(i)) then
1213 if(k >= mrad(i)) then
1214 dz = zl(i,k+1) - zl(i,k)
1215 tem = 0.5 * xlamde(i,k) * dz
1216 factor = 1. + tem
1217 qcdo(i,k,ntke)=((1.-tem)*qcdo(i,k+1,ntke)+tem*
1218 & (tke(i,k)+tke(i,k+1)))/factor
1219 endif
1220 endif
1221 enddo
1222 enddo
1223!
1224!----------------------------------------------------------------------
1226!
1227 do i=1,im
1228 ad(i,1) = 1.0
1229 f1(i,1) = tke(i,1)
1230 enddo
1231!
1232 do k = 1,km1
1233 do i=1,im
1234 dtodsd = dt2/del(i,k)
1235 dtodsu = dt2/del(i,k+1)
1236 dsig = prsl(i,k)-prsl(i,k+1)
1237 rdz = rdzt(i,k)
1238 tem1 = dsig * dkq(i,k) * rdz
1239 dsdz2 = tem1 * rdz
1240 au(i,k) = -dtodsd*dsdz2
1241 al(i,k) = -dtodsu*dsdz2
1242 ad(i,k) = ad(i,k)-au(i,k)
1243 ad(i,k+1)= 1.-al(i,k)
1244 tem2 = dsig * rdz
1245!
1246 if(pcnvflg(i) .and. k < kpbl(i)) then
1247 ptem = 0.5 * tem2 * xmf(i,k)
1248 ptem1 = dtodsd * ptem
1249 ptem2 = dtodsu * ptem
1250 tem = tke(i,k) + tke(i,k+1)
1251 ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke)
1252 f1(i,k) = f1(i,k)-(ptem-tem)*ptem1
1253 f1(i,k+1) = tke(i,k+1)+(ptem-tem)*ptem2
1254 else
1255 f1(i,k+1) = tke(i,k+1)
1256 endif
1257!
1258 if(scuflg(i)) then
1259 if(k >= mrad(i) .and. k < krad(i)) then
1260 ptem = 0.5 * tem2 * xmfd(i,k)
1261 ptem1 = dtodsd * ptem
1262 ptem2 = dtodsu * ptem
1263 tem = tke(i,k) + tke(i,k+1)
1264 ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke)
1265 f1(i,k) = f1(i,k) + (ptem - tem) * ptem1
1266 f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2
1267 endif
1268 endif
1269!
1270 enddo
1271 enddo
1272!
1274!
1275 call tridit(im,km,1,al,ad,au,f1,au,f1)
1276!
1277! recover tendency of tke
1278!
1279 do k = 1,km
1280 do i = 1,im
1281! f1(i,k) = max(f1(i,k), tkmin)
1282 qtend = (f1(i,k)-q1(i,k,ntke))*rdt
1283 rtg(i,k,ntke) = rtg(i,k,ntke)+qtend
1284 enddo
1285 enddo
1286!
1288!
1289 do i=1,im
1290 ad(i,1) = 1.
1291 f1(i,1) = t1(i,1) + dtdz1(i) * heat(i)
1292 f2(i,1) = q1(i,1,1) + dtdz1(i) * evap(i)
1293 enddo
1294 if(ntrac1 >= 2) then
1295 do kk = 2, ntrac1
1296 is = (kk-1) * km
1297 do i = 1, im
1298 f2(i,1+is) = q1(i,1,kk)
1299 enddo
1300 enddo
1301 endif
1302!
1303 do k = 1,km1
1304 do i = 1,im
1305 dtodsd = dt2/del(i,k)
1306 dtodsu = dt2/del(i,k+1)
1307 dsig = prsl(i,k)-prsl(i,k+1)
1308 rdz = rdzt(i,k)
1309 tem1 = dsig * dkt(i,k) * rdz
1310 dsdzt = tem1 * gocp
1311 dsdz2 = tem1 * rdz
1312 au(i,k) = -dtodsd*dsdz2
1313 al(i,k) = -dtodsu*dsdz2
1314 ad(i,k) = ad(i,k)-au(i,k)
1315 ad(i,k+1)= 1.-al(i,k)
1316 tem2 = dsig * rdz
1317!
1318 if(pcnvflg(i) .and. k < kpbl(i)) then
1319 ptem = 0.5 * tem2 * xmf(i,k)
1320 ptem1 = dtodsd * ptem
1321 ptem2 = dtodsu * ptem
1322 tem = t1(i,k) + t1(i,k+1)
1323 ptem = tcko(i,k) + tcko(i,k+1)
1324 f1(i,k) = f1(i,k)+dtodsd*dsdzt-(ptem-tem)*ptem1
1325 f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+(ptem-tem)*ptem2
1326 tem = q1(i,k,1) + q1(i,k+1,1)
1327 ptem = qcko(i,k,1) + qcko(i,k+1,1)
1328 f2(i,k) = f2(i,k) - (ptem - tem) * ptem1
1329 f2(i,k+1) = q1(i,k+1,1) + (ptem - tem) * ptem2
1330 else
1331 f1(i,k) = f1(i,k)+dtodsd*dsdzt
1332 f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
1333 f2(i,k+1) = q1(i,k+1,1)
1334 endif
1335!
1336 if(scuflg(i)) then
1337 if(k >= mrad(i) .and. k < krad(i)) then
1338 ptem = 0.5 * tem2 * xmfd(i,k)
1339 ptem1 = dtodsd * ptem
1340 ptem2 = dtodsu * ptem
1341 ptem = tcdo(i,k) + tcdo(i,k+1)
1342 tem = t1(i,k) + t1(i,k+1)
1343 f1(i,k) = f1(i,k) + (ptem - tem) * ptem1
1344 f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2
1345 tem = q1(i,k,1) + q1(i,k+1,1)
1346 ptem = qcdo(i,k,1) + qcdo(i,k+1,1)
1347 f2(i,k) = f2(i,k) + (ptem - tem) * ptem1
1348 f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2
1349 endif
1350 endif
1351 enddo
1352 enddo
1353!
1354 if(ntrac1 >= 2) then
1355 do kk = 2, ntrac1
1356 is = (kk-1) * km
1357 do k = 1, km1
1358 do i = 1, im
1359 if(pcnvflg(i) .and. k < kpbl(i)) then
1360 dtodsd = dt2/del(i,k)
1361 dtodsu = dt2/del(i,k+1)
1362 dsig = prsl(i,k)-prsl(i,k+1)
1363 tem = dsig * rdzt(i,k)
1364 ptem = 0.5 * tem * xmf(i,k)
1365 ptem1 = dtodsd * ptem
1366 ptem2 = dtodsu * ptem
1367 tem1 = qcko(i,k,kk) + qcko(i,k+1,kk)
1368 tem2 = q1(i,k,kk) + q1(i,k+1,kk)
1369 f2(i,k+is) = f2(i,k+is) - (tem1 - tem2) * ptem1
1370 f2(i,k+1+is)= q1(i,k+1,kk) + (tem1 - tem2) * ptem2
1371 else
1372 f2(i,k+1+is) = q1(i,k+1,kk)
1373 endif
1374!
1375 if(scuflg(i)) then
1376 if(k >= mrad(i) .and. k < krad(i)) then
1377 dtodsd = dt2/del(i,k)
1378 dtodsu = dt2/del(i,k+1)
1379 dsig = prsl(i,k)-prsl(i,k+1)
1380 tem = dsig * rdzt(i,k)
1381 ptem = 0.5 * tem * xmfd(i,k)
1382 ptem1 = dtodsd * ptem
1383 ptem2 = dtodsu * ptem
1384 tem1 = qcdo(i,k,kk) + qcdo(i,k+1,kk)
1385 tem2 = q1(i,k,kk) + q1(i,k+1,kk)
1386 f2(i,k+is) = f2(i,k+is) + (tem1 - tem2) * ptem1
1387 f2(i,k+1+is)= f2(i,k+1+is) - (tem1 - tem2) * ptem2
1388 endif
1389 endif
1390!
1391 enddo
1392 enddo
1393 enddo
1394 endif
1395!
1397!
1398 call tridin(im,km,ntrac1,al,ad,au,f1,f2,au,f1,f2)
1399!
1400! recover tendencies of heat and moisture
1401!
1402 do k = 1,km
1403 do i = 1,im
1404 ttend = (f1(i,k)-t1(i,k))*rdt
1405 qtend = (f2(i,k)-q1(i,k,1))*rdt
1406 tdt(i,k) = tdt(i,k)+ttend
1407 rtg(i,k,1) = rtg(i,k,1)+qtend
1408 dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend
1409 dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend
1410 enddo
1411 enddo
1412 if (ldiag3d .and. .not. gen_tend) then
1413 idtend = dtidx(index_of_temperature,index_of_process_pbl)
1414 if(idtend>=1) then
1415 dtend(:,:,idtend) = dtend(:,:,idtend) + delt*rdt*(f1-t1)
1416 endif
1417 idtend = dtidx(100+ntqv,index_of_process_pbl)
1418 if(idtend>=1) then
1419 dtend(:,:,idtend) = dtend(:,:,idtend) +delt*rdt*(f2-q1(:,:,1))
1420 endif
1421 endif
1422!
1423 if(ntrac1 >= 2) then
1424 do kk = 2, ntrac1
1425 is = (kk-1) * km
1426 do k = 1, km
1427 do i = 1, im
1428 qtend = (f2(i,k+is)-q1(i,k,kk))*rdt
1429 rtg(i,k,kk) = rtg(i,k,kk)+qtend
1430 enddo
1431 enddo
1432 enddo
1433 endif
1434!
1436!
1437 if(dspheat) then
1438 do k = 1,km1
1439 do i = 1,im
1440! tem = min(diss(i,k), dspmax)
1441! ttend = tem / cp
1442 ttend = diss(i,k) / cp
1443 tdt(i,k) = tdt(i,k) + dspfac * ttend
1444 enddo
1445 enddo
1446 endif
1447!
1449!
1450 do i=1,im
1451 ad(i,1) = 1.0 + dtdz1(i) * stress(i) / spd1(i)
1452 f1(i,1) = u1(i,1)
1453 f2(i,1) = v1(i,1)
1454 enddo
1455!
1456 do k = 1,km1
1457 do i=1,im
1458 dtodsd = dt2/del(i,k)
1459 dtodsu = dt2/del(i,k+1)
1460 dsig = prsl(i,k)-prsl(i,k+1)
1461 rdz = rdzt(i,k)
1462 tem1 = dsig * dku(i,k) * rdz
1463 dsdz2 = tem1*rdz
1464 au(i,k) = -dtodsd*dsdz2
1465 al(i,k) = -dtodsu*dsdz2
1466 ad(i,k) = ad(i,k)-au(i,k)
1467 ad(i,k+1)= 1.-al(i,k)
1468 tem2 = dsig * rdz
1469!
1470 if(pcnvflg(i) .and. k < kpbl(i)) then
1471 ptem = 0.5 * tem2 * xmf(i,k)
1472 ptem1 = dtodsd * ptem
1473 ptem2 = dtodsu * ptem
1474 tem = u1(i,k) + u1(i,k+1)
1475 ptem = ucko(i,k) + ucko(i,k+1)
1476 f1(i,k) = f1(i,k) - (ptem - tem) * ptem1
1477 f1(i,k+1) = u1(i,k+1) + (ptem - tem) * ptem2
1478 tem = v1(i,k) + v1(i,k+1)
1479 ptem = vcko(i,k) + vcko(i,k+1)
1480 f2(i,k) = f2(i,k) - (ptem - tem) * ptem1
1481 f2(i,k+1) = v1(i,k+1) + (ptem - tem) * ptem2
1482 else
1483 f1(i,k+1) = u1(i,k+1)
1484 f2(i,k+1) = v1(i,k+1)
1485 endif
1486!
1487 if(scuflg(i)) then
1488 if(k >= mrad(i) .and. k < krad(i)) then
1489 ptem = 0.5 * tem2 * xmfd(i,k)
1490 ptem1 = dtodsd * ptem
1491 ptem2 = dtodsu * ptem
1492 tem = u1(i,k) + u1(i,k+1)
1493 ptem = ucdo(i,k) + ucdo(i,k+1)
1494 f1(i,k) = f1(i,k) + (ptem - tem) *ptem1
1495 f1(i,k+1) = f1(i,k+1) - (ptem - tem) *ptem2
1496 tem = v1(i,k) + v1(i,k+1)
1497 ptem = vcdo(i,k) + vcdo(i,k+1)
1498 f2(i,k) = f2(i,k) + (ptem - tem) * ptem1
1499 f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2
1500 endif
1501 endif
1502!
1503 enddo
1504 enddo
1505!
1507!
1508 call tridi2(im,km,al,ad,au,f1,f2,au,f1,f2)
1509!
1510! recover tendencies of momentum
1511!
1512 do k = 1,km
1513 do i = 1,im
1514 utend = (f1(i,k)-u1(i,k))*rdt
1515 vtend = (f2(i,k)-v1(i,k))*rdt
1516 du(i,k) = du(i,k)+utend
1517 dv(i,k) = dv(i,k)+vtend
1518 dusfc(i) = dusfc(i)+conw*del(i,k)*utend
1519 dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend
1520 enddo
1521 enddo
1522 if (ldiag3d .and. .not. gen_tend) then
1523 idtend=dtidx(index_of_x_wind,index_of_process_pbl)
1524 if(idtend>=1) then
1525 dtend(:,:,idtend) = dtend(:,:,idtend) + delt*rdt*(f1-u1)
1526 endif
1527 idtend=dtidx(index_of_y_wind,index_of_process_pbl)
1528 if(idtend>=1) then
1529 dtend(:,:,idtend) = dtend(:,:,idtend) + delt*rdt*(f2-v1)
1530 endif
1531 endif
1532!
1533!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1535!
1536 do i = 1, im
1537 hpbl(i) = hpblx(i)
1538 kpbl(i) = kpblx(i)
1539 enddo
1540!
1541!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1542 return
1543 end subroutine satmedmfvdif_run
1544
1545 end module satmedmfvdif
subroutine tridit(l, n, nt, cl, cm, cu, rt, au, at)
This subroutine solves tridiagonal problem for TKE.
Definition tridi.f:164
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
subroutine mfscu(im, ix, km, kmscu, ntcw, ntrac1, delt, cnvflg, zl, zm, q1, t1, u1, v1, plyr, pix, thlx, thvx, thlvx, gdx, thetae, radj, krad, mrad, radmin, buo, xmfd, tcdo, qcdo, ucdo, vcdo, xlamde)
This subroutine computes mass flux and downdraft parcel properties for stratocumulus-top-driven turbu...
Definition mfscu.f:16
subroutine mfpblt(im, ix, km, kmpbl, ntcw, ntrac1, delt, cnvflg, zl, zm, q1, t1, u1, v1, plyr, pix, thlx, thvx, gdx, hpbl, kpbl, vpert, buo, xmf, tcko, qcko, ucko, vcko, xlamue)
This subroutine computes mass flux and updraft parcel properties for thermals driven by surface heati...
Definition mfpblt.f:20
This module contains the subroutine that calculates mass flux and updraft parcel properties for therm...
Definition mfpblt.f:9