| blob | ac95882eb4babd2367a54bb5dcd2242e9a184363 |
1 !WRF:model_layer:physics
2 !
3 module module_bl_shinhong
4 !
5 contains
6 !
7 !-------------------------------------------------------------------------------
8 !
9 subroutine shinhong(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, &
10 rublten,rvblten,rthblten, &
11 rqvblten,rqcblten,rqiblten,flag_qi, &
12 cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, &
13 dz8w,psfc, &
14 znu,znw,p_top, &
15 znt,ust,hpbl,psim,psih, &
16 xland,hfx,qfx,wspd,br, &
17 dt,kpbl2d, &
18 exch_h, &
19 u10,v10, &
20 shinhong_tke_diag,tke_pbl,el_pbl,corf, &
21 dx,dy, &
22 ids,ide, jds,jde, kds,kde, &
23 ims,ime, jms,jme, kms,kme, &
24 its,ite, jts,jte, kts,kte, &
25 !optional
26 ctopo,ctopo2, &
27 wstar,delta, &
28 regime )
29 !-------------------------------------------------------------------------------
30 implicit none
31 !-------------------------------------------------------------------------------
32 !-- u3d 3d u-velocity interpolated to theta points (m/s)
33 !-- v3d 3d v-velocity interpolated to theta points (m/s)
34 !-- th3d 3d potential temperature (k)
35 !-- t3d temperature (k)
36 !-- qv3d 3d water vapor mixing ratio (kg/kg)
37 !-- qc3d 3d cloud mixing ratio (kg/kg)
38 !-- qi3d 3d ice mixing ratio (kg/kg)
39 ! (note: if P_QI<PARAM_FIRST_SCALAR this should be zero filled)
40 !-- p3d 3d pressure (pa)
41 !-- p3di 3d pressure (pa) at interface level
42 !-- pi3d 3d exner function (dimensionless)
43 !-- rublten u tendency due to
44 ! pbl parameterization (m/s/s)
45 !-- rvblten v tendency due to
46 ! pbl parameterization (m/s/s)
47 !-- rthblten theta tendency due to
48 ! pbl parameterization (K/s)
49 !-- rqvblten qv tendency due to
50 ! pbl parameterization (kg/kg/s)
51 !-- rqcblten qc tendency due to
52 ! pbl parameterization (kg/kg/s)
53 !-- rqiblten qi tendency due to
54 ! pbl parameterization (kg/kg/s)
55 !-- cp heat capacity at constant pressure for dry air (j/kg/k)
56 !-- g acceleration due to gravity (m/s^2)
57 !-- rovcp r/cp
58 !-- rd gas constant for dry air (j/kg/k)
59 !-- rovg r/g
60 !-- ep1 constant for virtual temperature (r_v/r_d - 1)
61 !-- ep2 constant for specific humidity calculation
62 !-- karman von karman constant
63 !-- xlv latent heat of vaporization (j/kg)
64 !-- rv gas constant for water vapor (j/kg/k)
65 !-- dz8w dz between full levels (m)
66 !-- psfc pressure at the surface (pa)
67 !-- znu eta values on half (mass) levels
68 !-- znw eta values on full (w) levels
69 !-- p_top pressure top of the model (pa)
70 !-- znt roughness length (m)
71 !-- ust u* in similarity theory (m/s)
72 !-- hpbl pbl height (m)
73 !-- psim similarity stability function for momentum
74 !-- psih similarity stability function for heat
75 !-- xland land mask (1 for land, 2 for water)
76 !-- hfx upward heat flux at the surface (w/m^2)
77 !-- qfx upward moisture flux at the surface (kg/m^2/s)
78 !-- wspd wind speed at lowest model level (m/s)
79 !-- br bulk richardson number in surface layer
80 !-- u10 u-wind speed at 10 m (m/s)
81 !-- v10 v-wind speed at 10 m (m/s)
82 !-- dt time step (s)
83 !-- ids start index for i in domain
84 !-- ide end index for i in domain
85 !-- jds start index for j in domain
86 !-- jde end index for j in domain
87 !-- kds start index for k in domain
88 !-- kde end index for k in domain
89 !-- ims start index for i in memory
90 !-- ime end index for i in memory
91 !-- jms start index for j in memory
92 !-- jme end index for j in memory
93 !-- kms start index for k in memory
94 !-- kme end index for k in memory
95 !-- its start index for i in tile
96 !-- ite end index for i in tile
97 !-- jts start index for j in tile
98 !-- jte end index for j in tile
99 !-- kts start index for k in tile
100 !-- kte end index for k in tile
101 !-------------------------------------------------------------------------------
102 !
103 integer,parameter :: ndiff = 3
104 real,parameter :: rcl = 1.0
105 !
106 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
107 ims,ime, jms,jme, kms,kme, &
108 its,ite, jts,jte, kts,kte
109 integer, intent(in ) :: shinhong_tke_diag
110 !
111 real, intent(in ) :: dt,cp,g,rovcp,rovg,rd,xlv,rv
112 real, intent(in ) :: ep1,ep2,karman
113 real, intent(in ) :: dx,dy
114 !
115 integer, dimension( ims:ime, jms:jme ) , &
116 intent(out ) :: kpbl2d
117 !
118 real, dimension( ims:ime, kms:kme, jms:jme ) , &
119 intent(in ) :: u3d, &
120 v3d
121 real, dimension( ims:ime, kms:kme, jms:jme ) , &
122 intent(in ) :: qv3d, &
123 qc3d, &
124 qi3d, &
125 p3d, &
126 pi3d, &
127 th3d, &
128 t3d, &
129 dz8w
130 real, dimension( ims:ime, kms:kme, jms:jme ) , &
131 intent(in ) :: p3di
132 !
133 real, dimension( ims:ime, kms:kme, jms:jme ) , &
134 intent(inout) :: rublten, &
135 rvblten, &
136 rthblten, &
137 rqvblten, &
138 rqcblten
139 real, dimension( ims:ime, kms:kme, jms:jme ) , &
140 intent(inout) :: exch_h
141 real, dimension( ims:ime, kms:kme, jms:jme ) , &
142 intent(inout) :: tke_pbl, &
143 el_pbl
144 !
145 real, dimension( ims:ime, jms:jme ) , &
146 intent(in ) :: xland, &
147 hfx, &
148 qfx, &
149 corf, &
150 br, &
151 psfc
152 real, dimension( ims:ime, jms:jme ) , &
153 intent(in ) :: &
154 psim, &
155 psih
156 !
157 real, dimension( ims:ime, jms:jme ) , &
158 intent(inout) :: u10, &
159 v10
160 real, dimension( ims:ime, jms:jme ) , &
161 intent(inout) :: znt, &
162 ust, &
163 hpbl, &
164 wspd
165 !
166 logical, intent(in) :: flag_qi
167 !
168 ! optional
169 !
170 real, dimension( ims:ime, kms:kme, jms:jme ) , &
171 intent(inout), optional :: rqiblten
172 !
173 real, dimension( ims:ime, jms:jme ) , &
174 intent(inout), optional :: wstar, &
175 delta
176 real, dimension( ims:ime, jms:jme ) , &
177 intent(inout), optional :: regime
178 !
179 real, dimension( ims:ime, jms:jme ) , &
180 intent(in ), optional :: ctopo, &
181 ctopo2
182 !
183 real, dimension( kms:kme ) , &
184 intent(in ), optional :: znu, &
185 znw
186 !
187 real, optional, intent(in ) :: p_top
188 !
189 ! local
190 !
191 integer :: i,j,k
192 real, dimension( its:ite, kts:kte*ndiff ) :: rqvbl2dt, &
193 qv2d
194 real, dimension( its:ite, kts:kte ) :: pdh
195 real, dimension( its:ite, kts:kte+1 ) :: pdhi
196 real, dimension( its:ite ) :: &
197 dusfc, &
198 dvsfc, &
199 dtsfc, &
200 dqsfc
201 !
202 qv2d(its:ite,:) = 0.0
203 !
204 do j = jts,jte
205 do k = kts,kte+1
206 do i = its,ite
207 if(k.le.kte)pdh(i,k) = p3d(i,k,j)
208 pdhi(i,k) = p3di(i,k,j)
209 enddo
210 enddo
211 do k = kts,kte
212 do i = its,ite
213 qv2d(i,k) = qv3d(i,k,j)
214 qv2d(i,k+kte) = qc3d(i,k,j)
215 if(present(rqiblten)) qv2d(i,k+kte+kte) = qi3d(i,k,j)
216 enddo
217 enddo
218 !
219 call shinhong2d(J=j,ux=u3d(ims,kms,j),vx=v3d(ims,kms,j) &
220 ,tx=t3d(ims,kms,j) &
221 ,qx=qv2d(its,kts) &
222 ,p2d=pdh(its,kts),p2di=pdhi(its,kts) &
223 ,pi2d=pi3d(ims,kms,j) &
224 ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) &
225 ,ttnp=rthblten(ims,kms,j),qtnp=rqvbl2dt(its,kts),ndiff=ndiff &
226 ,cp=cp,g=g,rovcp=rovcp,rd=rd,rovg=rovg &
227 ,xlv=xlv,rv=rv &
228 ,ep1=ep1,ep2=ep2,karman=karman &
229 ,dz8w2d=dz8w(ims,kms,j) &
230 ,psfcpa=psfc(ims,j),znt=znt(ims,j),ust=ust(ims,j) &
231 ,hpbl=hpbl(ims,j) &
232 ,regime=regime(ims,j),psim=psim(ims,j) &
233 ,psih=psih(ims,j),xland=xland(ims,j) &
234 ,hfx=hfx(ims,j),qfx=qfx(ims,j) &
235 ,wspd=wspd(ims,j),br=br(ims,j) &
236 ,dusfc=dusfc,dvsfc=dvsfc,dtsfc=dtsfc,dqsfc=dqsfc &
237 ,dt=dt,rcl=1.0,kpbl1d=kpbl2d(ims,j) &
238 ,exch_hx=exch_h(ims,kms,j) &
239 ,wstar=wstar(ims,j) &
240 ,delta=delta(ims,j) &
241 ,u10=u10(ims,j),v10=v10(ims,j) &
242 ,ctopo=ctopo(ims,j),ctopo2=ctopo2(ims,j) &
243 ,shinhong_tke_diag=shinhong_tke_diag &
244 ,tke=tke_pbl(ims,kms,j),el_pbl=el_pbl(ims,kms,j) &
245 ,corf=corf(ims,j) &
246 ,dx=dx,dy=dy &
247 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
248 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
249 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
250 !
251 do k = kts,kte
252 do i = its,ite
253 rthblten(i,k,j) = rthblten(i,k,j)/pi3d(i,k,j)
254 rqvblten(i,k,j) = rqvbl2dt(i,k)
255 rqcblten(i,k,j) = rqvbl2dt(i,k+kte)
256 if(present(rqiblten)) rqiblten(i,k,j) = rqvbl2dt(i,k+kte+kte)
257 enddo
258 enddo
259 !
260 enddo
261 !
262 end subroutine shinhong
263 !-------------------------------------------------------------------------------
264 !
265 !-------------------------------------------------------------------------------
266 subroutine shinhong2d(j,ux,vx,tx,qx,p2d,p2di,pi2d, &
267 utnp,vtnp,ttnp,qtnp,ndiff, &
268 cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, &
269 dz8w2d,psfcpa, &
270 znt,ust,hpbl,psim,psih, &
271 xland,hfx,qfx,wspd,br, &
272 dusfc,dvsfc,dtsfc,dqsfc, &
273 dt,rcl,kpbl1d, &
274 exch_hx, &
275 wstar,delta, &
276 shinhong_tke_diag, &
277 tke,el_pbl,corf, &
278 u10,v10, &
279 ctopo,ctopo2, &
280 dx,dy, &
281 ids,ide, jds,jde, kds,kde, &
282 ims,ime, jms,jme, kms,kme, &
283 its,ite, jts,jte, kts,kte, &
284 !optional
285 regime )
286 !-------------------------------------------------------------------------------
287 implicit none
288 !-------------------------------------------------------------------------------
289 !
290 ! the shinhongpbl (shin and hong 2015) is based on the les study of shin
291 ! and hong (2013). the major ingredients of the shinhongpbl are
292 ! 1) the prescribed nonlocal heat transport profile fit to the les and
293 ! 2) inclusion of explicit scale dependency functions for vertical
294 ! transport in convective pbl.
295 ! so, the shinhongpbl works at the gray zone resolution of convective pbl.
296 ! note that honnert et al. (2011) first suggested explicit scale dependency
297 ! function, and shin and hong (2013) further classified the function by
298 ! stability (u*/w*) in convective pbl and calculated the function for
299 ! nonlocal and local transport separately.
300 ! vertical mixing in the stable boundary layer and free atmosphere follows
301 ! hong (2010) and hong et al. (2006), same as the ysupbl scheme.
302 !
303 ! shinhongpbl:
304 ! coded and implemented by hyeyum hailey shin (ncar)
305 ! summer 2014
306 !
307 ! ysupbl:
308 ! coded by song-you hong (yonsei university) and implemented by
309 ! song-you hong (yonsei university) and jimy dudhia (ncar)
310 ! summer 2002
311 !
312 ! references:
313 ! shin and hong (2015) mon. wea. rev.
314 ! shin and hong (2013) j. atmos. sci.
315 ! honnert, masson, and couvreux (2011) j. atmos. sci.
316 ! hong (2010) quart. j. roy. met. soc
317 ! hong, noh, and dudhia (2006), mon. wea. rev.
318 !
319 !-------------------------------------------------------------------------------
320 !
321 real,parameter :: xkzminm = 0.1,xkzminh = 0.01
322 real,parameter :: xkzmax = 1000.,rimin = -100.
323 real,parameter :: rlam = 30.,prmin = 0.25,prmax = 4.
324 real,parameter :: brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4
325 real,parameter :: afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0
326 real,parameter :: phifac = 8.,sfcfrac = 0.1
327 real,parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001
328 real,parameter :: h1 = 0.33333333, h2 = 0.6666667
329 real,parameter :: ckz = 0.001,zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16.
330 real,parameter :: tmin=1.e-2
331 real,parameter :: gamcrt = 3.,gamcrq = 2.e-3
332 real,parameter :: xka = 2.4e-5
333 integer,parameter :: imvdif = 1
334 !
335 ! tunable parameters for tke
336 !
337 real,parameter :: epsq2l = 0.01,c_1 = 1.0,gamcre = 0.224
338 !
339 ! tunable parameters for prescribed nonlocal transport profile
340 !
341 real,parameter :: mltop = 1.0,sfcfracn1 = 0.075
342 real,parameter :: nlfrac = 0.7,enlfrac = -0.4
343 real,parameter :: a11 = 1.0,a12 = -1.15
344 real,parameter :: ezfac = 1.5
345 real,parameter :: cpent = -0.4,rigsmax = 100.
346 real,parameter :: entfmin = 1.0, entfmax = 5.0
347 !
348 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
349 ims,ime, jms,jme, kms,kme, &
350 its,ite, jts,jte, kts,kte, &
351 j,ndiff
352 integer, intent(in ) :: shinhong_tke_diag
353 !
354 real, intent(in ) :: dt,rcl,cp,g,rovcp,rovg,rd,xlv,rv
355 real, intent(in ) :: ep1,ep2,karman
356 real, intent(in ) :: dx,dy
357 !
358 integer, dimension( ims:ime ) , &
359 intent(out ) :: kpbl1d
360 !
361 real, dimension( ims:ime, kms:kme ) , &
362 intent(in ) :: dz8w2d, &
363 pi2d
364 real, dimension( ims:ime, kms:kme ) , &
365 intent(in ) :: ux, &
366 vx
367 real, dimension( ims:ime, kms:kme ) , &
368 intent(in ) :: tx
369 !
370 real, dimension( its:ite, kts:kte*ndiff ) , &
371 intent(in ) :: qx
372 real, dimension( its:ite, kts:kte+1 ) , &
373 intent(in ) :: p2di
374 real, dimension( its:ite, kts:kte ) , &
375 intent(in ) :: p2d
376 !
377 real, dimension( ims:ime, kms:kme ) , &
378 intent(inout) :: utnp, &
379 vtnp, &
380 ttnp
381 real, dimension( ims:ime, kms:kme ) , &
382 intent(inout) :: exch_hx
383 real, dimension( ims:ime, kms:kme ) , &
384 intent(inout) :: tke, &
385 el_pbl
386 real, dimension( its:ite, kts:kte*ndiff ) , &
387 intent(inout) :: qtnp
388 !
389 real, dimension( ims:ime ) , &
390 intent(in ) :: xland, &
391 hfx, &
392 qfx
393 real, dimension( ims:ime ) , &
394 intent(in ) :: br, &
395 psim, &
396 psih, &
397 psfcpa
398 real, dimension( ims:ime ) , &
399 intent(in ) :: corf
400 !
401 real, dimension( ims:ime ) , &
402 intent(inout) :: ust, &
403 hpbl, &
404 znt
405 real, dimension( ims:ime ) , &
406 intent(inout) :: wspd
407 real, dimension( ims:ime ) , &
408 intent(inout) :: u10, &
409 v10
410 !
411 real, dimension( ims:ime ) , &
412 optional , &
413 intent(in ) :: ctopo, &
414 ctopo2
415 real, dimension( ims:ime ) , &
416 optional , &
417 intent(inout) :: regime
418 real, dimension( its:ite ) , &
419 intent(out ) :: wstar, &
420 delta
421 !
422 ! local vars
423 !
424 integer :: n,i,k,l,ic,is,nwmass
425 integer :: klpbl, kqc, kqi
426 integer :: lmh,lmxl
427 !
428 real :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0
429 real :: ss,ri,qmean,tmean,alpha,chi,zk,rl2,dk,sri
430 real :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz
431 real :: utend,vtend,ttend,qtend
432 real :: dtstep,govrthv
433 real :: cont, conq, conw, conwrc
434 real :: delxy,pu1,pth1,pq1
435 real :: dex,hgame_c
436 real :: zfacdx
437 real :: amf1,amf2,bmf2,amf3,bmf3,amf4,bmf4,sflux0,snlflux0
438 real :: mlfrac,ezfrac,sfcfracn
439 real :: uwst,uwstx,csfac
440 real :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, &
441 dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr, &
442 prfac,prfac2,phim8z
443 real :: cenlfrac
444 !
445 integer, dimension( its:ite ) :: kpbl
446 real, dimension( its:ite ) :: hol
447 real, dimension( its:ite ) :: deltaoh
448 real, dimension( its:ite ) :: rigs, &
449 enlfrac2, &
450 cslen
451 real, dimension( its:ite ) :: &
452 rhox, &
453 govrth, &
454 zl1,thermal, &
455 wscale, &
456 hgamt,hgamq, &
457 brdn,brup, &
458 phim,phih, &
459 dusfc,dvsfc, &
460 dtsfc,dqsfc, &
461 prpbl, &
462 wspd1
463 real, dimension( its:ite ) :: &
464 ust3, &
465 wstar3, &
466 hgamu,hgamv, &
467 wm2, we, &
468 bfxpbl, &
469 hfxpbl,qfxpbl, &
470 ufxpbl,vfxpbl, &
471 dthvx
472 real, dimension( its:ite ) :: &
473 brcr, &
474 sflux, &
475 zol1, &
476 brcr_sbro
477 real, dimension( its:ite ) :: &
478 efxpbl, &
479 hpbl_cbl, &
480 epshol, &
481 ct
482 !
483 real, dimension( its:ite, kts:kte ) :: xkzm,xkzh, &
484 f1,f2, &
485 r1,r2, &
486 ad,au, &
487 cu, &
488 al, &
489 xkzq, &
490 zfac
491 real, dimension( its:ite, kts:kte ) :: &
492 thx,thvx, &
493 del, &
494 dza, &
495 dzq, &
496 xkzom, &
497 xkzoh, &
498 za
499 real, dimension( its:ite, kts:kte ) :: &
500 wscalek
501 real, dimension( its:ite, kts:kte ) :: &
502 xkzml,xkzhl, &
503 zfacent,entfac
504 real, dimension( its:ite, kts:kte ) :: &
505 mf, &
506 zfacmf, &
507 entfacmf
508 real, dimension( its:ite, kts:kte ) :: &
509 q2x, &
510 hgame2d, &
511 tflux_e, &
512 qflux_e, &
513 tvflux_e
514 real, dimension( its:ite, kts:kte+1 ) :: zq
515 real, dimension( its:ite, kts:kte, ndiff ) :: r3,f3
516 !
517 real, dimension( kts:kte ) :: &
518 uxk,vxk, &
519 txk,thxk,thvxk, &
520 q2xk, &
521 hgame
522 real, dimension( kts:kte ) :: &
523 ps1d,pb1d,eps1d,pt1d, &
524 xkze1d,eflx_l1d,eflx_nl1d, &
525 ptke1
526 real, dimension( kts+1:kte ) :: &
527 s2,gh,rig,el, &
528 akmk,akhk, &
529 mfk,ufxpblk,vfxpblk,qfxpblk
530 real, dimension( kts:kte+1 ) :: zqk
531 real, dimension( kts:kte*ndiff ) :: qxk
532 !
533 logical, dimension( its:ite ) :: pblflg, &
534 sfcflg, &
535 stable
536 logical, dimension( ndiff ) :: ifvmix
537 !
538 !-------------------------------------------------------------------------------
539 !
540 klpbl = kte
541 lmh = 1
542 lmxl = 1
543 !
544 cont=cp/g
545 conq=xlv/g
546 conw=1./g
547 conwrc = conw*sqrt(rcl)
548 conpr = bfac*karman*sfcfrac
549 !
550 ! k-start index for cloud and rain
551 !
552 kqc = 1 + kte
553 kqi = 1 + kte*2
554 nwmass = 3
555 ifvmix(:) = .true.
556 !
557 do k = kts,kte
558 do i = its,ite
559 thx(i,k) = tx(i,k)/pi2d(i,k)
560 enddo
561 enddo
562 !
563 do k = kts,kte
564 do i = its,ite
565 tvcon = (1.+ep1*qx(i,k))
566 thvx(i,k) = thx(i,k)*tvcon
567 enddo
568 enddo
569 !
570 do i = its,ite
571 tvcon = (1.+ep1*qx(i,1))
572 rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon)
573 govrth(i) = g/thx(i,1)
574 enddo
575 !
576 !-----compute the height of full- and half-sigma levels above ground
577 ! level, and the layer thicknesses.
578 !
579 do i = its,ite
580 zq(i,1) = 0.
581 enddo
582 !
583 do k = kts,kte
584 do i = its,ite
585 zq(i,k+1) = dz8w2d(i,k)+zq(i,k)
586 enddo
587 enddo
588 !
589 do k = kts,kte
590 do i = its,ite
591 za(i,k) = 0.5*(zq(i,k)+zq(i,k+1))
592 dzq(i,k) = zq(i,k+1)-zq(i,k)
593 del(i,k) = p2di(i,k)-p2di(i,k+1)
594 enddo
595 enddo
596 !
597 do i = its,ite
598 dza(i,1) = za(i,1)
599 enddo
600 !
601 do k = kts+1,kte
602 do i = its,ite
603 dza(i,k) = za(i,k)-za(i,k-1)
604 enddo
605 enddo
606 !
607 !
608 !-----initialize vertical tendencies and
609 !
610 utnp(its:ite,:) = 0.
611 vtnp(its:ite,:) = 0.
612 ttnp(its:ite,:) = 0.
613 qtnp(its:ite,:) = 0.
614 !
615 do i = its,ite
616 wspd1(i) = sqrt(ux(i,1)*ux(i,1)+vx(i,1)*vx(i,1))+1.e-9
617 enddo
618 !
619 !---- compute vertical diffusion
620 !
621 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
622 ! compute preliminary variables
623 !
624 dtstep = dt
625 dt2 = 2.*dtstep
626 rdt = 1./dt2
627 !
628 do i = its,ite
629 bfxpbl(i) = 0.0
630 hfxpbl(i) = 0.0
631 qfxpbl(i) = 0.0
632 ufxpbl(i) = 0.0
633 vfxpbl(i) = 0.0
634 hgamu(i) = 0.0
635 hgamv(i) = 0.0
636 delta(i) = 0.0
637 enddo
638 !
639 do i = its,ite
640 efxpbl(i) = 0.0
641 hpbl_cbl(i) = 0.0
642 epshol(i) = 0.0
643 ct(i) = 0.0
644 enddo
645 !
646 do i = its,ite
647 deltaoh(i) = 0.0
648 rigs(i) = 0.0
649 enlfrac2(i) = 0.0
650 cslen(i) = 0.0
651 enddo
652 !
653 do k = kts,klpbl
654 do i = its,ite
655 wscalek(i,k) = 0.0
656 enddo
657 enddo
658 !
659 do k = kts,klpbl
660 do i = its,ite
661 zfac(i,k) = 0.0
662 enddo
663 enddo
664 !
665 do k = kts,kte
666 do i = its,ite
667 q2x(i,k) = 2.*tke(i,k)
668 enddo
669 enddo
670 !
671 do k = kts,kte
672 do i = its,ite
673 el_pbl(i,k) = 0.0
674 hgame2d(i,k) = 0.0
675 tflux_e(i,k) = 0.0
676 qflux_e(i,k) = 0.0
677 tvflux_e(i,k) = 0.0
678 enddo
679 enddo
680 !
681 do k = kts,kte
682 do i = its,ite
683 mf(i,k) = 0.0
684 zfacmf(i,k) = 0.0
685 enddo
686 enddo
687 !
688 do k = kts,klpbl-1
689 do i = its,ite
690 xkzom(i,k) = xkzminm
691 xkzoh(i,k) = xkzminh
692 enddo
693 enddo
694 !
695 do i = its,ite
696 dusfc(i) = 0.
697 dvsfc(i) = 0.
698 dtsfc(i) = 0.
699 dqsfc(i) = 0.
700 enddo
701 !
702 do i = its,ite
703 hgamt(i) = 0.
704 hgamq(i) = 0.
705 wscale(i) = 0.
706 kpbl(i) = 1
707 hpbl(i) = zq(i,1)
708 hpbl_cbl(i) = zq(i,1)
709 zl1(i) = za(i,1)
710 thermal(i)= thvx(i,1)
711 pblflg(i) = .true.
712 sfcflg(i) = .true.
713 sflux(i) = hfx(i)/rhox(i)/cp + qfx(i)/rhox(i)*ep1*thx(i,1)
714 if(br(i).gt.0.0) sfcflg(i) = .false.
715 enddo
716 !
717 ! compute the first guess of pbl height
718 !
719 do i = its,ite
720 stable(i) = .false.
721 brup(i) = br(i)
722 brcr(i) = brcr_ub
723 enddo
724 !
725 do k = 2,klpbl
726 do i = its,ite
727 if(.not.stable(i))then
728 brdn(i) = brup(i)
729 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
730 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
731 kpbl(i) = k
732 stable(i) = brup(i).gt.brcr(i)
733 endif
734 enddo
735 enddo
736 !
737 do i = its,ite
738 k = kpbl(i)
739 if(brdn(i).ge.brcr(i))then
740 brint = 0.
741 elseif(brup(i).le.brcr(i))then
742 brint = 1.
743 else
744 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
745 endif
746 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
747 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
748 if(kpbl(i).le.1) pblflg(i) = .false.
749 enddo
750 !
751 do i = its,ite
752 fm = psim(i)
753 fh = psih(i)
754 zol1(i) = max(br(i)*fm*fm/fh,rimin)
755 if(sfcflg(i))then
756 zol1(i) = min(zol1(i),-zfmin)
757 else
758 zol1(i) = max(zol1(i),zfmin)
759 endif
760 hol1 = zol1(i)*hpbl(i)/zl1(i)*sfcfrac
761 epshol(i) = hol1
762 if(sfcflg(i))then
763 phim(i) = (1.-aphi16*hol1)**(-1./4.)
764 phih(i) = (1.-aphi16*hol1)**(-1./2.)
765 bfx0 = max(sflux(i),0.)
766 hfx0 = max(hfx(i)/rhox(i)/cp,0.)
767 qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.)
768 wstar3(i) = (govrth(i)*bfx0*hpbl(i))
769 wstar(i) = (wstar3(i))**h1
770 else
771 phim(i) = (1.+aphi5*hol1)
772 phih(i) = phim(i)
773 wstar(i) = 0.
774 wstar3(i) = 0.
775 endif
776 ust3(i) = ust(i)**3.
777 wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1
778 wscale(i) = min(wscale(i),ust(i)*aphi16)
779 wscale(i) = max(wscale(i),ust(i)/aphi5)
780 enddo
781 !
782 ! compute the surface variables for pbl height estimation
783 ! under unstable conditions
784 !
785 do i = its,ite
786 if(sfcflg(i).and.sflux(i).gt.0.0)then
787 gamfac = bfac/rhox(i)/wscale(i)
788 hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
789 hgamq(i) = min(gamfac*qfx(i),gamcrq)
790 vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac
791 thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
792 hgamt(i) = max(hgamt(i),0.0)
793 hgamq(i) = max(hgamq(i),0.0)
794 brint = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.)
795 hgamu(i) = brint*ux(i,1)
796 hgamv(i) = brint*vx(i,1)
797 else
798 pblflg(i) = .false.
799 endif
800 enddo
801 !
802 ! enhance the pbl height by considering the thermal
803 !
804 do i = its,ite
805 if(pblflg(i))then
806 kpbl(i) = 1
807 hpbl(i) = zq(i,1)
808 endif
809 enddo
810 !
811 do i = its,ite
812 if(pblflg(i))then
813 stable(i) = .false.
814 brup(i) = br(i)
815 brcr(i) = brcr_ub
816 endif
817 enddo
818 !
819 do k = 2,klpbl
820 do i = its,ite
821 if(.not.stable(i).and.pblflg(i))then
822 brdn(i) = brup(i)
823 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
824 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
825 kpbl(i) = k
826 stable(i) = brup(i).gt.brcr(i)
827 endif
828 enddo
829 enddo
830 !
831 do i = its,ite
832 if(pblflg(i)) then
833 k = kpbl(i)
834 if(brdn(i).ge.brcr(i))then
835 brint = 0.
836 elseif(brup(i).le.brcr(i))then
837 brint = 1.
838 else
839 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
840 endif
841 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
842 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
843 if(kpbl(i).le.1) pblflg(i) = .false.
844 if (wstar(i) .ne. 0) then
845 uwst = abs(ust(i)/wstar(i)-0.5)
846 uwstx = -80.*uwst+14.
847 csfac = 0.5*(tanh(uwstx)+3.)
848 else
849 csfac = 1
850 endif
851 cslen(i) = csfac*hpbl(i)
852 endif
853 enddo
854 !
855 ! stable boundary layer
856 !
857 do i = its,ite
858 hpbl_cbl(i) = hpbl(i)
859 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
860 brup(i) = br(i)
861 stable(i) = .false.
862 else
863 stable(i) = .true.
864 endif
865 enddo
866 !
867 do i = its,ite
868 if((.not.stable(i)).and.((xland(i)-1.5).ge.0))then
869 wspd10 = u10(i)*u10(i) + v10(i)*v10(i)
870 wspd10 = sqrt(wspd10)
871 ross = wspd10 / (cori*znt(i))
872 brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3)
873 endif
874 enddo
875 !
876 do i = its,ite
877 if(.not.stable(i))then
878 if((xland(i)-1.5).ge.0)then
879 brcr(i) = brcr_sbro(i)
880 else
881 brcr(i) = brcr_sb
882 endif
883 endif
884 enddo
885 !
886 do k = 2,klpbl
887 do i = its,ite
888 if(.not.stable(i))then
889 brdn(i) = brup(i)
890 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
891 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
892 kpbl(i) = k
893 stable(i) = brup(i).gt.brcr(i)
894 endif
895 enddo
896 enddo
897 !
898 do i = its,ite
899 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
900 k = kpbl(i)
901 if(brdn(i).ge.brcr(i))then
902 brint = 0.
903 elseif(brup(i).le.brcr(i))then
904 brint = 1.
905 else
906 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
907 endif
908 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
909 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
910 if(kpbl(i).le.1) pblflg(i) = .false.
911 endif
912 enddo
913 !
914 ! scale dependency for nonlocal momentum and moisture transport
915 !
916 delxy=sqrt(dx*dy)
917 !
918 do i = its,ite
919 pu1=pu(delxy,cslen(i))
920 pq1=pq(delxy,cslen(i))
921 if(pblflg(i)) then
922 hgamu(i) = hgamu(i)*pu1
923 hgamv(i) = hgamv(i)*pu1
924 hgamq(i) = hgamq(i)*pq1
925 endif
926 enddo
927 !
928 ! estimate the entrainment parameters
929 !
930 delxy=sqrt(dx*dy)
931 !
932 do i = its,ite
933 if(pblflg(i)) then
934 k = kpbl(i) - 1
935 prpbl(i) = 1.0
936 wm3 = wstar3(i) + 5. * ust3(i)
937 wm2(i) = wm3**h2
938 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
939 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin)
940 dthx = max(thx(i,k+1)-thx(i,k),tmin)
941 dqx = min(qx(i,k+1)-qx(i,k),0.0)
942 we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
943 hfxpbl(i) = we(i)*dthx
944 pq1=pq(delxy,cslen(i))
945 qfxpbl(i) = we(i)*dqx*pq1
946 !
947 pu1=pu(delxy,cslen(i))
948 dux = ux(i,k+1)-ux(i,k)
949 dvx = vx(i,k+1)-vx(i,k)
950 if(dux.gt.tmin) then
951 ufxpbl(i) = max(prpbl(i)*we(i)*dux*pu1,-ust(i)*ust(i))
952 elseif(dux.lt.-tmin) then
953 ufxpbl(i) = min(prpbl(i)*we(i)*dux*pu1,ust(i)*ust(i))
954 else
955 ufxpbl(i) = 0.0
956 endif
957 if(dvx.gt.tmin) then
958 vfxpbl(i) = max(prpbl(i)*we(i)*dvx*pu1,-ust(i)*ust(i))
959 elseif(dvx.lt.-tmin) then
960 vfxpbl(i) = min(prpbl(i)*we(i)*dvx*pu1,ust(i)*ust(i))
961 else
962 vfxpbl(i) = 0.0
963 endif
964 delb = govrth(i)*d3*hpbl(i)
965 delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.)
966 delb = govrth(i)*dthvx(i)
967 deltaoh(i) = d1*hpbl(i) + d2*wm2(i)/delb
968 deltaoh(i) = max(ezfac*deltaoh(i),hpbl(i)-za(i,kpbl(i)-1)-1.)
969 deltaoh(i) = min(deltaoh(i) ,hpbl(i))
970 if ((dux .ne. 0) .or. (dvx .ne. 0)) then
971 rigs(i) = govrth(i)*dthvx(i)*deltaoh(i)/(dux**2.+dvx**2.)
972 else
973 rigs(i) = rigsmax
974 endif
975 rigs(i) = max(min(rigs(i), rigsmax),rimin)
976 if((rigs(i).gt.0) .and. (abs(rigs(i)+cpent) .le. 1.e-6))then
977 cenlfrac = entfmax
978 else
979 cenlfrac = rigs(i)/(rigs(i)+cpent)
980 endif
981 cenlfrac = min(cenlfrac,entfmax)
982 enlfrac2(i) = max(wm3/wstar3(i)*cenlfrac, entfmin)
983 enlfrac2(i) = enlfrac2(i)*enlfrac
984 endif
985 enddo
986 !
987 do k = kts,klpbl
988 do i = its,ite
989 if(pblflg(i))then
990 entfacmf(i,k) = sqrt(((zq(i,k+1)-hpbl(i))/deltaoh(i))**2.)
991 endif
992 if(pblflg(i).and.k.ge.kpbl(i))then
993 entfac(i,k) = ((zq(i,k+1)-hpbl(i))/deltaoh(i))**2.
994 else
995 entfac(i,k) = 1.e30
996 endif
997 enddo
998 enddo
999 !
1000 ! compute diffusion coefficients below pbl
1001 !
1002 do k = kts,klpbl
1003 do i = its,ite
1004 if(k.lt.kpbl(i)) then
1005 zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.)
1006 zfacent(i,k) = (1.-zfac(i,k))**3.
1007 wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
1008 if(sfcflg(i)) then
1009 prfac = conpr
1010 prfac2 = 15.9*wstar3(i)/ust3(i)/(1.+4.*karman*wstar3(i)/ust3(i))
1011 prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2.
1012 else
1013 prfac = 0.
1014 prfac2 = 0.
1015 prnumfac = 0.
1016 phim8z = 1.+aphi5*zol1(i)*zq(i,k+1)/zl1(i)
1017 wscalek(i,k) = ust(i)/phim8z
1018 wscalek(i,k) = max(wscalek(i,k),0.001)
1019 endif
1020 prnum0 = (phih(i)/phim(i)+prfac)
1021 prnum0 = max(min(prnum0,prmax),prmin)
1022 xkzm(i,k) = wscalek(i,k)*karman*zq(i,k+1)*zfac(i,k)**pfac
1023 prnum = 1. + (prnum0-1.)*exp(prnumfac)
1024 xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac)
1025 prnum0 = prnum0/(1.+prfac2*karman*sfcfrac)
1026 prnum = 1. + (prnum0-1.)*exp(prnumfac)
1027 xkzh(i,k) = xkzm(i,k)/prnum
1028 xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
1029 xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
1030 xkzq(i,k) = xkzq(i,k)+xkzoh(i,k)
1031 xkzm(i,k) = min(xkzm(i,k),xkzmax)
1032 xkzh(i,k) = min(xkzh(i,k),xkzmax)
1033 xkzq(i,k) = min(xkzq(i,k),xkzmax)
1034 endif
1035 enddo
1036 enddo
1037 !
1038 ! compute diffusion coefficients over pbl (free atmosphere)
1039 !
1040 do k = kts,kte-1
1041 do i = its,ite
1042 if(k.ge.kpbl(i)) then
1043 ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k)) &
1044 +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k))) &
1045 /(dza(i,k+1)*dza(i,k+1))+1.e-9
1046 govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
1047 ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1))
1048 ! in cloud
1049 if(imvdif.eq.1.and.nwmass.ge.3)then
1050 if((qx(i,kqc+k-1)+qx(i,kqi+k-1)).gt.0.01e-3 &
1051 .and.(qx(i,kqc+k)+qx(i,kqi+k)).gt.0.01e-3) then
1052 qmean = 0.5*(qx(i,k)+qx(i,k+1))
1053 tmean = 0.5*(tx(i,k)+tx(i,k+1))
1054 alpha = xlv*qmean/rd/tmean
1055 chi = xlv*xlv*qmean/cp/rv/tmean/tmean
1056 ri = (1.+alpha)*(ri-g*g/ss/tmean/cp*((chi-alpha)/(1.+chi)))
1057 endif
1058 endif
1059 zk = karman*zq(i,k+1)
1060 rlamdz = min(max(0.1*dza(i,k+1),rlam),300.)
1061 rlamdz = min(dza(i,k+1),rlamdz)
1062 rl2 = (zk*rlamdz/(rlamdz+zk))**2
1063 dk = rl2*sqrt(ss)
1064 if(ri.lt.0.)then
1065 ! unstable regime
1066 ri = max(ri, rimin)
1067 sri = sqrt(-ri)
1068 xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
1069 xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
1070 else
1071 ! stable regime
1072 xkzh(i,k) = dk/(1+5.*ri)**2
1073 prnum = 1.0+2.1*ri
1074 prnum = min(prnum,prmax)
1075 xkzm(i,k) = xkzh(i,k)*prnum
1076 endif
1077 !
1078 xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
1079 xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
1080 xkzm(i,k) = min(xkzm(i,k),xkzmax)
1081 xkzh(i,k) = min(xkzh(i,k),xkzmax)
1082 xkzml(i,k) = xkzm(i,k)
1083 xkzhl(i,k) = xkzh(i,k)
1084 endif
1085 enddo
1086 enddo
1087 !
1088 ! prescribe nonlocal heat transport below pbl
1089 !
1090 do i = its,ite
1091 deltaoh(i) = deltaoh(i)/hpbl(i)
1092 enddo
1093 !
1094 delxy=sqrt(dx*dy)
1095 do i = its,ite
1096 mlfrac = mltop-deltaoh(i)
1097 ezfrac = mltop+deltaoh(i)
1098 zfacmf(i,1) = min(max((zq(i,2)/hpbl(i)),zfmin),1.)
1099 sfcfracn = max(sfcfracn1,zfacmf(i,1))
1100 !
1101 sflux0 = (a11+a12*sfcfracn)*sflux(i)
1102 snlflux0 = nlfrac*sflux0
1103 amf1 = snlflux0/sfcfracn
1104 if (pblflg(i)) then
1105 amf2 = -snlflux0/(mlfrac-sfcfracn)
1106 bmf2 = -mlfrac*amf2
1107 endif
1108 if ((deltaoh(i) .eq. 0) .and. (enlfrac2(i) .eq. 0)) then
1109 amf3 = 0.
1110 else
1111 amf3 = snlflux0*enlfrac2(i)/deltaoh(i)
1112 endif
1113 bmf3 = -amf3*mlfrac
1114 hfxpbl(i) = amf3+bmf3
1115 pth1=pthnl(delxy,cslen(i))
1116 hfxpbl(i) = hfxpbl(i)*pth1
1117 !
1118 do k = kts,klpbl
1119 zfacmf(i,k) = max((zq(i,k+1)/hpbl(i)),zfmin)
1120 if(pblflg(i).and.k.lt.kpbl(i)) then
1121 if(zfacmf(i,k).le.sfcfracn) then
1122 mf(i,k) = amf1*zfacmf(i,k)
1123 else if (zfacmf(i,k).le.mlfrac) then
1124 mf(i,k) = amf2*zfacmf(i,k)+bmf2
1125 endif
1126 mf(i,k) = mf(i,k)+hfxpbl(i)*exp(-entfacmf(i,k))
1127 mf(i,k) = mf(i,k)*pth1
1128 endif
1129 enddo
1130 enddo
1131 !
1132 ! compute tridiagonal matrix elements for heat
1133 !
1134 do k = kts,kte
1135 do i = its,ite
1136 au(i,k) = 0.
1137 al(i,k) = 0.
1138 ad(i,k) = 0.
1139 f1(i,k) = 0.
1140 enddo
1141 enddo
1142 !
1143 do i = its,ite
1144 ad(i,1) = 1.
1145 f1(i,1) = thx(i,1)-300.+hfx(i)/cont/del(i,1)*dt2
1146 enddo
1147 !
1148 delxy=sqrt(dx*dy)
1149 do k = kts,kte-1
1150 do i = its,ite
1151 dtodsd = dt2/del(i,k)
1152 dtodsu = dt2/del(i,k+1)
1153 dsig = p2d(i,k)-p2d(i,k+1)
1154 rdz = 1./dza(i,k+1)
1155 tem1 = dsig*xkzh(i,k)*rdz
1156 if(pblflg(i).and.k.lt.kpbl(i)) then
1157 dsdzt = tem1*(-mf(i,k)/xkzh(i,k))
1158 f1(i,k) = f1(i,k)+dtodsd*dsdzt
1159 f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt
1160 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1161 xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
1162 xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
1163 xkzh(i,k) = max(xkzh(i,k),xkzoh(i,k))
1164 xkzh(i,k) = min(xkzh(i,k),xkzmax)
1165 f1(i,k+1) = thx(i,k+1)-300.
1166 else
1167 f1(i,k+1) = thx(i,k+1)-300.
1168 endif
1169 tem1 = dsig*xkzh(i,k)*rdz
1170 dsdz2 = tem1*rdz
1171 au(i,k) = -dtodsd*dsdz2
1172 al(i,k) = -dtodsu*dsdz2
1173 !
1174 ! scale dependency for local heat transport
1175 !
1176 zfacdx=0.2*hpbl(i)/zq(i,k+1)
1177 delxy=sqrt(dx*dy)*max(zfacdx,1.0)
1178 pth1=pthl(delxy,hpbl(i))
1179 if(pblflg(i).and.k.lt.kpbl(i)) then
1180 au(i,k) = au(i,k)*pth1
1181 al(i,k) = al(i,k)*pth1
1182 endif
1183 ad(i,k) = ad(i,k)-au(i,k)
1184 ad(i,k+1) = 1.-al(i,k)
1185 exch_hx(i,k+1) = xkzh(i,k)
1186 enddo
1187 enddo
1188 !
1189 ! copies here to avoid duplicate input args for tridin
1190 !
1191 do k = kts,kte
1192 do i = its,ite
1193 cu(i,k) = au(i,k)
1194 r1(i,k) = f1(i,k)
1195 enddo
1196 enddo
1197 !
1198 call tridin_ysu(al,ad,cu,r1,au,f1,its,ite,kts,kte,1)
1199 !
1200 ! recover tendencies of heat
1201 !
1202 do k = kte,kts,-1
1203 do i = its,ite
1204 ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k)
1205 ttnp(i,k) = ttnp(i,k)+ttend
1206 dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)/pi2d(i,k)
1207 if(k.eq.kte) then
1208 tflux_e(i,k) = ttend*dz8w2d(i,k)
1209 else
1210 tflux_e(i,k) = tflux_e(i,k+1) + ttend*dz8w2d(i,k)
1211 endif
1212 enddo
1213 enddo
1214 !
1215 ! compute tridiagonal matrix elements for moisture, clouds, and gases
1216 !
1217 do k = kts,kte
1218 do i = its,ite
1219 au(i,k) = 0.
1220 al(i,k) = 0.
1221 ad(i,k) = 0.
1222 enddo
1223 enddo
1224 !
1225 do ic = 1,ndiff
1226 do i = its,ite
1227 do k = kts,kte
1228 f3(i,k,ic) = 0.
1229 enddo
1230 enddo
1231 enddo
1232 !
1233 do i = its,ite
1234 ad(i,1) = 1.
1235 f3(i,1,1) = qx(i,1)+qfx(i)*g/del(i,1)*dt2
1236 enddo
1237 !
1238 if(ndiff.ge.2) then
1239 do ic = 2,ndiff
1240 is = (ic-1) * kte
1241 do i = its,ite
1242 f3(i,1,ic) = qx(i,1+is)
1243 enddo
1244 enddo
1245 endif
1246 !
1247 do k = kts,kte-1
1248 do i = its,ite
1249 if(k.ge.kpbl(i)) then
1250 xkzq(i,k) = xkzh(i,k)
1251 endif
1252 enddo
1253 enddo
1254 !
1255 do k = kts,kte-1
1256 do i = its,ite
1257 dtodsd = dt2/del(i,k)
1258 dtodsu = dt2/del(i,k+1)
1259 dsig = p2d(i,k)-p2d(i,k+1)
1260 rdz = 1./dza(i,k+1)
1261 tem1 = dsig*xkzq(i,k)*rdz
1262 if(pblflg(i).and.k.lt.kpbl(i)) then
1263 dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzq(i,k))
1264 f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq
1265 f3(i,k+1,1) = qx(i,k+1)-dtodsu*dsdzq
1266 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1267 xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
1268 xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k))
1269 xkzq(i,k) = max(xkzq(i,k),xkzoh(i,k))
1270 xkzq(i,k) = min(xkzq(i,k),xkzmax)
1271 f3(i,k+1,1) = qx(i,k+1)
1272 else
1273 f3(i,k+1,1) = qx(i,k+1)
1274 endif
1275 tem1 = dsig*xkzq(i,k)*rdz
1276 dsdz2 = tem1*rdz
1277 au(i,k) = -dtodsd*dsdz2
1278 al(i,k) = -dtodsu*dsdz2
1279 !
1280 ! scale dependency for local moisture transport
1281 !
1282 zfacdx=0.2*hpbl(i)/zq(i,k+1)
1283 delxy=sqrt(dx*dy)*max(zfacdx,1.0)
1284 pq1=pq(delxy,hpbl(i))
1285 if(pblflg(i).and.k.lt.kpbl(i)) then
1286 au(i,k) = au(i,k)*pq1
1287 al(i,k) = al(i,k)*pq1
1288 endif
1289 ad(i,k) = ad(i,k)-au(i,k)
1290 ad(i,k+1) = 1.-al(i,k)
1291 ! exch_hx(i,k+1) = xkzh(i,k)
1292 enddo
1293 enddo
1294 !
1295 if(ndiff.ge.2) then
1296 do ic = 2,ndiff
1297 is = (ic-1) * kte
1298 do k = kts,kte-1
1299 do i = its,ite
1300 f3(i,k+1,ic) = qx(i,k+1+is)
1301 enddo
1302 enddo
1303 enddo
1304 endif
1305 !
1306 ! copies here to avoid duplicate input args for tridin
1307 !
1308 do k = kts,kte
1309 do i = its,ite
1310 cu(i,k) = au(i,k)
1311 enddo
1312 enddo
1313 !
1314 do ic = 1,ndiff
1315 do k = kts,kte
1316 do i = its,ite
1317 r3(i,k,ic) = f3(i,k,ic)
1318 enddo
1319 enddo
1320 enddo
1321 !
1322 ! solve tridiagonal problem for moisture, clouds, and gases
1323 !
1324 call tridin_ysu(al,ad,cu,r3,au,f3,its,ite,kts,kte,ndiff)
1325 !
1326 ! recover tendencies of heat and moisture
1327 !
1328 do k = kte,kts,-1
1329 do i = its,ite
1330 qtend = (f3(i,k,1)-qx(i,k))*rdt
1331 qtnp(i,k) = qtnp(i,k)+qtend
1332 dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k)
1333 if(k.eq.kte) then
1334 qflux_e(i,k) = qtend*dz8w2d(i,k)
1335 else
1336 qflux_e(i,k) = qflux_e(i,k+1) + qtend*dz8w2d(i,k)
1337 endif
1338 tvflux_e(i,k) = tflux_e(i,k) + qflux_e(i,k)*ep1*thx(i,k)
1339 enddo
1340 enddo
1341 !
1342 do k = kts,kte
1343 do i = its,ite
1344 if(pblflg(i).and.k.lt.kpbl(i)) then
1345 hgame_c=c_1*0.2*2.5*(g/thvx(i,k))*wstar(i)/(0.25*(q2x(i,k+1)+q2x(i,k)))
1346 hgame_c=min(hgame_c,gamcre)
1347 if(k.eq.kte)then
1348 hgame2d(i,k)=hgame_c*0.5*tvflux_e(i,k)*hpbl(i)
1349 hgame2d(i,k)=max(hgame2d(i,k),0.0)
1350 else
1351 hgame2d(i,k)=hgame_c*0.5*(tvflux_e(i,k)+tvflux_e(i,k+1))*hpbl(i)
1352 hgame2d(i,k)=max(hgame2d(i,k),0.0)
1353 endif
1354 endif
1355 enddo
1356 enddo
1357 !
1358 if(ndiff.ge.2) then
1359 do ic = 2,ndiff
1360 if(ifvmix(ic)) then
1361 is = (ic-1) * kte
1362 do k = kte,kts,-1
1363 do i = its,ite
1364 qtend = (f3(i,k,ic)-qx(i,k+is))*rdt
1365 qtnp(i,k+is) = qtnp(i,k+is)+qtend
1366 enddo
1367 enddo
1368 endif
1369 enddo
1370 endif
1371 !
1372 ! compute tridiagonal matrix elements for momentum
1373 !
1374 do i = its,ite
1375 do k = kts,kte
1376 au(i,k) = 0.
1377 al(i,k) = 0.
1378 ad(i,k) = 0.
1379 f1(i,k) = 0.
1380 f2(i,k) = 0.
1381 enddo
1382 enddo
1383 !
1384 do i = its,ite
1385 ! paj: ctopo=1 if topo_wind=0 (default)
1386 ! mchen add this line to make sure NMM can still work with YSU PBL
1387 ad(i,1) = 1.+ctopo(i)*ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 &
1388 *(wspd1(i)/wspd(i))**2
1389 f1(i,1) = ux(i,1)
1390 f2(i,1) = vx(i,1)
1391 enddo
1392 !
1393 delxy=sqrt(dx*dy)
1394 do k = kts,kte-1
1395 do i = its,ite
1396 dtodsd = dt2/del(i,k)
1397 dtodsu = dt2/del(i,k+1)
1398 dsig = p2d(i,k)-p2d(i,k+1)
1399 rdz = 1./dza(i,k+1)
1400 tem1 = dsig*xkzm(i,k)*rdz
1401 if(pblflg(i).and.k.lt.kpbl(i))then
1402 dsdzu = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k))
1403 dsdzv = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k))
1404 f1(i,k) = f1(i,k)+dtodsd*dsdzu
1405 f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu
1406 f2(i,k) = f2(i,k)+dtodsd*dsdzv
1407 f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv
1408 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1409 xkzm(i,k) = prpbl(i)*xkzh(i,k)
1410 xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
1411 xkzm(i,k) = max(xkzm(i,k),xkzom(i,k))
1412 xkzm(i,k) = min(xkzm(i,k),xkzmax)
1413 f1(i,k+1) = ux(i,k+1)
1414 f2(i,k+1) = vx(i,k+1)
1415 else
1416 f1(i,k+1) = ux(i,k+1)
1417 f2(i,k+1) = vx(i,k+1)
1418 endif
1419 tem1 = dsig*xkzm(i,k)*rdz
1420 dsdz2 = tem1*rdz
1421 au(i,k) = -dtodsd*dsdz2
1422 al(i,k) = -dtodsu*dsdz2
1423 !
1424 ! scale dependency for local momentum transport
1425 !
1426 zfacdx=0.2*hpbl(i)/zq(i,k+1)
1427 delxy=sqrt(dx*dy)*max(zfacdx,1.0)
1428 pu1=pu(delxy,hpbl(i))
1429 if(pblflg(i).and.k.lt.kpbl(i)) then
1430 au(i,k) = au(i,k)*pu1
1431 al(i,k) = al(i,k)*pu1
1432 endif
1433 ad(i,k) = ad(i,k)-au(i,k)
1434 ad(i,k+1) = 1.-al(i,k)
1435 enddo
1436 enddo
1437 !
1438 ! copies here to avoid duplicate input args for tridin
1439 !
1440 do k = kts,kte
1441 do i = its,ite
1442 cu(i,k) = au(i,k)
1443 r1(i,k) = f1(i,k)
1444 r2(i,k) = f2(i,k)
1445 enddo
1446 enddo
1447 !
1448 ! solve tridiagonal problem for momentum
1449 !
1450 call tridi1n(al,ad,cu,r1,r2,au,f1,f2,its,ite,kts,kte,1)
1451 !
1452 ! recover tendencies of momentum
1453 !
1454 do k = kte,kts,-1
1455 do i = its,ite
1456 utend = (f1(i,k)-ux(i,k))*rdt
1457 vtend = (f2(i,k)-vx(i,k))*rdt
1458 utnp(i,k) = utnp(i,k)+utend
1459 vtnp(i,k) = vtnp(i,k)+vtend
1460 dusfc(i) = dusfc(i) + utend*conwrc*del(i,k)
1461 dvsfc(i) = dvsfc(i) + vtend*conwrc*del(i,k)
1462 enddo
1463 enddo
1464 !
1465 do i = its,ite
1466 kpbl1d(i) = kpbl(i)
1467 enddo
1468 !
1469 ! paj: ctopo2=1 if topo_wind=0 (default)
1470 !
1471 do i = its,ite
1472 u10(i) = ctopo2(i)*u10(i)+(1-ctopo2(i))*ux(i,1)
1473 v10(i) = ctopo2(i)*v10(i)+(1-ctopo2(i))*vx(i,1)
1474 enddo
1475 !
1476 !---- calculate sgs tke which is consistent with shinhongpbl algorithm
1477 !
1478 if (shinhong_tke_diag.eq.1) then
1479 !
1480 tke_calculation: do i = its,ite
1481 do k = kts+1,kte
1482 s2(k) = 0.0
1483 gh(k) = 0.0
1484 rig(k) = 0.0
1485 el(k) = 0.0
1486 akmk(k) = 0.0
1487 akhk(k) = 0.0
1488 mfk(k) = 0.0
1489 ufxpblk(k) = 0.0
1490 vfxpblk(k) = 0.0
1491 qfxpblk(k) = 0.0
1492 enddo
1493 !
1494 do k = kts,kte
1495 uxk(k) = 0.0
1496 vxk(k) = 0.0
1497 txk(k) = 0.0
1498 thxk(k) = 0.0
1499 thvxk(k) = 0.0
1500 q2xk(k) = 0.0
1501 hgame(k) = 0.0
1502 ps1d(k) = 0.0
1503 pb1d(k) = 0.0
1504 eps1d(k) = 0.0
1505 pt1d(k) = 0.0
1506 xkze1d(k) = 0.0
1507 eflx_l1d(k) = 0.0
1508 eflx_nl1d(k) = 0.0
1509 ptke1(k) = 1.0
1510 enddo
1511 !
1512 do k = kts,kte+1
1513 zqk(k) = 0.0
1514 enddo
1515 !
1516 do k = kts,kte*ndiff
1517 qxk(k) = 0.0
1518 enddo
1519 !
1520 do k = kts,kte
1521 uxk(k) = ux(i,k)
1522 vxk(k) = vx(i,k)
1523 txk(k) = tx(i,k)
1524 thxk(k) = thx(i,k)
1525 thvxk(k) = thvx(i,k)
1526 q2xk(k) = q2x(i,k)
1527 hgame(k) = hgame2d(i,k)
1528 enddo
1529 !
1530 do k = kts,kte-1
1531 if(pblflg(i).and.k.le.kpbl(i)) then
1532 zfacdx = 0.2*hpbl(i)/za(i,k)
1533 delxy = sqrt(dx*dy)*max(zfacdx,1.0)
1534 ptke1(k+1) = ptke(delxy,hpbl(i))
1535 endif
1536 enddo
1537 !
1538 do k = kts,kte+1
1539 zqk(k) = zq(i,k)
1540 enddo
1541 !
1542 do k = kts,kte*ndiff
1543 qxk(k) = qx(i,k)
1544 enddo
1545 !
1546 do k = kts+1,kte
1547 akmk(k) = xkzm(i,k-1)
1548 akhk(k) = xkzh(i,k-1)
1549 mfk(k) = mf(i,k-1)/xkzh(i,k-1)
1550 ufxpblk(k) = ufxpbl(i)*zfacent(i,k-1)/xkzm(i,k-1)
1551 vfxpblk(k) = vfxpbl(i)*zfacent(i,k-1)/xkzm(i,k-1)
1552 qfxpblk(k) = qfxpbl(i)*zfacent(i,k-1)/xkzq(i,k-1)
1553 enddo
1554 !
1555 if(pblflg(i)) then
1556 k = kpbl(i) - 1
1557 dex = 0.25*(q2xk(k+2)-q2xk(k))
1558 efxpbl(i) = we(i)*dex
1559 endif
1560 !
1561 delxy=sqrt(dx*dy)
1562 !
1563 !---- find the mixing length
1564 !
1565 call mixlen(lmh,uxk,vxk,txk,thxk,qxk(kts),qxk(kte+1) &
1566 ,q2xk,zqk,ust(i),corf(i),epshol(i) &
1567 ,s2,gh,rig,el &
1568 ,hpbl(i),kpbl(i),lmxl,ct(i) &
1569 ,hgamu(i),hgamv(i),hgamq(i),pblflg(i) &
1570 ,mfk,ufxpblk,vfxpblk,qfxpblk &
1571 ,ep1,karman,cp &
1572 ,ids,ide,jds,jde,kds,kde &
1573 ,ims,ime,jms,jme,kms,kme &
1574 ,its,ite,jts,jte,kts,kte )
1575 !
1576 !---- solve for the production/dissipation of the turbulent kinetic energy
1577 !
1578 call prodq2(lmh,dt,ust(i),s2,rig,q2xk,el,zqk,akmk,akhk &
1579 ,uxk,vxk,thxk,thvxk &
1580 ,hgamu(i),hgamv(i),hgamq(i),delxy &
1581 ,hpbl(i),pblflg(i),kpbl(i) &
1582 ,mfk,ufxpblk,vfxpblk,qfxpblk &
1583 ,ep1 &
1584 ,ids,ide,jds,jde,kds,kde &
1585 ,ims,ime,jms,jme,kms,kme &
1586 ,its,ite,jts,jte,kts,kte )
1587 !
1588 !
1589 !---- carry out the vertical diffusion of turbulent kinetic energy
1590 !
1591 call vdifq(lmh,dt,q2xk,el,zqk &
1592 ,akhk,ptke1 &
1593 ,hgame,hpbl(i),pblflg(i),kpbl(i) &
1594 ,efxpbl(i) &
1595 ,ids,ide,jds,jde,kds,kde &
1596 ,ims,ime,jms,jme,kms,kme &
1597 ,its,ite,jts,jte,kts,kte )
1598 !
1599 !---- save the new tke and mixing length.
1600 !
1601 do k = kts,kte
1602 q2x(i,k) = amax1(q2xk(k),epsq2l)
1603 tke(i,k) = 0.5*q2x(i,k)
1604 if(k/=kts) el_pbl(i,k) = el(k) ! el is not defined at kte
1605 enddo
1606 !
1607 enddo tke_calculation
1608 endif
1609 !
1610 !---- end of tke calculation
1611 !
1612 !
1613 !---- end of vertical diffusion
1614 !
1615 end subroutine shinhong2d
1616 !-------------------------------------------------------------------------------
1617 !
1618 !-------------------------------------------------------------------------------
1619 subroutine tridi1n(cl,cm,cu,r1,r2,au,f1,f2,its,ite,kts,kte,nt)
1620 !-------------------------------------------------------------------------------
1621 implicit none
1622 !-------------------------------------------------------------------------------
1623 !
1624 integer, intent(in ) :: its,ite, kts,kte, nt
1625 !
1626 real, dimension( its:ite, kts+1:kte+1 ) , &
1627 intent(in ) :: cl
1628 !
1629 real, dimension( its:ite, kts:kte ) , &
1630 intent(in ) :: cm, &
1631 r1
1632 real, dimension( its:ite, kts:kte,nt ) , &
1633 intent(in ) :: r2
1634 !
1635 real, dimension( its:ite, kts:kte ) , &
1636 intent(inout) :: au, &
1637 cu, &
1638 f1
1639 real, dimension( its:ite, kts:kte,nt ) , &
1640 intent(inout) :: f2
1641 !
1642 real :: fk
1643 integer :: i,k,l,n,it
1644 !
1645 !-------------------------------------------------------------------------------
1646 !
1647 l = ite
1648 n = kte
1649 !
1650 do i = its,l
1651 fk = 1./cm(i,1)
1652 au(i,1) = fk*cu(i,1)
1653 f1(i,1) = fk*r1(i,1)
1654 enddo
1655 !
1656 do it = 1,nt
1657 do i = its,l
1658 fk = 1./cm(i,1)
1659 f2(i,1,it) = fk*r2(i,1,it)
1660 enddo
1661 enddo
1662 !
1663 do k = kts+1,n-1
1664 do i = its,l
1665 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1666 au(i,k) = fk*cu(i,k)
1667 f1(i,k) = fk*(r1(i,k)-cl(i,k)*f1(i,k-1))
1668 enddo
1669 enddo
1670 !
1671 do it = 1,nt
1672 do k = kts+1,n-1
1673 do i = its,l
1674 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1675 f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
1676 enddo
1677 enddo
1678 enddo
1679 !
1680 do i = its,l
1681 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1682 f1(i,n) = fk*(r1(i,n)-cl(i,n)*f1(i,n-1))
1683 enddo
1684 !
1685 do it = 1,nt
1686 do i = its,l
1687 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1688 f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
1689 enddo
1690 enddo
1691 !
1692 do k = n-1,kts,-1
1693 do i = its,l
1694 f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1)
1695 enddo
1696 enddo
1697 !
1698 do it = 1,nt
1699 do k = n-1,kts,-1
1700 do i = its,l
1701 f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1702 enddo
1703 enddo
1704 enddo
1705 !
1706 end subroutine tridi1n
1707 !-------------------------------------------------------------------------------
1708 !
1709 !-------------------------------------------------------------------------------
1710 subroutine tridin_ysu(cl,cm,cu,r2,au,f2,its,ite,kts,kte,nt)
1711 !-------------------------------------------------------------------------------
1712 implicit none
1713 !-------------------------------------------------------------------------------
1714 !
1715 integer, intent(in ) :: its,ite, kts,kte, nt
1716 !
1717 real, dimension( its:ite, kts+1:kte+1 ) , &
1718 intent(in ) :: cl
1719 !
1720 real, dimension( its:ite, kts:kte ) , &
1721 intent(in ) :: cm
1722 real, dimension( its:ite, kts:kte,nt ) , &
1723 intent(in ) :: r2
1724 !
1725 real, dimension( its:ite, kts:kte ) , &
1726 intent(inout) :: au, &
1727 cu
1728 real, dimension( its:ite, kts:kte,nt ) , &
1729 intent(inout) :: f2
1730 !
1731 real :: fk
1732 integer :: i,k,l,n,it
1733 !
1734 !-------------------------------------------------------------------------------
1735 !
1736 l = ite
1737 n = kte
1738 !
1739 do it = 1,nt
1740 do i = its,l
1741 fk = 1./cm(i,1)
1742 au(i,1) = fk*cu(i,1)
1743 f2(i,1,it) = fk*r2(i,1,it)
1744 enddo
1745 enddo
1746 !
1747 do it = 1,nt
1748 do k = kts+1,n-1
1749 do i = its,l
1750 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1751 au(i,k) = fk*cu(i,k)
1752 f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
1753 enddo
1754 enddo
1755 enddo
1756 !
1757 do it = 1,nt
1758 do i = its,l
1759 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1760 f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
1761 enddo
1762 enddo
1763 !
1764 do it = 1,nt
1765 do k = n-1,kts,-1
1766 do i = its,l
1767 f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1768 enddo
1769 enddo
1770 enddo
1771 !
1772 end subroutine tridin_ysu
1773 !-------------------------------------------------------------------------------
1774 !
1775 !-------------------------------------------------------------------------------
1776 subroutine mixlen(lmh,u,v,t,the,q,cwm,q2,z,ustar,corf,epshol, &
1777 s2,gh,ri,el,hpbl,lpbl,lmxl,ct, &
1778 hgamu,hgamv,hgamq,pblflg, &
1779 mf,ufxpbl,vfxpbl,qfxpbl, &
1780 p608,vkarman,cp, &
1781 ids,ide,jds,jde,kds,kde, &
1782 ims,ime,jms,jme,kms,kme, &
1783 its,ite,jts,jte,kts,kte)
1784 !-------------------------------------------------------------------------------
1785 implicit none
1786 !-------------------------------------------------------------------------------
1787 ! qnse model constants
1788 !-------------------------------------------------------------------------------
1789 real,parameter :: blckdr=0.0063,cn=0.75
1790 real,parameter :: eps1=1.e-12,epsl=0.32,epsru=1.e-7,epsrs=1.e-7
1791 real,parameter :: el0max=1000.,el0min=1.,elfc=0.23*0.5
1792 real,parameter :: alph=0.30,beta=1./273.,g=9.81,btg=beta*g
1793 real,parameter :: a1=0.659888514560862645,a2x=0.6574209922667784586
1794 real,parameter :: b1=11.87799326209552761,b2=7.226971804046074028
1795 real,parameter :: c1=0.000830955950095854396
1796 real,parameter :: adnh= 9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg
1797 real,parameter :: adnm=18.*a1*a1*a2x*(b2-3.*a2x)*btg
1798 real,parameter :: bdnh= 3.*a2x*(7.*a1+b2)*btg,bdnm= 6.*a1*a1
1799 !-------------------------------------------------------------------------------
1800 ! free term in the equilibrium equation for (l/q)**2
1801 !-------------------------------------------------------------------------------
1802 real,parameter :: aeqh=9.*a1*a2x*a2x*b1*btg*btg &
1803 +9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg
1804 real,parameter :: aeqm=3.*a1*a2x*b1*(3.*a2x+3.*b2*c1+18.*a1*c1-b2) &
1805 *btg+18.*a1*a1*a2x*(b2-3.*a2x)*btg
1806 !-------------------------------------------------------------------------------
1807 ! forbidden turbulence area
1808 !-------------------------------------------------------------------------------
1809 real,parameter :: requ=-aeqh/aeqm
1810 real,parameter :: epsgh=1.e-9,epsgm=requ*epsgh
1811 !-------------------------------------------------------------------------------
1812 ! near isotropy for shear turbulence, ww/q2 lower limit
1813 !-------------------------------------------------------------------------------
1814 real,parameter :: ubryl=(18.*requ*a1*a1*a2x*b2*c1*btg &
1815 +9.*a1*a2x*a2x*b2*btg*btg) &
1816 /(requ*adnm+adnh)
1817 real,parameter :: ubry=(1.+epsrs)*ubryl,ubry3=3.*ubry
1818 real,parameter :: aubh=27.*a1*a2x*a2x*b2*btg*btg-adnh*ubry3
1819 real,parameter :: aubm=54.*a1*a1*a2x*b2*c1*btg -adnm*ubry3
1820 real,parameter :: bubh=(9.*a1*a2x+3.*a2x*b2)*btg-bdnh*ubry3
1821 real,parameter :: bubm=18.*a1*a1*c1 -bdnm*ubry3
1822 real,parameter :: cubr=1.-ubry3,rcubr=1./cubr
1823 !-------------------------------------------------------------------------------
1824 ! k profile constants
1825 !-------------------------------------------------------------------------------
1826 real,parameter :: elcbl=0.77
1827 !-------------------------------------------------------------------------------
1828 !
1829 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
1830 ims,ime, jms,jme, kms,kme, &
1831 its,ite, jts,jte, kts,kte
1832 integer, intent(in ) :: lmh,lmxl,lpbl
1833 !
1834 real, intent(in ) :: p608,vkarman,cp
1835 real, intent(in ) :: hpbl,corf,ustar,hgamu,hgamv,hgamq
1836 real, intent(inout) :: ct,epshol
1837 !
1838 real, dimension( kts:kte ) , &
1839 intent(in ) :: cwm, &
1840 q, &
1841 q2, &
1842 t, &
1843 the, &
1844 u, &
1845 v
1846 !
1847 real, dimension( kts+1:kte ) , &
1848 intent(in ) :: mf, &
1849 ufxpbl, &
1850 vfxpbl, &
1851 qfxpbl
1852 !
1853 real, dimension( kts:kte+1 ) , &
1854 intent(in ) :: z
1855 !
1856 real, dimension( kts+1:kte ) , &
1857 intent(out ) :: el, &
1858 ri, &
1859 gh, &
1860 s2
1861 !
1862 logical,intent(in) :: pblflg
1863 !
1864 ! local vars
1865 !
1866 integer :: k,lpblm
1867 real :: suk,svk,elocp
1868 real :: a,aden,b,bden,aubr,bubr,blmx,el0,eloq2x,ghl,s2l, &
1869 qol2st,qol2un,qdzl,rdz,sq,srel,szq,tem,thm,vkrmz,rlambda, &
1870 rlb,rln,f
1871 real :: ckp
1872 real, dimension( kts:kte ) :: q1, &
1873 en2
1874 real, dimension( kts+1:kte ) :: dth, &
1875 elm, &
1876 rel
1877 !
1878 !-------------------------------------------------------------------------------
1879 !
1880 elocp=2.72e6/cp
1881 ct=0.
1882 !
1883 do k = kts,kte
1884 q1(k) = 0.
1885 enddo
1886 !
1887 do k = kts+1,kte
1888 dth(k) = the(k)-the(k-1)
1889 enddo
1890 !
1891 do k = kts+2,kte
1892 if(dth(k)>0..and.dth(k-1)<=0.)then
1893 dth(k)=dth(k)+ct
1894 exit
1895 endif
1896 enddo
1897 !
1898 ! compute local gradient richardson number
1899 !
1900 do k = kte,kts+1,-1
1901 rdz=2./(z(k+1)-z(k-1))
1902 s2l=((u(k)-u(k-1))**2+(v(k)-v(k-1))**2)*rdz*rdz ! s**2
1903 if(pblflg.and.k.le.lpbl)then
1904 suk=(u(k)-u(k-1))*rdz
1905 svk=(v(k)-v(k-1))*rdz
1906 s2l=(suk-hgamu/hpbl-ufxpbl(k))*suk+(svk-hgamv/hpbl-vfxpbl(k))*svk
1907 endif
1908 s2l=max(s2l,epsgm)
1909 s2(k)=s2l
1910 !
1911 tem=(t(k)+t(k-1))*0.5
1912 thm=(the(k)+the(k-1))*0.5
1913 a=thm*p608
1914 b=(elocp/tem-1.-p608)*thm
1915 ghl=(dth(k)*((q(k)+q(k-1)+cwm(k)+cwm(k-1))*(0.5*p608)+1.) &
1916 +(q(k)-q(k-1)+cwm(k)-cwm(k-1))*a &
1917 +(cwm(k)-cwm(k-1))*b)*rdz ! dtheta/dz
1918 if(pblflg.and.k.le.lpbl)then
1919 ghl=ghl-mf(k)-(hgamq/hpbl+qfxpbl(k))*a
1920 endif
1921 if(abs(ghl)<=epsgh)ghl=epsgh
1922 !
1923 en2(k)=ghl*g/thm ! n**2
1924 gh(k)=ghl
1925 ri(k)=en2(k)/s2l
1926 enddo
1927 !
1928 ! find maximum mixing lengths and the level of the pbl top
1929 !
1930 do k = kte,kts+1,-1
1931 s2l=s2(k)
1932 ghl=gh(k)
1933 if(ghl>=epsgh)then
1934 if(s2l/ghl<=requ)then
1935 elm(k)=epsl
1936 else
1937 aubr=(aubm*s2l+aubh*ghl)*ghl
1938 bubr= bubm*s2l+bubh*ghl
1939 qol2st=(-0.5*bubr+sqrt(bubr*bubr*0.25-aubr*cubr))*rcubr
1940 eloq2x=1./qol2st
1941 elm(k)=max(sqrt(eloq2x*q2(k)),epsl)
1942 endif
1943 else
1944 aden=(adnm*s2l+adnh*ghl)*ghl
1945 bden= bdnm*s2l+bdnh*ghl
1946 qol2un=-0.5*bden+sqrt(bden*bden*0.25-aden)
1947 eloq2x=1./(qol2un+epsru) ! repsr1/qol2un
1948 elm(k)=max(sqrt(eloq2x*q2(k)),epsl)
1949 endif
1950 enddo
1951 !
1952 do k = lpbl,lmh,-1
1953 q1(k)=sqrt(q2(k))
1954 enddo
1955 !
1956 szq=0.
1957 sq =0.
1958 do k = kte,kts+1,-1
1959 qdzl=(q1(k)+q1(k-1))*(z(k)-z(k-1))
1960 szq=(z(k)+z(k-1)-z(lmh)-z(lmh))*qdzl+szq
1961 sq=qdzl+sq
1962 enddo
1963 !
1964 ! computation of asymptotic l in blackadar formula
1965 !
1966 el0=min(alph*szq*0.5/sq,el0max)
1967 el0=max(el0 ,el0min)
1968 !
1969 ! above the pbl top
1970 !
1971 lpblm=min(lpbl+1,kte)
1972 do k = kte,lpblm,-1
1973 el(k)=(z(k+1)-z(k-1))*elfc
1974 rel(k)=el(k)/elm(k)
1975 enddo
1976 !
1977 ! inside the pbl
1978 !
1979 epshol=min(epshol,0.0)
1980 ckp=elcbl*((1.0-8.0*epshol)**(1./3.))
1981 if(lpbl>lmh)then
1982 do k = lpbl,lmh+1,-1
1983 vkrmz=(z(k)-z(lmh))*vkarman
1984 if(pblflg) then
1985 vkrmz=ckp*(z(k)-z(lmh))*vkarman
1986 el(k)=vkrmz/(vkrmz/el0+1.)
1987 else
1988 el(k)=vkrmz/(vkrmz/el0+1.)
1989 endif
1990 rel(k)=el(k)/elm(k)
1991 enddo
1992 endif
1993 !
1994 do k = lpbl-1,lmh+2,-1
1995 srel=min(((rel(k-1)+rel(k+1))*0.5+rel(k))*0.5,rel(k))
1996 el(k)=max(srel*elm(k),epsl)
1997 enddo
1998 !
1999 ! mixing length for the qnse model in stable case
2000 !
2001 f=max(corf,eps1)
2002 rlambda=f/(blckdr*ustar)
2003 do k = kte,kts+1,-1
2004 if(en2(k)>=0.0)then ! stable case
2005 vkrmz=(z(k)-z(lmh))*vkarman
2006 rlb=rlambda+1./vkrmz
2007 rln=sqrt(2.*en2(k)/q2(k))/cn
2008 el(k)=1./(rlb+rln)
2009 endif
2010 enddo
2011 !
2012 end subroutine mixlen
2013 !-------------------------------------------------------------------------------
2014 !
2015 !-------------------------------------------------------------------------------
2016 subroutine prodq2(lmh,dtturbl,ustar,s2,ri,q2,el,z,akm,akh, &
2017 uxk,vxk,thxk,thvxk, &
2018 hgamu,hgamv,hgamq,delxy, &
2019 hpbl,pblflg,kpbl, &
2020 mf,ufxpbl,vfxpbl,qfxpbl, &
2021 p608, &
2022 ids,ide,jds,jde,kds,kde, &
2023 ims,ime,jms,jme,kms,kme, &
2024 its,ite,jts,jte,kts,kte)
2025 !-------------------------------------------------------------------------------
2026 implicit none
2027 !-------------------------------------------------------------------------------
2028 !
2029 real,parameter :: epsq2l = 0.01,c0 = 0.55,ceps = 16.6,g = 9.81
2030 !
2031 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
2032 ims,ime, jms,jme, kms,kme, &
2033 its,ite, jts,jte, kts,kte
2034 integer, intent(in ) :: lmh,kpbl
2035 !
2036 real, intent(in ) :: p608,dtturbl,ustar
2037 real, intent(in ) :: hgamu,hgamv,hgamq,delxy,hpbl
2038 !
2039 logical, intent(in ) :: pblflg
2040 !
2041 real, dimension( kts:kte ) , &
2042 intent(in ) :: uxk, &
2043 vxk, &
2044 thxk, &
2045 thvxk
2046 real, dimension( kts+1:kte ) , &
2047 intent(in ) :: s2, &
2048 ri, &
2049 akm, &
2050 akh, &
2051 el, &
2052 mf, &
2053 ufxpbl, &
2054 vfxpbl, &
2055 qfxpbl
2056 !
2057 real, dimension( kts:kte+1 ) , &
2058 intent(in ) :: z
2059 !
2060 real, dimension( kts:kte ) , &
2061 intent(inout) :: q2
2062 !
2063 ! local vars
2064 !
2065 integer :: k
2066 !
2067 real :: s2l,q2l,deltaz,akml,akhl,en2,pr,bpr,dis,rc02
2068 real :: suk,svk,gthvk,govrthvk,pru,prv
2069 real :: thm,disel
2070 !
2071 !-------------------------------------------------------------------------------
2072 !
2073 rc02=2.0/(c0*c0)
2074 !
2075 ! start of production/dissipation loop
2076 !
2077 main_integration: do k = kts+1,kte
2078 deltaz=0.5*(z(k+1)-z(k-1))
2079 s2l=s2(k)
2080 q2l=q2(k)
2081 suk=(uxk(k)-uxk(k-1))/deltaz
2082 svk=(vxk(k)-vxk(k-1))/deltaz
2083 gthvk=(thvxk(k)-thvxk(k-1))/deltaz
2084 govrthvk=g/(0.5*(thvxk(k)+thvxk(k-1)))
2085 akml=akm(k)
2086 akhl=akh(k)
2087 en2=ri(k)*s2l !n**2
2088 thm=(thxk(k)+thxk(k-1))*0.5
2089 !
2090 ! turbulence production term
2091 !
2092 if(pblflg.and.k.le.kpbl)then
2093 pru=(akml*(suk-hgamu/hpbl-ufxpbl(k)))*suk
2094 prv=(akml*(svk-hgamv/hpbl-vfxpbl(k)))*svk
2095 else
2096 pru=akml*suk*suk
2097 prv=akml*svk*svk
2098 endif
2099 pr=pru+prv
2100 !
2101 ! buoyancy production
2102 !
2103 if(pblflg.and.k.le.kpbl)then
2104 bpr=(akhl*(gthvk-mf(k)-(hgamq/hpbl+qfxpbl(k))*p608*thm))*govrthvk
2105 else
2106 bpr=akhl*gthvk*govrthvk
2107 endif
2108 !
2109 ! dissipation
2110 !
2111 disel=min(delxy,ceps*el(k))
2112 dis=(q2l)**1.5/disel
2113 !
2114 q2l=q2l+2.0*(pr-bpr-dis)*dtturbl
2115 q2(k)=amax1(q2l,epsq2l)
2116 !
2117 ! end of production/dissipation loop
2118 !
2119 enddo main_integration
2120 !
2121 ! lower boundary condition for q2
2122 !
2123 q2(kts)=amax1(rc02*ustar*ustar,epsq2l)
2124 !
2125 end subroutine prodq2
2126 !-------------------------------------------------------------------------------
2127 !
2128 !-------------------------------------------------------------------------------
2129 subroutine vdifq(lmh,dtdif,q2,el,z, &
2130 akhk,ptke1, &
2131 hgame,hpbl,pblflg,kpbl, &
2132 efxpbl, &
2133 ids,ide,jds,jde,kds,kde, &
2134 ims,ime,jms,jme,kms,kme, &
2135 its,ite,jts,jte,kts,kte)
2136 !-------------------------------------------------------------------------------
2137 implicit none
2138 !-------------------------------------------------------------------------------
2139 !
2140 real,parameter :: c_k=1.0,esq=5.0
2141 !
2142 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
2143 ims,ime, jms,jme, kms,kme, &
2144 its,ite, jts,jte, kts,kte
2145 integer, intent(in ) :: lmh,kpbl
2146 !
2147 real, intent(in ) :: dtdif,hpbl,efxpbl
2148 !
2149 logical, intent(in ) :: pblflg
2150 !
2151 real, dimension( kts:kte ) , &
2152 intent(in ) :: hgame, &
2153 ptke1
2154 real, dimension( kts+1:kte ) , &
2155 intent(in ) :: el, &
2156 akhk
2157 real, dimension( kts:kte+1 ) , &
2158 intent(in ) :: z
2159 !
2160 real, dimension( kts:kte ) , &
2161 intent(inout) :: q2
2162 !
2163 ! local vars
2164 !
2165 integer :: k
2166 !
2167 real :: aden,akqs,bden,besh,besm,cden,cf,dtozs,ell,eloq2,eloq4
2168 real :: elqdz,esh,esm,esqhf,ghl,gml,q1l,rden,rdz
2169 real :: zak
2170 !
2171 real, dimension( kts+1:kte ) :: zfacentk
2172 real, dimension( kts+2:kte ) :: akq, &
2173 cm, &
2174 cr, &
2175 dtoz, &
2176 rsq2
2177 !
2178 !-------------------------------------------------------------------------------
2179 !
2180 ! vertical turbulent diffusion
2181 !
2182 esqhf=0.5*esq
2183 do k = kts+1,kte
2184 zak=0.5*(z(k)+z(k-1)) !zak of vdifq = za(k-1) of shinhong2d
2185 zfacentk(k)=(zak/hpbl)**3.0
2186 enddo
2187 !
2188 do k = kte,kts+2,-1
2189 dtoz(k)=(dtdif+dtdif)/(z(k+1)-z(k-1))
2190 akq(k)=c_k*(akhk(k)/(z(k+1)-z(k-1))+akhk(k-1)/(z(k)-z(k-2)))
2191 akq(k)=akq(k)*ptke1(k)
2192 cr(k)=-dtoz(k)*akq(k)
2193 enddo
2194 !
2195 akqs=c_k*akhk(kts+1)/(z(kts+2)-z(kts))
2196 akqs=akqs*ptke1(kts+1)
2197 cm(kte)=dtoz(kte)*akq(kte)+1.
2198 rsq2(kte)=q2(kte)
2199 !
2200 do k = kte-1,kts+2,-1
2201 cf=-dtoz(k)*akq(k+1)/cm(k+1)
2202 cm(k)=-cr(k+1)*cf+(akq(k+1)+akq(k))*dtoz(k)+1.
2203 rsq2(k)=-rsq2(k+1)*cf+q2(k)
2204 if(pblflg.and.k.lt.kpbl) then
2205 rsq2(k)=rsq2(k)-dtoz(k)*(2.0*hgame(k)/hpbl)*akq(k+1)*(z(k+1)-z(k)) &
2206 +dtoz(k)*(2.0*hgame(k-1)/hpbl)*akq(k)*(z(k)-z(k-1))
2207 rsq2(k)=rsq2(k)-dtoz(k)*2.0*efxpbl*zfacentk(k+1) &
2208 +dtoz(k)*2.0*efxpbl*zfacentk(k)
2209 endif
2210 enddo
2211 !
2212 dtozs=(dtdif+dtdif)/(z(kts+2)-z(kts))
2213 cf=-dtozs*akq(lmh+2)/cm(lmh+2)
2214 !
2215 if(pblflg.and.((lmh+1).lt.kpbl)) then
2216 q2(lmh+1)=(dtozs*akqs*q2(lmh)-rsq2(lmh+2)*cf+q2(lmh+1) &
2217 -dtozs*(2.0*hgame(lmh+1)/hpbl)*akq(lmh+2)*(z(lmh+2)-z(lmh+1)) &
2218 +dtozs*(2.0*hgame(lmh)/hpbl)*akqs*(z(lmh+1)-z(lmh)))
2219 q2(lmh+1)=q2(lmh+1)-dtozs*2.0*efxpbl*zfacentk(lmh+2) &
2220 +dtozs*2.0*efxpbl*zfacentk(lmh+1)
2221 q2(lmh+1)=q2(lmh+1)/((akq(lmh+2)+akqs)*dtozs-cr(lmh+2)*cf+1.)
2222 else
2223 q2(lmh+1)=(dtozs*akqs*q2(lmh)-rsq2(lmh+2)*cf+q2(lmh+1)) &
2224 /((akq(lmh+2)+akqs)*dtozs-cr(lmh+2)*cf+1.)
2225 endif
2226 !
2227 do k = lmh+2,kte
2228 q2(k)=(-cr(k)*q2(k-1)+rsq2(k))/cm(k)
2229 enddo
2230 !
2231 end subroutine vdifq
2232 !-------------------------------------------------------------------------------
2233 !
2234 !-------------------------------------------------------------------------------
2235 subroutine shinhonginit(rublten,rvblten,rthblten,rqvblten, &
2236 rqcblten,rqiblten, &
2237 tke_pbl, &
2238 p_qi,p_first_scalar, &
2239 restart, allowed_to_read, &
2240 ids, ide, jds, jde, kds, kde, &
2241 ims, ime, jms, jme, kms, kme, &
2242 its, ite, jts, jte, kts, kte )
2243 !-------------------------------------------------------------------------------
2244 implicit none
2245 !-------------------------------------------------------------------------------
2246 !
2247 real,parameter :: epsq2l = 0.01
2248 logical , intent(in) :: restart, allowed_to_read
2249 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
2250 ims, ime, jms, jme, kms, kme, &
2251 its, ite, jts, jte, kts, kte
2252 integer , intent(in) :: p_qi,p_first_scalar
2253 real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: &
2254 rublten, &
2255 rvblten, &
2256 rthblten, &
2257 rqvblten, &
2258 rqcblten, &
2259 rqiblten
2260 real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: &
2261 tke_pbl
2262 integer :: i, j, k, itf, jtf, ktf
2263 !
2264 jtf = min0(jte,jde-1)
2265 ktf = min0(kte,kde-1)
2266 itf = min0(ite,ide-1)
2267 !
2268 if(.not.restart)then
2269 do j = jts,jtf
2270 do k = kts,ktf
2271 do i = its,itf
2272 rublten(i,k,j) = 0.
2273 rvblten(i,k,j) = 0.
2274 rthblten(i,k,j) = 0.
2275 rqvblten(i,k,j) = 0.
2276 rqcblten(i,k,j) = 0.
2277 tke_pbl(i,k,j) = epsq2l/2.
2278 enddo
2279 enddo
2280 enddo
2281 endif
2282 !
2283 if (p_qi .ge. p_first_scalar .and. .not.restart) then
2284 do j = jts,jtf
2285 do k = kts,ktf
2286 do i = its,itf
2287 rqiblten(i,k,j) = 0.
2288 enddo
2289 enddo
2290 enddo
2291 endif
2292 !
2293 end subroutine shinhonginit
2294 !-------------------------------------------------------------------------------
2295 !
2296 !-------------------------------------------------------------------------------
2297 function pu(d,h)
2298 !-------------------------------------------------------------------------------
2299 implicit none
2300 !-------------------------------------------------------------------------------
2301 real :: pu
2302 real,parameter :: pmin = 0.0,pmax = 1.0
2303 real,parameter :: a1 = 1.0, a2 = 0.070, a3 = 1.0, a4 = 0.142, a5 = 0.071
2304 real,parameter :: b1 = 2.0, b2 = 0.6666667
2305 real :: d,h,doh,num,den
2306 !
2307 if (h .ne. 0) then
2308 doh=d/h
2309 num=a1*(doh)**b1+a2*(doh)**b2
2310 den=a3*(doh)**b1+a4*(doh)**b2+a5
2311 pu=num/den
2312 else
2313 pu=1.
2314 endif
2315 pu=max(pu,pmin)
2316 pu=min(pu,pmax)
2317 !
2318 return
2319 end function
2320 !-------------------------------------------------------------------------------
2321 !
2322 !-------------------------------------------------------------------------------
2323 function pq(d,h)
2324 !-------------------------------------------------------------------------------
2325 implicit none
2326 !-------------------------------------------------------------------------------
2327 real :: pq
2328 real,parameter :: pmin = 0.0,pmax = 1.0
2329 real,parameter :: a1 = 1.0, a2 = -0.098, a3 = 1.0, a4 = 0.106, a5 = 0.5
2330 real,parameter :: b1 = 2.0
2331 real :: d,h,doh,num,den
2332 !
2333 if (h .ne. 0) then
2334 doh=d/h
2335 num=a1*(doh)**b1+a2
2336 den=a3*(doh)**b1+a4
2337 pq=a5*num/den+(1.-a5)
2338 else
2339 pq=1.
2340 endif
2341 pq=max(pq,pmin)
2342 pq=min(pq,pmax)
2343 !
2344 return
2345 end function
2346 !-------------------------------------------------------------------------------
2347 !
2348 !-------------------------------------------------------------------------------
2349 function pthnl(d,h)
2350 !-------------------------------------------------------------------------------
2351 implicit none
2352 !-------------------------------------------------------------------------------
2353 real :: pthnl
2354 real,parameter :: pmin = 0.0,pmax = 1.0
2355 real,parameter :: a1 = 1.000, a2 = 0.936, a3 = -1.110, &
2356 a4 = 1.000, a5 = 0.312, a6 = 0.329, a7 = 0.243
2357 real,parameter :: b1 = 2.0, b2 = 0.875
2358 real :: d,h,doh,num,den
2359 !
2360 if (h .ne. 0) then
2361 doh=d/h
2362 num=a1*(doh)**b1+a2*(doh)**b2+a3
2363 den=a4*(doh)**b1+a5*(doh)**b2+a6
2364 pthnl=a7*num/den+(1.-a7)
2365 else
2366 pthnl=1.
2367 endif
2368 pthnl=max(pthnl,pmin)
2369 pthnl=min(pthnl,pmax)
2370 !
2371 return
2372 end function
2373 !-------------------------------------------------------------------------------
2374 !
2375 !-------------------------------------------------------------------------------
2376 function pthl(d,h)
2377 !-------------------------------------------------------------------------------
2378 implicit none
2379 !-------------------------------------------------------------------------------
2380 real :: pthl
2381 real,parameter :: pmin = 0.0,pmax = 1.0
2382 real,parameter :: a1 = 1.000, a2 = 0.870, a3 = -0.913, &
2383 a4 = 1.000, a5 = 0.153, a6 = 0.278, a7 = 0.280
2384 real,parameter :: b1 = 2.0, b2 = 0.5
2385 real :: d,h,doh,num,den
2386 !
2387 if (h .ne. 0) then
2388 doh=d/h
2389 num=a1*(doh)**b1+a2*(doh)**b2+a3
2390 den=a4*(doh)**b1+a5*(doh)**b2+a6
2391 pthl=a7*num/den+(1.-a7)
2392 else
2393 pthl=1.
2394 endif
2395 pthl=max(pthl,pmin)
2396 pthl=min(pthl,pmax)
2397 !
2398 return
2399 end function
2400 !-------------------------------------------------------------------------------
2401 !
2402 !-------------------------------------------------------------------------------
2403 function ptke(d,h)
2404 !-------------------------------------------------------------------------------
2405 implicit none
2406 !-------------------------------------------------------------------------------
2407 real :: ptke
2408 real,parameter :: pmin = 0.0,pmax = 1.0
2409 real,parameter :: a1 = 1.000, a2 = 0.070, &
2410 a3 = 1.000, a4 = 0.142, a5 = 0.071
2411 real,parameter :: b1 = 2.0, b2 = 0.6666667
2412 real :: d,h,doh,num,den
2413 !
2414 if (h .ne. 0) then
2415 doh=d/h
2416 num=a1*(doh)**b1+a2*(doh)**b2
2417 den=a3*(doh)**b1+a4*(doh)**b2+a5
2418 ptke=num/den
2419 else
2420 ptke=1.
2421 endif
2422 ptke=max(ptke,pmin)
2423 ptke=min(ptke,pmax)
2424 !
2425 return
2426 end function
2427 !-------------------------------------------------------------------------------
2428 !
2429 !-------------------------------------------------------------------------------
2430 end module module_bl_shinhong
2431 !-------------------------------------------------------------------------------