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)
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,
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
77 real(kind=kind_phys),
dimension(:,:),
intent(inout) :: du, dv,
79 real(kind=kind_phys),
dimension(:,:,:),
intent(inout) :: rtg
81 real(kind=kind_phys),
dimension(:,:,:),
intent(inout),
optional ::&
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,
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
95 character(len=*),
intent(out) :: errmsg
96 integer,
intent(out) :: errflg
100 integer,
parameter :: kp = kind_phys
101 integer i,is,k,kk,km1,kmpbl,kp1, ntloc
103 logical pblflg(im), sfcflg(im), flg(im)
105 real(kind=kind_phys),
dimension(im) :: phih, phim
106 &, rbdn, rbup, sflux, z0, crb, zol, thermal
109 real(kind=kind_phys),
dimension(im,km) :: theta, thvx, zl, a1, ad
111 real(kind=kind_phys),
dimension(im,km-1):: xkzo, xkzmo, al, au
114 real(kind=kind_phys) zi(im,km+1), a2(im,km*(ntrac+1))
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
121 real(kind=kind_phys),
parameter :: one=1.0_kp, zero=0.0_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
129 &, prmin=0.25_kp, prmax=4.0_kp, vk=0.4_kp,
131 real(kind=kind_phys) :: gravi, cont, conq, gocp, go2
160 zi(i,k) = phii(i,k) * gravi
161 zl(i,k) = phil(i,k) * gravi
162 dt2odel(i,k) = dt2 / del(i,k)
166 zi(i,km+1) = phii(i,km+1) * gravi
171 rdzt(i,k) = one / (zl(i,k+1) - zl(i,k))
178 tx1(i) = one / prsi(i,1)
185 if (k <= kinver(i))
then
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
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)
210 z0(i) = 0.01_kp * zorl(i)
215 if(rbsoil(i) > zero) sfcflg(i) = .false.
228 tx1(i) = tx1(i) + max(q1(i,k,ntcw+kk-1), qlmin)
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))
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))
251 thermal(i) = thvx(i,1)
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),
264 if (.not.flg(i))
then
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)
270 flg(i) = rbup(i) > crb(i)
277 if (rbdn(i) >= crb(i))
then
279 elseif(rbup(i) <= crb(i))
then
282 rbint = (crb(i)-rbdn(i)) / (rbup(i)-rbdn(i))
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
295 zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
297 zol(i) = min(zol(i),-zfmin)
299 zol(i) = max(zol(i),zfmin)
301 zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1)
305 tem = one / max(one - aphi16*zol1, 1.0e-8_kp)
307 phim(i) = sqrt(phih(i))
309 phim(i) = one + aphi5*zol1
325 if (.not.flg(i))
then
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)
331 flg(i) = rbup(i) > crb(i)
338 if(rbdn(i) >= crb(i))
then
340 elseif(rbup(i) <= crb(i))
then
343 rbint = (crb(i)-rbdn(i)) / (rbup(i)-rbdn(i))
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
356 tem = phih(i)/phim(i) + cfac*vk*sfcfrac
358 tem = phih(i)/phim(i)
360 prnum(i,1) = min(prmin,max(prmax,tem))
364 if(zol(i) > zolcr)
then
374 if(k >= kpbl(i))
then
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)
384 prnum(i,kp1) = min(one + 2.1_kp*ri, prmax)
387 prnum(i,kp1) = prnum(i,1)
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))
402 a1(i,1) = t1(i,1) + beta(i) * heat(i)
403 a2(i,1) = q1(i,1,1) + beta(i) * evap(i)
414 a2(i,1+is) = q1(i,1,k)
423 dtodsd = dt2odel(i,k)
424 dtodsu = dt2odel(i,kp1)
425 dsig = prsl(i,k)-prsl(i,kp1)
427 tem1 = dsig * dkt(i,k) * rdz
429 au(i,k) = -dtodsd*dsdz2
430 al(i,k) = -dtodsu*dsdz2
432 ad(i,k) = ad(i,k)-au(i,k)
433 ad(i,kp1) = one - al(i,k)
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)
450 a2(i,kp1+is) = q1(i,kp1,kk)
459 call tridin(im,km,ntloc,al,ad,au,a1,a2,au,a1,a2)
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
474 if(ldiag3d .and. .not. gen_tend)
then
475 idtend = dtidx(index_of_temperature,index_of_process_pbl)
477 dtend(:,:,idtend) = dtend(:,:,idtend) + (a1-t1)
479 idtend = dtidx(ntqv+100,index_of_process_pbl)
481 dtend(:,:,idtend) = dtend(:,:,idtend) + a2-q1(:,:,1)
485 dtsfc(i) = dtsfc(i) * cont
486 dqsfc(i) = dqsfc(i) * conq
495 qtend = (a2(i,k+is)-q1(i,k,kk))*rdt
496 rtg(i,k,kk) = rtg(i,k,kk) + qtend
501 if(ldiag3d .and. ntoz>0 .and. .not. gen_tend)
then
502 idtend=dtidx(100+ntoz,index_of_process_pbl)
508 qtend = (a2(i,k+is)-q1(i,k,kk))
509 dtend(i,k,idtend) = dtend(i,k,idtend) + qtend
519 ad(i,1) = one + beta(i) * stress(i) / spd1(i)
527 dtodsd = dt2odel(i,k)
528 dtodsu = dt2odel(i,kp1)
529 dsig = prsl(i,k)-prsl(i,kp1)
531 tem1 = dsig*dku(i,k)*rdz
533 au(i,k) = -dtodsd*dsdz2
534 al(i,k) = -dtodsu*dsdz2
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)
544 call tridi2(im,km,al,ad,au,a1,a2,au,a1,a2)
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
559 if (ldiag3d .and. .not. gen_tend)
then
560 idtend = dtidx(index_of_x_wind,index_of_process_pbl)
562 dtend(:,:,idtend) = dtend(:,:,idtend) + a1-u1
564 idtend = dtidx(index_of_y_wind,index_of_process_pbl)
566 dtend(:,:,idtend) = dtend(:,:,idtend) + a1-v1
576 a1(i,1) = q1(i,1,ntke)
582 dtodsd = dt2odel(i,k)
583 dtodsu = dt2odel(i,kp1)
584 dsig = prsl(i,k)-prsl(i,kp1)
586 tem1 = dsig*dku(i,k)*(rdz+rdz)
588 au(i,k) = -dtodsd*dsdz2
589 al(i,k) = -dtodsu*dsdz2
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)
597 call tridi1(im,km,al,ad,au,a1,au,a1)
601 qtend = (a1(i,k)-q1(i,k,ntke))*rdt
602 rtg(i,k,ntke) = rtg(i,k,ntke) + qtend