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