CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
moninshoc.f
1
3
5 module moninshoc
6
7 use mfpbl_mod
8 use tridi_mod
9
10 contains
11
12 subroutine moninshoc_init (do_shoc, errmsg, errflg)
13
14 implicit none
15 logical, intent(in) :: do_shoc
16 character(len=*), intent(out) :: errmsg
17 integer, intent(out) :: errflg
18
19 ! Initialize CCPP error handling variables
20 errmsg = ''
21 errflg = 0
22
23 ! Consistency checks
24 if (.not. do_shoc) then
25 errflg = 1
26 write(errmsg,'(*(a))') 'Logic error: do_shoc = .false.'
27 return
28 end if
29
30 end subroutine moninshoc_init
31
32!!!!! ========================================================== !!!!!
33! subroutine 'moninshoc' computes pbl height and applies vertical diffusion
34! using the coefficient provided by the SHOC scheme (from previous step)
35! 2015-05-04 - Shrinivas Moorthi - original version based on monin
36! 2018-03-21 - Shrinivas Moorthi - fixed a bug related to tke vertical diffusion
37! and gneralized the tke location in tracer array
38! 2018-03-23 - Shrinivas Moorthi - used twice the momentum diffusion coefficient
39! for tke as in Deardorff (1980) - added tridi1
40!
44 subroutine moninshoc_run (im,km,ntrac,ntcw,ncnd,dv,du,tau,rtg,
45 & u1,v1,t1,q1,tkh,prnum,ntke,
46 & psk,rbsoil,zorl,u10m,v10m,fm,fh,
47 & tsea,heat,evap,stress,spd1,kpbl,
48 & prsi,del,prsl,prslk,phii,phil,delt,
49 & dusfc,dvsfc,dtsfc,dqsfc,dkt,hpbl,
50 & kinver,xkzm_m,xkzm_h,xkzm_s,xkzminv,
51 & grav,rd,cp,hvap,fv,ntoz,dtend,dtidx,
52 & index_of_temperature,index_of_x_wind,
53 & index_of_y_wind,index_of_process_pbl,
54 & gen_tend,ldiag3d,ntqv,errmsg,errflg)
55!
56 use machine , only : kind_phys
57 use funcphys , only : fpvs
58
59 implicit none
60!
61! arguments
62!
63 integer, intent(in) :: im,
64 & km, ntrac, ntcw, ncnd, ntke, ntoz
65 integer, dimension(:), intent(in) :: kinver
66 real(kind=kind_phys), intent(in) :: delt,
67 & xkzm_m, xkzm_h, xkzm_s, xkzminv
68 real(kind=kind_phys), intent(in) :: grav,
69 & rd, cp, hvap, fv
70 real(kind=kind_phys), dimension(:), intent(in) :: psk,
71 & rbsoil, zorl, u10m, v10m, fm, fh, tsea, heat, evap, stress, spd1
72 real(kind=kind_phys), dimension(:,:), intent(in) :: u1, v1,
73 & t1, tkh, del, prsl, phil, prslk
74 real(kind=kind_phys), dimension(:,:), intent(in) :: prsi, phii
75 real(kind=kind_phys), dimension(:,:,:), intent(in) :: q1
76
77 real(kind=kind_phys), dimension(:,:), intent(inout) :: du, dv,
78 & tau
79 real(kind=kind_phys), dimension(:,:,:), intent(inout) :: rtg
80
81 real(kind=kind_phys), dimension(:,:,:), intent(inout), optional ::&
82 & dtend
83 integer, dimension(:,:), intent(in) :: dtidx
84 integer, intent(in) :: index_of_temperature, index_of_x_wind,
85 & index_of_y_wind, index_of_process_pbl, ntqv
86 logical, intent(in) :: ldiag3d,
87 & gen_tend
88
89 integer, dimension(:), intent(out) :: kpbl
90 real(kind=kind_phys), dimension(:), intent(out) :: dusfc,
91 & dvsfc, dtsfc, dqsfc, hpbl
92 real(kind=kind_phys), dimension(:,:), intent(out) :: prnum
93 real(kind=kind_phys), dimension(:,:), intent(out) :: dkt
94
95 character(len=*), intent(out) :: errmsg
96 integer, intent(out) :: errflg
97!
98! locals
99!
100 integer, parameter :: kp = kind_phys
101 integer i,is,k,kk,km1,kmpbl,kp1, ntloc
102!
103 logical pblflg(im), sfcflg(im), flg(im)
104
105 real(kind=kind_phys), dimension(im) :: phih, phim
106 &, rbdn, rbup, sflux, z0, crb, zol, thermal
107 &, beta, tx1
108!
109 real(kind=kind_phys), dimension(im,km) :: theta, thvx, zl, a1, ad
110 &, dt2odel
111 real(kind=kind_phys), dimension(im,km-1):: xkzo, xkzmo, al, au
112 &, dku, rdzt
113!
114 real(kind=kind_phys) zi(im,km+1), a2(im,km*(ntrac+1))
115!
116 real(kind=kind_phys) dsdz2, dsdzq, dsdzt, dsig, dt2, rdt
117 &, dtodsd, dtodsu, rdz, tem, tem1
118 &, ttend, utend, vtend, qtend
119 &, spdk2, rbint, ri, zol1, robn, bvf2
120!
121 real(kind=kind_phys), parameter :: one=1.0_kp, zero=0.0_kp
122 &, zolcr=0.2_kp,
123 & zolcru=-0.5_kp, rimin=-100.0_kp, sfcfrac=0.1_kp,
124 & crbcon=0.25_kp, crbmin=0.15_kp, crbmax=0.35_kp,
125 & qmin=1.0e-8_kp, zfmin=1.0d-8, qlmin=1.0e-12_kp,
126 & aphi5=5.0_kp, aphi16=16.0_kp, f0=1.0e-4_kp
127 &, dkmin=zero, dkmax=1000.0_kp
128! &, dkmin=zero, dkmax=1000., xkzminv=0.3
129 &, prmin=0.25_kp, prmax=4.0_kp, vk=0.4_kp,
130 & cfac=6.5_kp
131 real(kind=kind_phys) :: gravi, cont, conq, gocp, go2
132 integer :: idtend
133
134 gravi = one / grav
135 cont = cp * gravi
136 conq = hvap * gravi
137 gocp = grav / cp
138 go2 = grav * 0.5_kp
139
140! Initialize CCPP error handling variables
141 errmsg = ''
142 errflg = 0
143!
144! Set intent(out) variables
145 dkt = zero
146!
147!-----------------------------------------------------------------------
148!
149! compute preliminary variables
150!
151 dt2 = delt
152 rdt = one / dt2
153 km1 = km - 1
154 kmpbl = km / 2
155!
156 rtg = zero
157!
158 do k=1,km
159 do i=1,im
160 zi(i,k) = phii(i,k) * gravi
161 zl(i,k) = phil(i,k) * gravi
162 dt2odel(i,k) = dt2 / del(i,k)
163 enddo
164 enddo
165 do i=1,im
166 zi(i,km+1) = phii(i,km+1) * gravi
167 enddo
168!
169 do k = 1,km1
170 do i=1,im
171 rdzt(i,k) = one / (zl(i,k+1) - zl(i,k))
172 prnum(i,k) = one
173 enddo
174 enddo
175! Setup backgrond diffision
176 do i=1,im
177 prnum(i,km) = one
178 tx1(i) = one / prsi(i,1)
179 enddo
180 do k = 1,km1
181 do i=1,im
182 xkzo(i,k) = zero
183 xkzmo(i,k) = zero
184! if (k < kinver(i)) then
185 if (k <= kinver(i)) then
186! vertical background diffusivity for heat and momentum
187 tem1 = one - prsi(i,k+1) * tx1(i)
188 tem1 = min(one, exp(-tem1 * tem1 * 10.0_kp))
189 xkzo(i,k) = xkzm_h * tem1
190 xkzmo(i,k) = xkzm_m * tem1
191 endif
192 enddo
193 enddo
194!
195! diffusivity in the inversion layer is set to be xkzminv (m^2/s)
196!
197 do k = 1,kmpbl
198 do i=1,im
199 if(zi(i,k+1) > 250.0_kp) then
200 tem1 = (t1(i,k+1)-t1(i,k)) * rdzt(i,k)
201 if(tem1 > 1.0e-5_kp) then
202 xkzo(i,k) = min(xkzo(i,k),xkzminv)
203 endif
204 endif
205 enddo
206 enddo
207!
208!
209 do i = 1,im
210 z0(i) = 0.01_kp * zorl(i)
211 kpbl(i) = 1
212 hpbl(i) = zi(i,1)
213 pblflg(i) = .true.
214 sfcflg(i) = .true.
215 if(rbsoil(i) > zero) sfcflg(i) = .false.
216 dusfc(i) = zero
217 dvsfc(i) = zero
218 dtsfc(i) = zero
219 dqsfc(i) = zero
220 enddo
221!
222 do k = 1,km
223 do i=1,im
224 tx1(i) = zero
225 enddo
226 do kk=1,ncnd
227 do i=1,im
228 tx1(i) = tx1(i) + max(q1(i,k,ntcw+kk-1), qlmin)
229 enddo
230 enddo
231 do i = 1,im
232 theta(i,k) = t1(i,k) * psk(i) / prslk(i,k)
233 thvx(i,k) = theta(i,k)*(one+fv*max(q1(i,k,1),qmin)-tx1(i))
234 enddo
235 enddo
236!
237 do i = 1,im
238 sflux(i) = heat(i) + evap(i)*fv*theta(i,1)
239 if (.not.sfcflg(i) .or. sflux(i) <= zero) pblflg(i)=.false.
240 beta(i) = dt2 / (zi(i,2)-zi(i,1))
241 enddo
242!
243! compute the pbl height
244!
245! write(0,*)' IN moninbl u10=',u10m(1:5),' v10=',v10m(1:5)
246 do i=1,im
247 flg(i) = .false.
248 rbup(i) = rbsoil(i)
249!
250 if (pblflg(i)) then
251 thermal(i) = thvx(i,1)
252 crb(i) = crbcon
253 else
254 thermal(i) = tsea(i)*(one+fv*max(q1(i,1,1),qmin))
255 tem = max(one, sqrt(u10m(i)*u10m(i) + v10m(i)*v10m(i)))
256 robn = tem / (f0 * z0(i))
257 tem1 = 1.0e-7_kp * robn
258 crb(i) = max(min(0.16_kp * (tem1 ** (-0.18_kp)), crbmax),
259 & crbmin)
260 endif
261 enddo
262 do k = 1, kmpbl
263 do i = 1, im
264 if (.not.flg(i)) then
265 rbdn(i) = rbup(i)
266 spdk2 = max((u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)), one)
267 rbup(i) = (thvx(i,k)-thermal(i))*phil(i,k)
268 & / (thvx(i,1)*spdk2)
269 kpbl(i) = k
270 flg(i) = rbup(i) > crb(i)
271 endif
272 enddo
273 enddo
274 do i = 1,im
275 if(kpbl(i) > 1) then
276 k = kpbl(i)
277 if (rbdn(i) >= crb(i)) then
278 rbint = zero
279 elseif(rbup(i) <= crb(i)) then
280 rbint = one
281 else
282 rbint = (crb(i)-rbdn(i)) / (rbup(i)-rbdn(i))
283 endif
284 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
285 if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1
286 else
287 hpbl(i) = zl(i,1)
288 kpbl(i) = 1
289 endif
290 enddo
291!
292! compute similarity parameters
293!
294 do i=1,im
295 zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
296 if (sfcflg(i)) then
297 zol(i) = min(zol(i),-zfmin)
298 else
299 zol(i) = max(zol(i),zfmin)
300 endif
301 zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1)
302 if (sfcflg(i)) then
303! phim(i) = (1.-aphi16*zol1)**(-1./4.)
304! phih(i) = (1.-aphi16*zol1)**(-1./2.)
305 tem = one / max(one - aphi16*zol1, 1.0e-8_kp)
306 phih(i) = sqrt(tem)
307 phim(i) = sqrt(phih(i))
308 else
309 phim(i) = one + aphi5*zol1
310 phih(i) = phim(i)
311 endif
312 enddo
313!
314! enhance the pbl height by considering the thermal excess
315!
316 do i=1,im
317 flg(i) = .true.
318 if (pblflg(i)) then
319 flg(i) = .false.
320 rbup(i) = rbsoil(i)
321 endif
322 enddo
323 do k = 2, kmpbl
324 do i = 1, im
325 if (.not.flg(i)) then
326 rbdn(i) = rbup(i)
327 spdk2 = max((u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)), one)
328 rbup(i) = (thvx(i,k)-thermal(i)) * phil(i,k)
329 & / (thvx(i,1)*spdk2)
330 kpbl(i) = k
331 flg(i) = rbup(i) > crb(i)
332 endif
333 enddo
334 enddo
335 do i = 1,im
336 if (pblflg(i)) then
337 k = kpbl(i)
338 if(rbdn(i) >= crb(i)) then
339 rbint = zero
340 elseif(rbup(i) <= crb(i)) then
341 rbint = one
342 else
343 rbint = (crb(i)-rbdn(i)) / (rbup(i)-rbdn(i))
344 endif
345 if (k > 1) then
346 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
347 if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1
348 if(kpbl(i) <= 1) then
349 pblflg(i) = .false.
350 endif
351 else
352 pblflg(i) = .false.
353 endif
354 endif
355 if (pblflg(i)) then
356 tem = phih(i)/phim(i) + cfac*vk*sfcfrac
357 else
358 tem = phih(i)/phim(i)
359 endif
360 prnum(i,1) = min(prmin,max(prmax,tem))
361 enddo
362!
363 do i = 1, im
364 if(zol(i) > zolcr) then
365 kpbl(i) = 1
366 endif
367 enddo
368!
369! compute Prandtl number above boundary layer
370!
371 do k = 1, km1
372 kp1 = k + 1
373 do i=1,im
374 if(k >= kpbl(i)) then
375 rdz = rdzt(i,k)
376 tem = u1(i,k) - u1(i,kp1)
377 tem1 = v1(i,k) - v1(i,kp1)
378 tem = (tem*tem + tem1*tem1) * rdz * rdz
379 bvf2 = go2*(thvx(i,kp1)-thvx(i,k))*rdz / (t1(i,k)+t1(i,kp1))
380 ri = max(bvf2/tem,rimin)
381 if(ri < zero) then ! unstable regime
382 prnum(i,kp1) = one
383 else
384 prnum(i,kp1) = min(one + 2.1_kp*ri, prmax)
385 endif
386 elseif (k > 1) then
387 prnum(i,kp1) = prnum(i,1)
388 endif
389!
390! prnum(i,kp1) = one
391 prnum(i,kp1) = max(prmin, min(prmax, prnum(i,kp1)))
392 tem = tkh(i,kp1) * prnum(i,kp1)
393 dku(i,k) = max(min(tem+xkzmo(i,k), dkmax), xkzmo(i,k))
394 dkt(i,k) = max(min(tkh(i,kp1)+xkzo(i,k), dkmax), xkzo(i,k))
395 enddo
396 enddo
397!
398! compute tridiagonal matrix elements for heat and moisture
399!
400 do i=1,im
401 ad(i,1) = one
402 a1(i,1) = t1(i,1) + beta(i) * heat(i)
403 a2(i,1) = q1(i,1,1) + beta(i) * evap(i)
404 enddo
405
406 ntloc = 1
407 if(ntrac > 1) then
408 is = 0
409 do k = 2, ntrac
410 if (k /= ntke) then
411 ntloc = ntloc + 1
412 is = is + km
413 do i = 1, im
414 a2(i,1+is) = q1(i,1,k)
415 enddo
416 endif
417 enddo
418 endif
419!
420 do k = 1,km1
421 kp1 = k + 1
422 do i = 1,im
423 dtodsd = dt2odel(i,k)
424 dtodsu = dt2odel(i,kp1)
425 dsig = prsl(i,k)-prsl(i,kp1)
426 rdz = rdzt(i,k)
427 tem1 = dsig * dkt(i,k) * rdz
428 dsdz2 = tem1 * rdz
429 au(i,k) = -dtodsd*dsdz2
430 al(i,k) = -dtodsu*dsdz2
431!
432 ad(i,k) = ad(i,k)-au(i,k)
433 ad(i,kp1) = one - al(i,k)
434 dsdzt = tem1 * gocp
435 a1(i,k) = a1(i,k) + dtodsd*dsdzt
436 a1(i,kp1) = t1(i,kp1) - dtodsu*dsdzt
437 a2(i,kp1) = q1(i,kp1,1)
438!
439 enddo
440 enddo
441!
442 if(ntrac > 1) then
443 is = 0
444 do kk = 2, ntrac
445 if (kk /= ntke) then
446 is = is + km
447 do k = 1, km1
448 kp1 = k + 1
449 do i = 1, im
450 a2(i,kp1+is) = q1(i,kp1,kk)
451 enddo
452 enddo
453 endif
454 enddo
455 endif
456!
457! solve tridiagonal problem for heat, moisture and tracers
458!
459 call tridin(im,km,ntloc,al,ad,au,a1,a2,au,a1,a2)
460
461!
462! recover tendencies of heat and moisture
463!
464 do k = 1,km
465 do i = 1,im
466 ttend = (a1(i,k)-t1(i,k)) * rdt
467 qtend = (a2(i,k)-q1(i,k,1)) * rdt
468 tau(i,k) = tau(i,k) + ttend
469 rtg(i,k,1) = rtg(i,k,1) + qtend
470 dtsfc(i) = dtsfc(i) + del(i,k)*ttend
471 dqsfc(i) = dqsfc(i) + del(i,k)*qtend
472 enddo
473 enddo
474 if(ldiag3d .and. .not. gen_tend) then
475 idtend = dtidx(index_of_temperature,index_of_process_pbl)
476 if(idtend>=1) then
477 dtend(:,:,idtend) = dtend(:,:,idtend) + (a1-t1)
478 endif
479 idtend = dtidx(ntqv+100,index_of_process_pbl)
480 if(idtend>=1) then
481 dtend(:,:,idtend) = dtend(:,:,idtend) + a2-q1(:,:,1)
482 endif
483 endif
484 do i = 1,im
485 dtsfc(i) = dtsfc(i) * cont
486 dqsfc(i) = dqsfc(i) * conq
487 enddo
488 if(ntrac > 1) then
489 is = 0
490 do kk = 2, ntrac
491 if (kk /= ntke) then
492 is = is + km
493 do k = 1, km
494 do i = 1, im
495 qtend = (a2(i,k+is)-q1(i,k,kk))*rdt
496 rtg(i,k,kk) = rtg(i,k,kk) + qtend
497 enddo
498 enddo
499 endif
500 enddo
501 if(ldiag3d .and. ntoz>0 .and. .not. gen_tend) then
502 idtend=dtidx(100+ntoz,index_of_process_pbl)
503 if(idtend>=1) then
504 kk = ntoz
505 is = (kk-1) * km
506 do k = 1, km
507 do i = 1, im
508 qtend = (a2(i,k+is)-q1(i,k,kk))
509 dtend(i,k,idtend) = dtend(i,k,idtend) + qtend
510 enddo
511 enddo
512 endif
513 endif
514 endif
515!
516! compute tridiagonal matrix elements for momentum
517!
518 do i=1,im
519 ad(i,1) = one + beta(i) * stress(i) / spd1(i)
520 a1(i,1) = u1(i,1)
521 a2(i,1) = v1(i,1)
522 enddo
523!
524 do k = 1,km1
525 kp1 = k + 1
526 do i=1,im
527 dtodsd = dt2odel(i,k)
528 dtodsu = dt2odel(i,kp1)
529 dsig = prsl(i,k)-prsl(i,kp1)
530 rdz = rdzt(i,k)
531 tem1 = dsig*dku(i,k)*rdz
532 dsdz2 = tem1 * rdz
533 au(i,k) = -dtodsd*dsdz2
534 al(i,k) = -dtodsu*dsdz2
535!
536 ad(i,k) = ad(i,k) - au(i,k)
537 ad(i,kp1) = one - al(i,k)
538 a1(i,kp1) = u1(i,kp1)
539 a2(i,kp1) = v1(i,kp1)
540!
541 enddo
542 enddo
543
544 call tridi2(im,km,al,ad,au,a1,a2,au,a1,a2)
545!
546! recover tendencies of momentum
547!
548 do k = 1,km
549 do i = 1,im
550 utend = (a1(i,k)-u1(i,k))*rdt
551 vtend = (a2(i,k)-v1(i,k))*rdt
552 du(i,k) = du(i,k) + utend
553 dv(i,k) = dv(i,k) + vtend
554 tem = del(i,k) * gravi
555 dusfc(i) = dusfc(i) + tem * utend
556 dvsfc(i) = dvsfc(i) + tem * vtend
557 enddo
558 enddo
559 if (ldiag3d .and. .not. gen_tend) then
560 idtend = dtidx(index_of_x_wind,index_of_process_pbl)
561 if(idtend>=1) then
562 dtend(:,:,idtend) = dtend(:,:,idtend) + a1-u1
563 endif
564 idtend = dtidx(index_of_y_wind,index_of_process_pbl)
565 if(idtend>=1) then
566 dtend(:,:,idtend) = dtend(:,:,idtend) + a1-v1
567 endif
568 endif
569!
570 if (ntke > 0) then ! solve tridiagonal problem for momentum and tke
571!
572! compute tridiagonal matrix elements for tke
573!
574 do i=1,im
575 ad(i,1) = one
576 a1(i,1) = q1(i,1,ntke)
577 enddo
578!
579 do k = 1,km1
580 kp1 = k + 1
581 do i=1,im
582 dtodsd = dt2odel(i,k)
583 dtodsu = dt2odel(i,kp1)
584 dsig = prsl(i,k)-prsl(i,kp1)
585 rdz = rdzt(i,k)
586 tem1 = dsig*dku(i,k)*(rdz+rdz)
587 dsdz2 = tem1 * rdz
588 au(i,k) = -dtodsd*dsdz2
589 al(i,k) = -dtodsu*dsdz2
590!
591 ad(i,k) = ad(i,k) - au(i,k)
592 ad(i,kp1) = one - al(i,k)
593 a1(i,kp1) = q1(i,kp1,ntke)
594 enddo
595 enddo
596
597 call tridi1(im,km,al,ad,au,a1,au,a1)
598!
599 do k = 1, km ! recover tendencies of tke
600 do i = 1, im
601 qtend = (a1(i,k)-q1(i,k,ntke))*rdt
602 rtg(i,k,ntke) = rtg(i,k,ntke) + qtend
603 enddo
604 enddo
605 endif
606!
607 return
608 end subroutine moninshoc_run
609
610 end module moninshoc
subroutine tridi1(l, n, cl, cm, cu, r1, au, a1)
Routine to solve the tridiagonal system to calculate temperature and moisture at .
Definition tridi.f:11
subroutine tridin(l, n, nt, cl, cm, cu, r1, r2, au, a1, a2)
Routine to solve the tridiagonal system to calculate u- and v-momentum at .
Definition tridi.f:95
subroutine tridi2(l, n, cl, cm, cu, r1, r2, au, a1, a2)
This subroutine ..
Definition tridi.f:50
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module contains the CCPP-compliant SHOC scheme.
Definition moninshoc.f:5