255 cp , cpv , cliq , cice , &
256 p608 , ep_2 , ep_3 , gtr , &
257 grav , g_inv , karman , p1000mb , &
258 rcp , r_d , r_v , rk , &
259 rvovrd , svp1 , svp2 , svp3 , &
260 xlf , xlv , xls , xlscp , &
261 xlvcp , tv0 , tv1 , tref , &
262 zero , half , one , two , &
263 onethird , twothirds , tkmin , t0c , &
272 real(kind_phys),
parameter :: cphm_st=5.0, cphm_unst=16.0, &
273 cphh_st=5.0, cphh_unst=16.0
276 real(kind_phys),
parameter :: &
285 &a1 = b1*( 1.0-3.0*g1 )/6.0, &
287 &c1 = g1 -1.0/( 3.0*a1*2.88449914061481660), &
288 &a2 = a1*( g1-c1 )/( g1*pr ), &
289 &g2 = b2/b1*( 1.0-c3 ) +2.0*a1/b1*( 3.0-2.0*c2 )
291 real(kind_phys),
parameter :: &
294 &e1c = 3.0*a2*b2*cc3, &
295 &e2c = 9.0*a1*a2*cc2, &
296 &e3c = 9.0*a2*a2*cc2*( 1.0-c5 ), &
297 &e4c = 12.0*a1*a2*cc2, &
302 real(kind_phys),
parameter :: qmin=0.0, zmax=1.0, sqfac=3.0
306 real(kind_phys),
parameter :: qkemin=1.e-4
307 real(kind_phys),
parameter :: tliq = 269.
310 real(kind_phys),
parameter :: rr2=0.7071068, rrp=0.3989423
320 real(kind_phys),
parameter :: ckmod=1.
325 real(kind_phys),
parameter :: scaleaware=1.
329 integer,
parameter :: bl_mynn_topdown = 0
331 integer,
parameter :: bl_mynn_edmf_dd = 0
334 integer,
parameter :: dheat_opt = 1
337 logical,
parameter :: env_subs = .false.
341 integer,
parameter :: bl_mynn_stfunc = 1
344 logical,
parameter :: debug_code = .false.
345 integer,
parameter :: idbg = 23
348 integer :: mynn_level
362 &initflag,restart,cycling, &
364 &u,v,w,th,sqv3d,sqc3d,sqi3d, &
366 &qnwfa,qnifa,qnbca,ozone, &
369 &ust,ch,hfx,qfx,rmol,wspd, &
370 &uoce,voce, & !ocean current
373 &nchem,kdvel,ndvel, & !Smoke/Chem variables
374 &chem3d,vdep,smoke_dbg, &
375 &frp,emis_ant_no, & ! JLS/RAR to adjust exchange coeffs
376 &mix_chem,enh_mix,rrfs_sd, & ! end smoke/chem variables
378 &rublten,rvblten,rthblten, &
379 &rqvblten,rqcblten,rqiblten, &
380 &rqncblten,rqniblten,rqsblten, &
381 &rqnwfablten,rqnifablten, &
382 &rqnbcablten,dozone, &
386 &dqke,qwt,qshear,qbuoy,qdiss, &
387 &qc_bl,qi_bl,cldfra_bl, &
388 &bl_mynn_tkeadvect, &
391 &bl_mynn_mixlength, &
397 &bl_mynn_mixscalars, &
399 &bl_mynn_cloudmix,bl_mynn_mixqt, &
400 &edmf_a,edmf_w,edmf_qt, &
401 &edmf_thl,edmf_ent,edmf_qc, &
402 &sub_thl3D,sub_sqv3D, &
403 &det_thl3D,det_sqv3D, &
404 &maxwidth,maxMF,ztop_plume, &
406 &spp_pbl,pattern_spp_pbl, &
408 &FLAG_QC,FLAG_QI,FLAG_QNC, &
410 &FLAG_QNWFA,FLAG_QNIFA, &
411 &FLAG_QNBCA,FLAG_OZONE, &
412 &IDS,IDE,JDS,JDE,KDS,KDE, &
413 &IMS,IME,JMS,JME,KMS,KME, &
414 &ITS,ITE,JTS,JTE,KTS,KTE )
418 integer,
intent(in) :: initflag
420 logical,
intent(in) :: restart,cycling
421 integer,
intent(in) :: tke_budget
422 integer,
intent(in) :: bl_mynn_cloudpdf
423 integer,
intent(in) :: bl_mynn_mixlength
424 integer,
intent(in) :: bl_mynn_edmf
425 logical,
intent(in) :: bl_mynn_tkeadvect
426 integer,
intent(in) :: bl_mynn_edmf_mom
427 integer,
intent(in) :: bl_mynn_edmf_tke
428 integer,
intent(in) :: bl_mynn_mixscalars
429 integer,
intent(in) :: bl_mynn_output
430 integer,
intent(in) :: bl_mynn_cloudmix
431 integer,
intent(in) :: bl_mynn_mixqt
432 integer,
intent(in) :: icloud_bl
433 real(kind_phys),
intent(in) :: closure
435 logical,
intent(in) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC,&
436 FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA, &
439 logical,
intent(in) :: mix_chem,enh_mix,rrfs_sd,smoke_dbg
441 integer,
intent(in) :: &
442 & IDS,IDE,JDS,JDE,KDS,KDE &
443 &,IMS,IME,JMS,JME,KMS,KME &
444 &,ITS,ITE,JTS,JTE,KTS,KTE
446#ifdef HARDCODE_VERTICAL
448# define kte HARDCODE_VERTICAL
457 real(kind_phys),
intent(in) :: delt
458 real(kind_phys),
dimension(:),
intent(in) :: dx
459 real(kind_phys),
dimension(:,:),
intent(in) :: dz, &
460 &u,v,w,th,sqv3D,p,exner,rho,T3D
461 real(kind_phys),
dimension(:,:),
intent(in) :: &
462 &sqc3D,sqi3D,sqs3D,qni,qnc,qnwfa,qnifa,qnbca
463 real(kind_phys),
dimension(:,:),
intent(in):: ozone
464 real(kind_phys),
dimension(:),
intent(in):: ust, &
466 real(kind_phys),
dimension(:,:),
intent(inout) :: &
468 real(kind_phys),
dimension(:,:),
intent(inout) :: &
470 real(kind_phys),
dimension(:,:),
intent(inout) :: &
471 &rublten,rvblten,rthblten,rqvblten,rqcblten, &
472 &rqiblten,rqsblten,rqniblten,rqncblten, &
473 &rqnwfablten,rqnifablten,rqnbcablten
474 real(kind_phys),
dimension(:,:),
intent(inout) :: dozone
475 real(kind_phys),
dimension(:,:),
intent(in) :: rthraten
477 real(kind_phys),
dimension(:,:),
intent(out) :: exch_h,exch_m
478 real(kind_phys),
dimension(:),
intent(in) :: xland, &
479 &ts,znt,hfx,qfx,uoce,voce
482 real(kind_phys),
dimension(:,:),
intent(inout),
optional :: &
483 & edmf_a,edmf_w,edmf_qt,edmf_thl,edmf_ent,edmf_qc, &
484 & sub_thl3D,sub_sqv3D,det_thl3D,det_sqv3D
489 real(kind_phys),
dimension(:),
intent(inout) :: Pblh
490 real(kind_phys),
dimension(:),
intent(inout) :: rmol
492 real(kind_phys),
dimension(ims:ime) :: psig_bl,psig_shcu
494 integer,
dimension(:),
intent(INOUT) :: &
496 integer,
dimension(:),
intent(INOUT) :: &
499 real(kind_phys),
dimension(:),
intent(out) :: &
500 &maxmf,maxwidth,ztop_plume
502 real(kind_phys),
dimension(:,:),
intent(inout) :: el_pbl
504 real(kind_phys),
dimension(:,:),
intent(inout),
optional :: &
505 &qWT,qSHEAR,qBUOY,qDISS,dqke
508 real(kind_phys),
dimension(kts:kte) :: &
509 &qwt1,qshear1,qbuoy1,qdiss1,dqke1,diss_heat
511 real(kind_phys),
dimension(:,:),
intent(out) :: Sh3D,Sm3D
513 real(kind_phys),
dimension(:,:),
intent(inout) :: &
514 &qc_bl,qi_bl,cldfra_bl
515 real(kind_phys),
dimension(kts:kte) :: qc_bl1D,qi_bl1D, &
516 &cldfra_bl1D,qc_bl1D_old,qi_bl1D_old,cldfra_bl1D_old
519 integer,
intent(IN ) :: nchem, kdvel, ndvel
520 real(kind_phys),
dimension(:,:,:),
intent(INOUT),
optional :: chem3d
521 real(kind_phys),
dimension(:,:),
intent(IN),
optional :: vdep
522 real(kind_phys),
dimension(:),
intent(IN),
optional :: frp
523 real(kind_phys),
dimension(:),
intent(IN) :: EMIS_ANT_NO
525 real(kind_phys),
dimension(kts:kte ,nchem) :: chem1
526 real(kind_phys),
dimension(kts:kte+1,nchem) :: s_awchem1
527 real(kind_phys),
dimension(ndvel) :: vd1
531 integer :: ITF,JTF,KTF, IMD,JMD
532 integer :: i,j,k,kproblem
533 real(kind_phys),
dimension(kts:kte) :: &
534 &thl,tl,qv1,qc1,qi1,qs1,sqw, &
535 &el, dfm, dfh, dfq, tcd, qcd, pdk, pdt, pdq, pdc, &
537 real(kind_phys),
dimension(kts:kte) :: &
538 &thetav,sh,sm,u1,v1,w1,p1, &
539 &ex1,dz1,th1,tk1,rho1,qke1,tsq1,qsq1,cov1, &
541 &du1,dv1,dth1,dqv1,dqc1,dqi1,dqs1,ozone1, &
542 &k_m1,k_h1,qni1,dqni1,qnc1,dqnc1,qnwfa1,qnifa1, &
543 &qnbca1,dqnwfa1,dqnifa1,dqnbca1,dozone1
546 real(kind_phys),
dimension(kts:kte) :: &
547 &dth1mf,dqv1mf,dqc1mf,du1mf,dv1mf
548 real(kind_phys),
dimension(kts:kte) :: &
549 &edmf_a1,edmf_w1,edmf_qt1,edmf_thl1, &
551 real(kind_phys),
dimension(kts:kte) :: &
552 &edmf_a_dd1,edmf_w_dd1,edmf_qt_dd1,edmf_thl_dd1, &
553 &edmf_ent_dd1,edmf_qc_dd1
554 real(kind_phys),
dimension(kts:kte) :: &
555 &sub_thl,sub_sqv,sub_u,sub_v, &
556 &det_thl,det_sqv,det_sqc,det_u,det_v
557 real(kind_phys),
dimension(kts:kte+1) :: &
558 &s_aw1,s_awthl1,s_awqt1, &
559 &s_awqv1,s_awqc1,s_awu1,s_awv1,s_awqke1, &
560 &s_awqnc1,s_awqni1,s_awqnwfa1,s_awqnifa1, &
562 real(kind_phys),
dimension(kts:kte+1) :: &
563 &sd_aw1,sd_awthl1,sd_awqt1, &
564 &sd_awqv1,sd_awqc1,sd_awu1,sd_awv1,sd_awqke1
566 real(kind_phys),
dimension(kts:kte+1) :: zw
567 real(kind_phys) :: cpm,sqcg,flt,fltv,flq,flqv,flqc, &
568 &pmz,phh,exnerg,zet,phi_m, &
569 &afk,abk,ts_decay, qc_bl2, qi_bl2, &
573 real(kind_phys),
dimension(ITS:ITE) :: maxKHtopdown
574 real(kind_phys),
dimension(kts:kte) :: KHtopdown,TKEprodTD
576 logical :: INITIALIZE_QKE,problem
579 integer,
intent(IN) :: spp_pbl
580 real(kind_phys),
dimension(:,:),
intent(IN),
optional :: pattern_spp_pbl
581 real(kind_phys),
dimension(KTS:KTE) :: rstoch_col
585 real(kind_phys) :: delt2
592 wsp = sqrt(u(i,k)**2 + v(i,k)**2)
593 if (abs(hfx(i)) > 1200. .or. abs(qfx(i)) > 0.001 .or. &
594 wsp > 200. .or. t3d(i,k) > 360. .or. t3d(i,k) < 160. .or. &
595 sqv3d(i,k)< 0.0 .or. sqc3d(i,k)< 0.0 )
then
598 print*,
"Incoming problem at: i=",i,
" k=1"
599 print*,
" QFX=",qfx(i),
" HFX=",hfx(i)
600 print*,
" wsp=",wsp,
" T=",t3d(i,k)
601 print*,
" qv=",sqv3d(i,k),
" qc=",sqc3d(i,k)
602 print*,
" u*=",ust(i),
" wspd=",wspd(i)
603 print*,
" xland=",xland(i),
" ts=",ts(i)
604 print*,
" z/L=",0.5*dz(i,1)*rmol(i),
" ps=",ps(i)
605 print*,
" znt=",znt(i),
" dx=",dx(i)
609 print*,
"===tk:",t3d(i,max(kproblem-3,1):min(kproblem+3,kte))
610 print*,
"===qv:",sqv3d(i,max(kproblem-3,1):min(kproblem+3,kte))
611 print*,
"===qc:",sqc3d(i,max(kproblem-3,1):min(kproblem+3,kte))
612 print*,
"===qi:",sqi3d(i,max(kproblem-3,1):min(kproblem+3,kte))
613 print*,
"====u:",u(i,max(kproblem-3,1):min(kproblem+3,kte))
614 print*,
"====v:",v(i,max(kproblem-3,1):min(kproblem+3,kte))
628 IF (bl_mynn_output > 0)
THEN
629 edmf_a(its:ite,kts:kte)=0.
630 edmf_w(its:ite,kts:kte)=0.
631 edmf_qt(its:ite,kts:kte)=0.
632 edmf_thl(its:ite,kts:kte)=0.
633 edmf_ent(its:ite,kts:kte)=0.
634 edmf_qc(its:ite,kts:kte)=0.
635 sub_thl3d(its:ite,kts:kte)=0.
636 sub_sqv3d(its:ite,kts:kte)=0.
637 det_thl3d(its:ite,kts:kte)=0.
638 det_sqv3d(its:ite,kts:kte)=0.
647 ktop_plume(its:ite)=0
648 ztop_plume(its:ite)=0.
651 maxkhtopdown(its:ite)=0.
659 IF (initflag > 0 .and. .not.restart)
THEN
662 IF ( (restart .or. cycling))
THEN
663 IF (maxval(qke(its:ite,kts)) < 0.0002)
THEN
664 initialize_qke = .true.
667 initialize_qke = .false.
671 initialize_qke = .true.
675 if (.not.restart .or. .not.cycling)
THEN
676 sh3d(its:ite,kts:kte)=0.
677 sm3d(its:ite,kts:kte)=0.
678 el_pbl(its:ite,kts:kte)=0.
679 tsq(its:ite,kts:kte)=0.
680 qsq(its:ite,kts:kte)=0.
681 cov(its:ite,kts:kte)=0.
682 cldfra_bl(its:ite,kts:kte)=0.
683 qc_bl(its:ite,kts:kte)=0.
684 qke(its:ite,kts:kte)=0.
688 cldfra_bl1d(kts:kte)=0.0
698 qc_bl1d_old(kts:kte)=0.0
699 cldfra_bl1d_old(kts:kte)=0.0
702 edmf_qc1(kts:kte)=0.0
703 edmf_a_dd1(kts:kte)=0.0
704 edmf_w_dd1(kts:kte)=0.0
705 edmf_qc_dd1(kts:kte)=0.0
717 IF (tke_budget .eq. 1)
THEN
740 if (icloud_bl > 0)
then
741 cldfra_bl1d(:)=cldfra_bl(i,:)
742 qc_bl1d(:)=qc_bl(i,:)
743 qi_bl1d(:)=qi_bl(i,:)
757 thetav(k)=th(i,k)*(1.+p608*sqv(k))
759 sqw(k)=sqv(k)+sqc(k)+sqi(k)
760 thl(k)=th1(k) - xlvcp/ex1(k)*sqc(k) &
761 & - xlscp/ex1(k)*(sqi(k))
770 zw(k)=zw(k-1)+dz(i,k-1)
772 IF (initialize_qke)
THEN
776 qke1(k)=5.*ust(i) * max((ust(i)*700. - zw(k))/(max(ust(i),0.01)*700.), 0.01)
787 rstoch_col(k)=pattern_spp_pbl(i,k)
794 zw(kte+1)=zw(kte)+dz(i,kte)
797 CALL get_pblh(kts,kte,pblh(i),thetav,&
798 & qke1,zw,dz1,xland(i),kpbl(i))
802 IF (scaleaware > 0.)
THEN
803 CALL scale_aware(dx(i),pblh(i),psig_bl(i),psig_shcu(i))
818 &pblh(i), th1, thetav, sh, sm, &
820 &el, qke1, tsq1, qsq1, cov1, &
821 &psig_bl(i), cldfra_bl1d, &
822 &bl_mynn_mixlength, &
825 &spp_pbl,rstoch_col )
827 IF (.not.restart)
THEN
839 IF (bl_mynn_tkeadvect)
THEN
863 IF (bl_mynn_tkeadvect)
THEN
869 if (tke_budget .eq. 1)
then
882 if (icloud_bl > 0)
then
883 cldfra_bl1d(:)=cldfra_bl(i,:)
884 qc_bl1d(:) =qc_bl(i,:)
885 qi_bl1d(:) =qi_bl(i,:)
886 cldfra_bl1d_old(:)=cldfra_bl(i,:)
887 qc_bl1d_old(:)=qc_bl(i,:)
888 qi_bl1d_old(:)=qi_bl(i,:)
897 dz1(kts:kte) =dz(i,kts:kte)
898 u1(kts:kte) =u(i,kts:kte)
899 v1(kts:kte) =v(i,kts:kte)
900 w1(kts:kte) =w(i,kts:kte)
901 th1(kts:kte) =th(i,kts:kte)
902 tk1(kts:kte) =t3d(i,kts:kte)
903 p1(kts:kte) =p(i,kts:kte)
904 ex1(kts:kte) =exner(i,kts:kte)
905 rho1(kts:kte) =rho(i,kts:kte)
906 sqv(kts:kte) =sqv3d(i,kts:kte)
907 sqc(kts:kte) =sqc3d(i,kts:kte)
908 qv1(kts:kte) =sqv(kts:kte)/(1.-sqv(kts:kte))
909 qc1(kts:kte) =sqc(kts:kte)/(1.-sqv(kts:kte))
910 qi1(kts:kte) =sqi(kts:kte)/(1.-sqv(kts:kte))
911 qs1(kts:kte) =sqs(kts:kte)/(1.-sqv(kts:kte))
922 qni1(kts:kte)=qni(i,kts:kte)
927 qnc1(kts:kte)=qnc(i,kts:kte)
931 IF (flag_qnwfa )
THEN
932 qnwfa1(kts:kte)=qnwfa(i,kts:kte)
936 IF (flag_qnifa )
THEN
937 qnifa1(kts:kte)=qnifa(i,kts:kte)
941 IF (flag_qnbca )
THEN
942 qnbca1(kts:kte)=qnbca(i,kts:kte)
946 IF (flag_ozone )
THEN
947 ozone1(kts:kte)=ozone(i,kts:kte)
951 el(kts:kte) =el_pbl(i,kts:kte)
952 qke1(kts:kte)=qke(i,kts:kte)
953 sh(kts:kte) =sh3d(i,kts:kte)
954 sm(kts:kte) =sm3d(i,kts:kte)
955 tsq1(kts:kte)=tsq(i,kts:kte)
956 qsq1(kts:kte)=qsq(i,kts:kte)
957 cov1(kts:kte)=cov(i,kts:kte)
959 rstoch_col(kts:kte)=pattern_spp_pbl(i,kts:kte)
961 rstoch_col(kts:kte)=0.0
1006 zw(k)=zw(k-1)+dz(i,k-1)
1009 sqw(k)= sqv(k)+sqc(k)+sqi(k)
1010 thl(k)= th1(k) - xlvcp/ex1(k)*sqc(k) &
1011 & - xlscp/ex1(k)*(sqi(k))
1016 thetav(k)=th1(k)*(1.+p608*sqv(k))
1018 zw(kte+1)=zw(kte)+dz(i,kte)
1021 if ( mix_chem )
then
1023 vd1(ic) = vdep(i,ic)
1027 chem1(k,ic) = chem3d(i,k,ic)
1044 CALL get_pblh(kts,kte,pblh(i),thetav,&
1045 & qke1,zw,dz1,xland(i),kpbl(i))
1051 if (scaleaware > 0.)
then
1052 call scale_aware(dx(i),pblh(i),psig_bl(i),psig_shcu(i))
1059 cpm=cp*(1.+0.84*qv1(kts))
1060 exnerg=(ps(i)/p1000mb)**rcp
1069 flqv = qfx(i)/rho1(kts)
1071 th_sfc = ts(i)/ex1(kts)
1075 flt =hfx(i)/(rho1(kts)*cpm )-xlvcp*flqc/ex1(kts)
1076 fltv=flt + flqv*p608*th_sfc
1079 rmol(i) = -karman*gtr*fltv/max(ust(i)**3,1.0e-6)
1080 zet = 0.5*dz(i,kts)*rmol(i)
1081 zet = max(zet, -20.)
1084 if (bl_mynn_stfunc == 0)
then
1086 if ( zet >= 0.0 )
then
1087 pmz = 1.0 + (cphm_st-1.0) * zet
1088 phh = 1.0 + cphh_st * zet
1090 pmz = 1.0/ (1.0-cphm_unst*zet)**0.25 - zet
1091 phh = 1.0/sqrt(1.0-cphh_unst*zet)
1106 &dx(i),dz1,zw,xland(i), &
1107 &thl,sqw,sqv,sqc,sqi,sqs, &
1108 &p1,ex1,tsq1,qsq1,cov1, &
1109 &sh,el,bl_mynn_cloudpdf, &
1110 &qc_bl1d,qi_bl1d,cldfra_bl1d, &
1112 &vt, vq, th1, sgm, rmol(i), &
1113 &spp_pbl, rstoch_col )
1118 if (bl_mynn_topdown.eq.1)
then
1119 call topdown_cloudrad(kts,kte,dz1,zw, &
1120 &xland(i),kpbl(i),pblh(i), &
1121 &sqc,sqi,sqw,thl,th1,ex1,p1,rho1,thetav, &
1122 &cldfra_bl1d,rthraten(i,:), &
1123 &maxkhtopdown(i),khtopdown,tkeprodtd )
1125 maxkhtopdown(i) = 0.0
1126 khtopdown(kts:kte) = 0.0
1127 tkeprodtd(kts:kte) = 0.0
1130 if (bl_mynn_edmf > 0)
then
1133 &kts,kte,delt,zw,dz1,p1,rho1, &
1134 &bl_mynn_edmf_mom, &
1135 &bl_mynn_edmf_tke, &
1136 &bl_mynn_mixscalars, &
1137 &u1,v1,w1,th1,thl,thetav,tk1, &
1138 &sqw,sqv,sqc,qke1, &
1139 &qnc1,qni1,qnwfa1,qnifa1,qnbca1, &
1141 &ust(i),flt,fltv,flq,flqv, &
1142 &pblh(i),kpbl(i),dx(i), &
1147 &edmf_a1,edmf_w1,edmf_qt1, &
1148 &edmf_thl1,edmf_ent1,edmf_qc1, &
1150 &s_aw1,s_awthl1,s_awqt1, &
1152 &s_awu1,s_awv1,s_awqke1, &
1153 &s_awqnc1,s_awqni1, &
1154 &s_awqnwfa1,s_awqnifa1,s_awqnbca1, &
1157 &det_thl,det_sqv,det_sqc, &
1160 &nchem,chem1,s_awchem1, &
1162 &qc_bl1d,cldfra_bl1d, &
1163 &qc_bl1d_old,cldfra_bl1d_old, &
1165 &flag_qnc,flag_qni, &
1166 &flag_qnwfa,flag_qnifa,flag_qnbca, &
1168 &maxwidth(i),ktop_plume(i), &
1169 &maxmf(i),ztop_plume(i), &
1170 &spp_pbl,rstoch_col )
1173 if (bl_mynn_edmf_dd == 1)
then
1174 call ddmf_jpl(kts,kte,delt,zw,dz1,p1, &
1175 &u1,v1,th1,thl,thetav,tk1, &
1176 &sqw,sqv,sqc,rho1,ex1, &
1179 &edmf_a_dd1,edmf_w_dd1,edmf_qt_dd1, &
1180 &edmf_thl_dd1,edmf_ent_dd1, &
1182 &sd_aw1,sd_awthl1,sd_awqt1, &
1183 &sd_awqv1,sd_awqc1,sd_awu1,sd_awv1, &
1185 &qc_bl1d,cldfra_bl1d, &
1194 &kts,kte,xland(i),closure, &
1196 &u1, v1, thl, thetav, sqc, sqw, &
1197 &qke1, tsq1, qsq1, cov1, &
1199 &rmol(i), flt, fltv, flq, &
1205 &qwt1,qshear1,qbuoy1,qdiss1, &
1207 &psig_bl(i),psig_shcu(i), &
1208 &cldfra_bl1d,bl_mynn_mixlength, &
1211 &spp_pbl,rstoch_col )
1218 &ust(i), flt, flq, pmz, phh, &
1219 &el, dfq, rho1, pdk, pdt, pdq, pdc, &
1220 &qke1, tsq1, qsq1, cov1, &
1221 &s_aw1, s_awqke1, bl_mynn_edmf_tke, &
1222 &qwt1, qdiss1, tke_budget )
1224 if (dheat_opt > 0)
then
1227 diss_heat(k) = min(max(1.0*(qke1(k)**1.5)/(b1*max(0.5*(el(k)+el(k+1)),1.))/cp, 0.0),0.002)
1229 diss_heat(k) = diss_heat(k) * exp(-10000./max(p1(k),1.))
1233 diss_heat(1:kte) = 0.
1240 &u1, v1, th1, tk1, qv1, &
1241 &qc1, qi1, kzero, qnc1, qni1, &
1242 &ps(i), p1, ex1, thl, &
1243 &sqv, sqc, sqi, kzero, sqw, &
1244 &qnwfa1, qnifa1, qnbca1, ozone1, &
1245 &ust(i),flt,flq,flqv,flqc, &
1246 &wspd(i),uoce(i),voce(i), &
1247 &tsq1, qsq1, cov1, &
1250 &du1, dv1, dth1, dqv1, &
1251 &dqc1, dqi1, dqs1, dqnc1, dqni1, &
1252 &dqnwfa1, dqnifa1, dqnbca1, &
1256 &s_aw1,s_awthl1,s_awqt1, &
1257 &s_awqv1,s_awqc1,s_awu1,s_awv1, &
1258 &s_awqnc1,s_awqni1, &
1259 &s_awqnwfa1,s_awqnifa1,s_awqnbca1, &
1260 &sd_aw1,sd_awthl1,sd_awqt1, &
1261 &sd_awqv1,sd_awqc1, &
1265 &det_thl,det_sqv,det_sqc, &
1267 &flag_qc,flag_qi,flag_qnc, &
1268 &flag_qni,flag_qs, &
1269 &flag_qnwfa,flag_qnifa, &
1272 &bl_mynn_cloudmix, &
1275 &bl_mynn_edmf_mom, &
1276 &bl_mynn_mixscalars )
1279 if ( mix_chem )
then
1281 call mynn_mix_chem(kts,kte,i, &
1282 &delt, dz1, pblh(i), &
1283 &nchem, kdvel, ndvel, &
1291 &enh_mix, smoke_dbg )
1293 call mynn_mix_chem(kts,kte,i, &
1294 &delt, dz1, pblh(i), &
1295 &nchem, kdvel, ndvel, &
1303 &enh_mix, smoke_dbg )
1307 chem3d(i,k,ic) = max(1.e-12, chem1(k,ic))
1312 call retrieve_exchange_coeffs(kts,kte, &
1313 &dfm, dfh, dz1, k_m1, k_h1 )
1316 exch_m(i,kts:kte) =k_m1(kts:kte)
1317 exch_h(i,kts:kte) =k_h1(kts:kte)
1318 rublten(i,kts:kte) =du1(kts:kte)
1319 rvblten(i,kts:kte) =dv1(kts:kte)
1320 rthblten(i,kts:kte)=dth1(kts:kte)
1321 rqvblten(i,kts:kte)=dqv1(kts:kte)
1322 if (bl_mynn_cloudmix > 0)
then
1323 if (flag_qc) rqcblten(i,kts:kte)=dqc1(kts:kte)
1324 if (flag_qi) rqiblten(i,kts:kte)=dqi1(kts:kte)
1325 if (flag_qs) rqsblten(i,kts:kte)=dqs1(kts:kte)
1327 if (flag_qc) rqcblten(i,:)=0.
1328 if (flag_qi) rqiblten(i,:)=0.
1329 if (flag_qs) rqsblten(i,:)=0.
1331 if (bl_mynn_cloudmix > 0 .and. bl_mynn_mixscalars > 0)
then
1332 if (flag_qnc) rqncblten(i,kts:kte) =dqnc1(kts:kte)
1333 if (flag_qni) rqniblten(i,kts:kte) =dqni1(kts:kte)
1334 if (flag_qnwfa) rqnwfablten(i,kts:kte)=dqnwfa1(kts:kte)
1335 if (flag_qnifa) rqnifablten(i,kts:kte)=dqnifa1(kts:kte)
1336 if (flag_qnbca) rqnbcablten(i,kts:kte)=dqnbca1(kts:kte)
1338 if (flag_qnc) rqncblten(i,:) =0.
1339 if (flag_qni) rqniblten(i,:) =0.
1340 if (flag_qnwfa) rqnwfablten(i,:)=0.
1341 if (flag_qnifa) rqnifablten(i,:)=0.
1342 if (flag_qnbca) rqnbcablten(i,:)=0.
1344 dozone(i,kts:kte)=dozone1(kts:kte)
1345 if (icloud_bl > 0)
then
1346 qc_bl(i,kts:kte) =qc_bl1d(kts:kte)
1347 qi_bl(i,kts:kte) =qi_bl1d(kts:kte)
1348 cldfra_bl(i,kts:kte)=cldfra_bl1d(kts:kte)
1350 el_pbl(i,kts:kte)=el(kts:kte)
1351 qke(i,kts:kte) =qke1(kts:kte)
1352 tsq(i,kts:kte) =tsq1(kts:kte)
1353 qsq(i,kts:kte) =qsq1(kts:kte)
1354 cov(i,kts:kte) =cov1(kts:kte)
1355 sh3d(i,kts:kte) =sh(kts:kte)
1356 sm3d(i,kts:kte) =sm(kts:kte)
1358 if (tke_budget .eq. 1)
then
1362 qshear1(k) =4.*(ust(i)**3*phi_m/(karman*dz(i,k)))-qshear1(k+1)
1363 qbuoy1(k) =4.*(-ust(i)**3*zet/(karman*dz(i,k)))-qbuoy1(k+1)
1366 qshear(i,k)=0.5*(qshear1(k)+qshear1(k+1))
1367 qbuoy(i,k) =0.5*(qbuoy1(k)+qbuoy1(k+1))
1369 qdiss(i,k) =qdiss1(k)
1370 dqke(i,k) =(qke1(k)-dqke(i,k))*0.5/delt
1382 if (bl_mynn_output > 0)
then
1383 if (bl_mynn_edmf > 0)
then
1384 edmf_a(i,kts:kte) =edmf_a1(kts:kte)
1385 edmf_w(i,kts:kte) =edmf_w1(kts:kte)
1386 edmf_qt(i,kts:kte) =edmf_qt1(kts:kte)
1387 edmf_thl(i,kts:kte) =edmf_thl1(kts:kte)
1388 edmf_ent(i,kts:kte) =edmf_ent1(kts:kte)
1389 edmf_qc(i,kts:kte) =edmf_qc1(kts:kte)
1390 sub_thl3d(i,kts:kte)=sub_thl(kts:kte)
1391 sub_sqv3d(i,kts:kte)=sub_sqv(kts:kte)
1392 det_thl3d(i,kts:kte)=det_thl(kts:kte)
1393 det_sqv3d(i,kts:kte)=det_sqv(kts:kte)
1406 if ( debug_code .and. (i .eq. idbg))
THEN
1407 if ( abs(qfx(i))>.001)print*,&
1408 "SUSPICIOUS VALUES AT: i=",i,
" QFX=",qfx(i)
1409 if ( abs(hfx(i))>1100.)print*,&
1410 "SUSPICIOUS VALUES AT: i=",i,
" HFX=",hfx(i)
1412 IF ( sh(k) < 0. .OR. sh(k)> 200.)print*,&
1413 "SUSPICIOUS VALUES AT: i,k=",i,k,
" sh=",sh(k)
1414 IF ( abs(vt(k)) > 2.0 )print*,&
1415 "SUSPICIOUS VALUES AT: i,k=",i,k,
" vt=",vt(k)
1416 IF ( abs(vq(k)) > 7000.)print*,&
1417 "SUSPICIOUS VALUES AT: i,k=",i,k,
" vq=",vq(k)
1418 IF ( qke(i,k) < -1. .OR. qke(i,k)> 200.)print*,&
1419 "SUSPICIOUS VALUES AT: i,k=",i,k,
" qke=",qke(i,k)
1420 IF ( el_pbl(i,k) < 0. .OR. el_pbl(i,k)> 1500.)print*,&
1421 "SUSPICIOUS VALUES AT: i,k=",i,k,
" el_pbl=",el_pbl(i,k)
1422 IF ( exch_m(i,k) < 0. .OR. exch_m(i,k)> 2000.)print*,&
1423 "SUSPICIOUS VALUES AT: i,k=",i,k,
" exxch_m=",exch_m(i,k)
1424 IF (icloud_bl > 0)
then
1425 IF ( cldfra_bl(i,k) < 0.0 .OR. cldfra_bl(i,k)> 1.)
THEN
1426 print*,
"SUSPICIOUS VALUES: CLDFRA_BL=",cldfra_bl(i,k),
" qc_bl=",qc_bl(i,k)
1444 IF (bl_mynn_tkeadvect)
THEN
1449#ifdef HARDCODE_VERTICAL
1519! & ust, rmo, pmz, phh, flt, flq, &
1520 & zi, theta, thetav, sh, sm, &
1522 & Qke, Tsq, Qsq, Cov, Psig_bl, cldfra_bl1D, &
1523 & bl_mynn_mixlength, &
1524 & edmf_w1,edmf_a1, &
1526 & spp_pbl,rstoch_col)
1530 integer,
intent(in) :: kts,kte
1531 integer,
intent(in) :: bl_mynn_mixlength
1532 logical,
intent(in) :: INITIALIZE_QKE
1534 real(kind_phys),
intent(in) :: rmo, Psig_bl, xland
1535 real(kind_phys),
intent(in) :: dx, ust, zi
1536 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz
1537 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
1538 real(kind_phys),
dimension(kts:kte),
intent(in) :: u,v,thl,&
1539 &qw,cldfra_bl1D,edmf_w1,edmf_a1
1540 real(kind_phys),
dimension(kts:kte),
intent(out) :: tsq,qsq,cov
1541 real(kind_phys),
dimension(kts:kte),
intent(inout) :: el,qke
1542 real(kind_phys),
dimension(kts:kte) :: &
1543 &ql,pdk,pdt,pdq,pdc,dtl,dqw,dtv, &
1544 &gm,gh,sm,sh,qkw,vt,vq
1546 real(kind_phys):: phm,vkz,elq,elv,b1l,b2l,pmz=1.,phh=1., &
1547 &flt=0.,fltv=0.,flq=0.,tmpq
1548 real(kind_phys),
dimension(kts:kte) :: theta,thetav
1549 real(kind_phys),
dimension(kts:kte) :: rstoch_col
1562 & u, v, thl, thetav, qw, &
1564 & dtl, dqw, dtv, gm, gh, sm, sh )
1569 IF (initialize_qke)
THEN
1571 qke(kts) = 1.5 * ust**2 * ( b1*pmz )**(2.0/3.0)
1575 qke(k)=qke(kts)*max((ust*700. - zw(k))/(max(ust,0.01)*700.), 0.01)
1579 phm = phh*b2 / ( b1*pmz )**(1.0/3.0)
1580 tsq(kts) = phm*( flt/ust )**2
1581 qsq(kts) = phm*( flq/ust )**2
1582 cov(kts) = phm*( flt/ust )*( flq/ust )
1586 el(k) = vkz/( 1.0 + vkz/100.0 )
1604 & rmo, flt, fltv, flq, &
1610 & qkw,psig_bl,cldfra_bl1d, &
1611 & bl_mynn_mixlength, &
1616 pdk(k) = elq*( sm(k)*gm(k) + &
1618 pdt(k) = elq* sh(k)*dtl(k)**2
1619 pdq(k) = elq* sh(k)*dqw(k)**2
1620 pdc(k) = elq* sh(k)*dtl(k)*dqw(k)
1624 vkz = karman*0.5*dz(kts)
1625 elv = 0.5*( el(kts+1)+el(kts) ) / vkz
1626 IF (initialize_qke)
THEN
1628 qke(kts) = 1.0 * max(ust,0.02)**2 * ( b1*pmz*elv )**(2.0/3.0)
1631 phm = phh*b2 / ( b1*pmz/elv**2 )**(1.0/3.0)
1632 tsq(kts) = phm*( flt/ust )**2
1633 qsq(kts) = phm*( flq/ust )**2
1634 cov(kts) = phm*( flt/ust )*( flq/ust )
1637 b1l = b1*0.25*( el(k+1)+el(k) )
1640 tmpq=min(max(b1l*( pdk(k+1)+pdk(k) ),qkemin),125.)
1642 IF (initialize_qke)
THEN
1643 qke(k) = tmpq**twothirds
1646 IF ( qke(k) .LE. 0.0 )
THEN
1649 b2l = b2*( b1l/b1 ) / sqrt( qke(k) )
1652 tsq(k) = b2l*( pdt(k+1)+pdt(k) )
1653 qsq(k) = b2l*( pdq(k+1)+pdq(k) )
1654 cov(k) = b2l*( pdc(k+1)+pdc(k) )
1664 IF (initialize_qke)
THEN
1665 qke(kts)=0.5*(qke(kts)+qke(kts+1))
1721 & u, v, thl, thetav, qw, &
1723 & dtl, dqw, dtv, gm, gh, sm, sh )
1727 integer,
intent(in) :: kts,kte
1729#ifdef HARDCODE_VERTICAL
1731# define kte HARDCODE_VERTICAL
1734 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz
1735 real(kind_phys),
dimension(kts:kte),
intent(in) :: u,v, &
1736 &thl,qw,ql,vt,vq,thetav
1737 real(kind_phys),
dimension(kts:kte),
intent(out) :: &
1738 &dtl,dqw,dtv,gm,gh,sm,sh
1742 real(kind_phys):: rfc,f1,f2,rf1,rf2,smc,shc, &
1743 &ri1,ri2,ri3,ri4,duz,dtz,dqz,vtt,vqq,dtq,dzk, &
1746 real(kind_phys):: a2fac
1754 f1 = b1*( g1-c1 ) +3.0*a2*( 1.0 -c2 )*( 1.0-c5 ) &
1755 & +2.0*a1*( 3.0-2.0*c2 )
1756 f2 = b1*( g1+g2 ) -3.0*a1*( 1.0 -c2 )
1757 rf1 = b1*( g1-c1 )/f1
1760 shc = 3.0*a2*( g1+g2 )
1764 ri3 = 4.0*rf2*smc -2.0*ri2
1768 dzk = 0.5 *( dz(k)+dz(k-1) )
1769 afk = dz(k)/( dz(k)+dz(k-1) )
1771 duz = ( u(k)-u(k-1) )**2 +( v(k)-v(k-1) )**2
1773 dtz = ( thl(k)-thl(k-1) )/( dzk )
1774 dqz = ( qw(k)-qw(k-1) )/( dzk )
1776 vtt = 1.0 +vt(k)*abk +vt(k-1)*afk
1777 vqq = tv0 +vq(k)*abk +vq(k-1)*afk
1778 dtq = vtt*dtz +vqq*dqz
1793 ri = -gh(k)/max( duz, 1.0e-10 )
1796 IF (ckmod .eq. 1)
THEN
1797 a2fac = 1./(1. + max(ri,0.0))
1803 f1 = b1*( g1-c1 ) +3.0*a2*a2fac *( 1.0 -c2 )*( 1.0-c5 ) &
1804 & +2.0*a1*( 3.0-2.0*c2 )
1805 f2 = b1*( g1+g2 ) -3.0*a1*( 1.0 -c2 )
1806 rf1 = b1*( g1-c1 )/f1
1808 smc = a1 /(a2*a2fac)* f1/f2
1809 shc = 3.0*(a2*a2fac)*( g1+g2 )
1813 ri3 = 4.0*rf2*smc -2.0*ri2
1817 rf = min( ri1*( ri + ri2-sqrt(ri**2 - ri3*ri + ri4) ), rfc )
1819 sh(k) = shc*( rfc-rf )/( 1.0-rf )
1820 sm(k) = smc*( rf1-rf )/( rf2-rf ) * sh(k)
1825#ifdef HARDCODE_VERTICAL
1853 & rmo, flt, fltv, flq, &
1859 & Psig_bl, cldfra_bl1D, &
1860 & bl_mynn_mixlength, &
1865 integer,
intent(in) :: kts,kte
1867#ifdef HARDCODE_VERTICAL
1869# define kte HARDCODE_VERTICAL
1872 integer,
intent(in) :: bl_mynn_mixlength
1873 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz
1874 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
1875 real(kind_phys),
intent(in) :: rmo,flt,fltv,flq,Psig_bl,xland
1876 real(kind_phys),
intent(in) :: dx,zi
1877 real(kind_phys),
dimension(kts:kte),
intent(in) :: u1,v1, &
1878 &qke,vt,vq,cldfra_bl1D,edmf_w1,edmf_a1
1879 real(kind_phys),
dimension(kts:kte),
intent(out) :: qkw, el
1880 real(kind_phys),
dimension(kts:kte),
intent(in) :: dtv
1881 real(kind_phys):: elt,vsc
1882 real(kind_phys),
dimension(kts:kte),
intent(in) :: theta
1883 real(kind_phys),
dimension(kts:kte) :: qtke,elBLmin,elBLavg,thetaw
1884 real(kind_phys):: wt,wt2,zi2,h1,h2,hs,elBLmin0,elBLavg0,cldavg
1888 real(kind_phys):: cns, & !< for surface layer (els) in stable conditions
1889 alp1, & !< for turbulent length scale (elt)
1890 alp2, & !< for buoyancy length scale (elb)
1891 alp3, & !< for buoyancy enhancement factor of elb
1892 alp4, & !< for surface layer (els) in unstable conditions
1893 alp5, & !< for BouLac mixing length or above PBLH
1900 real(kind_phys),
parameter :: minzi = 300.
1901 real(kind_phys),
parameter :: maxdz = 750.
1904 real(kind_phys),
parameter :: mindz = 300.
1907 real(kind_phys),
parameter :: ZSLH = 100.
1908 real(kind_phys),
parameter :: CSL = 2.
1912 real(kind_phys):: afk,abk,zwk,zwk1,dzk,qdz,vflx,bv,tau_cloud, &
1913 & wstar,elb,els,elf,el_stab,el_mf,el_stab_mf,elb_mf, &
1914 & PBLH_PLUS_ENT,Uonset,Ugrid,wt_u,el_les
1915 real(kind_phys),
parameter :: ctau = 1000.
1920 SELECT CASE(bl_mynn_mixlength)
1932 zi2 = min(10000.,zw(kte-2))
1933 h1=max(0.3*zi2,mindz)
1937 qkw(kts) = sqrt(max(qke(kts), qkemin))
1939 afk = dz(k)/( dz(k)+dz(k-1) )
1941 qkw(k) = sqrt(max(qke(k)*abk+qke(k-1)*afk, qkemin))
1950 DO WHILE (zwk .LE. zi2+h1)
1951 dzk = 0.5*( dz(k)+dz(k-1) )
1952 qdz = max( qkw(k)-qmin, 0.03 )*dzk
1960 vflx = ( vt(kts)+1.0 )*flt +( vq(kts)+tv0 )*flq
1961 vsc = ( gtr*elt*max( vflx, 0.0 ) )**onethird
1971 IF ( dtv(k) .GT. 0.0 )
THEN
1972 bv = sqrt( gtr*dtv(k) )
1973 elb = alp2*qkw(k) / bv &
1974 & *( 1.0 + alp3/alp2*&
1975 &sqrt( vsc/( bv*elt ) ) )
1976 elf = alp2 * qkw(k)/bv
1984 IF ( rmo .GT. 0.0 )
THEN
1985 els = karman*zwk/(1.0+cns*min( zwk*rmo, zmax ))
1987 els = karman*zwk*( 1.0 - alp4* zwk*rmo )**0.2
1994 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
1996 el(k) = min(elb/( elb/elt+elb/els+1.0 ),elf)
2002 ugrid = sqrt(u1(kts)**2 + v1(kts)**2)
2004 wt_u = (1.0 - min(max(ugrid - uonset, 0.0)/30.0, 0.5))
2015 h1=max(0.3*zi2,300.)
2019 qtke(kts)=max(0.5*qke(kts), 0.5*qkemin)
2020 thetaw(kts)=theta(kts)
2021 qkw(kts) = sqrt(max(qke(kts), qkemin))
2024 afk = dz(k)/( dz(k)+dz(k-1) )
2026 qkw(k) = sqrt(max(qke(k)*abk+qke(k-1)*afk, qkemin))
2027 qtke(k) = 0.5*(qkw(k)**2)
2028 thetaw(k)= theta(k)*abk + theta(k-1)*afk
2037 DO WHILE (zwk .LE. zi2+h1)
2038 dzk = 0.5*( dz(k)+dz(k-1) )
2039 qdz = min(max( qkw(k)-qmin, 0.01 ), 30.0)*dzk
2046 elt = min( max( alp1*elt/vsc, 8.), 400.)
2050 vsc = ( gtr*elt*max( vflx, 0.0 ) )**onethird
2057 CALL boulac_length(kts,kte,zw,dz,qtke,thetaw,elblmin,elblavg)
2063 IF ( dtv(k) .GT. 0.0 )
THEN
2064 bv = max( sqrt( gtr*dtv(k) ), 0.0001)
2065 elb = max(alp2*qkw(k), &
2066 & alp6*edmf_a1(k-1)*edmf_w1(k-1)) / bv &
2067 & *( 1.0 + alp3*sqrt( vsc/(bv*elt) ) )
2069 elf = 1.0 * qkw(k)/bv
2070 elblavg(k) = max(elblavg(k), alp6*edmf_a1(k-1)*edmf_w1(k-1)/bv)
2077 IF ( rmo .GT. 0.0 )
THEN
2078 els = karman*zwk/(1.0+cns*min( zwk*rmo, zmax ))
2080 els = karman*zwk*( 1.0 - alp4* zwk*rmo )**0.2
2084 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
2091 el(k) = sqrt( els**2/(1. + (els**2/elt**2)))
2092 el(k) = min(el(k), elb)
2093 el(k) = min(el(k), elf)
2094 el(k) = el(k)*(1.-wt) + alp5*elblavg(k)*wt
2097 el(k) = el(k)*psig_bl
2103 uonset = 3.5 + dz(kts)*0.1
2104 ugrid = sqrt(u1(kts)**2 + v1(kts)**2)
2118 h1=max(0.3*zi2,300.)
2122 qtke(kts)=max(0.5*qke(kts), 0.5*qkemin)
2123 qkw(kts) = sqrt(max(qke(kts), qkemin))
2126 afk = dz(k)/( dz(k)+dz(k-1) )
2128 qkw(k) = sqrt(max(qke(k)*abk+qke(k-1)*afk, qkemin))
2129 qtke(k) = 0.5*qkw(k)**2
2136 pblh_plus_ent = max(zi+h1, 100.)
2139 DO WHILE (zwk .LE. pblh_plus_ent)
2140 dzk = 0.5*( dz(k)+dz(k-1) )
2141 qdz = min(max( qkw(k)-qmin, 0.03 ), 30.0)*dzk
2148 elt = min( max(alp1*elt/vsc, 10.), 400.)
2152 vsc = ( gtr*elt*max( vflx, 0.0 ) )**onethird
2160 dzk = 0.5*( dz(k)+dz(k-1) )
2161 cldavg = 0.5*(cldfra_bl1d(k-1)+cldfra_bl1d(k))
2164 IF ( dtv(k) .GT. 0.0 )
THEN
2166 bv = max( sqrt( gtr*dtv(k) ), 0.001)
2168 elb_mf = max(alp2*qkw(k), &
2169 & alp6*edmf_a1(k-1)*edmf_w1(k-1)) / bv &
2170 & *( 1.0 + alp3*sqrt( vsc/( bv*elt ) ) )
2171 elb = min(max(alp5*qkw(k), alp6*edmf_a1(k)*edmf_w1(k))/bv, zwk)
2174 wstar = 1.25*(gtr*zi*max(vflx,1.0e-4))**onethird
2175 tau_cloud = min(max(ctau * wstar/grav, 30.), 150.)
2177 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
2178 tau_cloud = tau_cloud*(1.-wt) + 50.*wt
2179 elf = min(max(tau_cloud*sqrt(min(qtke(k),40.)), &
2180 & alp6*edmf_a1(k)*edmf_w1(k)/bv), zwk)
2201 wstar = 1.25*(gtr*zi*max(vflx,1.0e-4))**onethird
2202 tau_cloud = min(max(ctau * wstar/grav, 50.), 200.)
2204 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
2206 tau_cloud = tau_cloud*(1.-wt) + max(100.,dzk*0.25)*wt
2208 elb = min(tau_cloud*sqrt(min(qtke(k),40.)), zwk)
2213 elf = elf/(1. + (elf/800.))
2214 elb_mf = max(elb_mf, 0.01)
2217 IF ( rmo .GT. 0.0 )
THEN
2218 els = karman*zwk/(1.0+cns*min( zwk*rmo, zmax ))
2220 els = karman*zwk*( 1.0 - alp4* zwk*rmo )**0.2
2224 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
2227 el(k) = sqrt( els**2/(1. + (els**2/elt**2) +(els**2/elb_mf**2)))
2228 el(k) = el(k)*(1.-wt) + elf*wt
2231 el_les= min(els/(1. + (els/12.)), elb_mf)
2232 el(k) = el(k)*psig_bl + (1.-psig_bl)*el_les
2239#ifdef HARDCODE_VERTICAL
2277 integer,
intent(in) :: k,kts,kte
2278 real(kind_phys),
dimension(kts:kte),
intent(in) :: qtke,dz,theta
2279 real(kind_phys),
intent(out) :: lb1,lb2
2280 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
2283 integer :: izz, found
2284 real(kind_phys):: dlu,dld
2285 real(kind_phys):: dzt, zup, beta, zup_inf, bbb, tl, zdo, zdo_sup, zzz
2292 dlu=zw(kte+1)-zw(k)-dz(k)*0.5
2299 if (k .lt. kte)
then
2302 DO WHILE (found .EQ. 0)
2304 if (izz .lt. kte)
then
2306 zup=zup-beta*theta(k)*dzt
2308 zup=zup+beta*(theta(izz+1)+theta(izz))*dzt*0.5
2311 if (qtke(k).lt.zup .and. qtke(k).ge.zup_inf)
then
2312 bbb=(theta(izz+1)-theta(izz))/dzt
2313 if (bbb .ne. 0.)
then
2315 tl=(-beta*(theta(izz)-theta(k)) + &
2316 & sqrt( max(0.,(beta*(theta(izz)-theta(k)))**2 + &
2317 & 2.*bbb*beta*(qtke(k)-zup_inf))))/bbb/beta
2319 if (theta(izz) .ne. theta(k))
then
2320 tl=(qtke(k)-zup_inf)/(beta*(theta(izz)-theta(k)))
2348 if (k .gt. kts)
then
2352 DO WHILE (found .EQ. 0)
2354 if (izz .gt. kts)
then
2356 zdo=zdo+beta*theta(k)*dzt
2358 zdo=zdo-beta*(theta(izz-1)+theta(izz))*dzt*0.5
2361 if (qtke(k).lt.zdo .and. qtke(k).ge.zdo_sup)
then
2362 bbb=(theta(izz)-theta(izz-1))/dzt
2363 if (bbb .ne. 0.)
then
2364 tl=(beta*(theta(izz)-theta(k))+ &
2365 & sqrt( max(0.,(beta*(theta(izz)-theta(k)))**2 + &
2366 & 2.*bbb*beta*(qtke(k)-zdo_sup))))/bbb/beta
2368 if (theta(izz) .ne. theta(k))
then
2369 tl=(qtke(k)-zdo_sup)/(beta*(theta(izz)-theta(k)))
2392 dld = min(dld,zw(k+1))
2395 dlu=max(0.1,min(dlu,1000.))
2396 dld=max(0.1,min(dld,1000.))
2400 if (k .eq. kte)
then
2427 integer,
intent(in) :: kts,kte
2428 real(kind_phys),
dimension(kts:kte),
intent(in) :: qtke,dz,theta
2429 real(kind_phys),
dimension(kts:kte),
intent(out):: lb1,lb2
2430 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
2433 integer :: iz, izz, found
2434 real(kind_phys),
dimension(kts:kte) :: dlu,dld
2435 real(kind_phys),
parameter :: Lmax=2000.
2436 real(kind_phys):: dzt, zup, beta, zup_inf, bbb, tl, zdo, zdo_sup, zzz
2446 dlu(iz)=zw(kte+1)-zw(iz)-dz(iz)*0.5
2453 if (iz .lt. kte)
then
2457 DO WHILE (found .EQ. 0)
2459 if (izz .lt. kte)
then
2461 zup=zup-beta*theta(iz)*dzt
2463 zup=zup+beta*(theta(izz+1)+theta(izz))*dzt*0.5
2466 if (qtke(iz).lt.zup .and. qtke(iz).ge.zup_inf)
then
2467 bbb=(theta(izz+1)-theta(izz))/dzt
2468 if (bbb .ne. 0.)
then
2470 tl=(-beta*(theta(izz)-theta(iz)) + &
2471 & sqrt( max(0.,(beta*(theta(izz)-theta(iz)))**2 + &
2472 & 2.*bbb*beta*(qtke(iz)-zup_inf))))/bbb/beta
2474 if (theta(izz) .ne. theta(iz))
then
2475 tl=(qtke(iz)-zup_inf)/(beta*(theta(izz)-theta(iz)))
2503 if (iz .gt. kts)
then
2507 DO WHILE (found .EQ. 0)
2509 if (izz .gt. kts)
then
2511 zdo=zdo+beta*theta(iz)*dzt
2513 zdo=zdo-beta*(theta(izz-1)+theta(izz))*dzt*0.5
2516 if (qtke(iz).lt.zdo .and. qtke(iz).ge.zdo_sup)
then
2517 bbb=(theta(izz)-theta(izz-1))/dzt
2518 if (bbb .ne. 0.)
then
2519 tl=(beta*(theta(izz)-theta(iz))+ &
2520 & sqrt( max(0.,(beta*(theta(izz)-theta(iz)))**2 + &
2521 & 2.*bbb*beta*(qtke(iz)-zdo_sup))))/bbb/beta
2523 if (theta(izz) .ne. theta(iz))
then
2524 tl=(qtke(iz)-zdo_sup)/(beta*(theta(izz)-theta(iz)))
2547 dld(iz) = min(dld(iz),zw(iz+1))
2548 lb1(iz) = min(dlu(iz),dld(iz))
2550 dlu(iz)=max(0.1,min(dlu(iz),1000.))
2551 dld(iz)=max(0.1,min(dld(iz),1000.))
2552 lb2(iz) = sqrt(dlu(iz)*dld(iz))
2556 lb1(iz) = lb1(iz)/(1. + (lb1(iz)/lmax))
2557 lb2(iz) = lb2(iz)/(1. + (lb2(iz)/lmax))
2559 if (iz .eq. kte)
then
2560 lb1(kte) = lb1(kte-1)
2561 lb2(kte) = lb2(kte-1)
2623 & u, v, thl, thetav, ql, qw, &
2624 & qke, tsq, qsq, cov, &
2626 & rmo, flt, fltv, flq, &
2630 & Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, &
2631 & qWT1D,qSHEAR1D,qBUOY1D,qDISS1D, &
2633 & Psig_bl,Psig_shcu,cldfra_bl1D, &
2634 & bl_mynn_mixlength, &
2635 & edmf_w1,edmf_a1, &
2637 & spp_pbl,rstoch_col )
2641 integer,
intent(in) :: kts,kte
2643#ifdef HARDCODE_VERTICAL
2645# define kte HARDCODE_VERTICAL
2648 integer,
intent(in) :: bl_mynn_mixlength,tke_budget
2649 real(kind_phys),
intent(in) :: closure
2650 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz
2651 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
2652 real(kind_phys),
intent(in) :: rmo,flt,fltv,flq, &
2653 &Psig_bl,Psig_shcu,xland,dx,zi
2654 real(kind_phys),
dimension(kts:kte),
intent(in) :: u,v,thl,thetav,qw, &
2655 &ql,vt,vq,qke,tsq,qsq,cov,cldfra_bl1D,edmf_w1,edmf_a1, &
2658 real(kind_phys),
dimension(kts:kte),
intent(out) :: dfm,dfh,dfq, &
2659 &pdk,pdt,pdq,pdc,tcd,qcd,el
2661 real(kind_phys),
dimension(kts:kte),
intent(inout) :: &
2662 qWT1D,qSHEAR1D,qBUOY1D,qDISS1D
2663 real(kind_phys):: q3sq_old,dlsq1,qWTP_old,qWTP_new
2664 real(kind_phys):: dudz,dvdz,dTdz,upwp,vpwp,Tpwp
2666 real(kind_phys),
dimension(kts:kte) :: qkw,dtl,dqw,dtv,gm,gh,sm,sh
2670 real(kind_phys):: e6c,dzk,afk,abk,vtt,vqq, &
2671 &cw25,clow,cupp,gamt,gamq,smd,gamv,elq,elh
2673 real(kind_phys):: cldavg
2674 real(kind_phys),
dimension(kts:kte),
intent(in) :: theta
2676 real(kind_phys):: a2fac, duz, ri
2678 real:: auh,aum,adh,adm,aeh,aem,Req,Rsl,Rsl2, &
2679 gmelq,sm20,sh20,sm25max,sh25max,sm25min,sh25min, &
2680 sm_pbl,sh_pbl,zi2,wt,slht,wtpr
2682 DOUBLE PRECISION q2sq, t2sq, r2sq, c2sq, elsq, gmel, ghel
2683 DOUBLE PRECISION q3sq, t3sq, r3sq, c3sq, dlsq, qdiv
2684 DOUBLE PRECISION e1, e2, e3, e4, enum, eden, wden
2687 integer,
intent(in) :: spp_pbl
2688 real(kind_phys),
dimension(kts:kte) :: rstoch_col
2689 real(kind_phys):: Prnum, shb
2690 real(kind_phys),
parameter :: Prlimit = 5.0
2707 & u, v, thl, thetav, qw, &
2709 & dtl, dqw, dtv, gm, gh, sm, sh )
2714 & rmo, flt, fltv, flq, &
2720 & qkw,psig_bl,cldfra_bl1d, &
2721 & bl_mynn_mixlength, &
2726 dzk = 0.5 *( dz(k)+dz(k-1) )
2727 afk = dz(k)/( dz(k)+dz(k-1) )
2731 q2sq = b1*elsq*( sm(k)*gm(k)+sh(k)*gh(k) )
2733 sh20 = max(sh(k), 1e-5)
2734 sm20 = max(sm(k), 1e-5)
2735 sh(k)= max(sh(k), 1e-5)
2738 duz = ( u(k)-u(k-1) )**2 +( v(k)-v(k-1) )**2
2741 ri = -gh(k)/max( duz, 1.0e-10 )
2742 IF (ckmod .eq. 1)
THEN
2743 a2fac = 1./(1. + max(ri,0.0))
2754 prnum = min(0.76 + 4.0*max(ri,0.0), prlimit)
2763 IF ( debug_code )
THEN
2764 IF (sh(k)<0.0 .OR. sm(k)<0.0)
THEN
2765 print*,
"MYNN; mym_turbulence 2.0; sh=",sh(k),
" k=",k
2766 print*,
" gm=",gm(k),
" gh=",gh(k),
" sm=",sm(k)
2767 print*,
" q2sq=",q2sq,
" q3sq=",q3sq,
" q3/q2=",q3sq/q2sq
2768 print*,
" qke=",qke(k),
" el=",el(k),
" ri=",ri
2769 print*,
" PBLH=",zi,
" u=",u(k),
" v=",v(k)
2778 IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k)
2780 IF ( q3sq .LT. q2sq )
THEN
2782 qdiv = sqrt( q3sq/q2sq )
2799 sh(k) = sh(k) * qdiv
2800 sm(k) = sm(k) * qdiv
2817 e1 = q3sq - e1c*ghel*a2fac * qdiv**2
2818 e2 = q3sq - e2c*ghel*a2fac * qdiv**2
2819 e3 = e1 + e3c*ghel*a2fac**2 * qdiv**2
2820 e4 = e1 - e4c*ghel*a2fac * qdiv**2
2821 eden = e2*e4 + e3*e5c*gmel * qdiv**2
2822 eden = max( eden, 1.0d-20 )
2833 e1 = q3sq - e1c*ghel*a2fac
2834 e2 = q3sq - e2c*ghel*a2fac
2835 e3 = e1 + e3c*ghel*a2fac**2
2836 e4 = e1 - e4c*ghel*a2fac
2837 eden = e2*e4 + e3*e5c*gmel
2838 eden = max( eden, 1.0d-20 )
2842 sm(k) = q3sq*a1*( e3-3.0*c1*e4 )/eden
2846 sh(k) = q3sq*(a2*a2fac)*( e2+3.0*c1*e5c*gmel )/eden
2856 gmelq = max(gmel/q3sq, 1e-8)
2864 IF ( debug_code )
THEN
2865 IF ((sh(k)<sh25min .OR. sm(k)<sm25min .OR. &
2866 sh(k)>sh25max .OR. sm(k)>sm25max) )
THEN
2867 print*,
"In mym_turbulence 2.5: k=",k
2868 print*,
" sm=",sm(k),
" sh=",sh(k)
2869 print*,
" ri=",ri,
" Pr=",sm(k)/max(sh(k),1e-8)
2870 print*,
" gm=",gm(k),
" gh=",gh(k)
2871 print*,
" q2sq=",q2sq,
" q3sq=",q3sq, q3sq/q2sq
2872 print*,
" qke=",qke(k),
" el=",el(k)
2873 print*,
" PBLH=",zi,
" u=",u(k),
" v=",v(k)
2874 print*,
" SMnum=",q3sq*a1*( e3-3.0*c1*e4),
" SMdenom=",eden
2875 print*,
" SHnum=",q3sq*(a2*a2fac)*( e2+3.0*c1*e5c*gmel ),&
2881 IF ( sh(k) > sh25max ) sh(k) = sh25max
2882 IF ( sh(k) < sh25min ) sh(k) = sh25min
2894 shb = max(sh(k), 0.002)
2895 sm(k) = min(sm(k), prlimit*shb)
2898 IF ( closure .GE. 3.0 )
THEN
2899 t2sq = qdiv*b2*elsq*sh(k)*dtl(k)**2
2900 r2sq = qdiv*b2*elsq*sh(k)*dqw(k)**2
2901 c2sq = qdiv*b2*elsq*sh(k)*dtl(k)*dqw(k)
2902 t3sq = max( tsq(k)*abk+tsq(k-1)*afk, 0.0 )
2903 r3sq = max( qsq(k)*abk+qsq(k-1)*afk, 0.0 )
2904 c3sq = cov(k)*abk+cov(k-1)*afk
2907 c3sq = sign( min( abs(c3sq), sqrt(t3sq*r3sq) ), c3sq )
2909 vtt = 1.0 +vt(k)*abk +vt(k-1)*afk
2910 vqq = tv0 +vq(k)*abk +vq(k-1)*afk
2912 t2sq = vtt*t2sq +vqq*c2sq
2913 r2sq = vtt*c2sq +vqq*r2sq
2914 c2sq = max( vtt*t2sq+vqq*r2sq, 0.0d0 )
2915 t3sq = vtt*t3sq +vqq*c3sq
2916 r3sq = vtt*c3sq +vqq*r3sq
2917 c3sq = max( vtt*t3sq+vqq*r3sq, 0.0d0 )
2919 cw25 = e1*( e2 + 3.0*c1*e5c*gmel*qdiv**2 )/( 3.0*eden )
2923 IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k)
2928 auh = 27.*a1*((a2*a2fac)**2)*b2*(gtr)**2
2929 aum = 54.*(a1**2)*(a2*a2fac)*b2*c1*(gtr)
2930 adh = 9.*a1*((a2*a2fac)**2)*(12.*a1 + 3.*b2)*(gtr)**2
2931 adm = 18.*(a1**2)*(a2*a2fac)*(b2 - 3.*(a2*a2fac))*(gtr)
2933 aeh = (9.*a1*((a2*a2fac)**2)*b1 +9.*a1*((a2*a2fac)**2)* &
2934 (12.*a1 + 3.*b2))*(gtr)
2935 aem = 3.*a1*(a2*a2fac)*b1*(3.*(a2*a2fac) + 3.*b2*c1 + &
2936 (18.*a1*c1 - b2)) + &
2937 (18.)*(a1**2)*(a2*a2fac)*(b2 - 3.*(a2*a2fac))
2940 rsl = (auh + aum*req)/(3.*adh + 3.*adm*req)
2955 e2 = q3sq - e2c*ghel*a2fac * qdiv**2
2956 e3 = q3sq + e3c*ghel*a2fac**2 * qdiv**2
2957 e4 = q3sq - e4c*ghel*a2fac * qdiv**2
2958 eden = e2*e4 + e3 *e5c*gmel * qdiv**2
2963 wden = cc3*gtr**2 * dlsq**2/elsq * qdiv**2 &
2964 & *( e2*e4c*a2fac - e3c*e5c*gmel*a2fac**2 * qdiv**2 )
2966 IF ( wden .NE. 0.0 )
THEN
2968 clow = q3sq*( 0.12-cw25 )*eden/wden
2969 cupp = q3sq*( 0.76-cw25 )*eden/wden
2973 IF ( wden .GT. 0.0 )
THEN
2974 c3sq = min( max( c3sq, c2sq+clow ), c2sq+cupp )
2976 c3sq = max( min( c3sq, c2sq+clow ), c2sq+cupp )
2980 e1 = e2 + e5c*gmel * qdiv**2
2981 eden = max( eden, 1.0d-20 )
2986 e6c = 3.0*(a2*a2fac)*cc3*gtr * dlsq/elsq
2991 IF ( t2sq .GE. 0.0 )
THEN
2992 enum = max( qdiv*e6c*( t3sq-t2sq ), 0.0d0 )
2994 enum = min( qdiv*e6c*( t3sq-t2sq ), 0.0d0 )
2996 gamt =-e1 *
enum /eden
3001 IF ( r2sq .GE. 0.0 )
THEN
3002 enum = max( qdiv*e6c*( r3sq-r2sq ), 0.0d0 )
3004 enum = min( qdiv*e6c*( r3sq-r2sq ), 0.0d0 )
3006 gamq =-e1 *
enum /eden
3011 enum = max( qdiv*e6c*( c3sq-c2sq ), 0.0d0)
3015 smd = dlsq*enum*gtr/eden * qdiv**2 * (e3c*a2fac**2 + &
3016 & e4c*a2fac)*a1/(a2*a2fac)
3018 gamv = e1 *enum*gtr/eden
3026 IF ( debug_code )
THEN
3027 IF (sh(k)<-0.3 .OR. sm(k)<-0.3 .OR. &
3028 qke(k) < -0.1 .or. abs(smd) .gt. 2.0)
THEN
3029 print*,
" MYNN; mym_turbulence3.0; sh=",sh(k),
" k=",k
3030 print*,
" gm=",gm(k),
" gh=",gh(k),
" sm=",sm(k)
3031 print*,
" q2sq=",q2sq,
" q3sq=",q3sq,
" q3/q2=",q3sq/q2sq
3032 print*,
" qke=",qke(k),
" el=",el(k),
" ri=",ri
3033 print*,
" PBLH=",zi,
" u=",u(k),
" v=",v(k)
3048 cldavg = 0.5*(cldfra_bl1d(k-1) + cldfra_bl1d(k))
3049 IF (edmf_a1(k) > 0.001 .OR. cldavg > 0.02)
THEN
3051 sm(k) = max(sm(k), 0.03*min(10.*edmf_a1(k)*edmf_w1(k),1.0) )
3052 sh(k) = max(sh(k), 0.03*min(10.*edmf_a1(k)*edmf_w1(k),1.0) )
3054 sm(k) = max(sm(k), 0.05*min(cldavg,1.0) )
3055 sh(k) = max(sh(k), 0.05*min(cldavg,1.0) )
3063 pdk(k) = elq*( sm(k)*gm(k) &
3064 & +sh(k)*gh(k)+gamv ) + &
3066 pdt(k) = elh*( sh(k)*dtl(k)+gamt )*dtl(k)
3067 pdq(k) = elh*( sh(k)*dqw(k)+gamq )*dqw(k)
3068 pdc(k) = elh*( sh(k)*dtl(k)+gamt ) &
3070 & + elh*( sh(k)*dqw(k)+gamq )*dtl(k)*0.5
3077 dfm(k) = elq*sm(k) / dzk
3078 dfh(k) = elq*sh(k) / dzk
3085 IF (tke_budget .eq. 1)
THEN
3101 qshear1d(k) = elq*sm(k)*gm(k)
3109 qbuoy1d(k) = elq*(sh(k)*gh(k)+gamv)+0.5*tkeprodtd(k)
3132 tcd(k) = ( tcd(k+1)-tcd(k) )/( dzk )
3133 qcd(k) = ( qcd(k+1)-qcd(k) )/( dzk )
3136 if (spp_pbl==1)
then
3138 dfm(k)= dfm(k) + dfm(k)* rstoch_col(k) * 1.5 * max(exp(-max(zw(k)-8000.,0.0)/2000.),0.001)
3139 dfh(k)= dfh(k) + dfh(k)* rstoch_col(k) * 1.5 * max(exp(-max(zw(k)-8000.,0.0)/2000.),0.001)
3144#ifdef HARDCODE_VERTICAL
3201 & ust, flt, flq, pmz, phh, &
3203 & pdk, pdt, pdq, pdc, &
3204 & qke, tsq, qsq, cov, &
3205 & s_aw,s_awqke,bl_mynn_edmf_tke, &
3206 & qWT1D, qDISS1D,tke_budget)
3209 integer,
intent(in) :: kts,kte
3211#ifdef HARDCODE_VERTICAL
3213# define kte HARDCODE_VERTICAL
3216 real(kind_phys),
intent(in) :: closure
3217 integer,
intent(in) :: bl_mynn_edmf_tke,tke_budget
3218 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz, dfq, el, rho
3219 real(kind_phys),
dimension(kts:kte),
intent(inout) :: pdk, pdt, pdq, pdc
3220 real(kind_phys),
intent(in) :: flt, flq, pmz, phh
3221 real(kind_phys),
intent(in) :: ust, delt
3222 real(kind_phys),
dimension(kts:kte),
intent(inout) :: qke,tsq, qsq, cov
3224 real(kind_phys),
dimension(kts:kte+1),
intent(inout) :: s_awqke,s_aw
3227 real(kind_phys),
dimension(kts:kte),
intent(out) :: qWT1D, qDISS1D
3228 real(kind_phys),
dimension(kts:kte) :: tke_up,dzinv
3232 real(kind_phys),
dimension(kts:kte) :: qkw, bp, rp, df3q
3233 real(kind_phys):: vkz,pdk1,phm,pdt1,pdq1,pdc1,b1l,b2l,onoff
3234 real(kind_phys),
dimension(kts:kte) :: dtz
3235 real(kind_phys),
dimension(kts:kte) :: a,b,c,d,x
3237 real(kind_phys),
dimension(kts:kte) :: rhoinv
3238 real(kind_phys),
dimension(kts:kte+1) :: rhoz,kqdz,kmdz
3241 IF (bl_mynn_edmf_tke == 0)
THEN
3248 vkz = karman*0.5*dz(kts)
3254 qkw(k) = sqrt( max( qke(k), 0.0 ) )
3255 df3q(k)=sqfac*dfq(k)
3263 rhoinv(kts)=1./rho(kts)
3264 kqdz(kts) =rhoz(kts)*df3q(kts)
3265 kmdz(kts) =rhoz(kts)*dfq(kts)
3267 rhoz(k) =(rho(k)*dz(k-1) + rho(k-1)*dz(k))/(dz(k-1)+dz(k))
3268 rhoz(k) = max(rhoz(k),1e-4)
3269 rhoinv(k)=1./max(rho(k),1e-4)
3270 kqdz(k) = rhoz(k)*df3q(k)
3271 kmdz(k) = rhoz(k)*dfq(k)
3273 rhoz(kte+1)=rhoz(kte)
3274 kqdz(kte+1)=rhoz(kte+1)*df3q(kte)
3275 kmdz(kte+1)=rhoz(kte+1)*dfq(kte)
3279 kqdz(k) = max(kqdz(k), 0.5* s_aw(k))
3280 kqdz(k) = max(kqdz(k), -0.5*(s_aw(k)-s_aw(k+1)))
3281 kmdz(k) = max(kmdz(k), 0.5* s_aw(k))
3282 kmdz(k) = max(kmdz(k), -0.5*(s_aw(k)-s_aw(k+1)))
3286 pdk1 = 2.0*ust**3*pmz/( vkz )
3287 phm = 2.0/ust *phh/( vkz )
3293 pdk(kts) = pdk1 - pdk(kts+1)
3298 pdt(kts) = pdt(kts+1)
3299 pdq(kts) = pdq(kts+1)
3300 pdc(kts) = pdc(kts+1)
3305 b1l = b1*0.5*( el(k+1)+el(k) )
3306 bp(k) = 2.*qkw(k) / b1l
3307 rp(k) = pdk(k+1) + pdk(k)
3328 a(k)= - dtz(k)*kqdz(k)*rhoinv(k) &
3329 & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff
3330 b(k)=1. + dtz(k)*(kqdz(k)+kqdz(k+1))*rhoinv(k) &
3331 & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff &
3333 c(k)= - dtz(k)*kqdz(k+1)*rhoinv(k) &
3334 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff
3335 d(k)=rp(k)*delt + qke(k) &
3336 & + dtz(k)*rhoinv(k)*(s_awqke(k)-s_awqke(k+1))*onoff
3358 CALL tridiag2(kte,a,b,c,d,x)
3362 qke(k)=max(x(k), qkemin)
3363 qke(k)=min(qke(k), 150.)
3368 IF (tke_budget .eq. 1)
THEN
3373 qwt1d(k)=dzinv(k)*( &
3374 & (kqdz(k+1)*(tke_up(k+1)-tke_up(k))-kqdz(k)*tke_up(k)) &
3375 & + 0.5*rhoinv(k)*(s_aw(k+1)*tke_up(k+1) &
3376 & + (s_aw(k+1)-s_aw(k))*tke_up(k) &
3377 & + (s_awqke(k)-s_awqke(k+1)))*onoff)
3379 qwt1d(k)=dzinv(k)*( &
3380 & (kqdz(k+1)*(tke_up(k+1)-tke_up(k))-kqdz(k)*(tke_up(k)-tke_up(k-1))) &
3381 & + 0.5*rhoinv(k)*(s_aw(k+1)*tke_up(k+1) &
3382 & + (s_aw(k+1)-s_aw(k))*tke_up(k) &
3383 & - s_aw(k)*tke_up(k-1) &
3384 & + (s_awqke(k)-s_awqke(k+1)))*onoff)
3387 qwt1d(k)=dzinv(k)*(-kqdz(k)*(tke_up(k)-tke_up(k-1)) &
3388 & + 0.5*rhoinv(k)*(-s_aw(k)*tke_up(k)-s_aw(k)*tke_up(k-1)+s_awqke(k))*onoff)
3394 IF ( closure > 2.5 )
THEN
3398 b2l = b2*0.5*( el(k+1)+el(k) )
3399 bp(k) = 2.*qkw(k) / b2l
3400 rp(k) = pdq(k+1) + pdq(k)
3411 a(k)= - dtz(k)*kmdz(k)*rhoinv(k)
3412 b(k)=1. + dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + bp(k)*delt
3413 c(k)= - dtz(k)*kmdz(k+1)*rhoinv(k)
3414 d(k)=rp(k)*delt + qsq(k)
3423 CALL tridiag2(kte,a,b,c,d,x)
3427 qsq(k)=max(x(k),1e-17)
3432 IF ( qkw(k) .LE. 0.0 )
THEN
3435 b2l = b2*0.25*( el(k+1)+el(k) )/qkw(k)
3437 qsq(k) = b2l*( pdq(k+1)+pdq(k) )
3443 IF ( closure .GE. 3.0 )
THEN
3450 b2l = b2*0.5*( el(k+1)+el(k) )
3451 bp(k) = 2.*qkw(k) / b2l
3452 rp(k) = pdt(k+1) + pdt(k)
3469 a(k)= - dtz(k)*kmdz(k)*rhoinv(k)
3470 b(k)=1. + dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + bp(k)*delt
3471 c(k)= - dtz(k)*kmdz(k+1)*rhoinv(k)
3472 d(k)=rp(k)*delt + tsq(k)
3488 CALL tridiag2(kte,a,b,c,d,x)
3498 b2l = b2*0.5*( el(k+1)+el(k) )
3499 bp(k) = 2.*qkw(k) / b2l
3500 rp(k) = pdc(k+1) + pdc(k)
3517 a(k)= - dtz(k)*kmdz(k)*rhoinv(k)
3518 b(k)=1. + dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + bp(k)*delt
3519 c(k)= - dtz(k)*kmdz(k+1)*rhoinv(k)
3520 d(k)=rp(k)*delt + cov(k)
3536 CALL tridiag2(kte,a,b,c,d,x)
3547 IF ( qkw(k) .LE. 0.0 )
THEN
3550 b2l = b2*0.25*( el(k+1)+el(k) )/qkw(k)
3553 tsq(k) = b2l*( pdt(k+1)+pdt(k) )
3554 cov(k) = b2l*( pdc(k+1)+pdc(k) )
3562#ifdef HARDCODE_VERTICAL
3604 & dx, dz, zw, xland, &
3605 & thl, qw, qv, qc, qi, qs, &
3608 & Sh, el, bl_mynn_cloudpdf, &
3609 & qc_bl1D, qi_bl1D, &
3612 & Vt, Vq, th, sgm, rmo, &
3613 & spp_pbl,rstoch_col )
3617 integer,
intent(in) :: kts,kte, bl_mynn_cloudpdf
3619#ifdef HARDCODE_VERTICAL
3621# define kte HARDCODE_VERTICAL
3624 real(kind_phys),
intent(in) :: HFX1,rmo,xland
3625 real(kind_phys),
intent(in) :: dx,pblh1
3626 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz
3627 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
3628 real(kind_phys),
dimension(kts:kte),
intent(in) :: p,exner,thl,qw, &
3629 &qv,qc,qi,qs,tsq,qsq,cov,th
3631 real(kind_phys),
dimension(kts:kte),
intent(inout) :: vt,vq,sgm
3633 real(kind_phys),
dimension(kts:kte) :: alp,a,bet,b,ql,q1,RH
3634 real(kind_phys),
dimension(kts:kte),
intent(out) :: qc_bl1D,qi_bl1D, &
3636 DOUBLE PRECISION :: t3sq, r3sq, c3sq
3638 real(kind_phys):: qsl,esat,qsat,dqsl,cld0,q1k,qlk,eq1,qll, &
3639 &q2p,pt,rac,qt,t,xl,rsl,cpm,Fng,qww,alpha,beta,bb, &
3640 &ls,wt,wt2,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, &
3641 &qmq,qsat_tk,q1_rh,rh_hack,dzm1,zsl,maxqc
3642 real(kind_phys),
parameter :: qpct_sfc=0.025
3643 real(kind_phys),
parameter :: qpct_pbl=0.030
3644 real(kind_phys),
parameter :: qpct_trp=0.040
3645 real(kind_phys),
parameter :: rhcrit =0.83
3646 real(kind_phys),
parameter :: rhmax =1.02
3649 real(kind_phys):: erf
3652 real:: dth,dtl,dqw,dzk,els
3653 real(kind_phys),
dimension(kts:kte),
intent(in) :: Sh,el
3656 real(kind_phys) :: zagl,damp,PBLH2
3657 real(kind_phys) :: cfmax
3660 real(kind_phys) :: theta1, theta2, ht1, ht2
3664 integer,
intent(in) :: spp_pbl
3665 real(kind_phys),
dimension(kts:kte) :: rstoch_col
3666 real(kind_phys) :: qw_pert
3673 DO k = kte-3, kts, -1
3676 ht1 = 44307.692 * (1.0 - (p(k)/101325.)**0.190)
3677 ht2 = 44307.692 * (1.0 - (p(k+2)/101325.)**0.190)
3678 if ( (((theta2-theta1)/(ht2-ht1)) .lt. 10./1500. ) .AND. &
3679 & (ht1.lt.19000.) .and. (ht1.gt.4000.) )
then
3683 k_tropo = max(kts+2, k+2)
3687 SELECT CASE(bl_mynn_cloudpdf)
3708 qsl=ep_2*esat/max(1.e-4,(p(k)-ep_3*esat))
3710 dqsl = qsl*ep_2*xlv/( r_d*t**2 )
3712 alp(k) = 1.0/( 1.0+dqsl*xlvcp )
3713 bet(k) = dqsl*exner(k)
3717 t3sq = max( tsq(k), 0.0 )
3718 r3sq = max( qsq(k), 0.0 )
3720 c3sq = sign( min( abs(c3sq), sqrt(t3sq*r3sq) ), c3sq )
3721 r3sq = r3sq +bet(k)**2*t3sq -2.0*bet(k)*c3sq
3725 sgm(k) = sqrt( max( r3sq, 1.0d-10 ))
3727 q1(k) = qmq / sgm(k)
3729 cldfra_bl1d(k) = 0.5*( 1.0+erf( q1(k)*rr2 ) )
3732 eq1 = rrp*exp( -0.5*q1k*q1k )
3733 qll = max( cldfra_bl1d(k)*q1k + eq1, 0.0 )
3735 ql(k) = alp(k)*sgm(k)*qll
3737 liq_frac = min(1.0, max(0.0,(t-240.0)/29.0))
3738 qc_bl1d(k) = liq_frac*ql(k)
3739 qi_bl1d(k) = (1.0 - liq_frac)*ql(k)
3742 q2p = xlvcp/exner(k)
3743 pt = thl(k) +q2p*ql(k)
3746 qt = 1.0 +p608*qw(k) -(1.+p608)*(qc_bl1d(k)+qi_bl1d(k))*cldfra_bl1d(k)
3747 rac = alp(k)*( cldfra_bl1d(k)-qll*eq1 )*( q2p*qt-(1.+p608)*pt )
3752 vt(k) = qt-1.0 -rac*bet(k)
3753 vq(k) = p608*pt-tv0 +rac
3765 qsl=ep_2*esat/max(1.e-4,(p(k)-ep_3*esat))
3767 dqsl = qsl*ep_2*xlv/( r_d*t**2 )
3769 alp(k) = 1.0/( 1.0+dqsl*xlvcp )
3770 bet(k) = dqsl*exner(k)
3772 if (k .eq. kts)
then
3777 dth = 0.5*(thl(k+1)+thl(k)) - 0.5*(thl(k)+thl(max(k-1,kts)))
3778 dqw = 0.5*(qw(k+1) + qw(k)) - 0.5*(qw(k) + qw(max(k-1,kts)))
3779 sgm(k) = sqrt( max( (alp(k)**2 * max(el(k)**2,0.1) * &
3780 b2 * max(sh(k),0.03))/4. * &
3781 (dqw/dzk - bet(k)*(dth/dzk ))**2 , 1.0e-10) )
3783 q1(k) = qmq / sgm(k)
3784 cldfra_bl1d(k) = 0.5*( 1.0+erf( q1(k)*rr2 ) )
3790 eq1 = rrp*exp( -0.5*q1k*q1k )
3791 qll = max( cldfra_bl1d(k)*q1k + eq1, 0.0 )
3793 ql(k) = alp(k)*sgm(k)*qll
3794 liq_frac = min(1.0, max(0.0,(t-240.0)/29.0))
3795 qc_bl1d(k) = liq_frac*ql(k)
3796 qi_bl1d(k) = (1.0 - liq_frac)*ql(k)
3799 q2p = xlvcp/exner(k)
3800 pt = thl(k) +q2p*ql(k)
3803 qt = 1.0 +p608*qw(k) -(1.+p608)*(qc_bl1d(k)+qi_bl1d(k))*cldfra_bl1d(k)
3804 rac = alp(k)*( cldfra_bl1d(k)-qll*eq1 )*( q2p*qt-(1.+p608)*pt )
3809 vt(k) = qt-1.0 -rac*bet(k)
3810 vq(k) = p608*pt-tv0 +rac
3818 pblh2=max(10._kind_phys,pblh1)
3822 zagl = zagl + 0.5*(dz(k) + dzm1)
3828 rh(k) = max(min(rhmax, qw(k)/max(1.e-10,qsat_tk)),0.001_kind_phys)
3831 dqsl = qsat_tk*ep_2*xlv/( r_d*t**2 )
3832 alp(k) = 1.0/( 1.0+dqsl*xlvcp )
3833 bet(k) = dqsl*exner(k)
3835 rsl = xl*qsat_tk / (r_v*t**2)
3837 cpm = cp + qw(k)*cpv
3838 a(k) = 1./(1. + xl*rsl/cpm)
3842 qw_pert= qw(k) + qw(k)*0.5*rstoch_col(k)*real(spp_pbl)
3845 qmq = qw_pert - qsat_tk
3849 r3sq = max( qsq(k), 0.0 )
3851 sgm(k) = sqrt( r3sq )
3853 sgm(k) = min( sgm(k), qsat_tk*0.666 )
3857 wt = max(500. - max(dz(k)-100.,0.0), 0.0_kind_phys)/500.
3858 sgm(k) = sgm(k) + sgm(k)*0.2*(1.0-wt)
3861 qpct = qpct_pbl*wt + qpct_trp*(1.0-wt)
3862 qpct = min(qpct, max(qpct_sfc, qpct_pbl*zagl/500.) )
3863 sgm(k) = max( sgm(k), qsat_tk*qpct )
3865 q1(k) = qmq / sgm(k)
3870 wt2 = min(max( zagl - pblh2, 0.0 )/300., 1.0)
3872 if ((qi(k)+qs(k))>1.e-9 .and. (zagl .gt. pblh2))
then
3873 rh_hack =min(rhmax, rhcrit + wt2*0.045*(9.0 + log10(qi(k)+qs(k))))
3874 rh(k) =max(rh(k), rh_hack)
3876 q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
3877 q1(k) =max(q1_rh, q1(k) )
3880 if (qc(k)>1.e-6 .and. (zagl .gt. pblh2))
then
3881 rh_hack =min(rhmax, rhcrit + wt2*0.08*(6.0 + log10(qc(k))))
3882 rh(k) =max(rh(k), rh_hack)
3884 q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
3885 q1(k) =max(q1_rh, q1(k) )
3896 cldfra_bl1d(k) = max(0., min(1., 0.5+0.36*atan(1.8*(q1(k)+0.2))))
3901 maxqc = max(qw(k) - qsat_tk, 0.0)
3903 ql_water = sgm(k)*exp(1.2*q1k-1.)
3904 ql_ice = sgm(k)*exp(1.2*q1k-1.)
3905 elseif (q1k > 2.)
then
3906 ql_water = min(sgm(k)*q1k, maxqc)
3909 ql_water = min(sgm(k)*(exp(-1.) + 0.66*q1k + 0.086*q1k**2), maxqc)
3910 ql_ice = sgm(k)*(exp(-1.) + 0.66*q1k + 0.086*q1k**2)
3918 if (cldfra_bl1d(k) < 0.001)
then
3921 cldfra_bl1d(k) = 0.0
3924 liq_frac = min(1.0, max(0.0, (t-tice)/(tliq-tice)))
3925 qc_bl1d(k) = liq_frac*ql_water
3926 qi_bl1d(k) = (1.0-liq_frac)*ql_ice
3930 if (k .ge. k_tropo)
then
3939 if ((xland-1.5).GE.0)
then
3956 if (q1k .ge. 1.0)
then
3958 elseif (q1k .ge. -1.7 .and. q1k .lt. 1.0)
then
3959 fng = exp(-0.4*(q1k-1.0))
3960 elseif (q1k .ge. -2.5 .and. q1k .lt. -1.7)
then
3961 fng = 3.0 + exp(-3.8*(q1k+1.7))
3963 fng = min(23.9 + exp(-1.6*(q1k+2.5)), 60._kind_phys)
3966 cfmax = min(cldfra_bl1d(k), 0.6_kind_phys)
3968 zsl = min(max(25., 0.1*pblh2), 100.)
3969 wt = min(zagl/zsl, 1.0)
3980 beta = (th(k)/t)*(xl/cp) - 1.61*th(k)
3981 vt(k) = qww - cfmax*beta*bb*fng - 1.
3982 vq(k) = alpha + cfmax*beta*a(k)*fng - tv0
3990 fac_damp = min(zagl * 0.0025, 1.0)
3993 cld_factor = 1.0 + fac_damp*min((max(0.0, ( rh(k) - 0.92 )) / 0.145)**2, 0.37)
3994 cldfra_bl1d(k) = min( 1., cld_factor*cldfra_bl1d(k) )
4000 IF (bl_mynn_cloudpdf .LT. 0)
THEN
4002 cldfra_bl1d(k) = 0.0
4016#ifdef HARDCODE_VERTICAL
4029 &u,v,th,tk,qv,qc,qi,qs,qnc,qni, &
4031 &thl,sqv,sqc,sqi,sqs,sqw, &
4032 &qnwfa,qnifa,qnbca,ozone, &
4033 &ust,flt,flq,flqv,flqc,wspd, &
4038 &Du,Dv,Dth,Dqv,Dqc,Dqi,Dqs,Dqnc,Dqni, &
4039 &Dqnwfa,Dqnifa,Dqnbca,Dozone, &
4041 &s_aw,s_awthl,s_awqt,s_awqv,s_awqc, &
4044 &s_awqnwfa,s_awqnifa,s_awqnbca, &
4045 &sd_aw,sd_awthl,sd_awqt,sd_awqv, &
4046 &sd_awqc,sd_awu,sd_awv, &
4049 &det_thl,det_sqv,det_sqc, &
4051 &FLAG_QC,FLAG_QI,FLAG_QNC,FLAG_QNI, &
4053 &FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA, &
4055 &bl_mynn_cloudmix, &
4058 &bl_mynn_edmf_mom, &
4059 &bl_mynn_mixscalars )
4062 integer,
intent(in) :: kts,kte,i
4064#ifdef HARDCODE_VERTICAL
4066# define kte HARDCODE_VERTICAL
4069 integer,
intent(in) :: bl_mynn_cloudmix,bl_mynn_mixqt, &
4070 bl_mynn_edmf,bl_mynn_edmf_mom, &
4072 logical,
intent(IN) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QS, &
4073 &FLAG_QNC,FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA
4082 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: s_aw, &
4083 &s_awthl,s_awqt,s_awqnc,s_awqni,s_awqv,s_awqc,s_awu,s_awv, &
4084 &s_awqnwfa,s_awqnifa,s_awqnbca, &
4085 &sd_aw,sd_awthl,sd_awqt,sd_awqv,sd_awqc,sd_awu,sd_awv
4087 real(kind_phys),
dimension(kts:kte),
intent(in) :: sub_thl,sub_sqv, &
4088 &sub_u,sub_v,det_thl,det_sqv,det_sqc,det_u,det_v
4089 real(kind_phys),
dimension(kts:kte),
intent(in) :: u,v,th,tk,qv,qc,qi,&
4090 &qs,qni,qnc,rho,p,exner,dfq,dz,tsq,qsq,cov,tcd,qcd, &
4091 &cldfra_bl1d,diss_heat
4092 real(kind_phys),
dimension(kts:kte),
intent(inout) :: thl,sqw,sqv,sqc,&
4093 &sqi,sqs,qnwfa,qnifa,qnbca,ozone,dfm,dfh
4094 real(kind_phys),
dimension(kts:kte),
intent(inout) :: du,dv,dth,dqv, &
4095 &dqc,dqi,dqs,dqni,dqnc,dqnwfa,dqnifa,dqnbca,dozone
4096 real(kind_phys),
intent(in) :: flt,flq,flqv,flqc,uoce,voce
4097 real(kind_phys),
intent(in) :: ust,delt,psfc,wspd
4099 real(kind_phys):: wsp,wsp2,tk2,th2
4107 real(kind_phys),
dimension(kts:kte) :: dtz,dfhc,dfmc,delp
4108 real(kind_phys),
dimension(kts:kte) :: sqv2,sqc2,sqi2,sqs2,sqw2, &
4109 &qni2,qnc2,qnwfa2,qnifa2,qnbca2,ozone2
4110 real(kind_phys),
dimension(kts:kte) :: zfac,plumeKh,rhoinv
4111 real(kind_phys),
dimension(kts:kte) :: a,b,c,d,x
4112 real(kind_phys),
dimension(kts:kte+1) :: rhoz, &
4114 real(kind_phys):: rhs,gfluxm,gfluxp,dztop,maxdfh,mindfh,maxcf,maxKh,zw
4115 real(kind_phys):: t,esat,qsl,onoff,kh,km,dzk,rhosfc
4116 real(kind_phys):: ustdrag,ustdiff,qvflux
4117 real(kind_phys):: th_new,portion_qc,portion_qi,condensate,qsat
4122 real(kind_phys),
parameter :: nonloc = 1.0
4124 dztop=.5*(dz(kte)+dz(kte-1))
4129 IF (bl_mynn_edmf_mom == 0)
THEN
4137 rhosfc = psfc/(r_d*(tk(kts)+p608*qv(kts)))
4138 dtz(kts) =delt/dz(kts)
4140 rhoinv(kts)=1./rho(kts)
4141 khdz(kts) =rhoz(kts)*dfh(kts)
4142 kmdz(kts) =rhoz(kts)*dfm(kts)
4143 delp(kts) = psfc - (p(kts+1)*dz(kts) + p(kts)*dz(kts+1))/(dz(kts)+dz(kts+1))
4146 rhoz(k) =(rho(k)*dz(k-1) + rho(k-1)*dz(k))/(dz(k-1)+dz(k))
4147 rhoz(k) = max(rhoz(k),1e-4)
4148 rhoinv(k)=1./max(rho(k),1e-4)
4149 dzk = 0.5 *( dz(k)+dz(k-1) )
4150 khdz(k) = rhoz(k)*dfh(k)
4151 kmdz(k) = rhoz(k)*dfm(k)
4154 delp(k) = (p(k)*dz(k-1) + p(k-1)*dz(k))/(dz(k)+dz(k-1)) - &
4155 (p(k+1)*dz(k) + p(k)*dz(k+1))/(dz(k)+dz(k+1))
4157 delp(kte) =delp(kte-1)
4158 rhoz(kte+1)=rhoz(kte)
4159 khdz(kte+1)=rhoz(kte+1)*dfh(kte)
4160 kmdz(kte+1)=rhoz(kte+1)*dfm(kte)
4164 khdz(k) = max(khdz(k), 0.5*s_aw(k))
4165 khdz(k) = max(khdz(k), -0.5*(s_aw(k)-s_aw(k+1)))
4166 kmdz(k) = max(kmdz(k), 0.5*s_aw(k))
4167 kmdz(k) = max(kmdz(k), -0.5*(s_aw(k)-s_aw(k+1)))
4170 ustdrag = min(ust*ust,0.99)/wspd
4171 ustdiff = min(ust*ust,0.01)/wspd
4181 a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
4182 b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) &
4183 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4184 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4185 c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
4186 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4187 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4188 d(k)=u(k) + dtz(k)*uoce*ust**2/wspd &
4189 & - dtz(k)*rhoinv(k)*s_awu(k+1)*onoff &
4190 & + dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff &
4191 & + sub_u(k)*delt + det_u(k)*delt
4194 a(k)= -dtz(k)*kmdz(k)*rhoinv(k) &
4195 & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff &
4196 & + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
4197 b(k)=1.+ dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) &
4198 & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff &
4199 & + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
4200 c(k)= - dtz(k)*kmdz(k+1)*rhoinv(k) &
4201 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4202 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4203 d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff &
4204 & - dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff &
4205 & + sub_u(k)*delt + det_u(k)*delt
4227 CALL tridiag2(kte,a,b,c,d,x)
4232 du(k)=(x(k)-u(k))/delt
4242 a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
4243 b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) &
4244 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4245 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4246 c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
4247 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4248 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4249 d(k)=v(k) + dtz(k)*voce*ust**2/wspd &
4250 & - dtz(k)*rhoinv(k)*s_awv(k+1)*onoff &
4251 & + dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff &
4252 & + sub_v(k)*delt + det_v(k)*delt
4255 a(k)= -dtz(k)*kmdz(k)*rhoinv(k) &
4256 & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff &
4257 & + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
4258 b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) &
4259 & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff &
4260 & + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
4261 c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
4262 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4263 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4264 d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff &
4265 & - dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff &
4266 & + sub_v(k)*delt + det_v(k)*delt
4288 CALL tridiag2(kte,a,b,c,d,x)
4293 dv(k)=(x(k)-v(k))/delt
4318 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4319 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4320 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4321 d(k)=thl(k) + dtz(k)*rhosfc*flt*rhoinv(k) + tcd(k)*delt &
4322 & - dtz(k)*rhoinv(k)*s_awthl(k+1) -dtz(k)*rhoinv(k)*sd_awthl(k+1) + &
4323 & diss_heat(k)*delt + sub_thl(k)*delt + det_thl(k)*delt
4326 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k) + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)
4327 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4328 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1)) + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))
4329 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4330 d(k)=thl(k) + tcd(k)*delt + &
4331 & dtz(k)*rhoinv(k)*(s_awthl(k)-s_awthl(k+1)) + dtz(k)*rhoinv(k)*(sd_awthl(k)-sd_awthl(k+1)) + &
4332 & diss_heat(k)*delt + &
4333 & sub_thl(k)*delt + det_thl(k)*delt
4356 CALL tridiag2(kte,a,b,c,d,x)
4364IF (bl_mynn_mixqt > 0)
THEN
4388 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4389 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4390 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4391 d(k)=sqw(k) + dtz(k)*rhosfc*flq*rhoinv(k) + qcd(k)*delt - dtz(k)*rhoinv(k)*s_awqt(k+1) - dtz(k)*rhoinv(k)*sd_awqt(k+1)
4394 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k) + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)
4395 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4396 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1)) + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))
4397 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4398 d(k)=sqw(k) + qcd(k)*delt + dtz(k)*rhoinv(k)*(s_awqt(k)-s_awqt(k+1)) + dtz(k)*rhoinv(k)*(sd_awqt(k)-sd_awqt(k+1))
4419 CALL tridiag2(kte,a,b,c,d,sqw2)
4429IF (bl_mynn_mixqt == 0)
THEN
4434 IF (bl_mynn_cloudmix > 0 .AND. flag_qc)
THEN
4453 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4454 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4455 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4456 d(k)=sqc(k) + dtz(k)*rhosfc*flqc*rhoinv(k) + qcd(k)*delt &
4457 & - dtz(k)*rhoinv(k)*s_awqc(k+1) - dtz(k)*rhoinv(k)*sd_awqc(k+1) + &
4461 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k) + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)
4462 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4463 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1)) + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))
4464 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4465 d(k)=sqc(k) + qcd(k)*delt + dtz(k)*rhoinv(k)*(s_awqc(k)-s_awqc(k+1)) + dtz(k)*rhoinv(k)*(sd_awqc(k)-sd_awqc(k+1)) + &
4476 CALL tridiag2(kte,a,b,c,d,sqc2)
4488IF (bl_mynn_mixqt == 0)
THEN
4512 if (qvflux < 0.0)
then
4514 qvflux = max(qvflux, (min(0.9*sqv(kts) - 1e-8, 0.0)/dtz(kts)))
4518 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4519 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4520 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4521 d(k)=sqv(k) + dtz(k)*rhosfc*qvflux*rhoinv(k) + qcd(k)*delt &
4522 & - dtz(k)*rhoinv(k)*s_awqv(k+1) - dtz(k)*rhoinv(k)*sd_awqv(k+1) + &
4523 & sub_sqv(k)*delt + det_sqv(k)*delt
4526 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k) + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)
4527 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4528 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1)) + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))
4529 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1) - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)
4530 d(k)=sqv(k) + qcd(k)*delt + dtz(k)*rhoinv(k)*(s_awqv(k)-s_awqv(k+1)) + dtz(k)*rhoinv(k)*(sd_awqv(k)-sd_awqv(k+1)) + &
4531 & sub_sqv(k)*delt + det_sqv(k)*delt
4554 CALL tridiag2(kte,a,b,c,d,sqv2)
4567IF (bl_mynn_cloudmix > 0 .AND. flag_qi)
THEN
4571 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4572 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k)
4573 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4577 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4578 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k)
4579 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4603 CALL tridiag2(kte,a,b,c,d,sqi2)
4617IF (bl_mynn_cloudmix > 0 .AND. .false.)
THEN
4621 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4622 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k)
4623 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4627 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4628 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k)
4629 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4640 CALL tridiag2(kte,a,b,c,d,sqs2)
4653IF (bl_mynn_cloudmix > 0 .AND. flag_qni .AND. &
4654 bl_mynn_mixscalars > 0)
THEN
4658 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4659 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4660 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4661 d(k)=qni(k) - dtz(k)*rhoinv(k)*s_awqni(k+1)*nonloc
4664 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4665 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4666 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4667 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4668 d(k)=qni(k) + dtz(k)*rhoinv(k)*(s_awqni(k)-s_awqni(k+1))*nonloc
4678 CALL tridiag2(kte,a,b,c,d,x)
4694 IF (bl_mynn_cloudmix > 0 .AND. flag_qnc .AND. &
4695 bl_mynn_mixscalars > 0)
THEN
4699 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4700 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4701 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4702 d(k)=qnc(k) - dtz(k)*rhoinv(k)*s_awqnc(k+1)*nonloc
4705 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4706 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4707 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4708 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4709 d(k)=qnc(k) + dtz(k)*rhoinv(k)*(s_awqnc(k)-s_awqnc(k+1))*nonloc
4719 CALL tridiag2(kte,a,b,c,d,x)
4734IF (bl_mynn_cloudmix > 0 .AND. flag_qnwfa .AND. &
4735 bl_mynn_mixscalars > 0)
THEN
4739 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4740 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) - &
4741 & 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4742 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4743 d(k)=qnwfa(k) - dtz(k)*rhoinv(k)*s_awqnwfa(k+1)*nonloc
4746 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4747 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) + &
4748 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4749 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4750 d(k)=qnwfa(k) + dtz(k)*rhoinv(k)*(s_awqnwfa(k)-s_awqnwfa(k+1))*nonloc
4760 CALL tridiag2(kte,a,b,c,d,x)
4776IF (bl_mynn_cloudmix > 0 .AND. flag_qnifa .AND. &
4777 bl_mynn_mixscalars > 0)
THEN
4781 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4782 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) - &
4783 & 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4784 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4785 d(k)=qnifa(k) - dtz(k)*rhoinv(k)*s_awqnifa(k+1)*nonloc
4788 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4789 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) + &
4790 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4791 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4792 d(k)=qnifa(k) + dtz(k)*rhoinv(k)*(s_awqnifa(k)-s_awqnifa(k+1))*nonloc
4802 CALL tridiag2(kte,a,b,c,d,x)
4818IF (bl_mynn_cloudmix > 0 .AND. flag_qnbca .AND. &
4819 bl_mynn_mixscalars > 0)
THEN
4823 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4824 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) - &
4825 & 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4826 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4827 d(k)=qnbca(k) - dtz(k)*rhoinv(k)*s_awqnbca(k+1)*nonloc
4830 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4831 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) + &
4832 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4833 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4834 d(k)=qnbca(k) + dtz(k)*rhoinv(k)*(s_awqnbca(k)-s_awqnbca(k+1))*nonloc
4844 CALL tridiag2(kte,a,b,c,d,x)
4864 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4865 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k)
4866 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4870 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4871 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k)
4872 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4883 CALL tridiag2(kte,a,b,c,d,x)
4888 dozone(k)=(x(k)-ozone(k))/delt
4896 IF (bl_mynn_mixqt > 0)
THEN
4899 th_new = thl(k) + xlvcp/exner(k)*sqc(k) &
4900 & + xlscp/exner(k)*sqi(k)
4910 IF (sqc(k) > 0.0 .or. sqi(k) > 0.0)
THEN
4911 sqv2(k) = min(sqw2(k),qsat)
4912 portion_qc = sqc(k)/(sqc(k) + sqi(k))
4913 portion_qi = sqi(k)/(sqc(k) + sqi(k))
4914 condensate = max(sqw2(k) - qsat, 0.0)
4915 sqc2(k) = condensate*portion_qc
4916 sqi2(k) = condensate*portion_qi
4930 dqv(k)=(sqv2(k) - sqv(k))/delt
4934 IF (bl_mynn_cloudmix > 0)
THEN
4941 dqc(k)=(sqc2(k) - sqc(k))/delt
4953 IF (flag_qnc .AND. bl_mynn_mixscalars > 0)
THEN
4955 dqnc(k) = (qnc2(k)-qnc(k))/delt
4969 dqi(k)=(sqi2(k) - sqi(k))/delt
4983 dqs(k)=(sqs2(k) - sqs(k))/delt
4994 IF (flag_qni .AND. bl_mynn_mixscalars > 0)
THEN
4996 dqni(k)=(qni2(k)-qni(k))/delt
5016 CALL moisture_check(kte, delt, delp, exner, &
5017 sqv2, sqc2, sqi2, sqs2, thl, &
5018 dqv, dqc, dqi, dqs, dth )
5024 IF(dozone(k)*delt + ozone(k) < 0.)
THEN
5025 dozone(k)=-ozone(k)*0.99/delt
5034 dth(k)=(thl(k) + xlvcp/exner(k)*sqc2(k) &
5035 & + xlscp/exner(k)*(sqi2(k)) &
5045 dth(k)=(thl(k)+xlvcp/exner(k)*sqc2(k) - th(k))/delt
5056 IF (flag_qnwfa .AND. flag_qnifa .AND. &
5057 bl_mynn_mixscalars > 0)
THEN
5062 dqnwfa(k)=(qnwfa2(k) - qnwfa(k))/delt
5066 dqnifa(k)=(qnifa2(k) - qnifa(k))/delt
5078 IF (flag_qnbca .AND. bl_mynn_mixscalars > 0)
THEN
5080 dqnbca(k)=(qnbca2(k) - qnbca(k))/delt
5095 if (debug_code)
then
5098 wsp = sqrt(u(k)**2 + v(k)**2)
5099 wsp2 = sqrt((u(k)+du(k)*delt)**2 + (v(k)+du(k)*delt)**2)
5100 th2 = th(k) + dth(k)*delt
5102 if (wsp2 > 200. .or. tk2 > 360. .or. tk2 < 160.)
then
5104 print*,
"Outgoing problem at: i=",i,
" k=",k
5105 print*,
" incoming wsp=",wsp,
" outgoing wsp=",wsp2
5106 print*,
" incoming T=",th(k)*exner(k),
"outgoing T:",tk2
5107 print*,
" du=",du(k)*delt,
" dv=",dv(k)*delt,
" dth=",dth(k)*delt
5108 print*,
" km=",kmdz(k)*dz(k),
" kh=",khdz(k)*dz(k)
5109 print*,
" u*=",ust,
" wspd=",wspd,
"rhosfc=",rhosfc
5110 print*,
" LH=",flq*rhosfc*1004.,
" HFX=",flt*rhosfc*1004.
5111 print*,
" drag term=",ust**2/wspd*dtz(k)*rhosfc/rho(kts)
5116 print*,
"==thl:",thl(max(kproblem-3,1):min(kproblem+3,kte))
5117 print*,
"===qv:",sqv2(max(kproblem-3,1):min(kproblem+3,kte))
5118 print*,
"===qc:",sqc2(max(kproblem-3,1):min(kproblem+3,kte))
5119 print*,
"===qi:",sqi2(max(kproblem-3,1):min(kproblem+3,kte))
5120 print*,
"====u:",u(max(kproblem-3,1):min(kproblem+3,kte))
5121 print*,
"====v:",v(max(kproblem-3,1):min(kproblem+3,kte))
5125#ifdef HARDCODE_VERTICAL
5133 SUBROUTINE moisture_check(kte, delt, dp, exner, &
5134 qv, qc, qi, qs, th, &
5135 dqv, dqc, dqi, dqs, dth )
5151 integer,
intent(in) :: kte
5152 real(kind_phys),
intent(in) :: delt
5153 real(kind_phys),
dimension(kte),
intent(in) :: dp, exner
5154 real(kind_phys),
dimension(kte),
intent(inout) :: qv, qc, qi, qs, th
5155 real(kind_phys),
dimension(kte),
intent(inout) :: dqv, dqc, dqi, dqs, dth
5157 real(kind_phys):: dqc2, dqi2, dqs2, dqv2, sum, aa, dum
5158 real(kind_phys),
parameter :: qvmin = 1e-20, &
5163 dqc2 = max(0.0, qcmin-qc(k))
5164 dqi2 = max(0.0, qimin-qi(k))
5165 dqs2 = max(0.0, qimin-qs(k))
5168 dqc(k) = dqc(k) + dqc2/delt
5169 dqi(k) = dqi(k) + dqi2/delt
5170 dqs(k) = dqs(k) + dqs2/delt
5171 dqv(k) = dqv(k) - (dqc2+dqi2+dqs2)/delt
5172 dth(k) = dth(k) + xlvcp/exner(k)*(dqc2/delt) + &
5173 xlscp/exner(k)*((dqi2+dqs2)/delt)
5175 qc(k) = qc(k) + dqc2
5176 qi(k) = qi(k) + dqi2
5177 qs(k) = qs(k) + dqs2
5178 qv(k) = qv(k) - dqc2 - dqi2 - dqs2
5179 th(k) = th(k) + xlvcp/exner(k)*dqc2 + &
5180 xlscp/exner(k)*(dqi2+dqs2)
5183 dqv2 = max(0.0, qvmin-qv(k))
5184 dqv(k) = dqv(k) + dqv2/delt
5185 qv(k) = qv(k) + dqv2
5187 qv(k-1) = qv(k-1) - dqv2*dp(k)/dp(k-1)
5188 dqv(k-1) = dqv(k-1) - dqv2*dp(k)/dp(k-1)/delt
5190 qv(k) = max(qv(k),qvmin)
5191 qc(k) = max(qc(k),qcmin)
5192 qi(k) = max(qi(k),qimin)
5193 qs(k) = max(qs(k),qimin)
5198 if( dqv2 .gt. 1.e-20 )
then
5201 if( qv(k) .gt. 2.0*qvmin ) sum = sum + qv(k)*dp(k)
5203 aa = dqv2*dp(1)/max(1.e-20,sum)
5204 if( aa .lt. 0.5 )
then
5206 if( qv(k) .gt. 2.0*qvmin )
then
5209 dqv(k) = dqv(k) - dum/delt
5220 END SUBROUTINE moisture_check
5224 SUBROUTINE mynn_mix_chem(kts,kte,i, &
5226 nchem, kdvel, ndvel, &
5232 emis_ant_no, frp, rrfs_sd, &
5233 enh_mix, smoke_dbg )
5236 integer,
intent(in) :: kts,kte,i
5237 real(kind_phys),
dimension(kts:kte),
intent(in) :: dfh,dz,tcd,qcd
5238 real(kind_phys),
dimension(kts:kte),
intent(inout) :: rho
5239 real(kind_phys),
intent(in) :: flt
5240 real(kind_phys),
intent(in) :: delt,pblh
5241 integer,
intent(in) :: nchem, kdvel, ndvel
5242 real(kind_phys),
dimension( kts:kte+1),
intent(in) :: s_aw
5243 real(kind_phys),
dimension( kts:kte, nchem ),
intent(inout) :: chem1
5244 real(kind_phys),
dimension( kts:kte+1,nchem),
intent(in) :: s_awchem
5245 real(kind_phys),
dimension( ndvel ),
intent(in) :: vd1
5246 real(kind_phys),
intent(in) :: emis_ant_no,frp
5247 logical,
intent(in) :: rrfs_sd,enh_mix,smoke_dbg
5250 real(kind_phys),
dimension(kts:kte) :: dtz
5251 real(kind_phys),
dimension(kts:kte) :: a,b,c,d,x
5252 real(kind_phys):: rhs,dztop
5253 real(kind_phys):: t,dzk
5254 real(kind_phys):: hght
5255 real(kind_phys):: khdz_old, khdz_back
5256 integer :: k,kk,kmaxfire
5259 integer,
SAVE :: icall
5261 real(kind_phys),
dimension(kts:kte) :: rhoinv
5262 real(kind_phys),
dimension(kts:kte+1) :: rhoz,khdz
5263 real(kind_phys),
parameter :: NO_threshold = 10.0
5264 real(kind_phys),
parameter :: frp_threshold = 10.0
5265 real(kind_phys),
parameter :: pblh_threshold = 100.0
5267 dztop=.5*(dz(kte)+dz(kte-1))
5276 rhoinv(kts)=1./rho(kts)
5277 khdz(kts) =rhoz(kts)*dfh(kts)
5280 rhoz(k) =(rho(k)*dz(k-1) + rho(k-1)*dz(k))/(dz(k-1)+dz(k))
5281 rhoz(k) = max(rhoz(k),1e-4)
5282 rhoinv(k)=1./max(rho(k),1e-4)
5283 dzk = 0.5 *( dz(k)+dz(k-1) )
5284 khdz(k) = rhoz(k)*dfh(k)
5286 rhoz(kte+1)=rhoz(kte)
5287 khdz(kte+1)=rhoz(kte+1)*dfh(kte)
5291 khdz(k) = max(khdz(k), 0.5*s_aw(k))
5292 khdz(k) = max(khdz(k), -0.5*(s_aw(k)-s_aw(k+1)))
5296 IF ( rrfs_sd .and. enh_mix )
THEN
5299 khdz_back = pblh * 0.15 / dz(k)
5301 IF ( pblh < pblh_threshold )
THEN
5302 IF ( emis_ant_no > no_threshold )
THEN
5303 khdz(k) = max(1.1*khdz(k),sqrt((emis_ant_no / no_threshold)) / dz(k) * rhoz(k))
5306 IF ( frp > frp_threshold )
THEN
5307 kmaxfire = ceiling(log(frp))
5308 khdz(k) = max(1.1*khdz(k), (1. - k/(kmaxfire*2.)) * ((log(frp))**2.- 2.*log(frp)) / dz(k)*rhoz(k))
5322 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
5323 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)
5324 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)
5326 & - dtz(k)*vd1(ic)*chem1(k,ic) &
5327 & - dtz(k)*rhoinv(k)*s_awchem(k+1,ic)
5330 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)
5331 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
5332 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))
5333 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)
5334 d(k)=chem1(k,ic) + dtz(k)*rhoinv(k)*(s_awchem(k,ic)-s_awchem(k+1,ic))
5341 d(kte)=chem1(kte,ic)
5343 CALL tridiag3(kte,a,b,c,d,x)
5345 IF ( smoke_dbg )
THEN
5346 print*,
'aerosol mixing ic,chem1,chem2(k,ic)',ic,(chem1(kts:kts+10,ic)),(x(kts:kts+10))
5347 print*,
'aerosol PBL mixing ic,vd1(ic)',ic,vd1(ic)
5355 END SUBROUTINE mynn_mix_chem
5359 SUBROUTINE retrieve_exchange_coeffs(kts,kte,&
5360 &dfm,dfh,dz,K_m,K_h)
5364 integer ,
intent(in) :: kts,kte
5366 real(kind_phys),
dimension(KtS:KtE),
intent(in) :: dz,dfm,dfh
5368 real(kind_phys),
dimension(KtS:KtE),
intent(out) :: K_m, K_h
5372 real(kind_phys):: dzk
5378 dzk = 0.5 *( dz(k)+dz(k-1) )
5383 END SUBROUTINE retrieve_exchange_coeffs
5387 SUBROUTINE tridiag(n,a,b,c,d)
5397 integer,
intent(in):: n
5398 real(kind_phys),
dimension(n),
intent(in) :: a,b
5399 real(kind_phys),
dimension(n),
intent(inout) :: c,d
5403 real(kind_phys),
dimension(n) :: q
5410 p=1./(b(i)+a(i)*q(i-1))
5412 d(i)=(d(i)-a(i)*d(i-1))*p
5416 d(i)=d(i)+q(i)*d(i+1)
5419 END SUBROUTINE tridiag
5423 subroutine tridiag2(n,a,b,c,d,x)
5432 integer,
intent(in) :: n
5433 real(kind_phys),
dimension(n),
intent(in) :: a,b,c,d
5434 real(kind_phys),
dimension(n),
intent(out):: x
5435 real(kind_phys),
dimension(n) :: cp,dp
5444 m = b(i)-cp(i-1)*a(i)
5446 dp(i) = (d(i)-dp(i-1)*a(i))/m
5452 x(i) = dp(i)-cp(i)*x(i+1)
5455 end subroutine tridiag2
5458 subroutine tridiag3(kte,a,b,c,d,x)
5473 integer,
intent(in) :: kte
5474 integer,
parameter :: kts=1
5475 real(kind_phys),
dimension(kte) :: a,b,c,d
5476 real(kind_phys),
dimension(kte),
intent(out) :: x
5483 d(in)=d(in)-c(in)*d(in+1)/b(in+1)
5484 b(in)=b(in)-c(in)*a(in+1)/b(in+1)
5488 d(in)=d(in)-a(in)*d(in-1)/b(in-1)
5496 end subroutine tridiag3
5518 SUBROUTINE get_pblh(KTS,KTE,zi,thetav1D,qke1D,zw1D,dz1D,landsea,kzi)
5537 integer,
intent(in) :: KTS,KTE
5539#ifdef HARDCODE_VERTICAL
5541# define kte HARDCODE_VERTICAL
5544 real(kind_phys),
intent(out) :: zi
5545 real(kind_phys),
intent(in) :: landsea
5546 real(kind_phys),
dimension(kts:kte),
intent(in) :: thetav1D, qke1D, dz1D
5547 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw1D
5549 real(kind_phys):: PBLH_TKE,qtke,qtkem1,wt,maxqke,TKEeps,minthv
5550 real(kind_phys):: delt_thv
5551 real(kind_phys),
parameter :: sbl_lim = 200.
5552 real(kind_phys),
parameter :: sbl_damp = 400.
5553 integer :: I,J,K,kthv,ktke,kzi
5562 DO WHILE (zw1d(k) .LE. 200.)
5564 IF (minthv > thetav1d(k))
then
5565 minthv = thetav1d(k)
5575 IF((landsea-1.5).GE.0)
THEN
5587 IF (thetav1d(k) .GE. (minthv + delt_thv))
THEN
5588 zi = zw1d(k) - dz1d(k-1)* &
5589 & min((thetav1d(k)-(minthv + delt_thv))/ &
5590 & max(thetav1d(k)-thetav1d(k-1),1e-6),1.0)
5593 IF (k .EQ. kte-1) zi = zw1d(kts+1)
5594 IF (zi .NE. 0.0)
exit
5603 maxqke = max(qke1d(kts),0.)
5607 tkeeps = max(tkeeps,0.02)
5614 qtke =max(qke1d(k)/2.,0.)
5615 qtkem1=max(qke1d(k-1)/2.,0.)
5616 IF (qtke .LE. tkeeps)
THEN
5617 pblh_tke = zw1d(k) - dz1d(k-1)* &
5618 & min((tkeeps-qtke)/max(qtkem1-qtke, 1e-6), 1.0)
5620 pblh_tke = max(pblh_tke,zw1d(kts+1))
5624 IF (k .EQ. kte-1) pblh_tke = zw1d(kts+1)
5625 IF (pblh_tke .NE. 0.)
exit
5634 pblh_tke = min(pblh_tke,zi+350.)
5635 pblh_tke = max(pblh_tke,max(zi-350.,10.))
5637 wt=.5*tanh((zi - sbl_lim)/sbl_damp) + .5
5638 IF (maxqke <= 0.05)
THEN
5642 zi=pblh_tke*(1.-wt) + zi*wt
5647 IF ( zw1d(k) >= zi)
THEN
5653#ifdef HARDCODE_VERTICAL
5681 & kts,kte,dt,zw,dz,p,rho, &
5685 & u,v,w,th,thl,thv,tk, &
5687 & qnc,qni,qnwfa,qnifa,qnbca, &
5688 & exner,vt,vq,sgm, &
5689 & ust,flt,fltv,flq,flqv, &
5690 & pblh,kpbl,dx,landsea,ts, &
5691 ! outputs - updraft properties
5693 & edmf_qt,edmf_thl, &
5694 & edmf_ent,edmf_qc, &
5695 ! outputs - variables needed for solver
5696 & s_aw,s_awthl,s_awqt, &
5698 & s_awu,s_awv,s_awqke, &
5699 & s_awqnc,s_awqni, &
5700 & s_awqnwfa,s_awqnifa, &
5702 & sub_thl,sub_sqv, &
5704 & det_thl,det_sqv,det_sqc, &
5707 & nchem,chem1,s_awchem, &
5709 ! in/outputs - subgrid scale clouds
5710 & qc_bl1d,cldfra_bl1d, &
5711 & qc_bl1D_old,cldfra_bl1D_old, &
5712 ! inputs - flags for moist arrays
5715 & F_QNWFA,F_QNIFA,F_QNBCA, &
5718 & maxwidth,ktop,maxmf,ztop, &
5719 ! inputs for stochastic perturbations
5720 & spp_pbl,rstoch_col )
5723 integer,
intent(in) :: KTS,KTE,KPBL,momentum_opt,tke_opt,scalar_opt
5725#ifdef HARDCODE_VERTICAL
5727# define kte HARDCODE_VERTICAL
5731 integer,
intent(in) :: spp_pbl
5732 real(kind_phys),
dimension(kts:kte) :: rstoch_col
5734 real(kind_phys),
dimension(kts:kte),
intent(in) :: &
5735 &U,V,W,TH,THL,TK,QT,QV,QC, &
5736 &exner,dz,THV,P,rho,qke,qnc,qni,qnwfa,qnifa,qnbca
5737 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
5738 real(kind_phys),
intent(in) :: flt,fltv,flq,flqv,Psig_shcu, &
5739 &landsea,ts,dx,dt,ust,pblh
5740 logical,
optional :: F_QC,F_QI,F_QNC,F_QNI,F_QNWFA,F_QNIFA,F_QNBCA
5743 real(kind_phys),
dimension(kts:kte),
intent(out) :: edmf_a,edmf_w, &
5744 & edmf_qt,edmf_thl,edmf_ent,edmf_qc
5746 real(kind_phys),
dimension(kts:kte) :: edmf_th
5748 integer,
intent(out) :: ktop
5749 real(kind_phys),
intent(out) :: maxmf,ztop,maxwidth
5751 real(kind_phys),
dimension(kts:kte+1) :: s_aw, &
5752 &s_awthl,s_awqt,s_awqv,s_awqc,s_awqnc,s_awqni, &
5753 &s_awqnwfa,s_awqnifa,s_awqnbca,s_awu,s_awv, &
5756 real(kind_phys),
dimension(kts:kte),
intent(inout) :: &
5757 &qc_bl1d,cldfra_bl1d,qc_bl1d_old,cldfra_bl1d_old
5759 integer,
parameter :: nup=8, debug_mf=0
5760 real(kind_phys) :: nup2
5765 real(kind_phys),
dimension(kts:kte+1,1:NUP) :: &
5766 &UPW,UPTHL,UPQT,UPQC,UPQV, &
5767 &UPA,UPU,UPV,UPTHV,UPQKE,UPQNC, &
5768 &UPQNI,UPQNWFA,UPQNIFA,UPQNBCA
5770 real(kind_phys),
dimension(kts:kte,1:NUP) :: ENT,ENTf
5771 integer,
dimension(kts:kte,1:NUP) :: ENTi
5774 real(kind_phys):: fltv2,wstar,qstar,thstar,sigmaW,sigmaQT, &
5775 &sigmaTH,z0,pwmin,pwmax,wmin,wmax,wlv,Psig_w,maxw,maxqc,wpbl
5776 real(kind_phys):: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,QNCn,QNIn, &
5777 & QNWFAn,QNIFAn,QNBCAn, &
5778 & Wn2,Wn,EntEXP,EntEXM,EntW,BCOEFF,THVkm1,THVk,Pk,rho_int
5781 real(kind_phys),
parameter :: &
5788 real(kind_phys),
parameter :: &
5793 real(kind_phys),
parameter :: Atot = 0.10
5794 real(kind_phys),
parameter :: lmax = 1000.
5795 real(kind_phys),
parameter :: lmin = 300.
5796 real(kind_phys),
parameter :: dlmin = 0.
5797 real(kind_phys) :: minwidth
5798 real(kind_phys) :: dl
5799 real(kind_phys),
parameter :: dcut = 1.2
5803 real(kind_phys):: cn,c,l,n,an2,hux,wspd_pbl,cloud_base,width_flx
5806 integer,
intent(in) :: nchem
5807 real(kind_phys),
dimension(:, :) :: chem1
5808 real(kind_phys),
dimension(kts:kte+1, nchem) :: s_awchem
5809 real(kind_phys),
dimension(nchem) :: chemn
5810 real(kind_phys),
dimension(kts:kte+1,1:NUP, nchem) :: UPCHEM
5812 real(kind_phys),
dimension(kts:kte+1, nchem) :: edmf_chem
5813 logical,
intent(in) :: mix_chem
5816 real(kind_phys):: ERF
5818 logical :: superadiabatic
5821 real(kind_phys),
dimension(kts:kte),
intent(inout) :: vt, vq, sgm
5822 real(kind_phys):: sigq,xl,rsl,cpm,a,qmq,mf_cf,Aup,Q1,diffqt,qsat_tk,&
5823 Fng,qww,alpha,beta,bb,f,pt,t,q2p,b9,satvp,rhgrid, &
5824 Ac_mf,Ac_strat,qc_mf
5825 real(kind_phys),
parameter :: cf_thresh = 0.5
5828 real(kind_phys),
dimension(kts:kte) :: exneri,dzi,rhoz
5829 real(kind_phys):: THp, QTp, QCp, QCs, esat, qsl
5830 real(kind_phys):: csigma,acfac,ac_wsp
5833 integer :: overshoot
5834 real(kind_phys):: bvf, Frz, dzp
5838 real(kind_phys):: adjustment, flx1
5839 real(kind_phys),
parameter :: fluxportion=0.75
5844 real(kind_phys),
dimension(kts:kte) :: sub_thl,sub_sqv,sub_u,sub_v, &
5845 det_thl,det_sqv,det_sqc,det_u,det_v, &
5846 envm_a,envm_w,envm_thl,envm_sqv,envm_sqc, &
5848 real(kind_phys),
dimension(kts:kte+1) :: envi_a,envi_w
5849 real(kind_phys):: temp,sublim,qc_ent,qv_ent,qt_ent,thl_ent,detrate, &
5850 detrateUV,oow,exc_fac,aratio,detturb,qc_grid,qc_sgs, &
5851 qc_plume,exc_heat,exc_moist,tk_int,tvs
5852 real(kind_phys),
parameter :: Cdet = 1./45.
5853 real(kind_phys),
parameter :: dzpmax = 300.
5858 real(kind_phys),
parameter :: Csub=0.25
5861 real(kind_phys),
parameter :: pgfac = 0.00
5862 real(kind_phys):: Uk,Ukm1,Vk,Vkm1,dxsa
5892 if ( mix_chem )
then
5893 upchem(kts:kte+1,1:nup,1:nchem)=0.0
5904 if ( mix_chem )
then
5905 edmf_chem(kts:kte+1,1:nchem) = 0.0
5922 if ( mix_chem )
then
5923 s_awchem(kts:kte+1,1:nchem) = 0.0
5943 if (zw(k) > pblh + 500.)
exit
5946 if (w(k) < 0.)wpbl = 2.*w(k)
5947 maxw = max(maxw,abs(wpbl))
5950 if (zw(k)<=50.)k50=k
5953 qc_sgs = max(qc(k), qc_bl1d(k))
5954 if (qc_sgs> 1e-5 .and. (cldfra_bl1d(k) .ge. 0.5) .and. cloud_base == 9000.0)
then
5955 cloud_base = 0.5*(zw(k)+zw(k+1))
5960 maxw = max(0.,maxw - 1.0)
5961 psig_w = max(0.0, 1.0 - maxw)
5962 psig_w = min(psig_w, psig_shcu)
5966 if(psig_w == 0.0 .and. fltv > 0.0) fltv2 = -1.*fltv
5970 superadiabatic = .false.
5971 if ((landsea-1.5).ge.0)
then
5976 tvs = ts*(1.0+p608*qv(kts))
5979 if ((thv(k)-tvs)/(0.5*dz(k)) < hux)
then
5980 superadiabatic = .true.
5982 superadiabatic = .false.
5986 if ((thv(k)-thv(k-1))/(0.5*(dz(k)+dz(k-1))) < hux)
then
5987 superadiabatic = .true.
5989 superadiabatic = .false.
6004 maxwidth = min(dx*dcut, lmax)
6006 maxwidth = min(maxwidth, 1.1_kind_phys*pblh)
6008 if ((landsea-1.5) .lt. 0)
then
6009 maxwidth = min(maxwidth, 0.5_kind_phys*cloud_base)
6011 maxwidth = min(maxwidth, 0.9_kind_phys*cloud_base)
6014 wspd_pbl=sqrt(max(u(kts)**2 + v(kts)**2, 0.01_kind_phys))
6017 if ((landsea-1.5).LT.0)
then
6018 width_flx = max(min(1000.*(0.6*tanh((fltv - 0.040)/0.04) + .5),1000._kind_phys), 0._kind_phys)
6020 width_flx = max(min(1000.*(0.6*tanh((fltv - 0.007)/0.02) + .5),1000._kind_phys), 0._kind_phys)
6022 maxwidth = min(maxwidth, width_flx)
6026 if (maxwidth .ge. (lmax - 1.0) .and. fltv .gt. 0.2)minwidth = lmin + dlmin*min((fltv-0.2)/0.3, 1._kind_phys)
6028 if (maxwidth .le. minwidth)
then
6039if ( fltv2 > 0.002 .AND. (maxwidth > minwidth) .AND. superadiabatic)
then
6044 dl = (maxwidth - minwidth)/real(nup-1,kind=kind_phys)
6047 l = minwidth + dl*real(i-1)
6048 cn = cn + l**d * (l*l)/(dx*dx) * dl
6053 if ((landsea-1.5).LT.0)
then
6054 acfac = .5*tanh((fltv2 - 0.02)/0.05) + .5
6056 acfac = .5*tanh((fltv2 - 0.015)/0.04) + .5
6062 if (wspd_pbl .le. 10.)
then
6065 ac_wsp = 1.0 - min((wspd_pbl - 10.0)/15., 1.0)
6067 acfac = acfac * ac_wsp
6073 l = minwidth + dl*real(i-1)
6075 upa(1,i) = n*l*l/(dx*dx) * dl
6077 upa(1,i) = upa(1,i)*acfac
6078 an2 = an2 + upa(1,i)
6087 wstar=max(1.e-2,(gtr*fltv2*pblh)**(onethird))
6088 qstar=max(flq,1.0e-5)/wstar
6091 if ((landsea-1.5) .ge. 0)
then
6100 if ((landsea-1.5).GE.0)
then
6109 exc_fac = exc_fac * ac_wsp
6112 sigmaw =csigma*wstar*(z0/pblh)**(onethird)*(1 - 0.8*z0/pblh)
6113 sigmaqt=csigma*qstar*(z0/pblh)**(onethird)
6114 sigmath=csigma*thstar*(z0/pblh)**(onethird)
6118 wmin=min(sigmaw*pwmin,0.1)
6119 wmax=min(sigmaw*pwmax,0.5)
6123 wlv=wmin+(wmax-wmin)/nup2*(i-1)
6126 upw(1,i)=wmin + real(i)/real(nup)*(wmax-wmin)
6127 upu(1,i)=(u(kts)*dz(kts+1)+u(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6128 upv(1,i)=(v(kts)*dz(kts+1)+v(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6132 exc_heat = exc_fac*upw(1,i)*sigmath/sigmaw
6133 upthv(1,i)=(thv(kts)*dz(kts+1)+thv(kts+1)*dz(kts))/(dz(kts)+dz(kts+1)) &
6135 upthl(1,i)=(thl(kts)*dz(kts+1)+thl(kts+1)*dz(kts))/(dz(kts)+dz(kts+1)) &
6139 exc_moist=exc_fac*upw(1,i)*sigmaqt/sigmaw
6140 upqt(1,i)=(qt(kts)*dz(kts+1)+qt(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))&
6143 upqke(1,i)=(qke(kts)*dz(kts+1)+qke(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6144 upqnc(1,i)=(qnc(kts)*dz(kts+1)+qnc(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6145 upqni(1,i)=(qni(kts)*dz(kts+1)+qni(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6146 upqnwfa(1,i)=(qnwfa(kts)*dz(kts+1)+qnwfa(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6147 upqnifa(1,i)=(qnifa(kts)*dz(kts+1)+qnifa(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6148 upqnbca(1,i)=(qnbca(kts)*dz(kts+1)+qnbca(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6151 if ( mix_chem )
then
6154 upchem(1,i,ic)=(chem1(kts,ic)*dz(kts+1)+chem1(kts+1,ic)*dz(kts))/(dz(kts)+dz(kts+1))
6160 envm_thl(kts:kte)=thl(kts:kte)
6161 envm_sqv(kts:kte)=qv(kts:kte)
6162 envm_sqc(kts:kte)=qc(kts:kte)
6163 envm_u(kts:kte)=u(kts:kte)
6164 envm_v(kts:kte)=v(kts:kte)
6166 rhoz(k) = (rho(k)*dz(k+1)+rho(k+1)*dz(k))/(dz(k+1)+dz(k))
6168 rhoz(kte) = rho(kte)
6171 dxsa = 1. - min(max((12000.0-dx)/(12000.0-3000.0), 0.), 1.)
6177 l = minwidth + dl*real(i-1)
6181 wmin = 0.3 + l*0.0005
6182 ent(k,i) = 0.33/(min(max(upw(k-1,i),wmin),0.9)*l)
6190 ent(k,i) = max(ent(k,i),0.0003)
6194 IF(zw(k) >= min(pblh+1500., 4000.))
THEN
6195 ent(k,i)=ent(k,i) + (zw(k)-min(pblh+1500.,4000.))*5.0e-6
6199 ent(k,i) = ent(k,i) * (1.0 - rstoch_col(k))
6201 ent(k,i) = min(ent(k,i),0.9/(zw(k+1)-zw(k)))
6204 uk =(u(k)*dz(k+1)+u(k+1)*dz(k))/(dz(k+1)+dz(k))
6205 ukm1=(u(k-1)*dz(k)+u(k)*dz(k-1))/(dz(k-1)+dz(k))
6206 vk =(v(k)*dz(k+1)+v(k+1)*dz(k))/(dz(k+1)+dz(k))
6207 vkm1=(v(k-1)*dz(k)+v(k)*dz(k-1))/(dz(k-1)+dz(k))
6210 entexp= ent(k,i)*(zw(k+1)-zw(k))
6211 entexm= entexp*0.3333
6212 qtn =upqt(k-1,i) *(1.-entexp) + qt(k)*entexp
6213 thln=upthl(k-1,i)*(1.-entexp) + thl(k)*entexp
6214 un =upu(k-1,i) *(1.-entexm) + u(k)*entexm + dxsa*pgfac*(uk - ukm1)
6215 vn =upv(k-1,i) *(1.-entexm) + v(k)*entexm + dxsa*pgfac*(vk - vkm1)
6216 qken=upqke(k-1,i)*(1.-entexp) + qke(k)*entexp
6217 qncn=upqnc(k-1,i)*(1.-entexp) + qnc(k)*entexp
6218 qnin=upqni(k-1,i)*(1.-entexp) + qni(k)*entexp
6219 qnwfan=upqnwfa(k-1,i)*(1.-entexp) + qnwfa(k)*entexp
6220 qnifan=upqnifa(k-1,i)*(1.-entexp) + qnifa(k)*entexp
6221 qnbcan=upqnbca(k-1,i)*(1.-entexp) + qnbca(k)*entexp
6237 if ( mix_chem )
then
6242 chemn(ic)=upchem(k-1,i,ic)*(1.-entexp) + chem1(k,ic)*entexp
6247 pk =(p(k)*dz(k+1)+p(k+1)*dz(k))/(dz(k+1)+dz(k))
6252 thvk =(thv(k)*dz(k+1)+thv(k+1)*dz(k))/(dz(k+1)+dz(k))
6253 thvkm1=(thv(k-1)*dz(k)+thv(k)*dz(k-1))/(dz(k-1)+dz(k))
6256 b=grav*(thvn/thvk - 1.0)
6271 IF (upw(k-1,i) < 0.2 )
THEN
6272 wn = upw(k-1,i) + (-2. * ent(k,i) * upw(k-1,i) + bcoeff*b / max(upw(k-1,i),0.2)) * min(zw(k)-zw(k-1), 250.)
6274 wn = upw(k-1,i) + (-2. * ent(k,i) * upw(k-1,i) + bcoeff*b / upw(k-1,i)) * min(zw(k)-zw(k-1), 250.)
6278 IF(wn > upw(k-1,i) + min(1.25*(zw(k)-zw(k-1))/200., 2.0) )
THEN
6279 wn = upw(k-1,i) + min(1.25*(zw(k)-zw(k-1))/200., 2.0)
6282 IF(wn < upw(k-1,i) - min(1.25*(zw(k)-zw(k-1))/200., 2.0) )
THEN
6283 wn = upw(k-1,i) - min(1.25*(zw(k)-zw(k-1))/200., 2.0)
6285 wn = min(max(wn,0.0), 3.0)
6289 IF (k==kts+1 .AND. wn == 0.)
THEN
6294 IF (debug_mf == 1)
THEN
6295 IF (wn .GE. 3.0)
THEN
6297 print *,
" **** SUSPICIOUSLY LARGE W:"
6298 print *,
' QCn:',qcn,
' ENT=',ent(k,i),
' Nup2=',nup2
6299 print *,
'pblh:',pblh,
' Wn:',wn,
' UPW(k-1)=',upw(k-1,i)
6300 print *,
'K=',k,
' B=',b,
' dz=',zw(k)-zw(k-1)
6305 IF (wn <= 0.0 .AND. overshoot == 0)
THEN
6307 IF ( thvk-thvkm1 .GT. 0.0 )
THEN
6308 bvf = sqrt( gtr*(thvk-thvkm1)/dz(k) )
6310 frz = upw(k-1,i)/(bvf*dz(k))
6312 dzp = dz(k)*max(min(frz,1.0),0.0)
6326 aratio = min(upa(k-1,i)/(1.-upa(k-1,i)), 0.5)
6328 oow = -0.060/max(1.0,(0.5*(wn+upw(k-1,i))))
6329 detrate = min(max(oow*(wn-upw(k-1,i))/dz(k), detturb), .0002)
6330 detrateuv= min(max(oow*(wn-upw(k-1,i))/dz(k), detturb), .0001)
6331 envm_thl(k)=envm_thl(k) + (0.5*(thl_ent + upthl(k-1,i)) - thl(k))*detrate*aratio*min(dzp,dzpmax)
6332 qv_ent = 0.5*(max(qt_ent-qc_ent,0.) + max(upqt(k-1,i)-upqc(k-1,i),0.))
6333 envm_sqv(k)=envm_sqv(k) + (qv_ent-qv(k))*detrate*aratio*min(dzp,dzpmax)
6334 IF (upqc(k-1,i) > 1e-8)
THEN
6335 IF (qc(k) > 1e-6)
THEN
6338 qc_grid = cldfra_bl1d(k)*qc_bl1d(k)
6340 envm_sqc(k)=envm_sqc(k) + max(upa(k-1,i)*0.5*(qcn + upqc(k-1,i)) - qc_grid, 0.0)*detrate*aratio*min(dzp,dzpmax)
6342 envm_u(k) =envm_u(k) + (0.5*(un + upu(k-1,i)) - u(k))*detrateuv*aratio*min(dzp,dzpmax)
6343 envm_v(k) =envm_v(k) + (0.5*(vn + upv(k-1,i)) - v(k))*detrateuv*aratio*min(dzp,dzpmax)
6361 IF ( mix_chem )
THEN
6363 upchem(k,i,ic) = chemn(ic)
6372 IF (debug_mf == 1)
THEN
6373 IF (maxval(upw(:,i)) > 10.0 .OR. minval(upa(:,i)) < 0.0 .OR. &
6374 maxval(upa(:,i)) > atot .OR. nup2 > 10)
THEN
6376 print *,
'flq:',flq,
' fltv:',fltv2,
' Nup2=',nup2
6377 print *,
'pblh:',pblh,
' wstar:',wstar,
' ktop=',ktop
6378 print *,
'sigmaW=',sigmaw,
' sigmaTH=',sigmath,
' sigmaQT=',sigmaqt
6383 print *,
'UPA:',upa(:,i)
6384 print *,
'UPW:',upw(:,i)
6385 print *,
'UPTHL:',upthl(:,i)
6386 print *,
'UPQT:',upqt(:,i)
6387 print *,
'ENT:',ent(:,i)
6408 s_aw(k+1) = s_aw(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*psig_w
6409 s_awthl(k+1)= s_awthl(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upthl(k,i)*psig_w
6410 s_awqt(k+1) = s_awqt(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqt(k,i)*psig_w
6415 qc_plume = upqc(k,i)
6419 s_awqc(k+1) = s_awqc(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*qc_plume*psig_w
6420 s_awqv(k+1) = s_awqt(k+1) - s_awqc(k+1)
6424 if (momentum_opt > 0)
then
6427 s_awu(k+1) = s_awu(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upu(k,i)*psig_w
6428 s_awv(k+1) = s_awv(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upv(k,i)*psig_w
6433 if (tke_opt > 0)
then
6436 s_awqke(k+1)= s_awqke(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqke(k,i)*psig_w
6441 if ( mix_chem )
then
6445 s_awchem(k+1,ic) = s_awchem(k+1,ic) + rhoz(k)*upa(k,i)*upw(k,i)*upchem(k,i,ic)*psig_w
6451 if (scalar_opt > 0)
then
6454 s_awqnc(k+1) = s_awqnc(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqnc(k,i)*psig_w
6455 s_awqni(k+1) = s_awqni(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqni(k,i)*psig_w
6456 s_awqnwfa(k+1)= s_awqnwfa(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqnwfa(k,i)*psig_w
6457 s_awqnifa(k+1)= s_awqnifa(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqnifa(k,i)*psig_w
6458 s_awqnbca(k+1)= s_awqnbca(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqnbca(k,i)*psig_w
6466 IF (s_aw(kts+1) /= 0.)
THEN
6467 dzi(kts) = 0.5*(dz(kts)+dz(kts+1))
6468 flx1 = max(s_aw(kts+1)*(th(kts)-th(kts+1))/dzi(kts),1.0e-5)
6477 IF (flx1 > fluxportion*flt/dz(kts) .AND. flx1>0.0)
THEN
6478 adjustment= fluxportion*flt/dz(kts)/flx1
6479 s_aw = s_aw*adjustment
6480 s_awthl = s_awthl*adjustment
6481 s_awqt = s_awqt*adjustment
6482 s_awqc = s_awqc*adjustment
6483 s_awqv = s_awqv*adjustment
6484 s_awqnc = s_awqnc*adjustment
6485 s_awqni = s_awqni*adjustment
6486 s_awqnwfa = s_awqnwfa*adjustment
6487 s_awqnifa = s_awqnifa*adjustment
6488 s_awqnbca = s_awqnbca*adjustment
6489 IF (momentum_opt > 0)
THEN
6490 s_awu = s_awu*adjustment
6491 s_awv = s_awv*adjustment
6493 IF (tke_opt > 0)
THEN
6494 s_awqke= s_awqke*adjustment
6496 IF ( mix_chem )
THEN
6497 s_awchem = s_awchem*adjustment
6499 upa = upa*adjustment
6507 edmf_a(k) =edmf_a(k) +upa(k,i)
6508 edmf_w(k) =edmf_w(k) +upa(k,i)*upw(k,i)
6509 edmf_qt(k) =edmf_qt(k) +upa(k,i)*upqt(k,i)
6510 edmf_thl(k)=edmf_thl(k)+upa(k,i)*upthl(k,i)
6511 edmf_ent(k)=edmf_ent(k)+upa(k,i)*ent(k,i)
6512 edmf_qc(k) =edmf_qc(k) +upa(k,i)*upqc(k,i)
6518 if (edmf_a(k)>0.)
then
6519 edmf_w(k)=edmf_w(k)/edmf_a(k)
6520 edmf_qt(k)=edmf_qt(k)/edmf_a(k)
6521 edmf_thl(k)=edmf_thl(k)/edmf_a(k)
6522 edmf_ent(k)=edmf_ent(k)/edmf_a(k)
6523 edmf_qc(k)=edmf_qc(k)/edmf_a(k)
6524 edmf_a(k)=edmf_a(k)*psig_w
6526 if(edmf_a(k)*edmf_w(k) > maxmf) maxmf = edmf_a(k)*edmf_w(k)
6531 if ( mix_chem )
then
6535 edmf_chem(k,ic) = edmf_chem(k,ic) + rhoz(k)*upa(k,i)*upchem(k,i,ic)
6540 if (edmf_a(k)>0.)
then
6542 edmf_chem(k,ic) = edmf_chem(k,ic)/edmf_a(k)
6556 envi_w(k) = onethird*(edmf_w(k-1)+edmf_w(k)+edmf_w(k+1))
6557 envi_a(k) = onethird*(edmf_a(k-1)+edmf_a(k)+edmf_a(k+1))*adjustment
6560 envi_w(kts) = edmf_w(kts)
6561 envi_a(kts) = edmf_a(kts)
6564 envi_a(kte) = edmf_a(kte)
6567 envi_a(kte+1) = edmf_a(kte)
6571 IF (envi_w(kts) > 0.9*dz(kts)/dt)
THEN
6572 sublim = 0.9*dz(kts)/dt/envi_w(kts)
6580 envi_w(k)=csub*sublim*envi_w(k)*temp/(1.-temp)
6584 dzi(kts) = 0.5*(dz(kts)+dz(kts+1))
6585 sub_thl(kts)= 0.5*envi_w(kts)*envi_a(kts)* &
6586 (rho(kts+1)*thl(kts+1)-rho(kts)*thl(kts))/dzi(kts)/rhoz(k)
6587 sub_sqv(kts)= 0.5*envi_w(kts)*envi_a(kts)* &
6588 (rho(kts+1)*qv(kts+1)-rho(kts)*qv(kts))/dzi(kts)/rhoz(k)
6590 dzi(k) = 0.5*(dz(k)+dz(k+1))
6591 sub_thl(k)= 0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * &
6592 (rho(k+1)*thl(k+1)-rho(k)*thl(k))/dzi(k)/rhoz(k)
6593 sub_sqv(k)= 0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * &
6594 (rho(k+1)*qv(k+1)-rho(k)*qv(k))/dzi(k)/rhoz(k)
6598 det_thl(k)=cdet*(envm_thl(k)-thl(k))*envi_a(k)*psig_w
6599 det_sqv(k)=cdet*(envm_sqv(k)-qv(k))*envi_a(k)*psig_w
6600 det_sqc(k)=cdet*(envm_sqc(k)-qc(k))*envi_a(k)*psig_w
6603 IF (momentum_opt > 0)
THEN
6604 sub_u(kts)=0.5*envi_w(kts)*envi_a(kts)* &
6605 (rho(kts+1)*u(kts+1)-rho(kts)*u(kts))/dzi(kts)/rhoz(k)
6606 sub_v(kts)=0.5*envi_w(kts)*envi_a(kts)* &
6607 (rho(kts+1)*v(kts+1)-rho(kts)*v(kts))/dzi(kts)/rhoz(k)
6609 sub_u(k)=0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * &
6610 (rho(k+1)*u(k+1)-rho(k)*u(k))/dzi(k)/rhoz(k)
6611 sub_v(k)=0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * &
6612 (rho(k+1)*v(k+1)-rho(k)*v(k))/dzi(k)/rhoz(k)
6616 det_u(k) = cdet*(envm_u(k)-u(k))*envi_a(k)*psig_w
6617 det_v(k) = cdet*(envm_v(k)-v(k))*envi_a(k)*psig_w
6626 exneri(k) = (exner(k)*dz(k+1)+exner(k+1)*dz(k))/(dz(k+1)+dz(k))
6627 edmf_th(k)= edmf_thl(k) + xlvcp/exneri(k)*edmf_qc(k)
6628 dzi(k) = 0.5*(dz(k)+dz(k+1))
6636 if(0.5*(edmf_qc(k)+edmf_qc(k-1))>0.0 .and. (cldfra_bl1d(k) < cf_thresh))
THEN
6638 aup = (edmf_a(k)*dzi(k-1)+edmf_a(k-1)*dzi(k))/(dzi(k-1)+dzi(k))
6639 thp = (edmf_th(k)*dzi(k-1)+edmf_th(k-1)*dzi(k))/(dzi(k-1)+dzi(k))
6640 qtp = (edmf_qt(k)*dzi(k-1)+edmf_qt(k-1)*dzi(k))/(dzi(k-1)+dzi(k))
6646 qsl=ep_2*esat/max(1.e-7,(p(k)-ep_3*esat))
6649 if (edmf_qc(k)>0.0 .and. edmf_qc(k-1)>0.0)
then
6650 qcp = (edmf_qc(k)*dzi(k-1)+edmf_qc(k-1)*dzi(k))/(dzi(k-1)+dzi(k))
6652 qcp = max(edmf_qc(k),edmf_qc(k-1))
6659 rsl = xl*qsat_tk / (r_v*tk(k)**2)
6661 cpm = cp + qt(k)*cpv
6662 a = 1./(1. + xl*rsl/cpm)
6665 q2p = xlvcp/exner(k)
6666 pt = thl(k) +q2p*qcp*aup
6675 beta = pt*xl/(tk(k)*cp) - 1.61*pt
6690 sigq = 10. * aup * (qtp - qt(k))
6692 sigq = max(sigq, qsat_tk*0.02 )
6693 sigq = min(sigq, qsat_tk*0.25 )
6695 qmq = a * (qt(k) - qsat_tk)
6698 if ((landsea-1.5).GE.0)
then
6702 mf_cf = min(max(0.5 + 0.36 * atan(1.55*q1),0.01),0.6)
6703 mf_cf = max(mf_cf, 1.2 * aup)
6704 mf_cf = min(mf_cf, 5.0 * aup)
6709 mf_cf = min(max(0.5 + 0.36 * atan(1.55*q1),0.01),0.6)
6710 mf_cf = max(mf_cf, 1.8 * aup)
6711 mf_cf = min(mf_cf, 5.0 * aup)
6725 if ((landsea-1.5).GE.0)
then
6726 if (qcp * aup > 5e-5)
then
6727 qc_bl1d(k) = 1.86 * (qcp * aup) - 2.2e-5
6729 qc_bl1d(k) = 1.18 * (qcp * aup)
6731 cldfra_bl1d(k) = mf_cf
6734 if (qcp * aup > 5e-5)
then
6735 qc_bl1d(k) = 1.86 * (qcp * aup) - 2.2e-5
6737 qc_bl1d(k) = 1.18 * (qcp * aup)
6739 cldfra_bl1d(k) = mf_cf
6754 if (q1 .ge. 1.0)
then
6756 elseif (q1 .ge. -1.7 .and. q1 .lt. 1.0)
then
6757 fng = exp(-0.4*(q1-1.0))
6758 elseif (q1 .ge. -2.5 .and. q1 .lt. -1.7)
then
6759 fng = 3.0 + exp(-3.8*(q1+1.7))
6761 fng = min(23.9 + exp(-1.6*(q1+2.5)), 60.)
6765 vt(k) = qww - (1.5*aup)*beta*bb*fng - 1.
6766 vq(k) = alpha + (1.5*aup)*beta*a*fng - tv0
6774 maxqc = maxval(edmf_qc(1:ktop))
6775 if ( maxqc < 1.e-8) maxmf = -1.0*maxmf
6781if (edmf_w(1) > 4.0)
then
6783 print *,
'flq:',flq,
' fltv:',fltv2
6784 print *,
'pblh:',pblh,
' wstar:',wstar
6785 print *,
'sigmaW=',sigmaw,
' sigmaTH=',sigmath,
' sigmaQT=',sigmaqt
6813 print *,
' edmf_a',edmf_a(1:14)
6814 print *,
' edmf_w',edmf_w(1:14)
6815 print *,
' edmf_qt:',edmf_qt(1:14)
6816 print *,
' edmf_thl:',edmf_thl(1:14)
6821#ifdef HARDCODE_VERTICAL
6834real(kind_phys),
intent(in) :: QT,THL,P,zagl
6835real(kind_phys),
intent(out) :: THV
6836real(kind_phys),
intent(inout):: QC
6839real(kind_phys):: diff,exn,t,th,qs,qcold
6853 exn=(p/p1000mb)**rcp
6856 t=exn*thl + xlvcp*qc
6859 qc=0.5*qc + 0.5*max((qt-qs),0.)
6860 if (abs(qc-qcold)<diff)
exit
6863 t=exn*thl + xlvcp*qc
6868 if(zagl < 100.)qc=0.
6871 thv=(thl+xlvcp*qc)*(1.+qt*(rvovrd-1.)-rvovrd*qc)
6891subroutine condensation_edmf_r(QT,THL,P,zagl,THV,QC)
6896real(kind_phys),
intent(in) :: QT,THV,P,zagl
6897real(kind_phys),
intent(out) :: THL, QC
6900real(kind_phys):: diff,exn,t,th,qs,qcold
6907 exn=(p/p1000mb)**rcp
6917 t = exn*thv/(1.+qt*(rvovrd-1.)-rvovrd*qc)
6920 if (abs(qc-qcold)<diff)
exit
6922 thl = (t - xlv/cp*qc)/exn
6924end subroutine condensation_edmf_r
6933SUBROUTINE ddmf_jpl(kts,kte,dt,zw,dz,p, &
6934 &u,v,th,thl,thv,tk,qt,qv,qc, &
6936 &ust,wthl,wqt,pblh,kpbl, &
6937 &edmf_a_dd,edmf_w_dd, edmf_qt_dd, &
6938 &edmf_thl_dd,edmf_ent_dd,edmf_qc_dd, &
6939 &sd_aw,sd_awthl,sd_awqt, &
6940 &sd_awqv,sd_awqc,sd_awu,sd_awv, &
6942 &qc_bl1d,cldfra_bl1d, &
6945 integer,
intent(in) :: KTS,KTE,KPBL
6946 real(kind_phys),
dimension(kts:kte),
intent(in) :: &
6947 U,V,TH,THL,TK,QT,QV,QC,THV,P,rho,exner,dz
6948 real(kind_phys),
dimension(kts:kte),
intent(in) :: rthraten
6950 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: ZW
6951 real(kind_phys),
intent(in) :: WTHL,WQT
6952 real(kind_phys),
intent(in) :: dt,ust,pblh
6954 real(kind_phys),
dimension(kts:kte),
intent(out) :: &
6955 edmf_a_dd,edmf_w_dd, &
6956 edmf_qt_dd,edmf_thl_dd, edmf_ent_dd,edmf_qc_dd
6959 real(kind_phys),
dimension(kts:kte+1) :: &
6960 sd_aw, sd_awthl, sd_awqt, sd_awu, &
6961 sd_awv, sd_awqc, sd_awqv, sd_awqke, sd_aw2
6963 real(kind_phys),
dimension(kts:kte),
intent(in) :: &
6964 qc_bl1d, cldfra_bl1d
6966 integer,
parameter:: ndown = 5
6968 integer,
dimension(1:NDOWN) :: DD_initK
6969 real(kind_phys),
dimension(1:NDOWN) :: randNum
6971 real(kind_phys),
dimension(kts:kte+1,1:NDOWN) :: &
6972 DOWNW,DOWNTHL,DOWNQT,DOWNQC,DOWNA,DOWNU,DOWNV,DOWNTHV
6975 real(kind_phys),
dimension(KTS+1:KTE+1,1:NDOWN) :: ENT,ENTf
6976 integer,
dimension(KTS+1:KTE+1,1:NDOWN) :: ENTi
6979 integer :: K,I,ki, kminrad, qlTop, p700_ind, qlBase
6980 real(kind_phys):: wthv,wstar,qstar,thstar,sigmaW,sigmaQT, &
6981 sigmaTH,z0,pwmin,pwmax,wmin,wmax,wlv,wtv,went,mindownw
6982 real(kind_phys):: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,Wn2,Wn, &
6983 THVk,Pk,EntEXP,EntW,beta_dm,EntExp_M,rho_int
6984 real(kind_phys):: jump_thetav, jump_qt, jump_thetal, &
6985 refTHL, refTHV, refQT
6987 real(kind_phys):: minrad,zminrad, radflux, F0, wst_rad, wst_dd
6989 real(kind_phys):: sigq,xl,rsl,cpm,a,mf_cf,diffqt, &
6990 Fng,qww,alpha,beta,bb,f,pt,t,q2p,b9,satvp,rhgrid
6993 real(kind_phys),
parameter :: &
6994 &Wa=1., Wb=1.5, Z00=100., BCOEFF=0.2
6996 real(kind_phys),
parameter :: &
7004 integer,
parameter :: debug_mf=0
7006 dl = (1000.-500.)/real(ndown)
7046 do k = max(3,kpbl-2),kpbl+3
7047 if (qc(k).gt. 1.e-6 .AND. cldfra_bl1d(k).gt.0.5)
then
7053 do k = qltop, kts, -1
7054 if (qc(k) .gt. 1e-6)
then
7058 qlbase = (qltop+qlbase)/2
7071 radflux = rthraten(k) * exner(k)
7072 radflux = radflux * cp / grav * ( p(k) - p(k+1) )
7073 if ( radflux < 0.0 ) f0 = abs(radflux) + f0
7082 adn = min( 0.05 + f0*0.001, 0.3)
7108 p700_ind = minloc(abs(p-70000),1)
7109 jump_thetav = thv(p700_ind) - thv(1) - (thv(p700_ind)-thv(qltop+3))/(zw(p700_ind)-zw(qltop+3))*(zw(p700_ind)-zw(qltop))
7110 jump_qt = qc(p700_ind) + qv(p700_ind) - qc(1) - qv(1)
7111 jump_thetal = thl(p700_ind) - thl(1) - (thl(p700_ind)-thl(qltop+3))/(zw(p700_ind)-zw(qltop+3))*(zw(p700_ind)-zw(qltop))
7118 wst_rad = ( grav * zw(qltop) * f0 / (refthl * rho(qltop) * cp) ) ** (0.333)
7119 wst_rad = max(wst_rad, 0.1)
7120 wstar = max(0.,(grav/thv(1)*wthv*pblh)**(onethird))
7121 went = thv(1) / ( grav * jump_thetav * zw(qltop) ) * &
7122 (0.15 * (wstar**3 + 5*ust**3) + 0.35 * wst_rad**3 )
7123 qstar = abs(went*jump_qt/wst_rad)
7124 thstar = f0/rho(qltop)/cp/wst_rad - went*jump_thetav/wst_rad
7126 wst_dd = (0.15 * (wstar**3 + 5*ust**3) + 0.35 * wst_rad**3 ) ** (0.333)
7128 print*,
"qstar=",qstar,
" thstar=",thstar,
" wst_dd=",wst_dd
7129 print*,
"F0=",f0,
" wst_rad=",wst_rad,
" jump_thv=",jump_thetav
7130 print*,
"entrainment velocity=",went
7133 sigmaqt = 40 * qstar
7134 sigmath = 1.0 * thstar
7143 wlv=wmin+(wmax-wmin)/real(ndown)*(i-1)
7144 wtv=wmin+(wmax-wmin)/real(ndown)*i
7151 downa(ki,i)=adn/real(ndown)
7152 downu(ki,i)=(u(ki-1)*dz(ki) + u(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7153 downv(ki,i)=(v(ki-1)*dz(ki) + v(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7160 refthl = (thl(ki-1)*dz(ki) + thl(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7161 refthv = (thv(ki-1)*dz(ki) + thv(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7162 refqt = (qt(ki-1)*dz(ki) + qt(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7165 downqc(ki,i) = (qc(ki-1)*dz(ki) + qc(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7166 downqt(ki,i) = refqt
7167 downthv(ki,i)= refthv + 0.01 *downw(ki,i)*sigmath/sigmaw
7168 downthl(ki,i)= refthl + 0.01 *downw(ki,i)*sigmath/sigmaw
7180 dp = 500. + dl*real(i)
7182 DO k=dd_initk(i)-1,kts+1,-1
7185 wmin = 0.3 + dp*0.0005
7186 ent(k,i) = 0.33/(min(max(-1.0*downw(k+1,i),wmin),0.9)*dp)
7191 entexp =ent(k,i)*dz(k)
7192 entexp_m=ent(k,i)*0.333*dz(k)
7194 qtn =downqt(k+1,i) *(1.-entexp) + qt(k)*entexp
7195 thln=downthl(k+1,i)*(1.-entexp) + thl(k)*entexp
7196 un =downu(k+1,i) *(1.-entexp) + u(k)*entexp_m
7197 vn =downv(k+1,i) *(1.-entexp) + v(k)*entexp_m
7206 pk =(p(k-1)*dz(k)+p(k)*dz(k-1))/(dz(k)+dz(k-1))
7209 thvk =(thv(k-1)*dz(k)+thv(k)*dz(k-1))/(dz(k)+dz(k-1))
7210 b=grav*(thvn/thvk - 1.0)
7221 mindownw = min(downw(k+1,i),-0.2)
7222 wn = downw(k+1,i) + (-2.*ent(k,i)*downw(k+1,i) - &
7223 bcoeff*b/mindownw)*min(dz(k), 250.)
7227 IF (wn < downw(k+1,i) - min(1.25*dz(k)/200., -2.0))
THEN
7228 wn = downw(k+1,i) - min(1.25*dz(k)/200., -2.0)
7231 IF (wn > downw(k+1,i) + min(1.25*dz(k)/200., 2.0))
THEN
7232 wn = downw(k+1,i) + min(1.25*dz(k)/200., 2.0)
7234 wn = max(min(wn,0.0), -3.0)
7243 IF (wn .lt. 0.)
THEN
7251 downa(k,i) = downa(k+1,i)
7254 if (dd_initk(i) - k .lt. 2)
then
7276 edmf_a_dd(k) =edmf_a_dd(k) +downa(k-1,i)
7277 edmf_w_dd(k) =edmf_w_dd(k) +downa(k-1,i)*downw(k-1,i)
7278 edmf_qt_dd(k) =edmf_qt_dd(k) +downa(k-1,i)*downqt(k-1,i)
7279 edmf_thl_dd(k)=edmf_thl_dd(k)+downa(k-1,i)*downthl(k-1,i)
7280 edmf_ent_dd(k)=edmf_ent_dd(k)+downa(k-1,i)*ent(k-1,i)
7281 edmf_qc_dd(k) =edmf_qc_dd(k) +downa(k-1,i)*downqc(k-1,i)
7284 IF (edmf_a_dd(k) >0.)
THEN
7285 edmf_w_dd(k) =edmf_w_dd(k) /edmf_a_dd(k)
7286 edmf_qt_dd(k) =edmf_qt_dd(k) /edmf_a_dd(k)
7287 edmf_thl_dd(k)=edmf_thl_dd(k)/edmf_a_dd(k)
7288 edmf_ent_dd(k)=edmf_ent_dd(k)/edmf_a_dd(k)
7289 edmf_qc_dd(k) =edmf_qc_dd(k) /edmf_a_dd(k)
7298 rho_int = (rho(k)*dz(k+1)+rho(k+1)*dz(k))/(dz(k+1)+dz(k))
7300 sd_aw(k) =sd_aw(k) +rho_int*downa(k,i)*downw(k,i)
7301 sd_awthl(k)=sd_awthl(k)+rho_int*downa(k,i)*downw(k,i)*downthl(k,i)
7302 sd_awqt(k) =sd_awqt(k) +rho_int*downa(k,i)*downw(k,i)*downqt(k,i)
7303 sd_awqc(k) =sd_awqc(k) +rho_int*downa(k,i)*downw(k,i)*downqc(k,i)
7304 sd_awu(k) =sd_awu(k) +rho_int*downa(k,i)*downw(k,i)*downu(k,i)
7305 sd_awv(k) =sd_awv(k) +rho_int*downa(k,i)*downw(k,i)*downv(k,i)
7307 sd_awqv(k) = sd_awqt(k) - sd_awqc(k)
7310END SUBROUTINE ddmf_jpl
7314SUBROUTINE scale_aware(dx,PBL1,Psig_bl,Psig_shcu)
7325 real(kind_phys),
intent(in) :: dx,pbl1
7326 real(kind_phys),
intent(out) :: Psig_bl,Psig_shcu
7327 real(kind_phys) :: dxdh
7331 dxdh=max(2.5*dx,10.)/min(pbl1,3000.)
7339 psig_bl= ((dxdh**2) + 0.106*(dxdh**0.667))/((dxdh**2) +0.066*(dxdh**0.667) + 0.071)
7342 dxdh=max(2.5*dx,10.)/min(pbl1+500.,3500.)
7369 psig_shcu= ((dxdh**2) + 0.145*(dxdh**0.667))/((dxdh**2) +0.172*(dxdh**0.667) + 0.170)
7379 If(psig_bl > 1.0) psig_bl=1.0
7380 If(psig_bl < 0.0) psig_bl=0.0
7382 If(psig_shcu > 1.0) psig_shcu=1.0
7383 If(psig_shcu < 0.0) psig_shcu=0.0
7385 END SUBROUTINE scale_aware
7399 real(kind_phys),
intent(in):: t
7402 real(kind_phys),
parameter:: j0= .611583699e03
7403 real(kind_phys),
parameter:: j1= .444606896e02
7404 real(kind_phys),
parameter:: j2= .143177157e01
7405 real(kind_phys),
parameter:: j3= .264224321e-1
7406 real(kind_phys),
parameter:: j4= .299291081e-3
7407 real(kind_phys),
parameter:: j5= .203154182e-5
7408 real(kind_phys),
parameter:: j6= .702620698e-8
7409 real(kind_phys),
parameter:: j7= .379534310e-11
7410 real(kind_phys),
parameter:: j8=-.321582393e-13
7412 real(kind_phys),
parameter:: k0= .609868993e03
7413 real(kind_phys),
parameter:: k1= .499320233e02
7414 real(kind_phys),
parameter:: k2= .184672631e01
7415 real(kind_phys),
parameter:: k3= .402737184e-1
7416 real(kind_phys),
parameter:: k4= .565392987e-3
7417 real(kind_phys),
parameter:: k5= .521693933e-5
7418 real(kind_phys),
parameter:: k6= .307839583e-7
7419 real(kind_phys),
parameter:: k7= .105785160e-9
7420 real(kind_phys),
parameter:: k8= .161444444e-12
7422 xc=max(-80.,t - t0c)
7427 IF (t .GE. (t0c-6.))
THEN
7428 esat_blend = j0+xc*(j1+xc*(j2+xc*(j3+xc*(j4+xc*(j5+xc*(j6+xc*(j7+xc*j8)))))))
7429 ELSE IF (t .LE. tice)
THEN
7430 esat_blend = k0+xc*(k1+xc*(k2+xc*(k3+xc*(k4+xc*(k5+xc*(k6+xc*(k7+xc*k8)))))))
7432 esl = j0+xc*(j1+xc*(j2+xc*(j3+xc*(j4+xc*(j5+xc*(j6+xc*(j7+xc*j8)))))))
7433 esi = k0+xc*(k1+xc*(k2+xc*(k3+xc*(k4+xc*(k5+xc*(k6+xc*(k7+xc*k8)))))))
7434 chi = ((t0c-6.) - t)/((t0c-6.) - tice)
7450 real(kind_phys),
intent(in):: t, p
7451 real(kind_phys)::
qsat_blend,xc,esl,esi,rslf,rsif,chi
7453 real(kind_phys),
parameter:: j0= .611583699e03
7454 real(kind_phys),
parameter:: j1= .444606896e02
7455 real(kind_phys),
parameter:: j2= .143177157e01
7456 real(kind_phys),
parameter:: j3= .264224321e-1
7457 real(kind_phys),
parameter:: j4= .299291081e-3
7458 real(kind_phys),
parameter:: j5= .203154182e-5
7459 real(kind_phys),
parameter:: j6= .702620698e-8
7460 real(kind_phys),
parameter:: j7= .379534310e-11
7461 real(kind_phys),
parameter:: j8=-.321582393e-13
7463 real(kind_phys),
parameter:: k0= .609868993e03
7464 real(kind_phys),
parameter:: k1= .499320233e02
7465 real(kind_phys),
parameter:: k2= .184672631e01
7466 real(kind_phys),
parameter:: k3= .402737184e-1
7467 real(kind_phys),
parameter:: k4= .565392987e-3
7468 real(kind_phys),
parameter:: k5= .521693933e-5
7469 real(kind_phys),
parameter:: k6= .307839583e-7
7470 real(kind_phys),
parameter:: k7= .105785160e-9
7471 real(kind_phys),
parameter:: k8= .161444444e-12
7473 xc=max(-80.,t - t0c)
7475 IF (t .GE. (t0c-6.))
THEN
7476 esl = j0+xc*(j1+xc*(j2+xc*(j3+xc*(j4+xc*(j5+xc*(j6+xc*(j7+xc*j8)))))))
7477 esl = min(esl, p*0.15)
7479 ELSE IF (t .LE. tice)
THEN
7480 esi = k0+xc*(k1+xc*(k2+xc*(k3+xc*(k4+xc*(k5+xc*(k6+xc*(k7+xc*k8)))))))
7481 esi = min(esi, p*0.15)
7484 esl = j0+xc*(j1+xc*(j2+xc*(j3+xc*(j4+xc*(j5+xc*(j6+xc*(j7+xc*j8)))))))
7485 esl = min(esl, p*0.15)
7486 esi = k0+xc*(k1+xc*(k2+xc*(k3+xc*(k4+xc*(k5+xc*(k6+xc*(k7+xc*k8)))))))
7487 esi = min(esi, p*0.15)
7488 rslf = 0.622*esl/max(p-esl, 1e-5)
7489 rsif = 0.622*esi/max(p-esi, 1e-5)
7491 chi = ((t0c-6.) - t)/((t0c-6.) - tice)
7508 real(kind_phys),
intent(in):: t
7509 real(kind_phys)::
xl_blend,xlvt,xlst,chi
7512 IF (t .GE. t0c)
THEN
7513 xl_blend = xlv + (cpv-cliq)*(t-t0c)
7514 ELSE IF (t .LE. tice)
THEN
7515 xl_blend = xls + (cpv-cice)*(t-t0c)
7517 xlvt = xlv + (cpv-cliq)*(t-t0c)
7518 xlst = xls + (cpv-cice)*(t-t0c)
7520 chi = (t0c - t)/(t0c - tice)
7521 xl_blend = (1.-chi)*xlvt + chi*xlst
7537 real(kind_phys),
intent(in):: zet
7538 real(kind_phys):: dummy_0,dummy_1,dummy_11,dummy_2,dummy_22,dummy_3,dummy_33,dummy_4,dummy_44,dummy_psi
7539 real(kind_phys),
parameter :: am_st=6.1, bm_st=2.5, rbm_st=1./bm_st
7540 real(kind_phys),
parameter :: ah_st=5.3, bh_st=1.1, rbh_st=1./bh_st
7541 real(kind_phys),
parameter :: am_unst=10., ah_unst=34.
7542 real(kind_phys):: phi_m,phim
7544 if ( zet >= 0.0 )
then
7545 dummy_0=1+zet**bm_st
7546 dummy_1=zet+dummy_0**(rbm_st)
7547 dummy_11=1+dummy_0**(rbm_st-1)*zet**(bm_st-1)
7548 dummy_2=(-am_st/dummy_1)*dummy_11
7549 phi_m = 1-zet*dummy_2
7551 dummy_0 = (1.0-cphm_unst*zet)**0.25
7553 dummy_psi = 2.*log(0.5*(1.+dummy_0))+log(0.5*(1.+dummy_0**2))-2.*atan(dummy_0)+1.570796
7555 dummy_0=(1.-am_unst*zet)
7556 dummy_1=dummy_0**0.333333
7557 dummy_11=-0.33333*am_unst*dummy_0**(-0.6666667)
7558 dummy_2 = 0.33333*(dummy_1**2.+dummy_1+1.)
7559 dummy_22 = 0.3333*dummy_11*(2.*dummy_1+1.)
7560 dummy_3 = 0.57735*(2.*dummy_1+1.)
7561 dummy_33 = 1.1547*dummy_11
7562 dummy_4 = 1.5*log(dummy_2)-1.73205*atan(dummy_3)+1.813799364
7563 dummy_44 = (1.5/dummy_2)*dummy_22-1.73205*dummy_33/(1.+dummy_3**2)
7566 dummy_1 = 1./(1.+dummy_0)
7568 dummy_2 = ((1-phi_m)/zet+dummy_11*dummy_4+dummy_0*dummy_44)*dummy_1
7569 dummy_22 = -dummy_11*(dummy_psi+dummy_0*dummy_4)*dummy_1**2
7571 phi_m = 1.-zet*(dummy_2+dummy_22)
7589 real(kind_phys),
intent(in):: zet
7590 real(kind_phys):: dummy_0,dummy_1,dummy_11,dummy_2,dummy_22,dummy_3,dummy_33,dummy_4,dummy_44,dummy_psi
7591 real(kind_phys),
parameter :: am_st=6.1, bm_st=2.5, rbm_st=1./bm_st
7592 real(kind_phys),
parameter :: ah_st=5.3, bh_st=1.1, rbh_st=1./bh_st
7593 real(kind_phys),
parameter :: am_unst=10., ah_unst=34.
7594 real(kind_phys):: phh,phih
7596 if ( zet >= 0.0 )
then
7597 dummy_0=1+zet**bh_st
7598 dummy_1=zet+dummy_0**(rbh_st)
7599 dummy_11=1+dummy_0**(rbh_st-1)*zet**(bh_st-1)
7600 dummy_2=(-ah_st/dummy_1)*dummy_11
7601 phih = 1-zet*dummy_2
7603 dummy_0 = (1.0-cphh_unst*zet)**0.5
7605 dummy_psi = 2.*log(0.5*(1.+dummy_0))
7607 dummy_0=(1.-ah_unst*zet)
7608 dummy_1=dummy_0**0.333333
7609 dummy_11=-0.33333*ah_unst*dummy_0**(-0.6666667)
7610 dummy_2 = 0.33333*(dummy_1**2.+dummy_1+1.)
7611 dummy_22 = 0.3333*dummy_11*(2.*dummy_1+1.)
7612 dummy_3 = 0.57735*(2.*dummy_1+1.)
7613 dummy_33 = 1.1547*dummy_11
7614 dummy_4 = 1.5*log(dummy_2)-1.73205*atan(dummy_3)+1.813799364
7615 dummy_44 = (1.5/dummy_2)*dummy_22-1.73205*dummy_33/(1.+dummy_3**2)
7618 dummy_1 = 1./(1.+dummy_0)
7620 dummy_2 = ((1-phh)/zet+dummy_11*dummy_4+dummy_0*dummy_44)*dummy_1
7621 dummy_22 = -dummy_11*(dummy_psi+dummy_0*dummy_4)*dummy_1**2
7623 phih = 1.-zet*(dummy_2+dummy_22)
7628 SUBROUTINE topdown_cloudrad(kts,kte,dz1,zw,xland,kpbl,PBLH, &
7629 &sqc,sqi,sqw,thl,th1,ex1,p1,rho1,thetav, &
7630 &cldfra_bl1D,rthraten, &
7631 &maxKHtopdown,KHtopdown,TKEprodTD )
7634 integer,
intent(in) :: kte,kts
7635 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz1,sqc,sqi,sqw,&
7636 thl,th1,ex1,p1,rho1,thetav,cldfra_bl1d
7637 real(kind_phys),
dimension(kts:kte),
intent(in) :: rthraten
7638 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
7639 real(kind_phys),
intent(in) :: pblh
7640 real(kind_phys),
intent(in) :: xland
7641 integer ,
intent(in) :: kpbl
7643 real(kind_phys),
intent(out) :: maxkhtopdown
7644 real(kind_phys),
dimension(kts:kte),
intent(out) :: khtopdown,tkeprodtd
7646 real(kind_phys),
dimension(kts:kte) :: zfac,wscalek2,zfacent
7647 real(kind_phys) :: bfx0,sflux,wm2,wm3,h1,h2,bfxpbl,dthvx,tmp1
7648 real(kind_phys) :: temps,templ,zl1,wstar3_2
7649 real(kind_phys) :: ent_eff,radsum,radflux,we,rcldb,rvls,minrad,zminrad
7650 real(kind_phys),
parameter :: pfac =2.0, zfmin = 0.01, phifac=8.0
7651 integer :: k,kk,kminrad
7658 khtopdown(kts:kte)=0.0
7659 tkeprodtd(kts:kte)=0.0
7663 DO kk = max(1,kpbl-2),kpbl+3
7664 if (sqc(kk).gt. 1.e-6 .OR. sqi(kk).gt. 1.e-6 .OR. &
7665 cldfra_bl1d(kk).gt.0.5)
then
7668 if (rthraten(kk) < minrad)
then
7671 zminrad=zw(kk) + 0.5*dz1(kk)
7675 IF (max(kminrad,kpbl) < 2)cloudflg = .false.
7678 k = max(kpbl-1, kminrad-1)
7684 rvls=100.*6.112*exp(17.67*(templ-273.16)/(templ-29.65))*(ep_2/p1(k+1))
7685 temps=templ + (sqw(k)-rvls)/(cp/xlv + ep_2*xlv*rvls/(r_d*templ**2))
7686 rvls=100.*6.112*exp(17.67*(temps-273.15)/(temps-29.65))*(ep_2/p1(k+1))
7687 rcldb=max(sqw(k)-rvls,0.)
7690 dthvx = (thl(k+2) + th1(k+2)*p608*sqw(k+2)) &
7691 - (thl(k) + th1(k) *p608*sqw(k))
7692 dthvx = max(dthvx,0.1)
7693 tmp1 = xlvcp * rcldb/(ex1(k)*dthvx)
7696 ent_eff = 0.2 + 0.2*8.*tmp1
7699 DO kk = max(1,kpbl-3),kpbl+3
7700 radflux=rthraten(kk)*ex1(kk)
7701 radflux=radflux*cp/grav*(p1(kk)-p1(kk+1))
7702 if (radflux < 0.0 ) radsum=abs(radflux)+radsum
7706 if ((xland-1.5).GE.0)
THEN
7707 radsum=min(radsum,90.0)
7708 bfx0 = max(radsum/rho1(k)/cp,0.)
7710 radsum=min(0.25*radsum,30.0)
7711 bfx0 = max(radsum/rho1(k)/cp - max(sflux,0.0),0.)
7715 wm3 = grav/thetav(k)*bfx0*min(pblh,1500.)
7717 bfxpbl = - ent_eff * bfx0
7718 dthvx = max(thetav(k+1)-thetav(k),0.1)
7719 we = max(bfxpbl/dthvx,-sqrt(wm3**h2))
7723 zfac(kk) = min(max((1.-(zw(kk+1)-zl1)/(zminrad-zl1)),zfmin),1.)
7724 zfacent(kk) = 10.*max((zminrad-zw(kk+1))/zminrad,0.0)*(1.-zfac(kk))**3
7727 wscalek2(kk) = (phifac*karman*wm3*(zfac(kk)))**h1
7729 khtopdown(kk) = wscalek2(kk)*karman*(zminrad-zw(kk+1))*(1.-zfac(kk))**3
7730 khtopdown(kk) = max(khtopdown(kk),0.0)
7735 tkeprodtd(kk)=2.*ent_eff*wm3/max(pblh,100.)*zfacent(kk)
7736 tkeprodtd(kk)= max(tkeprodtd(kk),0.0)
7739 maxkhtopdown=maxval(khtopdown(:))
7741 END SUBROUTINE topdown_cloudrad
subroutine mym_turbulence(kts, kte, xland, closure, dz, dx, zw, u, v, thl, thetav, ql, qw, qke, tsq, qsq, cov, vt, vq, rmo, flt, fltv, flq, zi, theta, sh, sm, el, dfm, dfh, dfq, tcd, qcd, pdk, pdt, pdq, pdc, qwt1d, qshear1d, qbuoy1d, qdiss1d, tke_budget, psig_bl, psig_shcu, cldfra_bl1d, bl_mynn_mixlength, edmf_w1, edmf_a1, tkeprodtd, spp_pbl, rstoch_col)
This subroutine calculates the vertical diffusivity coefficients and the production terms for the tur...
subroutine get_pblh(kts, kte, zi, thetav1d, qke1d, zw1d, dz1d, landsea, kzi)
This subroutine calculates hybrid diagnotic boundary-layer height (PBLH).
subroutine dmp_mf(kts, kte, dt, zw, dz, p, rho, momentum_opt, tke_opt, scalar_opt, u, v, w, th, thl, thv, tk, qt, qv, qc, qke, qnc, qni, qnwfa, qnifa, qnbca, exner, vt, vq, sgm, ust, flt, fltv, flq, flqv, pblh, kpbl, dx, landsea, ts, edmf_a, edmf_w, edmf_qt, edmf_thl, edmf_ent, edmf_qc, s_aw, s_awthl, s_awqt, s_awqv, s_awqc, s_awu, s_awv, s_awqke, s_awqnc, s_awqni, s_awqnwfa, s_awqnifa, s_awqnbca, sub_thl, sub_sqv, sub_u, sub_v, det_thl, det_sqv, det_sqc, det_u, det_v, nchem, chem1, s_awchem, mix_chem, qc_bl1d, cldfra_bl1d, qc_bl1d_old, cldfra_bl1d_old, f_qc, f_qi, f_qnc, f_qni, f_qnwfa, f_qnifa, f_qnbca, psig_shcu, maxwidth, ktop, maxmf, ztop, spp_pbl, rstoch_col)
This subroutine is the Dynamic Multi-Plume (DMP) Mass-Flux Scheme.
subroutine mym_length(kts, kte, xland, dz, dx, zw, rmo, flt, fltv, flq, vt, vq, u1, v1, qke, dtv, el, zi, theta, qkw, psig_bl, cldfra_bl1d, bl_mynn_mixlength, edmf_w1, edmf_a1)
This subroutine calculates the mixing lengths.
subroutine mym_level2(kts, kte, dz, u, v, thl, thetav, qw, ql, vt, vq, dtl, dqw, dtv, gm, gh, sm, sh)
This subroutine calculates the level 2, non-dimensional wind shear and vertical temperature gradient...
real(kind_phys) function qsat_blend(t, p)
This function extends function "esat" and returns a "blended" saturation mixing ratio....
subroutine mynn_bl_driver(initflag, restart, cycling, delt, dz, dx, znt, u, v, w, th, sqv3d, sqc3d, sqi3d, sqs3d, qnc, qni, qnwfa, qnifa, qnbca, ozone, p, exner, rho, t3d, xland, ts, qsfc, ps, ust, ch, hfx, qfx, rmol, wspd, uoce, voce, qke, qke_adv, sh3d, sm3d, nchem, kdvel, ndvel, chem3d, vdep, smoke_dbg, frp, emis_ant_no, mix_chem, enh_mix, rrfs_sd, tsq, qsq, cov, rublten, rvblten, rthblten, rqvblten, rqcblten, rqiblten, rqncblten, rqniblten, rqsblten, rqnwfablten, rqnifablten, rqnbcablten, dozone, exch_h, exch_m, pblh, kpbl, el_pbl, dqke, qwt, qshear, qbuoy, qdiss, qc_bl, qi_bl, cldfra_bl, bl_mynn_tkeadvect, tke_budget, bl_mynn_cloudpdf, bl_mynn_mixlength, icloud_bl, closure, bl_mynn_edmf, bl_mynn_edmf_mom, bl_mynn_edmf_tke, bl_mynn_mixscalars, bl_mynn_output, bl_mynn_cloudmix, bl_mynn_mixqt, edmf_a, edmf_w, edmf_qt, edmf_thl, edmf_ent, edmf_qc, sub_thl3d, sub_sqv3d, det_thl3d, det_sqv3d, maxwidth, maxmf, ztop_plume, ktop_plume, spp_pbl, pattern_spp_pbl, rthraten, flag_qc, flag_qi, flag_qnc, flag_qni, flag_qs, flag_qnwfa, flag_qnifa, flag_qnbca, flag_ozone, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte)
This subroutine is the MYNN-EDNF PBL driver routine,which encompassed the majority of the subroutines...
subroutine mym_condensation(kts, kte, dx, dz, zw, xland, thl, qw, qv, qc, qi, qs, p, exner, tsq, qsq, cov, sh, el, bl_mynn_cloudpdf, qc_bl1d, qi_bl1d, cldfra_bl1d, pblh1, hfx1, vt, vq, th, sgm, rmo, spp_pbl, rstoch_col)
This subroutine calculates the nonconvective component of the subgrid cloud fraction and mixing ratio...
subroutine boulac_length0(k, kts, kte, zw, dz, qtke, theta, lb1, lb2)
This subroutine was taken from the BouLac scheme in WRF-ARW and modified for integration into the MYN...
subroutine condensation_edmf(qt, thl, p, zagl, thv, qc)
This subroutine.
subroutine mym_initialize(kts, kte, xland, dz, dx, zw, u, v, thl, qw, zi, theta, thetav, sh, sm, ust, rmo, el, qke, tsq, qsq, cov, psig_bl, cldfra_bl1d, bl_mynn_mixlength, edmf_w1, edmf_a1, initialize_qke, spp_pbl, rstoch_col)
This subroutine initializes the mixing length, TKE, , , and .
real(kind_phys) function xl_blend(t)
This function interpolates the latent heats of vaporization and sublimation into a single,...
subroutine boulac_length(kts, kte, zw, dz, qtke, theta, lb1, lb2)
This subroutine was taken from the BouLac scheme in WRF-ARW and modified for integration into the MYN...
subroutine mynn_tendencies(kts, kte, i, delt, dz, rho, u, v, th, tk, qv, qc, qi, qs, qnc, qni, psfc, p, exner, thl, sqv, sqc, sqi, sqs, sqw, qnwfa, qnifa, qnbca, ozone, ust, flt, flq, flqv, flqc, wspd, uoce, voce, tsq, qsq, cov, tcd, qcd, dfm, dfh, dfq, du, dv, dth, dqv, dqc, dqi, dqs, dqnc, dqni, dqnwfa, dqnifa, dqnbca, dozone, diss_heat, s_aw, s_awthl, s_awqt, s_awqv, s_awqc, s_awu, s_awv, s_awqnc, s_awqni, s_awqnwfa, s_awqnifa, s_awqnbca, sd_aw, sd_awthl, sd_awqt, sd_awqv, sd_awqc, sd_awu, sd_awv, sub_thl, sub_sqv, sub_u, sub_v, det_thl, det_sqv, det_sqc, det_u, det_v, flag_qc, flag_qi, flag_qnc, flag_qni, flag_qs, flag_qnwfa, flag_qnifa, flag_qnbca, cldfra_bl1d, bl_mynn_cloudmix, bl_mynn_mixqt, bl_mynn_edmf, bl_mynn_edmf_mom, bl_mynn_mixscalars)
This subroutine solves for tendencies of U, V, , qv, qc, and qi.
real(kind_phys) function esat_blend(t)
subroutine mym_predict(kts, kte, closure, delt, dz, ust, flt, flq, pmz, phh, el, dfq, rho, pdk, pdt, pdq, pdc, qke, tsq, qsq, cov, s_aw, s_awqke, bl_mynn_edmf_tke, qwt1d, qdiss1d, tke_budget)
This subroutine predicts the turbulent quantities at the next step.
This module defines model-specific constants/parameters.
This module contains the entity of MYNN-EDMF PBL scheme.