34 SUBROUTINE myjpbl_wrapper_run( &
38 ntcw,ntiw,ntrw,ntsw,ntgl, &
39 ugrs, vgrs, tgrs, qgrs, &
40 prsl, prsi, phii, hprime1, &
41 prsik_1, prslk_1, prslki, tsfc, qsfc, &
42 phy_myj_qsfc, phy_myj_thz0, phy_myj_qz0, &
43 phy_myj_uz0, phy_myj_vz0, phy_myj_z0base, &
44 phy_myj_akhs, phy_myj_akms, &
45 phy_myj_chkqlm, phy_myj_elflx, &
46 phy_myj_a1u, phy_myj_a1t, phy_myj_a1q, &
47 pblh, kpbl, kinver, slmsk, &
48 garea, ustar, cm, ch, wind, &
49 snowd, zorl, evap, hflx, &
50 dudt, dvdt, dtdt, dqdt, &
51 dusfc,dvsfc,dtsfc,dqsfc, &
52 dkt,xkzm_m, xkzm_h,xkzm_s, gamt,gamq, &
53 con_cp,con_g,con_rd, &
54 me, lprnt, gen_tend, ldiag3d, dtend, dtidx, &
55 index_of_temperature, index_of_x_wind, &
56 index_of_y_wind, index_of_process_pbl, &
57 ntqv, errmsg, errflg )
85 real,
parameter :: xkgdx=25000.,xkzinv=0.15
89 character(len=*),
intent(out) :: errmsg
90 integer,
intent(out) :: errflg
92 real(kind=kind_phys),
intent(inout),
optional :: dtend(:,:,:)
93 integer,
intent(in) :: dtidx(:,:)
94 integer,
intent(in) :: index_of_temperature, index_of_x_wind, &
95 index_of_y_wind, index_of_process_pbl, ntqv
98 integer,
intent(in) :: im, levs
99 integer,
intent(in) :: kdt, me
100 integer,
intent(in) :: ntrac,ntke,ntcw,ntiw,ntrw,ntsw,ntgl
101 logical,
intent(in) :: restart,do_myjsfc,lprnt,ldiag3d,gen_tend
102 real(kind=kind_phys),
intent(in) :: con_cp, con_g, con_rd
103 real(kind=kind_phys),
intent(in) :: dt_phs, xkzm_m, xkzm_h, xkzm_s
106 real(kind=kind_phys),
dimension(:),
intent(in) :: &
107 prsik_1, prslk_1, prslki, slmsk, garea, &
108 snowd, evap, hflx, cm, ch, wind, hprime1
109 real(kind=kind_phys),
dimension(:),
intent(inout) :: &
110 zorl, ustar, tsfc, qsfc
111 real(kind=kind_phys),
dimension(:),
intent(inout),
optional :: &
112 phy_myj_thz0, phy_myj_z0base, phy_myj_chkqlm, &
113 phy_myj_akhs, phy_myj_akms, phy_myj_qz0, &
114 phy_myj_qsfc, phy_myj_elflx, phy_myj_a1u, &
115 phy_myj_a1t, phy_myj_a1q, phy_myj_uz0, &
117 real(kind=kind_phys),
dimension(:),
intent(out) :: &
118 pblh,dusfc,dvsfc,dtsfc,dqsfc,gamt,gamq
119 integer,
dimension(:),
intent(out) :: kpbl
120 integer,
dimension(:),
intent(in) :: kinver
123 real(kind=kind_phys),
dimension(:,:),
intent(in) :: &
125 real(kind=kind_phys),
dimension(:,:),
intent(in) :: &
126 ugrs, vgrs, tgrs, prsl
129 real(kind=kind_phys),
dimension(:,:),
intent(inout) :: &
131 real(kind=kind_phys),
dimension(:,:),
intent(out) :: &
135 real(kind=kind_phys),
dimension(:,:,:),
intent(inout) :: &
139 integer :: ntsd, k, k1, i, kx1
140 integer :: i_min, i_max, k_min, k_max
142 logical :: lprnt1,lprnt2
143 integer :: ict, ide, lm, me1
144 real(kind=kfpt) :: dt_myj, tem, tem1, tem2, ptem
145 integer,
dimension(im) :: kpbl_myj
146 real(kind=kfpt),
dimension(1:levs-1):: epsl
147 real(kind=kfpt),
dimension(1:levs):: epsq2
148 real(kind=kfpt),
dimension(im) :: &
149 xland, sice, snowd1, ht, stdh, tsk, &
150 ustar1,z0,pblh_myj, &
152 real(kind=kfpt),
dimension(im,levs) :: &
153 u_myj, v_myj, t_myj, q_myj, th_myj, &
154 cw, dz_myj, pmid, q2, exner, del
155 real(kind=kfpt),
dimension(im,levs+1) :: pint
156 real(kind=kfpt),
dimension(im,levs) :: &
157 rublten,rvblten,rthblten,rqvblten,rqcblten
158 real(kind=kfpt),
dimension(im,levs) :: el_myj
159 real(kind=kfpt),
dimension(im) :: &
160 dusfc1,dvsfc1,dtsfc1,dqsfc1
161 real(kind=kfpt),
dimension(im) :: thlm,qlm
162 real(kind=kfpt),
dimension(im,13) :: phy_f2d_myj
163 real(kind=kfpt),
dimension(im,levs) :: xcofh &
165 real(kind=kind_phys) :: g, r_d, g_inv, cappa
166 real(kind=kind_phys) :: thz0, qz0, a1u, a1t, a1q
167 real(kind=kind_phys) :: z0m, aa1u, aa1t, z1uov, z1tox
168 real(kind=kind_phys) :: tmax,tmin,t_myj1
169 real(kind=kind_phys),
dimension(im) :: &
171 sm,work3,wind1,work4 &
172 ,rho,qfc1,gdx,xkzm_hx,xkzm_mx,tx1, tx2
175 real(kind=kind_phys),
dimension(im,levs) :: dkt2
176 integer :: uidx, vidx, tidx, qidx
194 if(me.eq.0)print*,
'ntsd=', ntsd
202 cappa = con_rd/con_cp
205 work3(i)=prsik_1(i) / prslk_1(i)
207 if(sice(i) < 0.7)sice(i)=0
208 sm(i)=1.;
if(slmsk(i) > 0.5 ) sm(i)=0.
211 sfcz(i)=phii(i,1)*g_inv
212 work4(i)=(1.e5/prsi(i,1))**cappa
213 thsfc(i)=tsfc(i)*work4(i)
219 u_myj(i,k)=ugrs(i,k1)
220 v_myj(i,k)=vgrs(i,k1)
221 t_myj(i,k)=tgrs(i,k1)
222 q_myj(i,k)=qgrs(i,k1,1)
223 cw(i,k) =qgrs(i,k1,ntcw)
229 q2(i,k) =max(0.02,qgrs(i,k1,ntke)*2.)
234 pmid(i,k) =prsl(i,k1)
235 exner(i,k)=(prsl(i,k1)*1.e-5)**cappa
236 th_myj(i,k)=tgrs(i,k1)/exner(i,k)
242 pint(i,k) =prsi(i,k1)
247 gdx(i) = sqrt(garea(i))
252 tx1(i) = 1.0 / prsi(i,1)
254 if(gdx(i) >= xkgdx)
then
258 tem = 1. / (xkgdx - 5.)
259 tem1 = (xkzm_h - 0.01) * tem
260 tem2 = (xkzm_m - 0.01) * tem
262 xkzm_hx(i) = 0.01 + tem1 * ptem
263 xkzm_mx(i) = 0.01 + tem2 * ptem
270 if (k < kinver(i))
then
272 ptem = prsi(i,k+1) * tx1(i)
274 tem1 = tem1 * tem1 * 10.0
275 xkzo(i,k) = xkzm_hx(i) * min(1.0, exp(-tem1))
276 xkzo(i,k) = min(xkzo(i,k),xkzinv)
278 if (ptem >= xkzm_s)
then
279 xkzmo(i,k) = xkzm_mx(i)
282 if (k == kx1 .and. k > 1) tx2(i) = 1.0 / prsi(i,k)
283 tem1 = 1.0 - prsi(i,k+1) * tx2(i)
284 tem1 = tem1 * tem1 * 5.0
285 xkzmo(i,k) = xkzm_mx(i) * min(1.0, exp(-tem1))
286 xkzmo(i,k) = min(xkzmo(i,k),xkzinv)
296 xcofh(i,k1)=xkzo(i,k)
297 el_myj(i,k1)=xkzmo(i,k)
304 xkzmo(i,k)=el_myj(i,k)
310 epsl(k)=sqrt(epsq2(k)*0.5)
315 epsq2(levs)=epsq2(levs-1)
320 del(i,k) = prsi(i,k1) - prsi(i,k1+1)
321 dz_myj(i,k) = (phii(i,k1+1)-phii(i,k1)) * g_inv
326 wind1(i)=max(wind(i),1.0)
329 if(.not.do_myjsfc)
then
331 if(sm(i).gt.0.5.and.sice(i).le.0.5)
then
332 z0m=max(0.018*g_inv*ustar(i)*ustar(i),1.59e-5)
333 z1uov=0.35*30.*sqrt(sqrt(z0m*ustar(i)/1.5e-5))/ustar(i)
334 aa1u=cm(i)*wind1(i)*z1uov
337 aa1t=ch(i)*wind1(i)*z1tox
353 phy_myj_akhs(i) = ch(i)*wind1(i)*(1.+a1t)
354 phy_myj_akms(i) = cm(i)*wind1(i)*(1.+a1u)
355 phy_myj_uz0(i) = u_myj(i,levs)*a1u/(1.+a1u)
356 phy_myj_vz0(i) = v_myj(i,levs)*a1u/(1.+a1u)
357 phy_myj_z0base(i)= z0m
361 qz0=max(evap(i)/phy_myj_akhs(i)+q_myj(i,levs),1.e-9)
362 thz0=hflx(i)/phy_myj_akhs(i)+th_myj(i,levs)
370 phy_myj_thz0(i) = thz0
373 if(cw(i,levs).gt.1.e-9)
then
374 phy_myj_chkqlm(i)= 0.
376 phy_myj_chkqlm(i)= 1.
383 phy_myj_akhs(i)=phy_myj_akhs(i)*wind1(i)/wind(i)
384 phy_myj_akms(i)=phy_myj_akms(i)*wind1(i)/wind(i)
390 phy_myj_elflx(i) = evap(i)
391 qsfc(i)=max(evap(i)/(ch(i)*wind1(i))+q_myj(i,levs),1.e-9)
392 tsfc1(i)=(hflx(i)/(ch(i)*wind1(i))+th_myj(i,levs))/work4(i)
393 phy_myj_qsfc(i) = qsfc(i)
394 thz0 = phy_myj_thz0(i)
395 thlm(i)=th_myj(i,levs)
400 phy_myj_thz0(i)=0.5*(thz0+ &
401 hflx(i)/phy_myj_akhs(i)+th_myj(i,levs))
405 phy_myj_qz0(i) = 0.5*(qz0+ &
406 max(evap(i)/phy_myj_akhs(i)+q_myj(i,levs),1.e-9))
421 if (me.eq.0.and.ntsd.lt.2)
then
422 print*,
'Qingfu test starts PBL'
423 print*,
'ntsd,me,im,levs,ict=',ntsd,me,im,levs,ict
424 print*,
'dt_phs,sfcz,dz_myj=',dt_phs,sfcz(1),dz_myj(1,5)
425 print*,
'pmid,pint,th_myj=',pmid(1,5),pint(1,5),th_myj(1,5)
426 print*,
't_myj,exner,q_myj=',t_myj(1,5),exner(1,5),q_myj(1,5)
427 print*,
'cw,u_myj,v_myj=',cw(1,5),u_myj(1,5),v_myj(1,5)
428 print*,
'tsfc,xland,sice,snowd=',tsfc(1),xland(1),sice(1),snowd(1)
429 print*,
'ustar,z0,pblh,kpbl=',ustar(1),z0(1),pblh(1),kpbl(1)
430 print*,
'q2,xcofh=',q2(1,5),xcofh(1,5)
434 print*,
'thlm,thsfc=',thlm(i),thsfc(i)
439 if(t_myj(i,k).gt.390..or.t_myj(i,k).lt.110.)
then
440 print*,
'Qingfu test starts PBL',i,k,t_myj(i,k)
441 print*,
'ntsd,me,im,levs,ict=',ntsd,me,im,levs,ict
442 print*,
'dt_phs,sfcz,dz_myj=',dt_phs,sfcz(i),dz_myj(i,k)
443 print*,
'pmid,pint,th_myj=',pmid(i,k),pint(i,k),th_myj(i,k)
444 print*,
't_myj,exner,q_myj=',t_myj(i,k),exner(i,k),q_myj(i,k)
445 print*,
'cw,u_myj,v_myj=',cw(i,k),u_myj(i,k),v_myj(i,k)
446 print*,
'tsfc,xland,sice,snowd=',tsfc(i),xland(i),sice(i),snowd(i)
447 print*,
'ustar,z0,pblh,kpbl=',ustar(i),z0(i),pblh(i),kpbl(i)
448 print*,
'q2,xcofh=',q2(i,k),xcofh(i,k)
458 if(tmax.lt.t_myj(i,k1))
then
463 if(tmin.gt.t_myj(i,k1))
then
486 ht(i)=phii(i,1)*g_inv
490 phy_f2d_myj(i,1) = phy_myj_qsfc(i)
491 phy_f2d_myj(i,2) = phy_myj_thz0(i)
492 phy_f2d_myj(i,3) = phy_myj_qz0(i)
493 phy_f2d_myj(i,4) = phy_myj_uz0(i)
494 phy_f2d_myj(i,5) = phy_myj_vz0(i)
495 phy_f2d_myj(i,6) = phy_myj_z0base(i)
496 phy_f2d_myj(i,7) = phy_myj_akhs(i)
497 phy_f2d_myj(i,8) = phy_myj_akms(i)
498 phy_f2d_myj(i,9) = phy_myj_chkqlm(i)
499 phy_f2d_myj(i,10) = phy_myj_elflx(i)
500 phy_f2d_myj(i,11) = phy_myj_a1u(i)
501 phy_f2d_myj(i,12) = phy_myj_a1t(i)
502 phy_f2d_myj(i,13) = phy_myj_a1q(i)
522 if(ntsd.eq.0.or.restart)
then
523 if(.not.restart) xcofh=0.
530 call myjpbl(ntsd,me,dt_myj,epsl,epsq2,ht,stdh,dz_myj,del &
531 ,pmid,pint,th_myj,t_myj,exner,q_myj,cw,u_myj,v_myj &
532 ,tsk,phy_f2d_myj(1:im,1),phy_f2d_myj(1:im,9) &
533 ,phy_f2d_myj(1:im,2),phy_f2d_myj(1:im,3) &
534 ,phy_f2d_myj(1:im,4),phy_f2d_myj(1:im,5) &
536 ,q2,xcofh,ustar1,z0,el_myj,pblh_myj,kpbl_myj,ct &
537 ,phy_f2d_myj(1:im,7),phy_f2d_myj(1:im,8) &
538 ,phy_f2d_myj(1:im,10),mixht,thlm,qlm &
539 ,rublten,rvblten,rthblten,rqvblten,rqcblten &
540 ,dusfc1,dvsfc1,dtsfc1,dqsfc1,xkzo,xkzmo,ict &
552 kpbl(i)=levs-kpbl_myj(i)
554 phy_myj_qsfc(i) = phy_f2d_myj(i,1)
555 phy_myj_thz0(i) = phy_f2d_myj(i,2)
556 phy_myj_qz0(i) = phy_f2d_myj(i,3)
557 phy_myj_uz0(i) = phy_f2d_myj(i,4)
558 phy_myj_vz0(i) = phy_f2d_myj(i,5)
559 phy_myj_z0base(i) = phy_f2d_myj(i,6)
560 phy_myj_akhs(i) = phy_f2d_myj(i,7)
561 phy_myj_akms(i) = phy_f2d_myj(i,8)
562 phy_myj_chkqlm(i) = phy_f2d_myj(i,9)
563 phy_myj_elflx(i) = phy_f2d_myj(i,10)
564 phy_myj_a1u(i) = phy_f2d_myj(i,11)
565 phy_myj_a1t(i) = phy_f2d_myj(i,12)
566 phy_myj_a1q(i) = phy_f2d_myj(i,13)
583 qgrs(:,k,ntke)=q2(:,k1)*0.5
592 dudt(i,k)=dudt(i,k)+rublten(i,k1)
593 dvdt(i,k)=dvdt(i,k)+rvblten(i,k1)
594 dtdt(i,k)=dtdt(i,k)+rthblten(i,k1)*exner(i,k1)
595 dqdt(i,k,1)=dqdt(i,k,1)+rqvblten(i,k1)
596 dqdt(i,k,ntcw)=dqdt(i,k,ntcw)+rqcblten(i,k1)
599 if (ldiag3d .and. .not. gen_tend)
then
600 uidx = dtidx(index_of_x_wind,index_of_process_pbl)
601 vidx = dtidx(index_of_y_wind,index_of_process_pbl)
602 tidx = dtidx(index_of_temperature,index_of_process_pbl)
603 qidx = dtidx(ntqv+100,index_of_process_pbl)
608 if(uidx>=1) dtend(:,k,uidx)=dtend(:,k,uidx)+rublten(:,k1)*dt_phs
609 if(vidx>=1) dtend(:,k,vidx)=dtend(:,k,vidx)+rvblten(:,k1)*dt_phs
610 if(tidx>=1) dtend(:,k,tidx)=dtend(:,k,tidx)+rthblten(:,k1)*exner(:,k1)*dt_phs
611 if(qidx>=1) dtend(:,k,qidx)=dtend(:,k,qidx)+rqvblten(:,k1)*dt_phs
618 if(tsfc(i).gt.350.)
then
619 print*,
'21tsfc,tsfc1,hflx=',tsfc(i),tsfc1(i),hflx(i)
620 print*,
'21qsfc,evap=',qsfc(i),evap(i)
629 t_myj1=t_myj(i,k1)+dtdt(i,k)*dt_phs
630 if(tmax.lt.t_myj1)
then
636 if(tmin.gt.t_myj1)
then
646 if(me.eq.me1.and.tmin.lt.113.6.or.tmax.gt.350.)
then
648 print*,
'bad bad tmin,tmax=',tmin,tmax,i_min,k_min,i_max,k_max
651 print*,
'delt,t_myj=',k,dtdt(i,k)*dt_phs,tgrs(i,k)
654 print*,
'ide,levs,ntsd=',ide,lm,ntsd,dt_myj
655 print*,
'epsl,epsq2,ht,stdh,xland,sice,snowd1=', &
656 epsl(i),epsq2(i),ht(i),stdh(i),xland(i),sice(i),snowd1(i)
657 print*,
'phy_f2d_myj=', &
658 (phy_f2d_myj(i,k),k=1,13)
659 print*,
'tsk(i),ustar1,z0,pblh_myj,kpbl_myj=', &
660 tsk(i),ustar1(i),z0(i),pblh_myj(i),kpbl_myj(i)
661 print*,
'mixht=',mixht(i)
663 print*,
'u,v,t=',k,u_myj(i,k),v_myj(i,k), &
667 print*,
'q,th,dz_myj=',k,q_myj(i,k),th_myj(i,k),dz_myj(i,k)
670 print*,
'del,pmid,pint,=',k,del(i,k), &
671 pmid(i,k),pint(i,k+1)
674 print*,
'exner,cw,q2=',k,exner(i,k),cw(i,k), &
678 print*,
'xcofh,el_myj,dkt=',k,xcofh(i,k),el_myj(i,k),dkt(i,k)
692 t_myj1=t_myj(i,k1)+dtdt(i,k)*dt_phs
693 if(tmax.lt.t_myj1)
then
699 if(tmin.gt.t_myj1)
then
706 print*,
'2after me i_min,k_min,i_max,k_max=',me,i_min,k_min,i_max,k_max
707 print*,
'ntsd,tmin,tmax=',ntsd,tmin,tmax
708 print*,
'dtdt(i,j)=',dtdt(i_max,k_max)*dt_phs,t_myj(i_max,k_max)
716 t_myj1=ugrs(i,k)+dudt(i,k)*dt_phs
718 if(tmax.lt.t_myj1)
then
723 if(tmin.gt.t_myj1)
then
730 print*,
'3after i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
731 print*,
'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
732 print*,
'dudt(i,k)=',dudt(i_max,k_max)*dt_phs,ugrs(i_max,k_max)
734 if(tmax.gt.200.or.tmin.lt.-200)
then
735 print*,
'bad,bad,bad=',dudt(i_max,k_max)*dt_phs,ugrs(i_max,k_max)
737 print*,
'k,dudt*dt_phs,ugrs=',k,dudt(i_max,k)*dt_phs,ugrs(i_max,k)
747 t_myj1=vgrs(i,k)+dvdt(i,k)*dt_phs
749 if(tmax.lt.t_myj1)
then
754 if(tmin.gt.t_myj1)
then
761 print*,
'4after i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
762 print*,
'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
763 print*,
'dvdt(i,k)=',dvdt(i_max,k_max)*dt_phs,vgrs(i_max,k_max)
771 t_myj1=q_myj(i,k1)+dqdt(i,k,1)*dt_phs
773 if(tmax.lt.t_myj1)
then
778 if(tmin.gt.t_myj1)
then
785 print*,
'5after i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
786 print*,
'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
787 print*,
'dqdt(i,k)=',dqdt(i_max,k_max,1)*dt_phs,qgrs(i_max,k_max,1)
795 t_myj1=cw(i,k1)+dqdt(i,k,ntcw)*dt_phs
797 if(tmax.lt.t_myj1)
then
802 if(tmin.gt.t_myj1)
then
809 print*,
'6after i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
810 print*,
'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
811 print*,
'dqdt(i,k,ntcw)=',dqdt(i_max,k_max,ntcw)*dt_phs,qgrs(i_max,k_max,ntcw)