| blob | 1e8b9f4b5297f5af8b02c93ade0dc23abacc3a83 |
1 !WRF:model_layer:physics
2 !
3 !
4 !
5 !
6 !
7 !
8 !
9 !
10 module module_bl_ysu
11 contains
12 !
13 !
14 !-------------------------------------------------------------------
15 !
16 subroutine ysu(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, &
17 rublten,rvblten,rthblten, &
18 rqvblten,rqcblten,rqiblten,flag_qi, &
19 cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, &
20 dz8w,psfc, &
21 znu,znw,mut,p_top, &
22 znt,ust,hpbl,psim,psih, &
23 xland,hfx,qfx,gz1oz0,wspd,br, &
24 dt,kpbl2d, &
25 exch_h, &
26 u10,v10, &
27 ctopo,ctopo2, &
28 ids,ide, jds,jde, kds,kde, &
29 ims,ime, jms,jme, kms,kme, &
30 its,ite, jts,jte, kts,kte, &
31 !optional
32 regime )
33 !-------------------------------------------------------------------
34 implicit none
35 !-------------------------------------------------------------------
36 !-- u3d 3d u-velocity interpolated to theta points (m/s)
37 !-- v3d 3d v-velocity interpolated to theta points (m/s)
38 !-- th3d 3d potential temperature (k)
39 !-- t3d temperature (k)
40 !-- qv3d 3d water vapor mixing ratio (kg/kg)
41 !-- qc3d 3d cloud mixing ratio (kg/kg)
42 !-- qi3d 3d ice mixing ratio (kg/kg)
43 ! (note: if P_QI<PARAM_FIRST_SCALAR this should be zero filled)
44 !-- p3d 3d pressure (pa)
45 !-- p3di 3d pressure (pa) at interface level
46 !-- pi3d 3d exner function (dimensionless)
47 !-- rr3d 3d dry air density (kg/m^3)
48 !-- rublten u tendency due to
49 ! pbl parameterization (m/s/s)
50 !-- rvblten v tendency due to
51 ! pbl parameterization (m/s/s)
52 !-- rthblten theta tendency due to
53 ! pbl parameterization (K/s)
54 !-- rqvblten qv tendency due to
55 ! pbl parameterization (kg/kg/s)
56 !-- rqcblten qc tendency due to
57 ! pbl parameterization (kg/kg/s)
58 !-- rqiblten qi tendency due to
59 ! pbl parameterization (kg/kg/s)
60 !-- cp heat capacity at constant pressure for dry air (j/kg/k)
61 !-- g acceleration due to gravity (m/s^2)
62 !-- rovcp r/cp
63 !-- rd gas constant for dry air (j/kg/k)
64 !-- rovg r/g
65 !-- dz8w dz between full levels (m)
66 !-- xlv latent heat of vaporization (j/kg)
67 !-- rv gas constant for water vapor (j/kg/k)
68 !-- psfc pressure at the surface (pa)
69 !-- znt roughness length (m)
70 !-- ust u* in similarity theory (m/s)
71 !-- hpbl pbl height (m)
72 !-- psim similarity stability function for momentum
73 !-- psih similarity stability function for heat
74 !-- xland land mask (1 for land, 2 for water)
75 !-- hfx upward heat flux at the surface (w/m^2)
76 !-- qfx upward moisture flux at the surface (kg/m^2/s)
77 !-- gz1oz0 log(z/z0) where z0 is roughness length
78 !-- wspd wind speed at lowest model level (m/s)
79 !-- u10 u-wind speed at 10 m (m/s)
80 !-- v10 v-wind speed at 10 m (m/s)
81 !-- br bulk richardson number in surface layer
82 !-- dt time step (s)
83 !-- rvovrd r_v divided by r_d (dimensionless)
84 !-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
85 !-- ep2 constant for specific humidity calculation
86 !-- karman von karman constant
87 !-- ids start index for i in domain
88 !-- ide end index for i in domain
89 !-- jds start index for j in domain
90 !-- jde end index for j in domain
91 !-- kds start index for k in domain
92 !-- kde end index for k in domain
93 !-- ims start index for i in memory
94 !-- ime end index for i in memory
95 !-- jms start index for j in memory
96 !-- jme end index for j in memory
97 !-- kms start index for k in memory
98 !-- kme end index for k in memory
99 !-- its start index for i in tile
100 !-- ite end index for i in tile
101 !-- jts start index for j in tile
102 !-- jte end index for j in tile
103 !-- kts start index for k in tile
104 !-- kte end index for k in tile
105 !-------------------------------------------------------------------
106 !
107 integer,parameter :: ndiff = 3
108 real,parameter :: rcl = 1.0
109 !
110 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
111 ims,ime, jms,jme, kms,kme, &
112 its,ite, jts,jte, kts,kte
113 !
114 real, intent(in ) :: dt,cp,g,rovcp,rovg,rd,xlv,rv
115 !
116 real, intent(in ) :: ep1,ep2,karman
117 !
118 real, dimension( ims:ime, kms:kme, jms:jme ) , &
119 intent(in ) :: qv3d, &
120 qc3d, &
121 qi3d, &
122 p3d, &
123 pi3d, &
124 th3d, &
125 t3d, &
126 dz8w
127 real, dimension( ims:ime, kms:kme, jms:jme ) , &
128 intent(in ) :: p3di
129 !
130 real, dimension( ims:ime, kms:kme, jms:jme ) , &
131 intent(inout) :: rublten, &
132 rvblten, &
133 rthblten, &
134 rqvblten, &
135 rqcblten
136 !
137 real, dimension( ims:ime, kms:kme, jms:jme ) , &
138 intent(inout) :: exch_h
139 real, dimension( ims:ime, jms:jme ) , &
140 intent(inout) :: u10, &
141 v10
142 !
143 real, dimension( ims:ime, jms:jme ) , &
144 intent(in ) :: xland, &
145 hfx, &
146 qfx, &
147 br, &
148 psfc
149 real, dimension( ims:ime, jms:jme ) , &
150 intent(in ) :: &
151 psim, &
152 psih, &
153 gz1oz0
154 real, dimension( ims:ime, jms:jme ) , &
155 intent(inout) :: znt, &
156 ust, &
157 hpbl, &
158 wspd
159 !
160 real, dimension( ims:ime, kms:kme, jms:jme ) , &
161 intent(in ) :: u3d, &
162 v3d
163 !
164 integer, dimension( ims:ime, jms:jme ) , &
165 intent(out ) :: kpbl2d
166 logical, intent(in) :: flag_qi
167 !optional
168 real, dimension( ims:ime, jms:jme ) , &
169 optional , &
170 intent(inout) :: regime
171 !
172 real, dimension( ims:ime, kms:kme, jms:jme ) , &
173 optional , &
174 intent(inout) :: rqiblten
175 !
176 real, dimension( kms:kme ) , &
177 optional , &
178 intent(in ) :: znu, &
179 znw
180 !
181 real, dimension( ims:ime, jms:jme ) , &
182 optional , &
183 intent(in ) :: mut
184 !
185 real, optional, intent(in ) :: p_top
186 !
187 real, dimension( ims:ime, jms:jme ) , &
188 optional , &
189 intent(in ) :: ctopo, &
190 ctopo2
191 !local
192 integer :: i,j,k
193 real, dimension( its:ite, kts:kte*ndiff ) :: rqvbl2dt, &
194 qv2d
195 real, dimension( its:ite, kts:kte ) :: pdh
196 real, dimension( its:ite, kts:kte+1 ) :: pdhi
197 real, dimension( its:ite ) :: &
198 dusfc, &
199 dvsfc, &
200 dtsfc, &
201 dqsfc
202 !
203 qv2d(:,:) = 0.0
204 do j = jts,jte
205 if(present(mut))then
206 ! For ARW we will replace p and p8w with dry hydrostatic pressure
207 do k = kts,kte+1
208 do i = its,ite
209 if(k.le.kte)pdh(i,k) = mut(i,j)*znu(k) + p_top
210 pdhi(i,k) = mut(i,j)*znw(k) + p_top
211 enddo
212 enddo
213 else
214 do k = kts,kte+1
215 do i = its,ite
216 if(k.le.kte)pdh(i,k) = p3d(i,k,j)
217 pdhi(i,k) = p3di(i,k,j)
218 enddo
219 enddo
220 endif
221 do k = kts,kte
222 do i = its,ite
223 qv2d(i,k) = qv3d(i,k,j)
224 qv2d(i,k+kte) = qc3d(i,k,j)
225 if(present(rqiblten)) qv2d(i,k+kte+kte) = qi3d(i,k,j)
226 enddo
227 enddo
228 !
229 call ysu2d(J=j,ux=u3d(ims,kms,j),vx=v3d(ims,kms,j) &
230 ,tx=t3d(ims,kms,j) &
231 ,qx=qv2d(its,kts) &
232 ,p2d=pdh(its,kts),p2di=pdhi(its,kts) &
233 ,pi2d=pi3d(ims,kms,j) &
234 ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) &
235 ,ttnp=rthblten(ims,kms,j),qtnp=rqvbl2dt(its,kts),ndiff=ndiff &
236 ,cp=cp,g=g,rovcp=rovcp,rd=rd,rovg=rovg &
237 ,xlv=xlv,rv=rv &
238 ,ep1=ep1,ep2=ep2,karman=karman &
239 ,dz8w2d=dz8w(ims,kms,j) &
240 ,psfcpa=psfc(ims,j),znt=znt(ims,j),ust=ust(ims,j) &
241 ,hpbl=hpbl(ims,j) &
242 ,regime=regime(ims,j),psim=psim(ims,j) &
243 ,psih=psih(ims,j),xland=xland(ims,j) &
244 ,hfx=hfx(ims,j),qfx=qfx(ims,j) &
245 ,wspd=wspd(ims,j),br=br(ims,j) &
246 ,dusfc=dusfc,dvsfc=dvsfc,dtsfc=dtsfc,dqsfc=dqsfc &
247 ,dt=dt,rcl=1.0,kpbl1d=kpbl2d(ims,j) &
248 ,exch_hx=exch_h(ims,kms,j) &
249 ,u10=u10(ims,j),v10=v10(ims,j) &
250 #if ( NMM_CORE != 1 )
251 ,ctopo=ctopo(ims,j),ctopo2=ctopo2(ims,j) &
252 #endif
253 ,gz1oz0=gz1oz0(ims,j) &
254 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
255 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
256 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
257 !
258 do k = kts,kte
259 do i = its,ite
260 rthblten(i,k,j) = rthblten(i,k,j)/pi3d(i,k,j)
261 rqvblten(i,k,j) = rqvbl2dt(i,k)
262 rqcblten(i,k,j) = rqvbl2dt(i,k+kte)
263 if(present(rqiblten)) rqiblten(i,k,j) = rqvbl2dt(i,k+kte+kte)
264 enddo
265 enddo
266 enddo
267 !
268 end subroutine ysu
269 !
270 !-------------------------------------------------------------------
271 !
272 subroutine ysu2d(j,ux,vx,tx,qx,p2d,p2di,pi2d, &
273 utnp,vtnp,ttnp,qtnp,ndiff, &
274 cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, &
275 dz8w2d,psfcpa, &
276 znt,ust,hpbl,psim,psih, &
277 xland,hfx,qfx,wspd,br, &
278 dusfc,dvsfc,dtsfc,dqsfc, &
279 dt,rcl,kpbl1d, &
280 exch_hx, &
281 u10,v10, &
282 ctopo,ctopo2, &
283 gz1oz0, &
284 ids,ide, jds,jde, kds,kde, &
285 ims,ime, jms,jme, kms,kme, &
286 its,ite, jts,jte, kts,kte, &
287 !optional
288 regime )
289 !-------------------------------------------------------------------
290 implicit none
291 !-------------------------------------------------------------------
292 !
293 ! this code is a revised vertical diffusion package ("ysupbl")
294 ! with a nonlocal turbulent mixing in the pbl after "mrfpbl".
295 ! the ysupbl (hong et al. 2006) is based on the study of noh
296 ! et al.(2003) and accumulated realism of the behavior of the
297 ! troen and mahrt (1986) concept implemented by hong and pan(1996).
298 ! the major ingredient of the ysupbl is the inclusion of an explicit
299 ! treatment of the entrainment processes at the entrainment layer.
300 ! this routine uses an implicit approach for vertical flux
301 ! divergence and does not require "miter" timesteps.
302 ! it includes vertical diffusion in the stable atmosphere
303 ! and moist vertical diffusion in clouds.
304 !
305 ! mrfpbl:
306 ! coded by song-you hong (ncep), implemented by jimy dudhia (ncar)
307 ! fall 1996
308 !
309 ! ysupbl:
310 ! coded by song-you hong (yonsei university) and implemented by
311 ! song-you hong (yonsei university) and jimy dudhia (ncar)
312 ! summer 2002
313 !
314 ! further modifications :
315 ! an enhanced stable layer mixing, april 2008
316 ! ==> increase pbl height when sfc is stable (hong 2010)
317 ! pressure-level diffusion, april 2009
318 ! ==> negligible differences
319 ! implicit forcing for momentum with clean up, july 2009
320 ! ==> prevents model blowup when sfc layer is too low
321 ! increase of lamda, maximum (30, 0.1 x del z) feb 2010
322 ! ==> prevents model blowup when delz is extremely large
323 ! revised prandtl number at surface, peggy lemone, feb 2010
324 ! ==> increase kh, decrease mixing due to counter-gradient term
325 ! revised thermal, shin et al. mon. wea. rev. , songyou hong, aug 2011
326 ! ==> reduce the thermal strength when z1 < 0.1 h
327 ! revised prandtl number for free convection, dudhia, mar 2012
328 ! ==> pr0 = 1 + bke (=0.272) when newtral, kh is reduced
329 ! minimum kzo = 0.01, lo = min (30m,delz), hong, mar 2012
330 ! ==> weaker mixing when stable, and les resolution in vertical
331 !
332 ! references:
333 !
334 ! hong (2010) quart. j. roy. met. soc
335 ! hong, noh, and dudhia (2006), mon. wea. rev.
336 ! hong and pan (1996), mon. wea. rev.
337 ! noh, chun, hong, and raasch (2003), boundary layer met.
338 ! troen and mahrt (1986), boundary layer met.
339 !
340 !-------------------------------------------------------------------
341 !
342 real,parameter :: xkzmin = 0.01,xkzmax = 1000.,rimin = -100.
343 real,parameter :: rlam = 30.,prmin = 0.25,prmax = 4.
344 real,parameter :: brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4
345 real,parameter :: afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0
346 real,parameter :: phifac = 8.,sfcfrac = 0.1
347 real,parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001
348 real,parameter :: h1 = 0.33333335, h2 = 0.6666667
349 real,parameter :: ckz = 0.001,zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16.
350 real,parameter :: tmin=1.e-2
351 real,parameter :: gamcrt = 3.,gamcrq = 2.e-3
352 real,parameter :: xka = 2.4e-5
353 integer,parameter :: imvdif = 1
354 !
355 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
356 ims,ime, jms,jme, kms,kme, &
357 its,ite, jts,jte, kts,kte, &
358 j,ndiff
359 !
360 real, intent(in ) :: dt,rcl,cp,g,rovcp,rovg,rd,xlv,rv
361 !
362 real, intent(in ) :: ep1,ep2,karman
363 !
364 real, dimension( ims:ime, kms:kme ), &
365 intent(in) :: dz8w2d, &
366 pi2d
367 !
368 real, dimension( ims:ime, kms:kme ) , &
369 intent(in ) :: tx
370 real, dimension( its:ite, kts:kte*ndiff ) , &
371 intent(in ) :: qx
372 !
373 real, dimension( ims:ime, kms:kme ) , &
374 intent(inout) :: utnp, &
375 vtnp, &
376 ttnp
377 real, dimension( its:ite, kts:kte*ndiff ) , &
378 intent(inout) :: qtnp
379 !
380 real, dimension( its:ite, kts:kte+1 ) , &
381 intent(in ) :: p2di
382 !
383 real, dimension( its:ite, kts:kte ) , &
384 intent(in ) :: p2d
385 !
386 !
387 real, dimension( ims:ime ) , &
388 intent(inout) :: ust, &
389 hpbl, &
390 znt
391 real, dimension( ims:ime ) , &
392 intent(in ) :: xland, &
393 hfx, &
394 qfx
395 !
396 real, dimension( ims:ime ), intent(inout) :: wspd
397 real, dimension( ims:ime ), intent(in ) :: br
398 !
399 real, dimension( ims:ime ), intent(in ) :: psim, &
400 psih
401 real, dimension( ims:ime ), intent(in ) :: gz1oz0
402 !
403 real, dimension( ims:ime ), intent(in ) :: psfcpa
404 integer, dimension( ims:ime ), intent(out ) :: kpbl1d
405 !
406 real, dimension( ims:ime, kms:kme ) , &
407 intent(in ) :: ux, &
408 vx
409 real, dimension( ims:ime ) , &
410 optional , &
411 intent(in ) :: ctopo, &
412 ctopo2
413 real, dimension( ims:ime ) , &
414 optional , &
415 intent(inout) :: regime
416 !
417 ! local vars
418 !
419 real, dimension( its:ite ) :: hol
420 real, dimension( its:ite, kts:kte+1 ) :: zq
421 !
422 real, dimension( its:ite, kts:kte ) :: &
423 thx,thvx, &
424 del, &
425 dza, &
426 dzq, &
427 xkzo, &
428 za
429 !
430 real, dimension( its:ite ) :: &
431 rhox, &
432 govrth, &
433 zl1,thermal, &
434 wscale, &
435 hgamt,hgamq, &
436 brdn,brup, &
437 phim,phih, &
438 dusfc,dvsfc, &
439 dtsfc,dqsfc, &
440 prpbl, &
441 wspd1
442 !
443 real, dimension( its:ite, kts:kte ) :: xkzm,xkzh, &
444 f1,f2, &
445 r1,r2, &
446 ad,au, &
447 cu, &
448 al, &
449 xkzq, &
450 zfac
451 !
452 !jdf added exch_hx
453 real, dimension( ims:ime, kms:kme ) , &
454 intent(inout) :: exch_hx
455 !
456 real, dimension( ims:ime ) , &
457 intent(inout) :: u10, &
458 v10
459 real, dimension( its:ite ) :: &
460 brcr, &
461 sflux, &
462 brcr_sbro
463 !
464 real, dimension( its:ite, kts:kte, ndiff) :: r3,f3
465 integer, dimension( its:ite ) :: kpbl
466 !
467 logical, dimension( its:ite ) :: pblflg, &
468 sfcflg, &
469 stable
470 !
471 integer :: n,i,k,l,ic,is
472 integer :: klpbl, ktrace1, ktrace2, ktrace3
473 !
474 !
475 real :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0
476 real :: ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri
477 real :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz
478 real :: utend,vtend,ttend,qtend
479 real :: dtstep,govrthv
480 real :: cont, conq, conw, conwrc
481 !
482 real, dimension( its:ite, kts:kte ) :: wscalek
483 real, dimension( its:ite ) :: delta
484 real, dimension( its:ite, kts:kte ) :: xkzml,xkzhl, &
485 zfacent,entfac
486 real, dimension( its:ite ) :: ust3, &
487 wstar3,wstar, &
488 hgamu,hgamv, &
489 wm2, we, &
490 bfxpbl, &
491 hfxpbl,qfxpbl, &
492 ufxpbl,vfxpbl, &
493 dthvx
494 real :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, &
495 dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr, &
496 prfac,prfac2,phim8z
497 !
498 !----------------------------------------------------------------------
499 !
500 klpbl = kte
501 !
502 cont=cp/g
503 conq=xlv/g
504 conw=1./g
505 conwrc = conw*sqrt(rcl)
506 conpr = bfac*karman*sfcfrac
507 !
508 ! k-start index for tracer diffusion
509 !
510 ktrace1 = 0
511 ktrace2 = 0 + kte
512 ktrace3 = 0 + kte*2
513 !
514 do k = kts,kte
515 do i = its,ite
516 thx(i,k) = tx(i,k)/pi2d(i,k)
517 enddo
518 enddo
519 !
520 do k = kts,kte
521 do i = its,ite
522 tvcon = (1.+ep1*qx(i,k))
523 thvx(i,k) = thx(i,k)*tvcon
524 enddo
525 enddo
526 !
527 do i = its,ite
528 tvcon = (1.+ep1*qx(i,1))
529 rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon)
530 govrth(i) = g/thx(i,1)
531 enddo
532 !
533 !-----compute the height of full- and half-sigma levels above ground
534 ! level, and the layer thicknesses.
535 !
536 do i = its,ite
537 zq(i,1) = 0.
538 enddo
539 !
540 do k = kts,kte
541 do i = its,ite
542 zq(i,k+1) = dz8w2d(i,k)+zq(i,k)
543 enddo
544 enddo
545 !
546 do k = kts,kte
547 do i = its,ite
548 za(i,k) = 0.5*(zq(i,k)+zq(i,k+1))
549 dzq(i,k) = zq(i,k+1)-zq(i,k)
550 del(i,k) = p2di(i,k)-p2di(i,k+1)
551 enddo
552 enddo
553 !
554 do i = its,ite
555 dza(i,1) = za(i,1)
556 enddo
557 !
558 do k = kts+1,kte
559 do i = its,ite
560 dza(i,k) = za(i,k)-za(i,k-1)
561 enddo
562 enddo
563 !
564 !
565 !-----initialize vertical tendencies and
566 !
567 utnp(:,:) = 0.
568 vtnp(:,:) = 0.
569 ttnp(:,:) = 0.
570 qtnp(:,:) = 0.
571 !
572 do i = its,ite
573 wspd1(i) = sqrt(ux(i,1)*ux(i,1)+vx(i,1)*vx(i,1))+1.e-9
574 enddo
575 !
576 !---- compute vertical diffusion
577 !
578 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
579 ! compute preliminary variables
580 !
581 dtstep = dt
582 dt2 = 2.*dtstep
583 rdt = 1./dt2
584 !
585 do i = its,ite
586 bfxpbl(i) = 0.0
587 hfxpbl(i) = 0.0
588 qfxpbl(i) = 0.0
589 ufxpbl(i) = 0.0
590 vfxpbl(i) = 0.0
591 hgamu(i) = 0.0
592 hgamv(i) = 0.0
593 delta(i) = 0.0
594 enddo
595 !
596 do k = kts,klpbl
597 do i = its,ite
598 wscalek(i,k) = 0.0
599 enddo
600 enddo
601 !
602 do k = kts,klpbl
603 do i = its,ite
604 zfac(i,k) = 0.0
605 enddo
606 enddo
607 do k = kts,klpbl-1
608 do i = its,ite
609 xkzo(i,k) = ckz*dza(i,k+1)
610 enddo
611 enddo
612 !
613 do i = its,ite
614 dusfc(i) = 0.
615 dvsfc(i) = 0.
616 dtsfc(i) = 0.
617 dqsfc(i) = 0.
618 enddo
619 !
620 do i = its,ite
621 hgamt(i) = 0.
622 hgamq(i) = 0.
623 wscale(i) = 0.
624 kpbl(i) = 1
625 hpbl(i) = zq(i,1)
626 zl1(i) = za(i,1)
627 thermal(i)= thvx(i,1)
628 pblflg(i) = .true.
629 sfcflg(i) = .true.
630 sflux(i) = hfx(i)/rhox(i)/cp + qfx(i)/rhox(i)*ep1*thx(i,1)
631 if(br(i).gt.0.0) sfcflg(i) = .false.
632 enddo
633 !
634 ! compute the first guess of pbl height
635 !
636 do i = its,ite
637 stable(i) = .false.
638 brup(i) = br(i)
639 brcr(i) = brcr_ub
640 enddo
641 !
642 do k = 2,klpbl
643 do i = its,ite
644 if(.not.stable(i))then
645 brdn(i) = brup(i)
646 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
647 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
648 kpbl(i) = k
649 stable(i) = brup(i).gt.brcr(i)
650 endif
651 enddo
652 enddo
653 !
654 do i = its,ite
655 k = kpbl(i)
656 if(brdn(i).ge.brcr(i))then
657 brint = 0.
658 elseif(brup(i).le.brcr(i))then
659 brint = 1.
660 else
661 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
662 endif
663 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
664 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
665 if(kpbl(i).le.1) pblflg(i) = .false.
666 enddo
667 !
668 do i = its,ite
669 fm = gz1oz0(i)-psim(i)
670 fh = gz1oz0(i)-psih(i)
671 hol(i) = max(br(i)*fm*fm/fh,rimin)
672 if(sfcflg(i))then
673 hol(i) = min(hol(i),-zfmin)
674 else
675 hol(i) = max(hol(i),zfmin)
676 endif
677 hol1 = hol(i)*hpbl(i)/zl1(i)*sfcfrac
678 hol(i) = -hol(i)*hpbl(i)/zl1(i)
679 if(sfcflg(i))then
680 phim(i) = (1.-aphi16*hol1)**(-1./4.)
681 phih(i) = (1.-aphi16*hol1)**(-1./2.)
682 bfx0 = max(sflux(i),0.)
683 hfx0 = max(hfx(i)/rhox(i)/cp,0.)
684 qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.)
685 wstar3(i) = (govrth(i)*bfx0*hpbl(i))
686 wstar(i) = (wstar3(i))**h1
687 else
688 phim(i) = (1.+aphi5*hol1)
689 phih(i) = phim(i)
690 wstar(i) = 0.
691 wstar3(i) = 0.
692 endif
693 ust3(i) = ust(i)**3.
694 wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1
695 wscale(i) = min(wscale(i),ust(i)*aphi16)
696 wscale(i) = max(wscale(i),ust(i)/aphi5)
697 enddo
698 !
699 ! compute the surface variables for pbl height estimation
700 ! under unstable conditions
701 !
702 do i = its,ite
703 if(sfcflg(i).and.sflux(i).gt.0.0)then
704 gamfac = bfac/rhox(i)/wscale(i)
705 hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
706 hgamq(i) = min(gamfac*qfx(i),gamcrq)
707 vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac
708 thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
709 hgamt(i) = max(hgamt(i),0.0)
710 hgamq(i) = max(hgamq(i),0.0)
711 brint = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.)
712 hgamu(i) = brint*ux(i,1)
713 hgamv(i) = brint*vx(i,1)
714 else
715 pblflg(i) = .false.
716 endif
717 enddo
718 !
719 ! enhance the pbl height by considering the thermal
720 !
721 do i = its,ite
722 if(pblflg(i))then
723 kpbl(i) = 1
724 hpbl(i) = zq(i,1)
725 endif
726 enddo
727 !
728 do i = its,ite
729 if(pblflg(i))then
730 stable(i) = .false.
731 brup(i) = br(i)
732 brcr(i) = brcr_ub
733 endif
734 enddo
735 !
736 do k = 2,klpbl
737 do i = its,ite
738 if(.not.stable(i).and.pblflg(i))then
739 brdn(i) = brup(i)
740 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
741 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
742 kpbl(i) = k
743 stable(i) = brup(i).gt.brcr(i)
744 endif
745 enddo
746 enddo
747 !
748 do i = its,ite
749 if(pblflg(i)) then
750 k = kpbl(i)
751 if(brdn(i).ge.brcr(i))then
752 brint = 0.
753 elseif(brup(i).le.brcr(i))then
754 brint = 1.
755 else
756 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
757 endif
758 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
759 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
760 if(kpbl(i).le.1) pblflg(i) = .false.
761 endif
762 enddo
763 !
764 ! stable boundary layer
765 !
766 do i = its,ite
767 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
768 brup(i) = br(i)
769 stable(i) = .false.
770 else
771 stable(i) = .true.
772 endif
773 enddo
774 !
775 do i = its,ite
776 if((.not.stable(i)).and.((xland(i)-1.5).ge.0))then
777 wspd10 = u10(i)*u10(i) + v10(i)*v10(i)
778 wspd10 = sqrt(wspd10)
779 ross = wspd10 / (cori*znt(i))
780 brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3)
781 endif
782 enddo
783 !
784 do i = its,ite
785 if(.not.stable(i))then
786 if((xland(i)-1.5).ge.0)then
787 brcr(i) = brcr_sbro(i)
788 else
789 brcr(i) = brcr_sb
790 endif
791 endif
792 enddo
793 do k = 2,klpbl
794 do i = its,ite
795 if(.not.stable(i))then
796 brdn(i) = brup(i)
797 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
798 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
799 kpbl(i) = k
800 stable(i) = brup(i).gt.brcr(i)
801 endif
802 enddo
803 enddo
804 !
805 do i = its,ite
806 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
807 k = kpbl(i)
808 if(brdn(i).ge.brcr(i))then
809 brint = 0.
810 elseif(brup(i).le.brcr(i))then
811 brint = 1.
812 else
813 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
814 endif
815 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
816 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
817 if(kpbl(i).le.1) pblflg(i) = .false.
818 endif
819 enddo
820 !
821 ! estimate the entrainment parameters
822 !
823 do i = its,ite
824 if(pblflg(i)) then
825 k = kpbl(i) - 1
826 prpbl(i) = 1.0
827 wm3 = wstar3(i) + 5. * ust3(i)
828 wm2(i) = wm3**h2
829 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
830 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin)
831 dthx = max(thx(i,k+1)-thx(i,k),tmin)
832 dqx = min(qx(i,k+1)-qx(i,k),0.0)
833 we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
834 hfxpbl(i) = we(i)*dthx
835 qfxpbl(i) = we(i)*dqx
836 !
837 dux = ux(i,k+1)-ux(i,k)
838 dvx = vx(i,k+1)-vx(i,k)
839 if(dux.gt.tmin) then
840 ufxpbl(i) = max(prpbl(i)*we(i)*dux,-ust(i)*ust(i))
841 elseif(dux.lt.-tmin) then
842 ufxpbl(i) = min(prpbl(i)*we(i)*dux,ust(i)*ust(i))
843 else
844 ufxpbl(i) = 0.0
845 endif
846 if(dvx.gt.tmin) then
847 vfxpbl(i) = max(prpbl(i)*we(i)*dvx,-ust(i)*ust(i))
848 elseif(dvx.lt.-tmin) then
849 vfxpbl(i) = min(prpbl(i)*we(i)*dvx,ust(i)*ust(i))
850 else
851 vfxpbl(i) = 0.0
852 endif
853 delb = govrth(i)*d3*hpbl(i)
854 delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.)
855 endif
856 enddo
857 !
858 do k = kts,klpbl
859 do i = its,ite
860 if(pblflg(i).and.k.ge.kpbl(i))then
861 entfac(i,k) = ((zq(i,k+1)-hpbl(i))/delta(i))**2.
862 else
863 entfac(i,k) = 1.e30
864 endif
865 enddo
866 enddo
867 !
868 ! compute diffusion coefficients below pbl
869 !
870 do k = kts,klpbl
871 do i = its,ite
872 if(k.lt.kpbl(i)) then
873 zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.)
874 zfacent(i,k) = (1.-zfac(i,k))**3.
875 wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
876 if(sfcflg(i)) then
877 prfac = conpr
878 prfac2 = 15.9*wstar3(i)/ust3(i)/(1.+4.*karman*wstar3(i)/ust3(i))
879 prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2.
880 else
881 prfac = 0.
882 prfac2 = 0.
883 prnumfac = 0.
884 phim8z = 1.+(phim(i)-1.)*zq(i,k+1)/sfcfrac/hpbl(i)
885 wscalek(i,k) = ust(i)/phim8z
886 wscalek(i,k) = max(wscalek(i,k),ust(i)/aphi5)
887 endif
888 prnum0 = (phih(i)/phim(i)+prfac)
889 prnum0 = min(prnum0,prmax)
890 prnum0 = max(prnum0,prmin)
891 xkzm(i,k) = wscalek(i,k)*karman*zq(i,k+1)*zfac(i,k)**pfac
892 prnum = 1. + (prnum0-1.)*exp(prnumfac)
893 xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac)
894 prnum0 = prnum0/(1.+prfac2*karman*sfcfrac)
895 prnum = 1. + (prnum0-1.)*exp(prnumfac)
896 xkzh(i,k) = xkzm(i,k)/prnum
897 xkzm(i,k) = min(xkzm(i,k),xkzmax)
898 xkzm(i,k) = max(xkzm(i,k),xkzo(i,k))
899 xkzh(i,k) = min(xkzh(i,k),xkzmax)
900 xkzh(i,k) = max(xkzh(i,k),xkzo(i,k))
901 xkzq(i,k) = min(xkzq(i,k),xkzmax)
902 xkzq(i,k) = max(xkzq(i,k),xkzo(i,k))
903 endif
904 enddo
905 enddo
906 !
907 ! compute diffusion coefficients over pbl (free atmosphere)
908 !
909 do k = kts,kte-1
910 do i = its,ite
911 if(k.ge.kpbl(i)) then
912 ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k)) &
913 +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k))) &
914 /(dza(i,k+1)*dza(i,k+1))+1.e-9
915 govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
916 ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1))
917 if(imvdif.eq.1.and.ndiff.ge.3)then
918 if((qx(i,ktrace2+k)+qx(i,ktrace3+k)).gt.0.01e-3.and.(qx(i &
919 ,ktrace2+k+1)+qx(i,ktrace3+k+1)).gt.0.01e-3)then
920 ! in cloud
921 qmean = 0.5*(qx(i,k)+qx(i,k+1))
922 tmean = 0.5*(tx(i,k)+tx(i,k+1))
923 alph = xlv*qmean/rd/tmean
924 chi = xlv*xlv*qmean/cp/rv/tmean/tmean
925 ri = (1.+alph)*(ri-g*g/ss/tmean/cp*((chi-alph)/(1.+chi)))
926 endif
927 endif
928 zk = karman*zq(i,k+1)
929 rlamdz = min(max(0.1*dza(i,k+1),rlam),300.)
930 rlamdz = min(dza(i,k+1),rlamdz)
931 rl2 = (zk*rlamdz/(rlamdz+zk))**2
932 dk = rl2*sqrt(ss)
933 if(ri.lt.0.)then
934 ! unstable regime
935 sri = sqrt(-ri)
936 xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
937 xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
938 else
939 ! stable regime
940 xkzh(i,k) = dk/(1+5.*ri)**2
941 prnum = 1.0+2.1*ri
942 prnum = min(prnum,prmax)
943 xkzm(i,k) = xkzh(i,k)*prnum
944 endif
945 !
946 xkzm(i,k) = min(xkzm(i,k),xkzmax)
947 xkzm(i,k) = max(xkzm(i,k),xkzo(i,k))
948 xkzh(i,k) = min(xkzh(i,k),xkzmax)
949 xkzh(i,k) = max(xkzh(i,k),xkzo(i,k))
950 xkzml(i,k) = xkzm(i,k)
951 xkzhl(i,k) = xkzh(i,k)
952 endif
953 enddo
954 enddo
955 !
956 ! compute tridiagonal matrix elements for heat
957 !
958 do k = kts,kte
959 do i = its,ite
960 au(i,k) = 0.
961 al(i,k) = 0.
962 ad(i,k) = 0.
963 f1(i,k) = 0.
964 enddo
965 enddo
966 !
967 do i = its,ite
968 ad(i,1) = 1.
969 f1(i,1) = thx(i,1)-300.+hfx(i)/cont/del(i,1)*dt2
970 enddo
971 !
972 do k = kts,kte-1
973 do i = its,ite
974 dtodsd = dt2/del(i,k)
975 dtodsu = dt2/del(i,k+1)
976 dsig = p2d(i,k)-p2d(i,k+1)
977 rdz = 1./dza(i,k+1)
978 tem1 = dsig*xkzh(i,k)*rdz
979 if(pblflg(i).and.k.lt.kpbl(i)) then
980 dsdzt = tem1*(-hgamt(i)/hpbl(i)-hfxpbl(i)*zfacent(i,k)/xkzh(i,k))
981 f1(i,k) = f1(i,k)+dtodsd*dsdzt
982 f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt
983 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
984 xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
985 xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
986 xkzh(i,k) = min(xkzh(i,k),xkzmax)
987 xkzh(i,k) = max(xkzh(i,k),xkzo(i,k))
988 f1(i,k+1) = thx(i,k+1)-300.
989 else
990 f1(i,k+1) = thx(i,k+1)-300.
991 endif
992 tem1 = dsig*xkzh(i,k)*rdz
993 dsdz2 = tem1*rdz
994 au(i,k) = -dtodsd*dsdz2
995 al(i,k) = -dtodsu*dsdz2
996 ad(i,k) = ad(i,k)-au(i,k)
997 ad(i,k+1) = 1.-al(i,k)
998 exch_hx(i,k+1) = xkzh(i,k)
999 enddo
1000 enddo
1001 !
1002 ! copies here to avoid duplicate input args for tridin
1003 !
1004 do k = kts,kte
1005 do i = its,ite
1006 cu(i,k) = au(i,k)
1007 r1(i,k) = f1(i,k)
1008 enddo
1009 enddo
1010 !
1011 call tridin_ysu(al,ad,cu,r1,au,f1,its,ite,kts,kte,1)
1012 !
1013 ! recover tendencies of heat
1014 !
1015 do k = kte,kts,-1
1016 do i = its,ite
1017 ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k)
1018 ttnp(i,k) = ttnp(i,k)+ttend
1019 dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)/pi2d(i,k)
1020 enddo
1021 enddo
1022 !
1023 ! compute tridiagonal matrix elements for moisture, clouds, and tracers
1024 !
1025 do k = kts,kte
1026 do i = its,ite
1027 au(i,k) = 0.
1028 al(i,k) = 0.
1029 ad(i,k) = 0.
1030 enddo
1031 enddo
1032 !
1033 do ic = 1,ndiff
1034 do i = its,ite
1035 do k = kts,kte
1036 f3(i,k,ic) = 0.
1037 enddo
1038 enddo
1039 enddo
1040 !
1041 do i = its,ite
1042 ad(i,1) = 1.
1043 f3(i,1,1) = qx(i,1)+qfx(i)*g/del(i,1)*dt2
1044 enddo
1045 !
1046 if(ndiff.ge.2) then
1047 do ic = 2,ndiff
1048 is = (ic-1) * kte
1049 do i = its,ite
1050 f3(i,1,ic) = qx(i,1+is)
1051 enddo
1052 enddo
1053 endif
1054 !
1055 do k = kts,kte
1056 do i = its,ite
1057 if(k.ge.kpbl(i)) then
1058 xkzq(i,k) = xkzh(i,k)
1059 endif
1060 enddo
1061 enddo
1062 !
1063 do k = kts,kte-1
1064 do i = its,ite
1065 dtodsd = dt2/del(i,k)
1066 dtodsu = dt2/del(i,k+1)
1067 dsig = p2d(i,k)-p2d(i,k+1)
1068 rdz = 1./dza(i,k+1)
1069 tem1 = dsig*xkzq(i,k)*rdz
1070 if(pblflg(i).and.k.lt.kpbl(i)) then
1071 dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzq(i,k))
1072 f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq
1073 f3(i,k+1,1) = qx(i,k+1)-dtodsu*dsdzq
1074 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1075 xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
1076 xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k))
1077 xkzq(i,k) = min(xkzq(i,k),xkzmax)
1078 xkzq(i,k) = max(xkzq(i,k),xkzo(i,k))
1079 f3(i,k+1,1) = qx(i,k+1)
1080 else
1081 f3(i,k+1,1) = qx(i,k+1)
1082 endif
1083 tem1 = dsig*xkzq(i,k)*rdz
1084 dsdz2 = tem1*rdz
1085 au(i,k) = -dtodsd*dsdz2
1086 al(i,k) = -dtodsu*dsdz2
1087 ad(i,k) = ad(i,k)-au(i,k)
1088 ad(i,k+1) = 1.-al(i,k)
1089 ! exch_hx(i,k+1) = xkzh(i,k)
1090 enddo
1091 enddo
1092 !
1093 if(ndiff.ge.2) then
1094 do ic = 2,ndiff
1095 is = (ic-1) * kte
1096 do k = kts,kte-1
1097 do i = its,ite
1098 f3(i,k+1,ic) = qx(i,k+1+is)
1099 enddo
1100 enddo
1101 enddo
1102 endif
1103 !
1104 ! copies here to avoid duplicate input args for tridin
1105 !
1106 do k = kts,kte
1107 do i = its,ite
1108 cu(i,k) = au(i,k)
1109 enddo
1110 enddo
1111 !
1112 do ic = 1,ndiff
1113 do k = kts,kte
1114 do i = its,ite
1115 r3(i,k,ic) = f3(i,k,ic)
1116 enddo
1117 enddo
1118 enddo
1119 !
1120 ! solve tridiagonal problem for moisture, clouds, and tracers
1121 !
1122 call tridin_ysu(al,ad,cu,r3,au,f3,its,ite,kts,kte,ndiff)
1123 !
1124 ! recover tendencies of heat and moisture
1125 !
1126 do k = kte,kts,-1
1127 do i = its,ite
1128 qtend = (f3(i,k,1)-qx(i,k))*rdt
1129 qtnp(i,k) = qtnp(i,k)+qtend
1130 dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k)
1131 enddo
1132 enddo
1133 !
1134 if(ndiff.ge.2) then
1135 do ic = 2,ndiff
1136 is = (ic-1) * kte
1137 do k = kte,kts,-1
1138 do i = its,ite
1139 qtend = (f3(i,k,ic)-qx(i,k+is))*rdt
1140 qtnp(i,k+is) = qtnp(i,k+is)+qtend
1141 enddo
1142 enddo
1143 enddo
1144 endif
1145 !
1146 ! compute tridiagonal matrix elements for momentum
1147 !
1148 do i = its,ite
1149 do k = kts,kte
1150 au(i,k) = 0.
1151 al(i,k) = 0.
1152 ad(i,k) = 0.
1153 f1(i,k) = 0.
1154 f2(i,k) = 0.
1155 enddo
1156 enddo
1157 !
1158 do i = its,ite
1159 ! paj: ctopo=1 if topo_wind=0 (default)
1160 ! mchen add this line to make sure NMM can still work with YSU PBL
1161 if(present(ctopo)) then
1162 ad(i,1) = 1.+ctopo(i)*ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 &
1163 *(wspd1(i)/wspd(i))**2
1164 else
1165 ad(i,1) = 1.+ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 &
1166 *(wspd1(i)/wspd(i))**2
1167 endif
1168 f1(i,1) = ux(i,1)
1169 f2(i,1) = vx(i,1)
1170 enddo
1171 !
1172 do k = kts,kte-1
1173 do i = its,ite
1174 dtodsd = dt2/del(i,k)
1175 dtodsu = dt2/del(i,k+1)
1176 dsig = p2d(i,k)-p2d(i,k+1)
1177 rdz = 1./dza(i,k+1)
1178 tem1 = dsig*xkzm(i,k)*rdz
1179 if(pblflg(i).and.k.lt.kpbl(i))then
1180 dsdzu = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k))
1181 dsdzv = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k))
1182 f1(i,k) = f1(i,k)+dtodsd*dsdzu
1183 f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu
1184 f2(i,k) = f2(i,k)+dtodsd*dsdzv
1185 f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv
1186 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1187 xkzm(i,k) = prpbl(i)*xkzh(i,k)
1188 xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
1189 xkzm(i,k) = min(xkzm(i,k),xkzmax)
1190 xkzm(i,k) = max(xkzm(i,k),xkzo(i,k))
1191 f1(i,k+1) = ux(i,k+1)
1192 f2(i,k+1) = vx(i,k+1)
1193 else
1194 f1(i,k+1) = ux(i,k+1)
1195 f2(i,k+1) = vx(i,k+1)
1196 endif
1197 tem1 = dsig*xkzm(i,k)*rdz
1198 dsdz2 = tem1*rdz
1199 au(i,k) = -dtodsd*dsdz2
1200 al(i,k) = -dtodsu*dsdz2
1201 ad(i,k) = ad(i,k)-au(i,k)
1202 ad(i,k+1) = 1.-al(i,k)
1203 enddo
1204 enddo
1205 !
1206 ! copies here to avoid duplicate input args for tridin
1207 !
1208 do k = kts,kte
1209 do i = its,ite
1210 cu(i,k) = au(i,k)
1211 r1(i,k) = f1(i,k)
1212 r2(i,k) = f2(i,k)
1213 enddo
1214 enddo
1215 !
1216 ! solve tridiagonal problem for momentum
1217 !
1218 call tridi1n(al,ad,cu,r1,r2,au,f1,f2,its,ite,kts,kte,1)
1219 !
1220 ! recover tendencies of momentum
1221 !
1222 do k = kte,kts,-1
1223 do i = its,ite
1224 utend = (f1(i,k)-ux(i,k))*rdt
1225 vtend = (f2(i,k)-vx(i,k))*rdt
1226 utnp(i,k) = utnp(i,k)+utend
1227 vtnp(i,k) = vtnp(i,k)+vtend
1228 dusfc(i) = dusfc(i) + utend*conwrc*del(i,k)
1229 dvsfc(i) = dvsfc(i) + vtend*conwrc*del(i,k)
1230 enddo
1231 enddo
1232 !
1233 ! paj: ctopo2=1 if topo_wind=0 (default)
1234 !
1235 do i = its,ite
1236 if(present(ctopo).and.present(ctopo2)) then ! mchen for NMM
1237 u10(i) = ctopo2(i)*u10(i)+(1-ctopo2(i))*ux(i,1)
1238 v10(i) = ctopo2(i)*v10(i)+(1-ctopo2(i))*vx(i,1)
1239 endif !mchen
1240 enddo
1241 !
1242 !---- end of vertical diffusion
1243 !
1244 do i = its,ite
1245 kpbl1d(i) = kpbl(i)
1246 enddo
1247 !
1248 end subroutine ysu2d
1249 !
1250 subroutine tridi1n(cl,cm,cu,r1,r2,au,f1,f2,its,ite,kts,kte,nt)
1251 !----------------------------------------------------------------
1252 implicit none
1253 !----------------------------------------------------------------
1254 !
1255 integer, intent(in ) :: its,ite, kts,kte, nt
1256 !
1257 real, dimension( its:ite, kts+1:kte+1 ) , &
1258 intent(in ) :: cl
1259 !
1260 real, dimension( its:ite, kts:kte ) , &
1261 intent(in ) :: cm, &
1262 r1
1263 real, dimension( its:ite, kts:kte,nt ) , &
1264 intent(in ) :: r2
1265 !
1266 real, dimension( its:ite, kts:kte ) , &
1267 intent(inout) :: au, &
1268 cu, &
1269 f1
1270 real, dimension( its:ite, kts:kte,nt ) , &
1271 intent(inout) :: f2
1272 !
1273 real :: fk
1274 integer :: i,k,l,n,it
1275 !
1276 !----------------------------------------------------------------
1277 !
1278 l = ite
1279 n = kte
1280 !
1281 do i = its,l
1282 fk = 1./cm(i,1)
1283 au(i,1) = fk*cu(i,1)
1284 f1(i,1) = fk*r1(i,1)
1285 enddo
1286 do it = 1,nt
1287 do i = its,l
1288 fk = 1./cm(i,1)
1289 f2(i,1,it) = fk*r2(i,1,it)
1290 enddo
1291 enddo
1292 do k = kts+1,n-1
1293 do i = its,l
1294 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1295 au(i,k) = fk*cu(i,k)
1296 f1(i,k) = fk*(r1(i,k)-cl(i,k)*f1(i,k-1))
1297 enddo
1298 enddo
1299 do it = 1,nt
1300 do k = kts+1,n-1
1301 do i = its,l
1302 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1303 f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
1304 enddo
1305 enddo
1306 enddo
1307 do i = its,l
1308 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1309 f1(i,n) = fk*(r1(i,n)-cl(i,n)*f1(i,n-1))
1310 enddo
1311 do it = 1,nt
1312 do i = its,l
1313 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1314 f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
1315 enddo
1316 enddo
1317 do k = n-1,kts,-1
1318 do i = its,l
1319 f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1)
1320 enddo
1321 enddo
1322 do it = 1,nt
1323 do k = n-1,kts,-1
1324 do i = its,l
1325 f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1326 enddo
1327 enddo
1328 enddo
1329 !
1330 end subroutine tridi1n
1331 !
1332 subroutine tridin_ysu(cl,cm,cu,r2,au,f2,its,ite,kts,kte,nt)
1333 !----------------------------------------------------------------
1334 implicit none
1335 !----------------------------------------------------------------
1336 !
1337 integer, intent(in ) :: its,ite, kts,kte, nt
1338 !
1339 real, dimension( its:ite, kts+1:kte+1 ) , &
1340 intent(in ) :: cl
1341 !
1342 real, dimension( its:ite, kts:kte ) , &
1343 intent(in ) :: cm
1344 real, dimension( its:ite, kts:kte,nt ) , &
1345 intent(in ) :: r2
1346 !
1347 real, dimension( its:ite, kts:kte ) , &
1348 intent(inout) :: au, &
1349 cu
1350 real, dimension( its:ite, kts:kte,nt ) , &
1351 intent(inout) :: f2
1352 !
1353 real :: fk
1354 integer :: i,k,l,n,it
1355 !
1356 !----------------------------------------------------------------
1357 !
1358 l = ite
1359 n = kte
1360 !
1361 do it = 1,nt
1362 do i = its,l
1363 fk = 1./cm(i,1)
1364 au(i,1) = fk*cu(i,1)
1365 f2(i,1,it) = fk*r2(i,1,it)
1366 enddo
1367 enddo
1368 do it = 1,nt
1369 do k = kts+1,n-1
1370 do i = its,l
1371 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1372 au(i,k) = fk*cu(i,k)
1373 f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
1374 enddo
1375 enddo
1376 enddo
1377 do it = 1,nt
1378 do i = its,l
1379 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1380 f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
1381 enddo
1382 enddo
1383 do it = 1,nt
1384 do k = n-1,kts,-1
1385 do i = its,l
1386 f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1387 enddo
1388 enddo
1389 enddo
1390 !
1391 end subroutine tridin_ysu
1392 !
1393 subroutine ysuinit(rublten,rvblten,rthblten,rqvblten, &
1394 rqcblten,rqiblten,p_qi,p_first_scalar, &
1395 restart, allowed_to_read, &
1396 ids, ide, jds, jde, kds, kde, &
1397 ims, ime, jms, jme, kms, kme, &
1398 its, ite, jts, jte, kts, kte )
1399 !-------------------------------------------------------------------
1400 implicit none
1401 !-------------------------------------------------------------------
1402 !
1403 logical , intent(in) :: restart, allowed_to_read
1404 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
1405 ims, ime, jms, jme, kms, kme, &
1406 its, ite, jts, jte, kts, kte
1407 integer , intent(in) :: p_qi,p_first_scalar
1408 real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: &
1409 rublten, &
1410 rvblten, &
1411 rthblten, &
1412 rqvblten, &
1413 rqcblten, &
1414 rqiblten
1415 integer :: i, j, k, itf, jtf, ktf
1416 !
1417 jtf = min0(jte,jde-1)
1418 ktf = min0(kte,kde-1)
1419 itf = min0(ite,ide-1)
1420 !
1421 if(.not.restart)then
1422 do j = jts,jtf
1423 do k = kts,ktf
1424 do i = its,itf
1425 rublten(i,k,j) = 0.
1426 rvblten(i,k,j) = 0.
1427 rthblten(i,k,j) = 0.
1428 rqvblten(i,k,j) = 0.
1429 rqcblten(i,k,j) = 0.
1430 enddo
1431 enddo
1432 enddo
1433 endif
1434 !
1435 if (p_qi .ge. p_first_scalar .and. .not.restart) then
1436 do j = jts,jtf
1437 do k = kts,ktf
1438 do i = its,itf
1439 rqiblten(i,k,j) = 0.
1440 enddo
1441 enddo
1442 enddo
1443 endif
1444 !
1445 end subroutine ysuinit
1446 !-------------------------------------------------------------------
1447 end module module_bl_ysu