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),
optional :: &
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),
optional :: 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),
optional :: &
499 real(kind_phys),
dimension(:),
intent(out),
optional :: &
500 &maxmf,maxwidth,ztop_plume
502 real(kind_phys),
dimension(:,:),
intent(inout),
optional :: 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),
optional :: Sh3D,Sm3D
513 real(kind_phys),
dimension(:,:),
intent(inout),
optional :: &
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))
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
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),1.0e-10))
1939 afk = dz(k)/( dz(k)+dz(k-1) )
1941 qkw(k) = sqrt(max(qke(k)*abk+qke(k-1)*afk,1.0e-3))
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 ) )**(1.0/3.0)
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.01)
2020 thetaw(kts)=theta(kts)
2021 qkw(kts) = sqrt(max(qke(kts),1.0e-10))
2024 afk = dz(k)/( dz(k)+dz(k-1) )
2026 qkw(k) = sqrt(max(qke(k)*abk+qke(k-1)*afk,1.0e-3))
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.03 ), 30.0)*dzk
2046 elt = min( max( alp1*elt/vsc, 10.), 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.01)
2123 qkw(kts) = sqrt(max(qke(kts),1.0e-4))
2126 afk = dz(k)/( dz(k)+dz(k-1) )
2128 qkw(k) = sqrt(max(qke(k)*abk+qke(k-1)*afk,1.0e-3))
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
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
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
3684 k_tropo = max(kts+2, k+2)
3688 SELECT CASE(bl_mynn_cloudpdf)
3709 qsl=ep_2*esat/max(1.e-4,(p(k)-ep_3*esat))
3711 dqsl = qsl*ep_2*xlv/( r_d*t**2 )
3713 alp(k) = 1.0/( 1.0+dqsl*xlvcp )
3714 bet(k) = dqsl*exner(k)
3718 t3sq = max( tsq(k), 0.0 )
3719 r3sq = max( qsq(k), 0.0 )
3721 c3sq = sign( min( abs(c3sq), sqrt(t3sq*r3sq) ), c3sq )
3722 r3sq = r3sq +bet(k)**2*t3sq -2.0*bet(k)*c3sq
3726 sgm(k) = sqrt( max( r3sq, 1.0d-10 ))
3728 q1(k) = qmq / sgm(k)
3730 cldfra_bl1d(k) = 0.5*( 1.0+erf( q1(k)*rr2 ) )
3733 eq1 = rrp*exp( -0.5*q1k*q1k )
3734 qll = max( cldfra_bl1d(k)*q1k + eq1, 0.0 )
3736 ql(k) = alp(k)*sgm(k)*qll
3738 liq_frac = min(1.0, max(0.0,(t-240.0)/29.0))
3739 qc_bl1d(k) = liq_frac*ql(k)
3740 qi_bl1d(k) = (1.0 - liq_frac)*ql(k)
3743 q2p = xlvcp/exner(k)
3744 pt = thl(k) +q2p*ql(k)
3747 qt = 1.0 +p608*qw(k) -(1.+p608)*(qc_bl1d(k)+qi_bl1d(k))*cldfra_bl1d(k)
3748 rac = alp(k)*( cldfra_bl1d(k)-qll*eq1 )*( q2p*qt-(1.+p608)*pt )
3753 vt(k) = qt-1.0 -rac*bet(k)
3754 vq(k) = p608*pt-tv0 +rac
3766 qsl=ep_2*esat/max(1.e-4,(p(k)-ep_3*esat))
3768 dqsl = qsl*ep_2*xlv/( r_d*t**2 )
3770 alp(k) = 1.0/( 1.0+dqsl*xlvcp )
3771 bet(k) = dqsl*exner(k)
3773 if (k .eq. kts)
then
3778 dth = 0.5*(thl(k+1)+thl(k)) - 0.5*(thl(k)+thl(max(k-1,kts)))
3779 dqw = 0.5*(qw(k+1) + qw(k)) - 0.5*(qw(k) + qw(max(k-1,kts)))
3780 sgm(k) = sqrt( max( (alp(k)**2 * max(el(k)**2,0.1) * &
3781 b2 * max(sh(k),0.03))/4. * &
3782 (dqw/dzk - bet(k)*(dth/dzk ))**2 , 1.0e-10) )
3784 q1(k) = qmq / sgm(k)
3785 cldfra_bl1d(k) = 0.5*( 1.0+erf( q1(k)*rr2 ) )
3791 eq1 = rrp*exp( -0.5*q1k*q1k )
3792 qll = max( cldfra_bl1d(k)*q1k + eq1, 0.0 )
3794 ql(k) = alp(k)*sgm(k)*qll
3795 liq_frac = min(1.0, max(0.0,(t-240.0)/29.0))
3796 qc_bl1d(k) = liq_frac*ql(k)
3797 qi_bl1d(k) = (1.0 - liq_frac)*ql(k)
3800 q2p = xlvcp/exner(k)
3801 pt = thl(k) +q2p*ql(k)
3804 qt = 1.0 +p608*qw(k) -(1.+p608)*(qc_bl1d(k)+qi_bl1d(k))*cldfra_bl1d(k)
3805 rac = alp(k)*( cldfra_bl1d(k)-qll*eq1 )*( q2p*qt-(1.+p608)*pt )
3810 vt(k) = qt-1.0 -rac*bet(k)
3811 vq(k) = p608*pt-tv0 +rac
3819 pblh2=max(10._kind_phys,pblh1)
3823 zagl = zagl + 0.5*(dz(k) + dzm1)
3829 rh(k) = max(min(rhmax, qw(k)/max(1.e-10,qsat_tk)),0.001_kind_phys)
3832 dqsl = qsat_tk*ep_2*xlv/( r_d*t**2 )
3833 alp(k) = 1.0/( 1.0+dqsl*xlvcp )
3834 bet(k) = dqsl*exner(k)
3836 rsl = xl*qsat_tk / (r_v*t**2)
3838 cpm = cp + qw(k)*cpv
3839 a(k) = 1./(1. + xl*rsl/cpm)
3843 qw_pert= qw(k) + qw(k)*0.5*rstoch_col(k)*real(spp_pbl)
3846 qmq = qw_pert - qsat_tk
3850 r3sq = max( qsq(k), 0.0 )
3852 sgm(k) = sqrt( r3sq )
3854 sgm(k) = min( sgm(k), qsat_tk*0.666 )
3858 wt = max(500. - max(dz(k)-100.,0.0), 0.0_kind_phys)/500.
3859 sgm(k) = sgm(k) + sgm(k)*0.2*(1.0-wt)
3862 qpct = qpct_pbl*wt + qpct_trp*(1.0-wt)
3863 qpct = min(qpct, max(qpct_sfc, qpct_pbl*zagl/500.) )
3864 sgm(k) = max( sgm(k), qsat_tk*qpct )
3866 q1(k) = qmq / sgm(k)
3871 wt2 = min(max( zagl - pblh2, 0.0 )/300., 1.0)
3873 if ((qi(k)+qs(k))>1.e-9 .and. (zagl .gt. pblh2))
then
3874 rh_hack =min(rhmax, rhcrit + wt2*0.045*(9.0 + log10(qi(k)+qs(k))))
3875 rh(k) =max(rh(k), rh_hack)
3877 q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
3878 q1(k) =max(q1_rh, q1(k) )
3881 if (qc(k)>1.e-6 .and. (zagl .gt. pblh2))
then
3882 rh_hack =min(rhmax, rhcrit + wt2*0.08*(6.0 + log10(qc(k))))
3883 rh(k) =max(rh(k), rh_hack)
3885 q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
3886 q1(k) =max(q1_rh, q1(k) )
3897 cldfra_bl1d(k) = max(0., min(1., 0.5+0.36*atan(1.8*(q1(k)+0.2))))
3902 maxqc = max(qw(k) - qsat_tk, 0.0)
3904 ql_water = sgm(k)*exp(1.2*q1k-1.)
3905 ql_ice = sgm(k)*exp(1.2*q1k-1.)
3906 elseif (q1k > 2.)
then
3907 ql_water = min(sgm(k)*q1k, maxqc)
3910 ql_water = min(sgm(k)*(exp(-1.) + 0.66*q1k + 0.086*q1k**2), maxqc)
3911 ql_ice = sgm(k)*(exp(-1.) + 0.66*q1k + 0.086*q1k**2)
3919 if (cldfra_bl1d(k) < 0.001)
then
3922 cldfra_bl1d(k) = 0.0
3925 liq_frac = min(1.0, max(0.0, (t-tice)/(tliq-tice)))
3926 qc_bl1d(k) = liq_frac*ql_water
3927 qi_bl1d(k) = (1.0-liq_frac)*ql_ice
3931 if (k .ge. k_tropo)
then
3940 if ((xland-1.5).GE.0)
then
3957 if (q1k .ge. 1.0)
then
3959 elseif (q1k .ge. -1.7 .and. q1k .lt. 1.0)
then
3960 fng = exp(-0.4*(q1k-1.0))
3961 elseif (q1k .ge. -2.5 .and. q1k .lt. -1.7)
then
3962 fng = 3.0 + exp(-3.8*(q1k+1.7))
3964 fng = min(23.9 + exp(-1.6*(q1k+2.5)), 60._kind_phys)
3967 cfmax = min(cldfra_bl1d(k), 0.6_kind_phys)
3969 zsl = min(max(25., 0.1*pblh2), 100.)
3970 wt = min(zagl/zsl, 1.0)
3981 beta = (th(k)/t)*(xl/cp) - 1.61*th(k)
3982 vt(k) = qww - cfmax*beta*bb*fng - 1.
3983 vq(k) = alpha + cfmax*beta*a(k)*fng - tv0
3991 fac_damp = min(zagl * 0.0025, 1.0)
3994 cld_factor = 1.0 + fac_damp*min((max(0.0, ( rh(k) - 0.92 )) / 0.145)**2, 0.37)
3995 cldfra_bl1d(k) = min( 1., cld_factor*cldfra_bl1d(k) )
4001 IF (bl_mynn_cloudpdf .LT. 0)
THEN
4003 cldfra_bl1d(k) = 0.0
4017#ifdef HARDCODE_VERTICAL
4030 &u,v,th,tk,qv,qc,qi,qs,qnc,qni, &
4032 &thl,sqv,sqc,sqi,sqs,sqw, &
4033 &qnwfa,qnifa,qnbca,ozone, &
4034 &ust,flt,flq,flqv,flqc,wspd, &
4039 &Du,Dv,Dth,Dqv,Dqc,Dqi,Dqs,Dqnc,Dqni, &
4040 &Dqnwfa,Dqnifa,Dqnbca,Dozone, &
4042 &s_aw,s_awthl,s_awqt,s_awqv,s_awqc, &
4045 &s_awqnwfa,s_awqnifa,s_awqnbca, &
4046 &sd_aw,sd_awthl,sd_awqt,sd_awqv, &
4047 &sd_awqc,sd_awu,sd_awv, &
4050 &det_thl,det_sqv,det_sqc, &
4052 &FLAG_QC,FLAG_QI,FLAG_QNC,FLAG_QNI, &
4054 &FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA, &
4056 &bl_mynn_cloudmix, &
4059 &bl_mynn_edmf_mom, &
4060 &bl_mynn_mixscalars )
4063 integer,
intent(in) :: kts,kte,i
4065#ifdef HARDCODE_VERTICAL
4067# define kte HARDCODE_VERTICAL
4070 integer,
intent(in) :: bl_mynn_cloudmix,bl_mynn_mixqt, &
4071 bl_mynn_edmf,bl_mynn_edmf_mom, &
4073 logical,
intent(IN) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QS, &
4074 &FLAG_QNC,FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA
4083 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: s_aw, &
4084 &s_awthl,s_awqt,s_awqnc,s_awqni,s_awqv,s_awqc,s_awu,s_awv, &
4085 &s_awqnwfa,s_awqnifa,s_awqnbca, &
4086 &sd_aw,sd_awthl,sd_awqt,sd_awqv,sd_awqc,sd_awu,sd_awv
4088 real(kind_phys),
dimension(kts:kte),
intent(in) :: sub_thl,sub_sqv, &
4089 &sub_u,sub_v,det_thl,det_sqv,det_sqc,det_u,det_v
4090 real(kind_phys),
dimension(kts:kte),
intent(in) :: u,v,th,tk,qv,qc,qi,&
4091 &qs,qni,qnc,rho,p,exner,dfq,dz,tsq,qsq,cov,tcd,qcd, &
4092 &cldfra_bl1d,diss_heat
4093 real(kind_phys),
dimension(kts:kte),
intent(inout) :: thl,sqw,sqv,sqc,&
4094 &sqi,sqs,qnwfa,qnifa,qnbca,ozone,dfm,dfh
4095 real(kind_phys),
dimension(kts:kte),
intent(inout) :: du,dv,dth,dqv, &
4096 &dqc,dqi,dqs,dqni,dqnc,dqnwfa,dqnifa,dqnbca,dozone
4097 real(kind_phys),
intent(in) :: flt,flq,flqv,flqc,uoce,voce
4098 real(kind_phys),
intent(in) :: ust,delt,psfc,wspd
4100 real(kind_phys):: wsp,wsp2,tk2,th2
4108 real(kind_phys),
dimension(kts:kte) :: dtz,dfhc,dfmc,delp
4109 real(kind_phys),
dimension(kts:kte) :: sqv2,sqc2,sqi2,sqs2,sqw2, &
4110 &qni2,qnc2,qnwfa2,qnifa2,qnbca2,ozone2
4111 real(kind_phys),
dimension(kts:kte) :: zfac,plumeKh,rhoinv
4112 real(kind_phys),
dimension(kts:kte) :: a,b,c,d,x
4113 real(kind_phys),
dimension(kts:kte+1) :: rhoz, &
4115 real(kind_phys):: rhs,gfluxm,gfluxp,dztop,maxdfh,mindfh,maxcf,maxKh,zw
4116 real(kind_phys):: t,esat,qsl,onoff,kh,km,dzk,rhosfc
4117 real(kind_phys):: ustdrag,ustdiff,qvflux
4118 real(kind_phys):: th_new,portion_qc,portion_qi,condensate,qsat
4123 real(kind_phys),
parameter :: nonloc = 1.0
4125 dztop=.5*(dz(kte)+dz(kte-1))
4130 IF (bl_mynn_edmf_mom == 0)
THEN
4138 rhosfc = psfc/(r_d*(tk(kts)+p608*qv(kts)))
4139 dtz(kts) =delt/dz(kts)
4141 rhoinv(kts)=1./rho(kts)
4142 khdz(kts) =rhoz(kts)*dfh(kts)
4143 kmdz(kts) =rhoz(kts)*dfm(kts)
4144 delp(kts) = psfc - (p(kts+1)*dz(kts) + p(kts)*dz(kts+1))/(dz(kts)+dz(kts+1))
4147 rhoz(k) =(rho(k)*dz(k-1) + rho(k-1)*dz(k))/(dz(k-1)+dz(k))
4148 rhoz(k) = max(rhoz(k),1e-4)
4149 rhoinv(k)=1./max(rho(k),1e-4)
4150 dzk = 0.5 *( dz(k)+dz(k-1) )
4151 khdz(k) = rhoz(k)*dfh(k)
4152 kmdz(k) = rhoz(k)*dfm(k)
4155 delp(k) = (p(k)*dz(k-1) + p(k-1)*dz(k))/(dz(k)+dz(k-1)) - &
4156 (p(k+1)*dz(k) + p(k)*dz(k+1))/(dz(k)+dz(k+1))
4158 delp(kte) =delp(kte-1)
4159 rhoz(kte+1)=rhoz(kte)
4160 khdz(kte+1)=rhoz(kte+1)*dfh(kte)
4161 kmdz(kte+1)=rhoz(kte+1)*dfm(kte)
4165 khdz(k) = max(khdz(k), 0.5*s_aw(k))
4166 khdz(k) = max(khdz(k), -0.5*(s_aw(k)-s_aw(k+1)))
4167 kmdz(k) = max(kmdz(k), 0.5*s_aw(k))
4168 kmdz(k) = max(kmdz(k), -0.5*(s_aw(k)-s_aw(k+1)))
4171 ustdrag = min(ust*ust,0.99)/wspd
4172 ustdiff = min(ust*ust,0.01)/wspd
4182 a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
4183 b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) &
4184 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4185 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4186 c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
4187 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4188 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4189 d(k)=u(k) + dtz(k)*uoce*ust**2/wspd &
4190 & - dtz(k)*rhoinv(k)*s_awu(k+1)*onoff &
4191 & + dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff &
4192 & + sub_u(k)*delt + det_u(k)*delt
4195 a(k)= -dtz(k)*kmdz(k)*rhoinv(k) &
4196 & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff &
4197 & + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
4198 b(k)=1.+ dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) &
4199 & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff &
4200 & + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
4201 c(k)= - dtz(k)*kmdz(k+1)*rhoinv(k) &
4202 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4203 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4204 d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff &
4205 & - dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff &
4206 & + sub_u(k)*delt + det_u(k)*delt
4233 du(k)=(x(k)-u(k))/delt
4243 a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
4244 b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) &
4245 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4246 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4247 c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
4248 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4249 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4250 d(k)=v(k) + dtz(k)*voce*ust**2/wspd &
4251 & - dtz(k)*rhoinv(k)*s_awv(k+1)*onoff &
4252 & + dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff &
4253 & + sub_v(k)*delt + det_v(k)*delt
4256 a(k)= -dtz(k)*kmdz(k)*rhoinv(k) &
4257 & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff &
4258 & + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
4259 b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) &
4260 & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff &
4261 & + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
4262 c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
4263 & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
4264 & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
4265 d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff &
4266 & - dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff &
4267 & + sub_v(k)*delt + det_v(k)*delt
4294 dv(k)=(x(k)-v(k))/delt
4319 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4320 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)
4321 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)
4322 d(k)=thl(k) + dtz(k)*rhosfc*flt*rhoinv(k) + tcd(k)*delt &
4323 & - dtz(k)*rhoinv(k)*s_awthl(k+1) -dtz(k)*rhoinv(k)*sd_awthl(k+1) + &
4324 & diss_heat(k)*delt + sub_thl(k)*delt + det_thl(k)*delt
4327 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)
4328 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4329 & 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))
4330 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)
4331 d(k)=thl(k) + tcd(k)*delt + &
4332 & dtz(k)*rhoinv(k)*(s_awthl(k)-s_awthl(k+1)) + dtz(k)*rhoinv(k)*(sd_awthl(k)-sd_awthl(k+1)) + &
4333 & diss_heat(k)*delt + &
4334 & sub_thl(k)*delt + det_thl(k)*delt
4365IF (bl_mynn_mixqt > 0)
THEN
4389 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4390 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)
4391 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)
4392 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)
4395 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)
4396 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4397 & 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))
4398 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)
4399 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))
4430IF (bl_mynn_mixqt == 0)
THEN
4435 IF (bl_mynn_cloudmix > 0 .AND. flag_qc)
THEN
4454 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4455 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)
4456 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)
4457 d(k)=sqc(k) + dtz(k)*rhosfc*flqc*rhoinv(k) + qcd(k)*delt &
4458 & - dtz(k)*rhoinv(k)*s_awqc(k+1) - dtz(k)*rhoinv(k)*sd_awqc(k+1) + &
4462 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)
4463 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4464 & 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))
4465 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)
4466 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)) + &
4489IF (bl_mynn_mixqt == 0)
THEN
4513 if (qvflux < 0.0)
then
4515 qvflux = max(qvflux, (min(0.9*sqv(kts) - 1e-8, 0.0)/dtz(kts)))
4519 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4520 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)
4521 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)
4522 d(k)=sqv(k) + dtz(k)*rhosfc*qvflux*rhoinv(k) + qcd(k)*delt &
4523 & - dtz(k)*rhoinv(k)*s_awqv(k+1) - dtz(k)*rhoinv(k)*sd_awqv(k+1) + &
4524 & sub_sqv(k)*delt + det_sqv(k)*delt
4527 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)
4528 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4529 & 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))
4530 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)
4531 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)) + &
4532 & sub_sqv(k)*delt + det_sqv(k)*delt
4568IF (bl_mynn_cloudmix > 0 .AND. flag_qi)
THEN
4572 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4573 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k)
4574 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4578 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4579 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k)
4580 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4618IF (bl_mynn_cloudmix > 0 .AND. .false.)
THEN
4622 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4623 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k)
4624 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4628 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4629 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k)
4630 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4654IF (bl_mynn_cloudmix > 0 .AND. flag_qni .AND. &
4655 bl_mynn_mixscalars > 0)
THEN
4659 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4660 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4661 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4662 d(k)=qni(k) - dtz(k)*rhoinv(k)*s_awqni(k+1)*nonloc
4665 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4666 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4667 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4668 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4669 d(k)=qni(k) + dtz(k)*rhoinv(k)*(s_awqni(k)-s_awqni(k+1))*nonloc
4695 IF (bl_mynn_cloudmix > 0 .AND. flag_qnc .AND. &
4696 bl_mynn_mixscalars > 0)
THEN
4700 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4701 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4702 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4703 d(k)=qnc(k) - dtz(k)*rhoinv(k)*s_awqnc(k+1)*nonloc
4706 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4707 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k) + &
4708 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4709 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4710 d(k)=qnc(k) + dtz(k)*rhoinv(k)*(s_awqnc(k)-s_awqnc(k+1))*nonloc
4735IF (bl_mynn_cloudmix > 0 .AND. flag_qnwfa .AND. &
4736 bl_mynn_mixscalars > 0)
THEN
4740 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4741 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) - &
4742 & 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4743 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4744 d(k)=qnwfa(k) - dtz(k)*rhoinv(k)*s_awqnwfa(k+1)*nonloc
4747 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4748 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) + &
4749 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4750 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4751 d(k)=qnwfa(k) + dtz(k)*rhoinv(k)*(s_awqnwfa(k)-s_awqnwfa(k+1))*nonloc
4777IF (bl_mynn_cloudmix > 0 .AND. flag_qnifa .AND. &
4778 bl_mynn_mixscalars > 0)
THEN
4782 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4783 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) - &
4784 & 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4785 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4786 d(k)=qnifa(k) - dtz(k)*rhoinv(k)*s_awqnifa(k+1)*nonloc
4789 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4790 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) + &
4791 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4792 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4793 d(k)=qnifa(k) + dtz(k)*rhoinv(k)*(s_awqnifa(k)-s_awqnifa(k+1))*nonloc
4819IF (bl_mynn_cloudmix > 0 .AND. flag_qnbca .AND. &
4820 bl_mynn_mixscalars > 0)
THEN
4824 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4825 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) - &
4826 & 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4827 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4828 d(k)=qnbca(k) - dtz(k)*rhoinv(k)*s_awqnbca(k+1)*nonloc
4831 a(k)= -dtz(k)*khdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*nonloc
4832 b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))*rhoinv(k) + &
4833 & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*nonloc
4834 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*nonloc
4835 d(k)=qnbca(k) + dtz(k)*rhoinv(k)*(s_awqnbca(k)-s_awqnbca(k+1))*nonloc
4865 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4866 b(k)=1.+dtz(k)*(khdz(k+1)+khdz(k))*rhoinv(k)
4867 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4871 a(k)= -dtz(k)*khdz(k)*rhoinv(k)
4872 b(k)=1.+dtz(k)*(khdz(k)+khdz(k+1))*rhoinv(k)
4873 c(k)= -dtz(k)*khdz(k+1)*rhoinv(k)
4889 dozone(k)=(x(k)-ozone(k))/delt
4897 IF (bl_mynn_mixqt > 0)
THEN
4900 th_new = thl(k) + xlvcp/exner(k)*sqc(k) &
4901 & + xlscp/exner(k)*sqi(k)
4911 IF (sqc(k) > 0.0 .or. sqi(k) > 0.0)
THEN
4912 sqv2(k) = min(sqw2(k),qsat)
4913 portion_qc = sqc(k)/(sqc(k) + sqi(k))
4914 portion_qi = sqi(k)/(sqc(k) + sqi(k))
4915 condensate = max(sqw2(k) - qsat, 0.0)
4916 sqc2(k) = condensate*portion_qc
4917 sqi2(k) = condensate*portion_qi
4931 dqv(k)=(sqv2(k) - sqv(k))/delt
4935 IF (bl_mynn_cloudmix > 0)
THEN
4942 dqc(k)=(sqc2(k) - sqc(k))/delt
4954 IF (flag_qnc .AND. bl_mynn_mixscalars > 0)
THEN
4956 dqnc(k) = (qnc2(k)-qnc(k))/delt
4970 dqi(k)=(sqi2(k) - sqi(k))/delt
4984 dqs(k)=(sqs2(k) - sqs(k))/delt
4995 IF (flag_qni .AND. bl_mynn_mixscalars > 0)
THEN
4997 dqni(k)=(qni2(k)-qni(k))/delt
5017 CALL moisture_check(kte, delt, delp, exner, &
5018 sqv2, sqc2, sqi2, sqs2, thl, &
5019 dqv, dqc, dqi, dqs, dth )
5025 IF(dozone(k)*delt + ozone(k) < 0.)
THEN
5026 dozone(k)=-ozone(k)*0.99/delt
5035 dth(k)=(thl(k) + xlvcp/exner(k)*sqc2(k) &
5036 & + xlscp/exner(k)*(sqi2(k)+sqs(k)) &
5046 dth(k)=(thl(k)+xlvcp/exner(k)*sqc2(k) - th(k))/delt
5057 IF (flag_qnwfa .AND. flag_qnifa .AND. &
5058 bl_mynn_mixscalars > 0)
THEN
5063 dqnwfa(k)=(qnwfa2(k) - qnwfa(k))/delt
5067 dqnifa(k)=(qnifa2(k) - qnifa(k))/delt
5079 IF (flag_qnbca .AND. bl_mynn_mixscalars > 0)
THEN
5081 dqnbca(k)=(qnbca2(k) - qnbca(k))/delt
5096 if (debug_code)
then
5099 wsp = sqrt(u(k)**2 + v(k)**2)
5100 wsp2 = sqrt((u(k)+du(k)*delt)**2 + (v(k)+du(k)*delt)**2)
5101 th2 = th(k) + dth(k)*delt
5103 if (wsp2 > 200. .or. tk2 > 360. .or. tk2 < 160.)
then
5105 print*,
"Outgoing problem at: i=",i,
" k=",k
5106 print*,
" incoming wsp=",wsp,
" outgoing wsp=",wsp2
5107 print*,
" incoming T=",th(k)*exner(k),
"outgoing T:",tk2
5108 print*,
" du=",du(k)*delt,
" dv=",dv(k)*delt,
" dth=",dth(k)*delt
5109 print*,
" km=",kmdz(k)*dz(k),
" kh=",khdz(k)*dz(k)
5110 print*,
" u*=",ust,
" wspd=",wspd,
"rhosfc=",rhosfc
5111 print*,
" LH=",flq*rhosfc*1004.,
" HFX=",flt*rhosfc*1004.
5112 print*,
" drag term=",ust**2/wspd*dtz(k)*rhosfc/rho(kts)
5117 print*,
"==thl:",thl(max(kproblem-3,1):min(kproblem+3,kte))
5118 print*,
"===qv:",sqv2(max(kproblem-3,1):min(kproblem+3,kte))
5119 print*,
"===qc:",sqc2(max(kproblem-3,1):min(kproblem+3,kte))
5120 print*,
"===qi:",sqi2(max(kproblem-3,1):min(kproblem+3,kte))
5121 print*,
"====u:",u(max(kproblem-3,1):min(kproblem+3,kte))
5122 print*,
"====v:",v(max(kproblem-3,1):min(kproblem+3,kte))
5126#ifdef HARDCODE_VERTICAL
5682 & kts,kte,dt,zw,dz,p,rho, &
5686 & u,v,w,th,thl,thv,tk, &
5688 & qnc,qni,qnwfa,qnifa,qnbca, &
5689 & exner,vt,vq,sgm, &
5690 & ust,flt,fltv,flq,flqv, &
5691 & pblh,kpbl,dx,landsea,ts, &
5692 ! outputs - updraft properties
5694 & edmf_qt,edmf_thl, &
5695 & edmf_ent,edmf_qc, &
5696 ! outputs - variables needed for solver
5697 & s_aw,s_awthl,s_awqt, &
5699 & s_awu,s_awv,s_awqke, &
5700 & s_awqnc,s_awqni, &
5701 & s_awqnwfa,s_awqnifa, &
5703 & sub_thl,sub_sqv, &
5705 & det_thl,det_sqv,det_sqc, &
5708 & nchem,chem1,s_awchem, &
5710 ! in/outputs - subgrid scale clouds
5711 & qc_bl1d,cldfra_bl1d, &
5712 & qc_bl1D_old,cldfra_bl1D_old, &
5713 ! inputs - flags for moist arrays
5716 & F_QNWFA,F_QNIFA,F_QNBCA, &
5719 & maxwidth,ktop,maxmf,ztop, &
5720 ! inputs for stochastic perturbations
5721 & spp_pbl,rstoch_col )
5724 integer,
intent(in) :: KTS,KTE,KPBL,momentum_opt,tke_opt,scalar_opt
5726#ifdef HARDCODE_VERTICAL
5728# define kte HARDCODE_VERTICAL
5732 integer,
intent(in) :: spp_pbl
5733 real(kind_phys),
dimension(kts:kte) :: rstoch_col
5735 real(kind_phys),
dimension(kts:kte),
intent(in) :: &
5736 &U,V,W,TH,THL,TK,QT,QV,QC, &
5737 &exner,dz,THV,P,rho,qke,qnc,qni,qnwfa,qnifa,qnbca
5738 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
5739 real(kind_phys),
intent(in) :: flt,fltv,flq,flqv,Psig_shcu, &
5740 &landsea,ts,dx,dt,ust,pblh
5741 logical,
optional :: F_QC,F_QI,F_QNC,F_QNI,F_QNWFA,F_QNIFA,F_QNBCA
5744 real(kind_phys),
dimension(kts:kte),
intent(out) :: edmf_a,edmf_w, &
5745 & edmf_qt,edmf_thl,edmf_ent,edmf_qc
5747 real(kind_phys),
dimension(kts:kte) :: edmf_th
5749 integer,
intent(out) :: ktop
5750 real(kind_phys),
intent(out) :: maxmf,ztop,maxwidth
5752 real(kind_phys),
dimension(kts:kte+1) :: s_aw, &
5753 &s_awthl,s_awqt,s_awqv,s_awqc,s_awqnc,s_awqni, &
5754 &s_awqnwfa,s_awqnifa,s_awqnbca,s_awu,s_awv, &
5757 real(kind_phys),
dimension(kts:kte),
intent(inout) :: &
5758 &qc_bl1d,cldfra_bl1d,qc_bl1d_old,cldfra_bl1d_old
5760 integer,
parameter :: nup=8, debug_mf=0
5761 real(kind_phys) :: nup2
5766 real(kind_phys),
dimension(kts:kte+1,1:NUP) :: &
5767 &UPW,UPTHL,UPQT,UPQC,UPQV, &
5768 &UPA,UPU,UPV,UPTHV,UPQKE,UPQNC, &
5769 &UPQNI,UPQNWFA,UPQNIFA,UPQNBCA
5771 real(kind_phys),
dimension(kts:kte,1:NUP) :: ENT,ENTf
5772 integer,
dimension(kts:kte,1:NUP) :: ENTi
5775 real(kind_phys):: fltv2,wstar,qstar,thstar,sigmaW,sigmaQT, &
5776 &sigmaTH,z0,pwmin,pwmax,wmin,wmax,wlv,Psig_w,maxw,maxqc,wpbl
5777 real(kind_phys):: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,QNCn,QNIn, &
5778 & QNWFAn,QNIFAn,QNBCAn, &
5779 & Wn2,Wn,EntEXP,EntEXM,EntW,BCOEFF,THVkm1,THVk,Pk,rho_int
5782 real(kind_phys),
parameter :: &
5789 real(kind_phys),
parameter :: &
5794 real(kind_phys),
parameter :: Atot = 0.10
5795 real(kind_phys),
parameter :: lmax = 1000.
5796 real(kind_phys),
parameter :: lmin = 300.
5797 real(kind_phys),
parameter :: dlmin = 0.
5798 real(kind_phys) :: minwidth
5799 real(kind_phys) :: dl
5800 real(kind_phys),
parameter :: dcut = 1.2
5804 real(kind_phys):: cn,c,l,n,an2,hux,wspd_pbl,cloud_base,width_flx
5807 integer,
intent(in) :: nchem
5808 real(kind_phys),
dimension(:, :) :: chem1
5809 real(kind_phys),
dimension(kts:kte+1, nchem) :: s_awchem
5810 real(kind_phys),
dimension(nchem) :: chemn
5811 real(kind_phys),
dimension(kts:kte+1,1:NUP, nchem) :: UPCHEM
5813 real(kind_phys),
dimension(kts:kte+1, nchem) :: edmf_chem
5814 logical,
intent(in) :: mix_chem
5817 real(kind_phys):: ERF
5819 logical :: superadiabatic
5822 real(kind_phys),
dimension(kts:kte),
intent(inout) :: vt, vq, sgm
5823 real(kind_phys):: sigq,xl,rsl,cpm,a,qmq,mf_cf,Aup,Q1,diffqt,qsat_tk,&
5824 Fng,qww,alpha,beta,bb,f,pt,t,q2p,b9,satvp,rhgrid, &
5825 Ac_mf,Ac_strat,qc_mf
5826 real(kind_phys),
parameter :: cf_thresh = 0.5
5829 real(kind_phys),
dimension(kts:kte) :: exneri,dzi,rhoz
5830 real(kind_phys):: THp, QTp, QCp, QCs, esat, qsl
5831 real(kind_phys):: csigma,acfac,ac_wsp
5834 integer :: overshoot
5835 real(kind_phys):: bvf, Frz, dzp
5839 real(kind_phys):: adjustment, flx1
5840 real(kind_phys),
parameter :: fluxportion=0.75
5845 real(kind_phys),
dimension(kts:kte) :: sub_thl,sub_sqv,sub_u,sub_v, &
5846 det_thl,det_sqv,det_sqc,det_u,det_v, &
5847 envm_a,envm_w,envm_thl,envm_sqv,envm_sqc, &
5849 real(kind_phys),
dimension(kts:kte+1) :: envi_a,envi_w
5850 real(kind_phys):: temp,sublim,qc_ent,qv_ent,qt_ent,thl_ent,detrate, &
5851 detrateUV,oow,exc_fac,aratio,detturb,qc_grid,qc_sgs, &
5852 qc_plume,exc_heat,exc_moist,tk_int,tvs
5853 real(kind_phys),
parameter :: Cdet = 1./45.
5854 real(kind_phys),
parameter :: dzpmax = 300.
5859 real(kind_phys),
parameter :: Csub=0.25
5862 real(kind_phys),
parameter :: pgfac = 0.00
5863 real(kind_phys):: Uk,Ukm1,Vk,Vkm1,dxsa
5893 if ( mix_chem )
then
5894 upchem(kts:kte+1,1:nup,1:nchem)=0.0
5905 if ( mix_chem )
then
5906 edmf_chem(kts:kte+1,1:nchem) = 0.0
5923 if ( mix_chem )
then
5924 s_awchem(kts:kte+1,1:nchem) = 0.0
5944 if (zw(k) > pblh + 500.)
exit
5947 if (w(k) < 0.)wpbl = 2.*w(k)
5948 maxw = max(maxw,abs(wpbl))
5951 if (zw(k)<=50.)k50=k
5954 qc_sgs = max(qc(k), qc_bl1d(k))
5955 if (qc_sgs> 1e-5 .and. (cldfra_bl1d(k) .ge. 0.5) .and. cloud_base == 9000.0)
then
5956 cloud_base = 0.5*(zw(k)+zw(k+1))
5961 maxw = max(0.,maxw - 1.0)
5962 psig_w = max(0.0, 1.0 - maxw)
5963 psig_w = min(psig_w, psig_shcu)
5967 if(psig_w == 0.0 .and. fltv > 0.0) fltv2 = -1.*fltv
5971 superadiabatic = .false.
5972 if ((landsea-1.5).ge.0)
then
5977 tvs = ts*(1.0+p608*qv(kts))
5980 if ((thv(k)-tvs)/(0.5*dz(k)) < hux)
then
5981 superadiabatic = .true.
5983 superadiabatic = .false.
5987 if ((thv(k)-thv(k-1))/(0.5*(dz(k)+dz(k-1))) < hux)
then
5988 superadiabatic = .true.
5990 superadiabatic = .false.
6005 maxwidth = min(dx*dcut, lmax)
6007 maxwidth = min(maxwidth, 1.1_kind_phys*pblh)
6009 if ((landsea-1.5) .lt. 0)
then
6010 maxwidth = min(maxwidth, 0.5_kind_phys*cloud_base)
6012 maxwidth = min(maxwidth, 0.9_kind_phys*cloud_base)
6015 wspd_pbl=sqrt(max(u(kts)**2 + v(kts)**2, 0.01_kind_phys))
6018 if ((landsea-1.5).LT.0)
then
6019 width_flx = max(min(1000.*(0.6*tanh((fltv - 0.040)/0.04) + .5),1000._kind_phys), 0._kind_phys)
6021 width_flx = max(min(1000.*(0.6*tanh((fltv - 0.007)/0.02) + .5),1000._kind_phys), 0._kind_phys)
6023 maxwidth = min(maxwidth, width_flx)
6027 if (maxwidth .ge. (lmax - 1.0) .and. fltv .gt. 0.2)minwidth = lmin + dlmin*min((fltv-0.2)/0.3, 1._kind_phys)
6029 if (maxwidth .le. minwidth)
then
6040if ( fltv2 > 0.002 .AND. (maxwidth > minwidth) .AND. superadiabatic)
then
6045 dl = (maxwidth - minwidth)/real(nup-1,kind=kind_phys)
6048 l = minwidth + dl*real(i-1)
6049 cn = cn + l**d * (l*l)/(dx*dx) * dl
6054 if ((landsea-1.5).LT.0)
then
6055 acfac = .5*tanh((fltv2 - 0.02)/0.05) + .5
6057 acfac = .5*tanh((fltv2 - 0.01)/0.03) + .5
6063 if (wspd_pbl .le. 10.)
then
6066 ac_wsp = 1.0 - min((wspd_pbl - 10.0)/15., 1.0)
6068 acfac = acfac * ac_wsp
6074 l = minwidth + dl*real(i-1)
6076 upa(1,i) = n*l*l/(dx*dx) * dl
6078 upa(1,i) = upa(1,i)*acfac
6079 an2 = an2 + upa(1,i)
6088 wstar=max(1.e-2,(gtr*fltv2*pblh)**(onethird))
6089 qstar=max(flq,1.0e-5)/wstar
6092 if ((landsea-1.5) .ge. 0)
then
6101 if ((landsea-1.5).GE.0)
then
6110 exc_fac = exc_fac * ac_wsp
6113 sigmaw =csigma*wstar*(z0/pblh)**(onethird)*(1 - 0.8*z0/pblh)
6114 sigmaqt=csigma*qstar*(z0/pblh)**(onethird)
6115 sigmath=csigma*thstar*(z0/pblh)**(onethird)
6119 wmin=min(sigmaw*pwmin,0.1)
6120 wmax=min(sigmaw*pwmax,0.5)
6124 wlv=wmin+(wmax-wmin)/nup2*(i-1)
6127 upw(1,i)=wmin + real(i)/real(nup)*(wmax-wmin)
6128 upu(1,i)=(u(kts)*dz(kts+1)+u(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6129 upv(1,i)=(v(kts)*dz(kts+1)+v(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6133 exc_heat = exc_fac*upw(1,i)*sigmath/sigmaw
6134 upthv(1,i)=(thv(kts)*dz(kts+1)+thv(kts+1)*dz(kts))/(dz(kts)+dz(kts+1)) &
6136 upthl(1,i)=(thl(kts)*dz(kts+1)+thl(kts+1)*dz(kts))/(dz(kts)+dz(kts+1)) &
6140 exc_moist=exc_fac*upw(1,i)*sigmaqt/sigmaw
6141 upqt(1,i)=(qt(kts)*dz(kts+1)+qt(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))&
6144 upqke(1,i)=(qke(kts)*dz(kts+1)+qke(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6145 upqnc(1,i)=(qnc(kts)*dz(kts+1)+qnc(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6146 upqni(1,i)=(qni(kts)*dz(kts+1)+qni(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6147 upqnwfa(1,i)=(qnwfa(kts)*dz(kts+1)+qnwfa(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6148 upqnifa(1,i)=(qnifa(kts)*dz(kts+1)+qnifa(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6149 upqnbca(1,i)=(qnbca(kts)*dz(kts+1)+qnbca(kts+1)*dz(kts))/(dz(kts)+dz(kts+1))
6152 if ( mix_chem )
then
6155 upchem(1,i,ic)=(chem1(kts,ic)*dz(kts+1)+chem1(kts+1,ic)*dz(kts))/(dz(kts)+dz(kts+1))
6161 envm_thl(kts:kte)=thl(kts:kte)
6162 envm_sqv(kts:kte)=qv(kts:kte)
6163 envm_sqc(kts:kte)=qc(kts:kte)
6164 envm_u(kts:kte)=u(kts:kte)
6165 envm_v(kts:kte)=v(kts:kte)
6167 rhoz(k) = (rho(k)*dz(k+1)+rho(k+1)*dz(k))/(dz(k+1)+dz(k))
6169 rhoz(kte) = rho(kte)
6172 dxsa = 1. - min(max((12000.0-dx)/(12000.0-3000.0), 0.), 1.)
6178 l = minwidth + dl*real(i-1)
6182 wmin = 0.3 + l*0.0005
6183 ent(k,i) = 0.33/(min(max(upw(k-1,i),wmin),0.9)*l)
6191 ent(k,i) = max(ent(k,i),0.0003)
6195 IF(zw(k) >= min(pblh+1500., 4000.))
THEN
6196 ent(k,i)=ent(k,i) + (zw(k)-min(pblh+1500.,4000.))*5.0e-6
6200 ent(k,i) = ent(k,i) * (1.0 - rstoch_col(k))
6202 ent(k,i) = min(ent(k,i),0.9/(zw(k+1)-zw(k)))
6205 uk =(u(k)*dz(k+1)+u(k+1)*dz(k))/(dz(k+1)+dz(k))
6206 ukm1=(u(k-1)*dz(k)+u(k)*dz(k-1))/(dz(k-1)+dz(k))
6207 vk =(v(k)*dz(k+1)+v(k+1)*dz(k))/(dz(k+1)+dz(k))
6208 vkm1=(v(k-1)*dz(k)+v(k)*dz(k-1))/(dz(k-1)+dz(k))
6211 entexp= ent(k,i)*(zw(k+1)-zw(k))
6212 entexm= entexp*0.3333
6213 qtn =upqt(k-1,i) *(1.-entexp) + qt(k)*entexp
6214 thln=upthl(k-1,i)*(1.-entexp) + thl(k)*entexp
6215 un =upu(k-1,i) *(1.-entexm) + u(k)*entexm + dxsa*pgfac*(uk - ukm1)
6216 vn =upv(k-1,i) *(1.-entexm) + v(k)*entexm + dxsa*pgfac*(vk - vkm1)
6217 qken=upqke(k-1,i)*(1.-entexp) + qke(k)*entexp
6218 qncn=upqnc(k-1,i)*(1.-entexp) + qnc(k)*entexp
6219 qnin=upqni(k-1,i)*(1.-entexp) + qni(k)*entexp
6220 qnwfan=upqnwfa(k-1,i)*(1.-entexp) + qnwfa(k)*entexp
6221 qnifan=upqnifa(k-1,i)*(1.-entexp) + qnifa(k)*entexp
6222 qnbcan=upqnbca(k-1,i)*(1.-entexp) + qnbca(k)*entexp
6238 if ( mix_chem )
then
6243 chemn(ic)=upchem(k-1,i,ic)*(1.-entexp) + chem1(k,ic)*entexp
6248 pk =(p(k)*dz(k+1)+p(k+1)*dz(k))/(dz(k+1)+dz(k))
6253 thvk =(thv(k)*dz(k+1)+thv(k+1)*dz(k))/(dz(k+1)+dz(k))
6254 thvkm1=(thv(k-1)*dz(k)+thv(k)*dz(k-1))/(dz(k-1)+dz(k))
6257 b=grav*(thvn/thvk - 1.0)
6272 IF (upw(k-1,i) < 0.2 )
THEN
6273 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.)
6275 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.)
6279 IF(wn > upw(k-1,i) + min(1.25*(zw(k)-zw(k-1))/200., 2.0) )
THEN
6280 wn = upw(k-1,i) + min(1.25*(zw(k)-zw(k-1))/200., 2.0)
6283 IF(wn < upw(k-1,i) - min(1.25*(zw(k)-zw(k-1))/200., 2.0) )
THEN
6284 wn = upw(k-1,i) - min(1.25*(zw(k)-zw(k-1))/200., 2.0)
6286 wn = min(max(wn,0.0), 3.0)
6290 IF (k==kts+1 .AND. wn == 0.)
THEN
6295 IF (debug_mf == 1)
THEN
6296 IF (wn .GE. 3.0)
THEN
6298 print *,
" **** SUSPICIOUSLY LARGE W:"
6299 print *,
' QCn:',qcn,
' ENT=',ent(k,i),
' Nup2=',nup2
6300 print *,
'pblh:',pblh,
' Wn:',wn,
' UPW(k-1)=',upw(k-1,i)
6301 print *,
'K=',k,
' B=',b,
' dz=',zw(k)-zw(k-1)
6306 IF (wn <= 0.0 .AND. overshoot == 0)
THEN
6308 IF ( thvk-thvkm1 .GT. 0.0 )
THEN
6309 bvf = sqrt( gtr*(thvk-thvkm1)/dz(k) )
6311 frz = upw(k-1,i)/(bvf*dz(k))
6313 dzp = dz(k)*max(min(frz,1.0),0.0)
6327 aratio = min(upa(k-1,i)/(1.-upa(k-1,i)), 0.5)
6329 oow = -0.060/max(1.0,(0.5*(wn+upw(k-1,i))))
6330 detrate = min(max(oow*(wn-upw(k-1,i))/dz(k), detturb), .0002)
6331 detrateuv= min(max(oow*(wn-upw(k-1,i))/dz(k), detturb), .0001)
6332 envm_thl(k)=envm_thl(k) + (0.5*(thl_ent + upthl(k-1,i)) - thl(k))*detrate*aratio*min(dzp,dzpmax)
6333 qv_ent = 0.5*(max(qt_ent-qc_ent,0.) + max(upqt(k-1,i)-upqc(k-1,i),0.))
6334 envm_sqv(k)=envm_sqv(k) + (qv_ent-qv(k))*detrate*aratio*min(dzp,dzpmax)
6335 IF (upqc(k-1,i) > 1e-8)
THEN
6336 IF (qc(k) > 1e-6)
THEN
6339 qc_grid = cldfra_bl1d(k)*qc_bl1d(k)
6341 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)
6343 envm_u(k) =envm_u(k) + (0.5*(un + upu(k-1,i)) - u(k))*detrateuv*aratio*min(dzp,dzpmax)
6344 envm_v(k) =envm_v(k) + (0.5*(vn + upv(k-1,i)) - v(k))*detrateuv*aratio*min(dzp,dzpmax)
6362 IF ( mix_chem )
THEN
6364 upchem(k,i,ic) = chemn(ic)
6373 IF (debug_mf == 1)
THEN
6374 IF (maxval(upw(:,i)) > 10.0 .OR. minval(upa(:,i)) < 0.0 .OR. &
6375 maxval(upa(:,i)) > atot .OR. nup2 > 10)
THEN
6377 print *,
'flq:',flq,
' fltv:',fltv2,
' Nup2=',nup2
6378 print *,
'pblh:',pblh,
' wstar:',wstar,
' ktop=',ktop
6379 print *,
'sigmaW=',sigmaw,
' sigmaTH=',sigmath,
' sigmaQT=',sigmaqt
6384 print *,
'UPA:',upa(:,i)
6385 print *,
'UPW:',upw(:,i)
6386 print *,
'UPTHL:',upthl(:,i)
6387 print *,
'UPQT:',upqt(:,i)
6388 print *,
'ENT:',ent(:,i)
6409 s_aw(k+1) = s_aw(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*psig_w
6410 s_awthl(k+1)= s_awthl(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upthl(k,i)*psig_w
6411 s_awqt(k+1) = s_awqt(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqt(k,i)*psig_w
6416 qc_plume = upqc(k,i)
6420 s_awqc(k+1) = s_awqc(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*qc_plume*psig_w
6421 s_awqv(k+1) = s_awqt(k+1) - s_awqc(k+1)
6425 if (momentum_opt > 0)
then
6428 s_awu(k+1) = s_awu(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upu(k,i)*psig_w
6429 s_awv(k+1) = s_awv(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upv(k,i)*psig_w
6434 if (tke_opt > 0)
then
6437 s_awqke(k+1)= s_awqke(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqke(k,i)*psig_w
6442 if ( mix_chem )
then
6446 s_awchem(k+1,ic) = s_awchem(k+1,ic) + rhoz(k)*upa(k,i)*upw(k,i)*upchem(k,i,ic)*psig_w
6452 if (scalar_opt > 0)
then
6455 s_awqnc(k+1) = s_awqnc(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqnc(k,i)*psig_w
6456 s_awqni(k+1) = s_awqni(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqni(k,i)*psig_w
6457 s_awqnwfa(k+1)= s_awqnwfa(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqnwfa(k,i)*psig_w
6458 s_awqnifa(k+1)= s_awqnifa(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqnifa(k,i)*psig_w
6459 s_awqnbca(k+1)= s_awqnbca(k+1) + rhoz(k)*upa(k,i)*upw(k,i)*upqnbca(k,i)*psig_w
6467 IF (s_aw(kts+1) /= 0.)
THEN
6468 dzi(kts) = 0.5*(dz(kts)+dz(kts+1))
6469 flx1 = max(s_aw(kts+1)*(th(kts)-th(kts+1))/dzi(kts),1.0e-5)
6478 IF (flx1 > fluxportion*flt/dz(kts) .AND. flx1>0.0)
THEN
6479 adjustment= fluxportion*flt/dz(kts)/flx1
6480 s_aw = s_aw*adjustment
6481 s_awthl = s_awthl*adjustment
6482 s_awqt = s_awqt*adjustment
6483 s_awqc = s_awqc*adjustment
6484 s_awqv = s_awqv*adjustment
6485 s_awqnc = s_awqnc*adjustment
6486 s_awqni = s_awqni*adjustment
6487 s_awqnwfa = s_awqnwfa*adjustment
6488 s_awqnifa = s_awqnifa*adjustment
6489 s_awqnbca = s_awqnbca*adjustment
6490 IF (momentum_opt > 0)
THEN
6491 s_awu = s_awu*adjustment
6492 s_awv = s_awv*adjustment
6494 IF (tke_opt > 0)
THEN
6495 s_awqke= s_awqke*adjustment
6497 IF ( mix_chem )
THEN
6498 s_awchem = s_awchem*adjustment
6500 upa = upa*adjustment
6508 edmf_a(k) =edmf_a(k) +upa(k,i)
6509 edmf_w(k) =edmf_w(k) +rhoz(k)*upa(k,i)*upw(k,i)
6510 edmf_qt(k) =edmf_qt(k) +rhoz(k)*upa(k,i)*upqt(k,i)
6511 edmf_thl(k)=edmf_thl(k)+rhoz(k)*upa(k,i)*upthl(k,i)
6512 edmf_ent(k)=edmf_ent(k)+rhoz(k)*upa(k,i)*ent(k,i)
6513 edmf_qc(k) =edmf_qc(k) +rhoz(k)*upa(k,i)*upqc(k,i)
6519 if (edmf_a(k)>0.)
then
6520 edmf_w(k)=edmf_w(k)/edmf_a(k)
6521 edmf_qt(k)=edmf_qt(k)/edmf_a(k)
6522 edmf_thl(k)=edmf_thl(k)/edmf_a(k)
6523 edmf_ent(k)=edmf_ent(k)/edmf_a(k)
6524 edmf_qc(k)=edmf_qc(k)/edmf_a(k)
6525 edmf_a(k)=edmf_a(k)*psig_w
6527 if(edmf_a(k)*edmf_w(k) > maxmf) maxmf = edmf_a(k)*edmf_w(k)
6532 if ( mix_chem )
then
6536 edmf_chem(k,ic) = edmf_chem(k,ic) + rhoz(k)*upa(k,i)*upchem(k,i,ic)
6541 if (edmf_a(k)>0.)
then
6543 edmf_chem(k,ic) = edmf_chem(k,ic)/edmf_a(k)
6557 envi_w(k) = onethird*(edmf_w(k-1)+edmf_w(k)+edmf_w(k+1))
6558 envi_a(k) = onethird*(edmf_a(k-1)+edmf_a(k)+edmf_a(k+1))*adjustment
6561 envi_w(kts) = edmf_w(kts)
6562 envi_a(kts) = edmf_a(kts)
6565 envi_a(kte) = edmf_a(kte)
6568 envi_a(kte+1) = edmf_a(kte)
6572 IF (envi_w(kts) > 0.9*dz(kts)/dt)
THEN
6573 sublim = 0.9*dz(kts)/dt/envi_w(kts)
6581 envi_w(k)=csub*sublim*envi_w(k)*temp/(1.-temp)
6585 dzi(kts) = 0.5*(dz(kts)+dz(kts+1))
6586 sub_thl(kts)= 0.5*envi_w(kts)*envi_a(kts)* &
6587 (rho(kts+1)*thl(kts+1)-rho(kts)*thl(kts))/dzi(kts)/rhoz(k)
6588 sub_sqv(kts)= 0.5*envi_w(kts)*envi_a(kts)* &
6589 (rho(kts+1)*qv(kts+1)-rho(kts)*qv(kts))/dzi(kts)/rhoz(k)
6591 dzi(k) = 0.5*(dz(k)+dz(k+1))
6592 sub_thl(k)= 0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * &
6593 (rho(k+1)*thl(k+1)-rho(k)*thl(k))/dzi(k)/rhoz(k)
6594 sub_sqv(k)= 0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * &
6595 (rho(k+1)*qv(k+1)-rho(k)*qv(k))/dzi(k)/rhoz(k)
6599 det_thl(k)=cdet*(envm_thl(k)-thl(k))*envi_a(k)*psig_w
6600 det_sqv(k)=cdet*(envm_sqv(k)-qv(k))*envi_a(k)*psig_w
6601 det_sqc(k)=cdet*(envm_sqc(k)-qc(k))*envi_a(k)*psig_w
6604 IF (momentum_opt > 0)
THEN
6605 sub_u(kts)=0.5*envi_w(kts)*envi_a(kts)* &
6606 (rho(kts+1)*u(kts+1)-rho(kts)*u(kts))/dzi(kts)/rhoz(k)
6607 sub_v(kts)=0.5*envi_w(kts)*envi_a(kts)* &
6608 (rho(kts+1)*v(kts+1)-rho(kts)*v(kts))/dzi(kts)/rhoz(k)
6610 sub_u(k)=0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * &
6611 (rho(k+1)*u(k+1)-rho(k)*u(k))/dzi(k)/rhoz(k)
6612 sub_v(k)=0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * &
6613 (rho(k+1)*v(k+1)-rho(k)*v(k))/dzi(k)/rhoz(k)
6617 det_u(k) = cdet*(envm_u(k)-u(k))*envi_a(k)*psig_w
6618 det_v(k) = cdet*(envm_v(k)-v(k))*envi_a(k)*psig_w
6627 exneri(k) = (exner(k)*dz(k+1)+exner(k+1)*dz(k))/(dz(k+1)+dz(k))
6628 edmf_th(k)= edmf_thl(k) + xlvcp/exneri(k)*edmf_qc(k)
6629 dzi(k) = 0.5*(dz(k)+dz(k+1))
6637 if(0.5*(edmf_qc(k)+edmf_qc(k-1))>0.0 .and. (cldfra_bl1d(k) < cf_thresh))
THEN
6639 aup = (edmf_a(k)*dzi(k-1)+edmf_a(k-1)*dzi(k))/(dzi(k-1)+dzi(k))
6640 thp = (edmf_th(k)*dzi(k-1)+edmf_th(k-1)*dzi(k))/(dzi(k-1)+dzi(k))
6641 qtp = (edmf_qt(k)*dzi(k-1)+edmf_qt(k-1)*dzi(k))/(dzi(k-1)+dzi(k))
6647 qsl=ep_2*esat/max(1.e-7,(p(k)-ep_3*esat))
6650 if (edmf_qc(k)>0.0 .and. edmf_qc(k-1)>0.0)
then
6651 qcp = (edmf_qc(k)*dzi(k-1)+edmf_qc(k-1)*dzi(k))/(dzi(k-1)+dzi(k))
6653 qcp = max(edmf_qc(k),edmf_qc(k-1))
6660 rsl = xl*qsat_tk / (r_v*tk(k)**2)
6662 cpm = cp + qt(k)*cpv
6663 a = 1./(1. + xl*rsl/cpm)
6666 q2p = xlvcp/exner(k)
6667 pt = thl(k) +q2p*qcp*aup
6676 beta = pt*xl/(tk(k)*cp) - 1.61*pt
6691 sigq = 10. * aup * (qtp - qt(k))
6693 sigq = max(sigq, qsat_tk*0.02 )
6694 sigq = min(sigq, qsat_tk*0.25 )
6696 qmq = a * (qt(k) - qsat_tk)
6699 if ((landsea-1.5).GE.0)
then
6703 mf_cf = min(max(0.5 + 0.36 * atan(1.55*q1),0.01),0.6)
6704 mf_cf = max(mf_cf, 1.2 * aup)
6705 mf_cf = min(mf_cf, 5.0 * aup)
6710 mf_cf = min(max(0.5 + 0.36 * atan(1.55*q1),0.01),0.6)
6711 mf_cf = max(mf_cf, 1.8 * aup)
6712 mf_cf = min(mf_cf, 5.0 * aup)
6726 if ((landsea-1.5).GE.0)
then
6727 if (qcp * aup > 5e-5)
then
6728 qc_bl1d(k) = 1.86 * (qcp * aup) - 2.2e-5
6730 qc_bl1d(k) = 1.18 * (qcp * aup)
6732 cldfra_bl1d(k) = mf_cf
6735 if (qcp * aup > 5e-5)
then
6736 qc_bl1d(k) = 1.86 * (qcp * aup) - 2.2e-5
6738 qc_bl1d(k) = 1.18 * (qcp * aup)
6740 cldfra_bl1d(k) = mf_cf
6755 if (q1 .ge. 1.0)
then
6757 elseif (q1 .ge. -1.7 .and. q1 .lt. 1.0)
then
6758 fng = exp(-0.4*(q1-1.0))
6759 elseif (q1 .ge. -2.5 .and. q1 .lt. -1.7)
then
6760 fng = 3.0 + exp(-3.8*(q1+1.7))
6762 fng = min(23.9 + exp(-1.6*(q1+2.5)), 60.)
6766 vt(k) = qww - (1.5*aup)*beta*bb*fng - 1.
6767 vq(k) = alpha + (1.5*aup)*beta*a*fng - tv0
6775 maxqc = maxval(edmf_qc(1:ktop))
6776 if ( maxqc < 1.e-8) maxmf = -1.0*maxmf
6782if (edmf_w(1) > 4.0)
then
6784 print *,
'flq:',flq,
' fltv:',fltv2
6785 print *,
'pblh:',pblh,
' wstar:',wstar
6786 print *,
'sigmaW=',sigmaw,
' sigmaTH=',sigmath,
' sigmaQT=',sigmaqt
6814 print *,
' edmf_a',edmf_a(1:14)
6815 print *,
' edmf_w',edmf_w(1:14)
6816 print *,
' edmf_qt:',edmf_qt(1:14)
6817 print *,
' edmf_thl:',edmf_thl(1:14)
6822#ifdef HARDCODE_VERTICAL
6934SUBROUTINE ddmf_jpl(kts,kte,dt,zw,dz,p, &
6935 &u,v,th,thl,thv,tk,qt,qv,qc, &
6937 &ust,wthl,wqt,pblh,kpbl, &
6938 &edmf_a_dd,edmf_w_dd, edmf_qt_dd, &
6939 &edmf_thl_dd,edmf_ent_dd,edmf_qc_dd, &
6940 &sd_aw,sd_awthl,sd_awqt, &
6941 &sd_awqv,sd_awqc,sd_awu,sd_awv, &
6943 &qc_bl1d,cldfra_bl1d, &
6946 integer,
intent(in) :: KTS,KTE,KPBL
6947 real(kind_phys),
dimension(kts:kte),
intent(in) :: &
6948 U,V,TH,THL,TK,QT,QV,QC,THV,P,rho,exner,dz
6949 real(kind_phys),
dimension(kts:kte),
intent(in) :: rthraten
6951 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: ZW
6952 real(kind_phys),
intent(in) :: WTHL,WQT
6953 real(kind_phys),
intent(in) :: dt,ust,pblh
6955 real(kind_phys),
dimension(kts:kte),
intent(out) :: &
6956 edmf_a_dd,edmf_w_dd, &
6957 edmf_qt_dd,edmf_thl_dd, edmf_ent_dd,edmf_qc_dd
6960 real(kind_phys),
dimension(kts:kte+1) :: &
6961 sd_aw, sd_awthl, sd_awqt, sd_awu, &
6962 sd_awv, sd_awqc, sd_awqv, sd_awqke, sd_aw2
6964 real(kind_phys),
dimension(kts:kte),
intent(in) :: &
6965 qc_bl1d, cldfra_bl1d
6967 integer,
parameter:: ndown = 5
6969 integer,
dimension(1:NDOWN) :: DD_initK
6970 real(kind_phys),
dimension(1:NDOWN) :: randNum
6972 real(kind_phys),
dimension(kts:kte+1,1:NDOWN) :: &
6973 DOWNW,DOWNTHL,DOWNQT,DOWNQC,DOWNA,DOWNU,DOWNV,DOWNTHV
6976 real(kind_phys),
dimension(KTS+1:KTE+1,1:NDOWN) :: ENT,ENTf
6977 integer,
dimension(KTS+1:KTE+1,1:NDOWN) :: ENTi
6980 integer :: K,I,ki, kminrad, qlTop, p700_ind, qlBase
6981 real(kind_phys):: wthv,wstar,qstar,thstar,sigmaW,sigmaQT, &
6982 sigmaTH,z0,pwmin,pwmax,wmin,wmax,wlv,wtv,went,mindownw
6983 real(kind_phys):: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,Wn2,Wn, &
6984 THVk,Pk,EntEXP,EntW,beta_dm,EntExp_M,rho_int
6985 real(kind_phys):: jump_thetav, jump_qt, jump_thetal, &
6986 refTHL, refTHV, refQT
6988 real(kind_phys):: minrad,zminrad, radflux, F0, wst_rad, wst_dd
6990 real(kind_phys):: sigq,xl,rsl,cpm,a,mf_cf,diffqt, &
6991 Fng,qww,alpha,beta,bb,f,pt,t,q2p,b9,satvp,rhgrid
6994 real(kind_phys),
parameter :: &
6995 &Wa=1., Wb=1.5, Z00=100., BCOEFF=0.2
6997 real(kind_phys),
parameter :: &
7005 integer,
parameter :: debug_mf=0
7007 dl = (1000.-500.)/real(ndown)
7047 do k = max(3,kpbl-2),kpbl+3
7048 if (qc(k).gt. 1.e-6 .AND. cldfra_bl1d(k).gt.0.5)
then
7054 do k = qltop, kts, -1
7055 if (qc(k) .gt. 1e-6)
then
7059 qlbase = (qltop+qlbase)/2
7072 radflux = rthraten(k) * exner(k)
7073 radflux = radflux * cp / grav * ( p(k) - p(k+1) )
7074 if ( radflux < 0.0 ) f0 = abs(radflux) + f0
7083 adn = min( 0.05 + f0*0.001, 0.3)
7109 p700_ind = minloc(abs(p-70000),1)
7110 jump_thetav = thv(p700_ind) - thv(1) - (thv(p700_ind)-thv(qltop+3))/(zw(p700_ind)-zw(qltop+3))*(zw(p700_ind)-zw(qltop))
7111 jump_qt = qc(p700_ind) + qv(p700_ind) - qc(1) - qv(1)
7112 jump_thetal = thl(p700_ind) - thl(1) - (thl(p700_ind)-thl(qltop+3))/(zw(p700_ind)-zw(qltop+3))*(zw(p700_ind)-zw(qltop))
7119 wst_rad = ( grav * zw(qltop) * f0 / (refthl * rho(qltop) * cp) ) ** (0.333)
7120 wst_rad = max(wst_rad, 0.1)
7121 wstar = max(0.,(grav/thv(1)*wthv*pblh)**(onethird))
7122 went = thv(1) / ( grav * jump_thetav * zw(qltop) ) * &
7123 (0.15 * (wstar**3 + 5*ust**3) + 0.35 * wst_rad**3 )
7124 qstar = abs(went*jump_qt/wst_rad)
7125 thstar = f0/rho(qltop)/cp/wst_rad - went*jump_thetav/wst_rad
7127 wst_dd = (0.15 * (wstar**3 + 5*ust**3) + 0.35 * wst_rad**3 ) ** (0.333)
7129 print*,
"qstar=",qstar,
" thstar=",thstar,
" wst_dd=",wst_dd
7130 print*,
"F0=",f0,
" wst_rad=",wst_rad,
" jump_thv=",jump_thetav
7131 print*,
"entrainment velocity=",went
7134 sigmaqt = 40 * qstar
7135 sigmath = 1.0 * thstar
7144 wlv=wmin+(wmax-wmin)/real(ndown)*(i-1)
7145 wtv=wmin+(wmax-wmin)/real(ndown)*i
7152 downa(ki,i)=adn/real(ndown)
7153 downu(ki,i)=(u(ki-1)*dz(ki) + u(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7154 downv(ki,i)=(v(ki-1)*dz(ki) + v(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7161 refthl = (thl(ki-1)*dz(ki) + thl(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7162 refthv = (thv(ki-1)*dz(ki) + thv(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7163 refqt = (qt(ki-1)*dz(ki) + qt(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7166 downqc(ki,i) = (qc(ki-1)*dz(ki) + qc(ki)*dz(ki-1)) /(dz(ki)+dz(ki-1))
7167 downqt(ki,i) = refqt
7168 downthv(ki,i)= refthv + 0.01 *downw(ki,i)*sigmath/sigmaw
7169 downthl(ki,i)= refthl + 0.01 *downw(ki,i)*sigmath/sigmaw
7181 dp = 500. + dl*real(i)
7183 DO k=dd_initk(i)-1,kts+1,-1
7186 wmin = 0.3 + dp*0.0005
7187 ent(k,i) = 0.33/(min(max(-1.0*downw(k+1,i),wmin),0.9)*dp)
7192 entexp =ent(k,i)*dz(k)
7193 entexp_m=ent(k,i)*0.333*dz(k)
7195 qtn =downqt(k+1,i) *(1.-entexp) + qt(k)*entexp
7196 thln=downthl(k+1,i)*(1.-entexp) + thl(k)*entexp
7197 un =downu(k+1,i) *(1.-entexp) + u(k)*entexp_m
7198 vn =downv(k+1,i) *(1.-entexp) + v(k)*entexp_m
7207 pk =(p(k-1)*dz(k)+p(k)*dz(k-1))/(dz(k)+dz(k-1))
7210 thvk =(thv(k-1)*dz(k)+thv(k)*dz(k-1))/(dz(k)+dz(k-1))
7211 b=grav*(thvn/thvk - 1.0)
7222 mindownw = min(downw(k+1,i),-0.2)
7223 wn = downw(k+1,i) + (-2.*ent(k,i)*downw(k+1,i) - &
7224 bcoeff*b/mindownw)*min(dz(k), 250.)
7228 IF (wn < downw(k+1,i) - min(1.25*dz(k)/200., -2.0))
THEN
7229 wn = downw(k+1,i) - min(1.25*dz(k)/200., -2.0)
7232 IF (wn > downw(k+1,i) + min(1.25*dz(k)/200., 2.0))
THEN
7233 wn = downw(k+1,i) + min(1.25*dz(k)/200., 2.0)
7235 wn = max(min(wn,0.0), -3.0)
7244 IF (wn .lt. 0.)
THEN
7252 downa(k,i) = downa(k+1,i)
7255 if (dd_initk(i) - k .lt. 2)
then
7277 edmf_a_dd(k) =edmf_a_dd(k) +downa(k-1,i)
7278 edmf_w_dd(k) =edmf_w_dd(k) +downa(k-1,i)*downw(k-1,i)
7279 edmf_qt_dd(k) =edmf_qt_dd(k) +downa(k-1,i)*downqt(k-1,i)
7280 edmf_thl_dd(k)=edmf_thl_dd(k)+downa(k-1,i)*downthl(k-1,i)
7281 edmf_ent_dd(k)=edmf_ent_dd(k)+downa(k-1,i)*ent(k-1,i)
7282 edmf_qc_dd(k) =edmf_qc_dd(k) +downa(k-1,i)*downqc(k-1,i)
7285 IF (edmf_a_dd(k) >0.)
THEN
7286 edmf_w_dd(k) =edmf_w_dd(k) /edmf_a_dd(k)
7287 edmf_qt_dd(k) =edmf_qt_dd(k) /edmf_a_dd(k)
7288 edmf_thl_dd(k)=edmf_thl_dd(k)/edmf_a_dd(k)
7289 edmf_ent_dd(k)=edmf_ent_dd(k)/edmf_a_dd(k)
7290 edmf_qc_dd(k) =edmf_qc_dd(k) /edmf_a_dd(k)
7299 rho_int = (rho(k)*dz(k+1)+rho(k+1)*dz(k))/(dz(k+1)+dz(k))
7301 sd_aw(k) =sd_aw(k) +rho_int*downa(k,i)*downw(k,i)
7302 sd_awthl(k)=sd_awthl(k)+rho_int*downa(k,i)*downw(k,i)*downthl(k,i)
7303 sd_awqt(k) =sd_awqt(k) +rho_int*downa(k,i)*downw(k,i)*downqt(k,i)
7304 sd_awqc(k) =sd_awqc(k) +rho_int*downa(k,i)*downw(k,i)*downqc(k,i)
7305 sd_awu(k) =sd_awu(k) +rho_int*downa(k,i)*downw(k,i)*downu(k,i)
7306 sd_awv(k) =sd_awv(k) +rho_int*downa(k,i)*downw(k,i)*downv(k,i)
7308 sd_awqv(k) = sd_awqt(k) - sd_awqc(k)
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...