361 &initflag,restart,cycling, &
363 &u,v,w,th,sqv3d,sqc3d,sqi3d, &
365 &qnwfa,qnifa,qnbca,ozone, &
368 &ust,ch,hfx,qfx,rmol,wspd, &
369 &uoce,voce, & !ocean current
372 &nchem,kdvel,ndvel, & !Smoke/Chem variables
373 &chem3d,vdep,smoke_dbg, &
374 &frp,emis_ant_no, & ! JLS/RAR to adjust exchange coeffs
375 &mix_chem,enh_mix,rrfs_sd, & ! end smoke/chem variables
377 &rublten,rvblten,rthblten, &
378 &rqvblten,rqcblten,rqiblten, &
379 &rqncblten,rqniblten,rqsblten, &
380 &rqnwfablten,rqnifablten, &
381 &rqnbcablten,dozone, &
385 &dqke,qwt,qshear,qbuoy,qdiss, &
386 &qc_bl,qi_bl,cldfra_bl, &
387 &bl_mynn_tkeadvect, &
390 &bl_mynn_mixlength, &
396 &bl_mynn_mixscalars, &
398 &bl_mynn_cloudmix,bl_mynn_mixqt, &
399 &edmf_a,edmf_w,edmf_qt, &
400 &edmf_thl,edmf_ent,edmf_qc, &
401 &sub_thl3D,sub_sqv3D, &
402 &det_thl3D,det_sqv3D, &
403 &maxwidth,maxMF,ztop_plume, &
405 &spp_pbl,pattern_spp_pbl, &
407 &FLAG_QC,FLAG_QI,FLAG_QNC, &
409 &FLAG_QNWFA,FLAG_QNIFA, &
410 &FLAG_QNBCA,FLAG_OZONE, &
411 &IDS,IDE,JDS,JDE,KDS,KDE, &
412 &IMS,IME,JMS,JME,KMS,KME, &
413 &ITS,ITE,JTS,JTE,KTS,KTE )
417 integer,
intent(in) :: initflag
419 logical,
intent(in) :: restart,cycling
420 integer,
intent(in) :: tke_budget
421 integer,
intent(in) :: bl_mynn_cloudpdf
422 integer,
intent(in) :: bl_mynn_mixlength
423 integer,
intent(in) :: bl_mynn_edmf
424 logical,
intent(in) :: bl_mynn_tkeadvect
425 integer,
intent(in) :: bl_mynn_edmf_mom
426 integer,
intent(in) :: bl_mynn_edmf_tke
427 integer,
intent(in) :: bl_mynn_mixscalars
428 integer,
intent(in) :: bl_mynn_output
429 integer,
intent(in) :: bl_mynn_cloudmix
430 integer,
intent(in) :: bl_mynn_mixqt
431 integer,
intent(in) :: icloud_bl
432 real(kind_phys),
intent(in) :: closure
434 logical,
intent(in) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC,&
435 FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA, &
438 logical,
intent(in) :: mix_chem,enh_mix,rrfs_sd,smoke_dbg
440 integer,
intent(in) :: &
441 & IDS,IDE,JDS,JDE,KDS,KDE &
442 &,IMS,IME,JMS,JME,KMS,KME &
443 &,ITS,ITE,JTS,JTE,KTS,KTE
445#ifdef HARDCODE_VERTICAL
447# define kte HARDCODE_VERTICAL
456 real(kind_phys),
intent(in) :: delt
457 real(kind_phys),
dimension(:),
intent(in) :: dx
458 real(kind_phys),
dimension(:,:),
intent(in) :: dz, &
459 &u,v,w,th,sqv3D,p,exner,rho,T3D
460 real(kind_phys),
dimension(:,:),
intent(in) :: &
461 &sqc3D,sqi3D,sqs3D,qni,qnc,qnwfa,qnifa,qnbca
462 real(kind_phys),
dimension(:,:),
intent(in):: ozone
463 real(kind_phys),
dimension(:),
intent(in):: ust, &
465 real(kind_phys),
dimension(:,:),
intent(inout),
optional :: &
467 real(kind_phys),
dimension(:,:),
intent(inout) :: &
469 real(kind_phys),
dimension(:,:),
intent(inout) :: &
470 &rublten,rvblten,rthblten,rqvblten,rqcblten, &
471 &rqiblten,rqsblten,rqniblten,rqncblten, &
472 &rqnwfablten,rqnifablten,rqnbcablten
473 real(kind_phys),
dimension(:,:),
intent(inout) :: dozone
474 real(kind_phys),
dimension(:,:),
intent(in) :: rthraten
476 real(kind_phys),
dimension(:,:),
intent(out),
optional :: exch_h,exch_m
477 real(kind_phys),
dimension(:),
intent(in) :: xland, &
478 &ts,znt,hfx,qfx,uoce,voce
481 real(kind_phys),
dimension(:,:),
intent(inout),
optional :: &
482 & edmf_a,edmf_w,edmf_qt,edmf_thl,edmf_ent,edmf_qc, &
483 & sub_thl3D,sub_sqv3D,det_thl3D,det_sqv3D
488 real(kind_phys),
dimension(:),
intent(inout) :: Pblh
489 real(kind_phys),
dimension(:),
intent(inout) :: rmol
491 real(kind_phys),
dimension(ims:ime) :: psig_bl,psig_shcu
493 integer,
dimension(:),
intent(INOUT) :: &
495 integer,
dimension(:),
intent(INOUT),
optional :: &
498 real(kind_phys),
dimension(:),
intent(out),
optional :: &
499 &maxmf,maxwidth,ztop_plume
501 real(kind_phys),
dimension(:,:),
intent(inout),
optional :: el_pbl
503 real(kind_phys),
dimension(:,:),
intent(inout),
optional :: &
504 &qWT,qSHEAR,qBUOY,qDISS,dqke
507 real(kind_phys),
dimension(kts:kte) :: &
508 &qwt1,qshear1,qbuoy1,qdiss1,dqke1,diss_heat
510 real(kind_phys),
dimension(:,:),
intent(out),
optional :: Sh3D,Sm3D
512 real(kind_phys),
dimension(:,:),
intent(inout),
optional :: &
513 &qc_bl,qi_bl,cldfra_bl
514 real(kind_phys),
dimension(kts:kte) :: qc_bl1D,qi_bl1D, &
515 &cldfra_bl1D,qc_bl1D_old,qi_bl1D_old,cldfra_bl1D_old
518 integer,
intent(IN ) :: nchem, kdvel, ndvel
519 real(kind_phys),
dimension(:,:,:),
intent(INOUT),
optional :: chem3d
520 real(kind_phys),
dimension(:,:),
intent(IN),
optional :: vdep
521 real(kind_phys),
dimension(:),
intent(IN),
optional :: frp
522 real(kind_phys),
dimension(:),
intent(IN) :: EMIS_ANT_NO
524 real(kind_phys),
dimension(kts:kte ,nchem) :: chem1
525 real(kind_phys),
dimension(kts:kte+1,nchem) :: s_awchem1
526 real(kind_phys),
dimension(ndvel) :: vd1
530 integer :: ITF,JTF,KTF, IMD,JMD
531 integer :: i,j,k,kproblem
532 real(kind_phys),
dimension(kts:kte) :: &
533 &thl,tl,qv1,qc1,qi1,qs1,sqw, &
534 &el, dfm, dfh, dfq, tcd, qcd, pdk, pdt, pdq, pdc, &
536 real(kind_phys),
dimension(kts:kte) :: &
537 &thetav,sh,sm,u1,v1,w1,p1, &
538 &ex1,dz1,th1,tk1,rho1,qke1,tsq1,qsq1,cov1, &
540 &du1,dv1,dth1,dqv1,dqc1,dqi1,dqs1,ozone1, &
541 &k_m1,k_h1,qni1,dqni1,qnc1,dqnc1,qnwfa1,qnifa1, &
542 &qnbca1,dqnwfa1,dqnifa1,dqnbca1,dozone1
545 real(kind_phys),
dimension(kts:kte) :: &
546 &dth1mf,dqv1mf,dqc1mf,du1mf,dv1mf
547 real(kind_phys),
dimension(kts:kte) :: &
548 &edmf_a1,edmf_w1,edmf_qt1,edmf_thl1, &
550 real(kind_phys),
dimension(kts:kte) :: &
551 &edmf_a_dd1,edmf_w_dd1,edmf_qt_dd1,edmf_thl_dd1, &
552 &edmf_ent_dd1,edmf_qc_dd1
553 real(kind_phys),
dimension(kts:kte) :: &
554 &sub_thl,sub_sqv,sub_u,sub_v, &
555 &det_thl,det_sqv,det_sqc,det_u,det_v
556 real(kind_phys),
dimension(kts:kte+1) :: &
557 &s_aw1,s_awthl1,s_awqt1, &
558 &s_awqv1,s_awqc1,s_awu1,s_awv1,s_awqke1, &
559 &s_awqnc1,s_awqni1,s_awqnwfa1,s_awqnifa1, &
561 real(kind_phys),
dimension(kts:kte+1) :: &
562 &sd_aw1,sd_awthl1,sd_awqt1, &
563 &sd_awqv1,sd_awqc1,sd_awu1,sd_awv1,sd_awqke1
565 real(kind_phys),
dimension(kts:kte+1) :: zw
566 real(kind_phys) :: cpm,sqcg,flt,fltv,flq,flqv,flqc, &
567 &pmz,phh,exnerg,zet,phi_m, &
568 &afk,abk,ts_decay, qc_bl2, qi_bl2, &
572 real(kind_phys),
dimension(ITS:ITE) :: maxKHtopdown
573 real(kind_phys),
dimension(kts:kte) :: KHtopdown,TKEprodTD
575 logical :: INITIALIZE_QKE,problem
578 integer,
intent(IN) :: spp_pbl
579 real(kind_phys),
dimension(:,:),
intent(IN),
optional :: pattern_spp_pbl
580 real(kind_phys),
dimension(KTS:KTE) :: rstoch_col
584 real(kind_phys) :: delt2
591 wsp = sqrt(u(i,k)**2 + v(i,k)**2)
592 if (abs(hfx(i)) > 1200. .or. abs(qfx(i)) > 0.001 .or. &
593 wsp > 200. .or. t3d(i,k) > 360. .or. t3d(i,k) < 160. .or. &
594 sqv3d(i,k)< 0.0 .or. sqc3d(i,k)< 0.0 )
then
597 print*,
"Incoming problem at: i=",i,
" k=1"
598 print*,
" QFX=",qfx(i),
" HFX=",hfx(i)
599 print*,
" wsp=",wsp,
" T=",t3d(i,k)
600 print*,
" qv=",sqv3d(i,k),
" qc=",sqc3d(i,k)
601 print*,
" u*=",ust(i),
" wspd=",wspd(i)
602 print*,
" xland=",xland(i),
" ts=",ts(i)
603 print*,
" z/L=",0.5*dz(i,1)*rmol(i),
" ps=",ps(i)
604 print*,
" znt=",znt(i),
" dx=",dx(i)
608 print*,
"===tk:",t3d(i,max(kproblem-3,1):min(kproblem+3,kte))
609 print*,
"===qv:",sqv3d(i,max(kproblem-3,1):min(kproblem+3,kte))
610 print*,
"===qc:",sqc3d(i,max(kproblem-3,1):min(kproblem+3,kte))
611 print*,
"===qi:",sqi3d(i,max(kproblem-3,1):min(kproblem+3,kte))
612 print*,
"====u:",u(i,max(kproblem-3,1):min(kproblem+3,kte))
613 print*,
"====v:",v(i,max(kproblem-3,1):min(kproblem+3,kte))
627 IF (bl_mynn_output > 0)
THEN
628 edmf_a(its:ite,kts:kte)=0.
629 edmf_w(its:ite,kts:kte)=0.
630 edmf_qt(its:ite,kts:kte)=0.
631 edmf_thl(its:ite,kts:kte)=0.
632 edmf_ent(its:ite,kts:kte)=0.
633 edmf_qc(its:ite,kts:kte)=0.
634 sub_thl3d(its:ite,kts:kte)=0.
635 sub_sqv3d(its:ite,kts:kte)=0.
636 det_thl3d(its:ite,kts:kte)=0.
637 det_sqv3d(its:ite,kts:kte)=0.
646 ktop_plume(its:ite)=0
647 ztop_plume(its:ite)=0.
650 maxkhtopdown(its:ite)=0.
658 IF (initflag > 0 .and. .not.restart)
THEN
661 IF ( (restart .or. cycling))
THEN
662 IF (maxval(qke(its:ite,kts)) < 0.0002)
THEN
663 initialize_qke = .true.
666 initialize_qke = .false.
670 initialize_qke = .true.
674 if (.not.restart .or. .not.cycling)
THEN
675 sh3d(its:ite,kts:kte)=0.
676 sm3d(its:ite,kts:kte)=0.
677 el_pbl(its:ite,kts:kte)=0.
678 tsq(its:ite,kts:kte)=0.
679 qsq(its:ite,kts:kte)=0.
680 cov(its:ite,kts:kte)=0.
681 cldfra_bl(its:ite,kts:kte)=0.
682 qc_bl(its:ite,kts:kte)=0.
683 qke(its:ite,kts:kte)=0.
687 cldfra_bl1d(kts:kte)=0.0
697 qc_bl1d_old(kts:kte)=0.0
698 cldfra_bl1d_old(kts:kte)=0.0
701 edmf_qc1(kts:kte)=0.0
702 edmf_a_dd1(kts:kte)=0.0
703 edmf_w_dd1(kts:kte)=0.0
704 edmf_qc_dd1(kts:kte)=0.0
716 IF (tke_budget .eq. 1)
THEN
739 if (icloud_bl > 0)
then
740 cldfra_bl1d(:)=cldfra_bl(i,:)
741 qc_bl1d(:)=qc_bl(i,:)
742 qi_bl1d(:)=qi_bl(i,:)
756 thetav(k)=th(i,k)*(1.+p608*sqv(k))
758 sqw(k)=sqv(k)+sqc(k)+sqi(k)
759 thl(k)=th1(k) - xlvcp/ex1(k)*sqc(k) &
760 & - xlscp/ex1(k)*(sqi(k))
769 zw(k)=zw(k-1)+dz(i,k-1)
771 IF (initialize_qke)
THEN
775 qke1(k)=5.*ust(i) * max((ust(i)*700. - zw(k))/(max(ust(i),0.01)*700.), 0.01)
786 rstoch_col(k)=pattern_spp_pbl(i,k)
793 zw(kte+1)=zw(kte)+dz(i,kte)
796 CALL get_pblh(kts,kte,pblh(i),thetav,&
797 & qke1,zw,dz1,xland(i),kpbl(i))
801 IF (scaleaware > 0.)
THEN
802 CALL scale_aware(dx(i),pblh(i),psig_bl(i),psig_shcu(i))
817 &pblh(i), th1, thetav, sh, sm, &
819 &el, qke1, tsq1, qsq1, cov1, &
820 &psig_bl(i), cldfra_bl1d, &
821 &bl_mynn_mixlength, &
824 &spp_pbl,rstoch_col )
826 IF (.not.restart)
THEN
838 IF (bl_mynn_tkeadvect)
THEN
862 IF (bl_mynn_tkeadvect)
THEN
868 if (tke_budget .eq. 1)
then
881 if (icloud_bl > 0)
then
882 cldfra_bl1d(:)=cldfra_bl(i,:)
883 qc_bl1d(:) =qc_bl(i,:)
884 qi_bl1d(:) =qi_bl(i,:)
885 cldfra_bl1d_old(:)=cldfra_bl(i,:)
886 qc_bl1d_old(:)=qc_bl(i,:)
887 qi_bl1d_old(:)=qi_bl(i,:)
896 dz1(kts:kte) =dz(i,kts:kte)
897 u1(kts:kte) =u(i,kts:kte)
898 v1(kts:kte) =v(i,kts:kte)
899 w1(kts:kte) =w(i,kts:kte)
900 th1(kts:kte) =th(i,kts:kte)
901 tk1(kts:kte) =t3d(i,kts:kte)
902 p1(kts:kte) =p(i,kts:kte)
903 ex1(kts:kte) =exner(i,kts:kte)
904 rho1(kts:kte) =rho(i,kts:kte)
905 sqv(kts:kte) =sqv3d(i,kts:kte)
906 sqc(kts:kte) =sqc3d(i,kts:kte)
907 qv1(kts:kte) =sqv(kts:kte)/(1.-sqv(kts:kte))
908 qc1(kts:kte) =sqc(kts:kte)/(1.-sqv(kts:kte))
909 qi1(kts:kte) =sqi(kts:kte)/(1.-sqv(kts:kte))
910 qs1(kts:kte) =sqs(kts:kte)/(1.-sqv(kts:kte))
921 qni1(kts:kte)=qni(i,kts:kte)
926 qnc1(kts:kte)=qnc(i,kts:kte)
930 IF (flag_qnwfa )
THEN
931 qnwfa1(kts:kte)=qnwfa(i,kts:kte)
935 IF (flag_qnifa )
THEN
936 qnifa1(kts:kte)=qnifa(i,kts:kte)
940 IF (flag_qnbca )
THEN
941 qnbca1(kts:kte)=qnbca(i,kts:kte)
945 IF (flag_ozone )
THEN
946 ozone1(kts:kte)=ozone(i,kts:kte)
950 el(kts:kte) =el_pbl(i,kts:kte)
951 qke1(kts:kte)=qke(i,kts:kte)
952 sh(kts:kte) =sh3d(i,kts:kte)
953 sm(kts:kte) =sm3d(i,kts:kte)
954 tsq1(kts:kte)=tsq(i,kts:kte)
955 qsq1(kts:kte)=qsq(i,kts:kte)
956 cov1(kts:kte)=cov(i,kts:kte)
958 rstoch_col(kts:kte)=pattern_spp_pbl(i,kts:kte)
960 rstoch_col(kts:kte)=0.0
1005 zw(k)=zw(k-1)+dz(i,k-1)
1008 sqw(k)= sqv(k)+sqc(k)+sqi(k)
1009 thl(k)= th1(k) - xlvcp/ex1(k)*sqc(k) &
1010 & - xlscp/ex1(k)*(sqi(k))
1015 thetav(k)=th1(k)*(1.+p608*sqv(k))
1017 zw(kte+1)=zw(kte)+dz(i,kte)
1020 if ( mix_chem )
then
1022 vd1(ic) = vdep(i,ic)
1026 chem1(k,ic) = chem3d(i,k,ic)
1043 CALL get_pblh(kts,kte,pblh(i),thetav,&
1044 & qke1,zw,dz1,xland(i),kpbl(i))
1050 if (scaleaware > 0.)
then
1051 call scale_aware(dx(i),pblh(i),psig_bl(i),psig_shcu(i))
1058 cpm=cp*(1.+0.84*qv1(kts))
1059 exnerg=(ps(i)/p1000mb)**rcp
1068 flqv = qfx(i)/rho1(kts)
1070 th_sfc = ts(i)/ex1(kts)
1074 flt =hfx(i)/(rho1(kts)*cpm )-xlvcp*flqc/ex1(kts)
1075 fltv=flt + flqv*p608*th_sfc
1078 rmol(i) = -karman*gtr*fltv/max(ust(i)**3,1.0e-6)
1079 zet = 0.5*dz(i,kts)*rmol(i)
1080 zet = max(zet, -20.)
1083 if (bl_mynn_stfunc == 0)
then
1085 if ( zet >= 0.0 )
then
1086 pmz = 1.0 + (cphm_st-1.0) * zet
1087 phh = 1.0 + cphh_st * zet
1089 pmz = 1.0/ (1.0-cphm_unst*zet)**0.25 - zet
1090 phh = 1.0/sqrt(1.0-cphh_unst*zet)
1105 &dx(i),dz1,zw,xland(i), &
1106 &thl,sqw,sqv,sqc,sqi,sqs, &
1107 &p1,ex1,tsq1,qsq1,cov1, &
1108 &sh,el,bl_mynn_cloudpdf, &
1109 &qc_bl1d,qi_bl1d,cldfra_bl1d, &
1111 &vt, vq, th1, sgm, rmol(i), &
1112 &spp_pbl, rstoch_col )
1117 if (bl_mynn_topdown.eq.1)
then
1118 call topdown_cloudrad(kts,kte,dz1,zw, &
1119 &xland(i),kpbl(i),pblh(i), &
1120 &sqc,sqi,sqw,thl,th1,ex1,p1,rho1,thetav, &
1121 &cldfra_bl1d,rthraten(i,:), &
1122 &maxkhtopdown(i),khtopdown,tkeprodtd )
1124 maxkhtopdown(i) = 0.0
1125 khtopdown(kts:kte) = 0.0
1126 tkeprodtd(kts:kte) = 0.0
1129 if (bl_mynn_edmf > 0)
then
1132 &kts,kte,delt,zw,dz1,p1,rho1, &
1133 &bl_mynn_edmf_mom, &
1134 &bl_mynn_edmf_tke, &
1135 &bl_mynn_mixscalars, &
1136 &u1,v1,w1,th1,thl,thetav,tk1, &
1137 &sqw,sqv,sqc,qke1, &
1138 &qnc1,qni1,qnwfa1,qnifa1,qnbca1, &
1140 &ust(i),flt,fltv,flq,flqv, &
1141 &pblh(i),kpbl(i),dx(i), &
1146 &edmf_a1,edmf_w1,edmf_qt1, &
1147 &edmf_thl1,edmf_ent1,edmf_qc1, &
1149 &s_aw1,s_awthl1,s_awqt1, &
1151 &s_awu1,s_awv1,s_awqke1, &
1152 &s_awqnc1,s_awqni1, &
1153 &s_awqnwfa1,s_awqnifa1,s_awqnbca1, &
1156 &det_thl,det_sqv,det_sqc, &
1159 &nchem,chem1,s_awchem1, &
1161 &qc_bl1d,cldfra_bl1d, &
1162 &qc_bl1d_old,cldfra_bl1d_old, &
1164 &flag_qnc,flag_qni, &
1165 &flag_qnwfa,flag_qnifa,flag_qnbca, &
1167 &maxwidth(i),ktop_plume(i), &
1168 &maxmf(i),ztop_plume(i), &
1169 &spp_pbl,rstoch_col )
1172 if (bl_mynn_edmf_dd == 1)
then
1173 call ddmf_jpl(kts,kte,delt,zw,dz1,p1, &
1174 &u1,v1,th1,thl,thetav,tk1, &
1175 &sqw,sqv,sqc,rho1,ex1, &
1178 &edmf_a_dd1,edmf_w_dd1,edmf_qt_dd1, &
1179 &edmf_thl_dd1,edmf_ent_dd1, &
1181 &sd_aw1,sd_awthl1,sd_awqt1, &
1182 &sd_awqv1,sd_awqc1,sd_awu1,sd_awv1, &
1184 &qc_bl1d,cldfra_bl1d, &
1193 &kts,kte,xland(i),closure, &
1195 &u1, v1, thl, thetav, sqc, sqw, &
1196 &qke1, tsq1, qsq1, cov1, &
1198 &rmol(i), flt, fltv, flq, &
1204 &qwt1,qshear1,qbuoy1,qdiss1, &
1206 &psig_bl(i),psig_shcu(i), &
1207 &cldfra_bl1d,bl_mynn_mixlength, &
1210 &spp_pbl,rstoch_col )
1217 &ust(i), flt, flq, pmz, phh, &
1218 &el, dfq, rho1, pdk, pdt, pdq, pdc, &
1219 &qke1, tsq1, qsq1, cov1, &
1220 &s_aw1, s_awqke1, bl_mynn_edmf_tke, &
1221 &qwt1, qdiss1, tke_budget )
1223 if (dheat_opt > 0)
then
1226 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)
1228 diss_heat(k) = diss_heat(k) * exp(-10000./max(p1(k),1.))
1232 diss_heat(1:kte) = 0.
1239 &u1, v1, th1, tk1, qv1, &
1240 &qc1, qi1, kzero, qnc1, qni1, &
1241 &ps(i), p1, ex1, thl, &
1242 &sqv, sqc, sqi, kzero, sqw, &
1243 &qnwfa1, qnifa1, qnbca1, ozone1, &
1244 &ust(i),flt,flq,flqv,flqc, &
1245 &wspd(i),uoce(i),voce(i), &
1246 &tsq1, qsq1, cov1, &
1249 &du1, dv1, dth1, dqv1, &
1250 &dqc1, dqi1, dqs1, dqnc1, dqni1, &
1251 &dqnwfa1, dqnifa1, dqnbca1, &
1255 &s_aw1,s_awthl1,s_awqt1, &
1256 &s_awqv1,s_awqc1,s_awu1,s_awv1, &
1257 &s_awqnc1,s_awqni1, &
1258 &s_awqnwfa1,s_awqnifa1,s_awqnbca1, &
1259 &sd_aw1,sd_awthl1,sd_awqt1, &
1260 &sd_awqv1,sd_awqc1, &
1264 &det_thl,det_sqv,det_sqc, &
1266 &flag_qc,flag_qi,flag_qnc, &
1267 &flag_qni,flag_qs, &
1268 &flag_qnwfa,flag_qnifa, &
1271 &bl_mynn_cloudmix, &
1274 &bl_mynn_edmf_mom, &
1275 &bl_mynn_mixscalars )
1278 if ( mix_chem )
then
1280 call mynn_mix_chem(kts,kte,i, &
1281 &delt, dz1, pblh(i), &
1282 &nchem, kdvel, ndvel, &
1290 &enh_mix, smoke_dbg )
1292 call mynn_mix_chem(kts,kte,i, &
1293 &delt, dz1, pblh(i), &
1294 &nchem, kdvel, ndvel, &
1302 &enh_mix, smoke_dbg )
1306 chem3d(i,k,ic) = max(1.e-12, chem1(k,ic))
1312 &dfm, dfh, dz1, k_m1, k_h1 )
1315 exch_m(i,kts:kte) =k_m1(kts:kte)
1316 exch_h(i,kts:kte) =k_h1(kts:kte)
1317 rublten(i,kts:kte) =du1(kts:kte)
1318 rvblten(i,kts:kte) =dv1(kts:kte)
1319 rthblten(i,kts:kte)=dth1(kts:kte)
1320 rqvblten(i,kts:kte)=dqv1(kts:kte)
1321 if (bl_mynn_cloudmix > 0)
then
1322 if (flag_qc) rqcblten(i,kts:kte)=dqc1(kts:kte)
1323 if (flag_qi) rqiblten(i,kts:kte)=dqi1(kts:kte)
1324 if (flag_qs) rqsblten(i,kts:kte)=dqs1(kts:kte)
1326 if (flag_qc) rqcblten(i,:)=0.
1327 if (flag_qi) rqiblten(i,:)=0.
1328 if (flag_qs) rqsblten(i,:)=0.
1330 if (bl_mynn_cloudmix > 0 .and. bl_mynn_mixscalars > 0)
then
1331 if (flag_qnc) rqncblten(i,kts:kte) =dqnc1(kts:kte)
1332 if (flag_qni) rqniblten(i,kts:kte) =dqni1(kts:kte)
1333 if (flag_qnwfa) rqnwfablten(i,kts:kte)=dqnwfa1(kts:kte)
1334 if (flag_qnifa) rqnifablten(i,kts:kte)=dqnifa1(kts:kte)
1335 if (flag_qnbca) rqnbcablten(i,kts:kte)=dqnbca1(kts:kte)
1337 if (flag_qnc) rqncblten(i,:) =0.
1338 if (flag_qni) rqniblten(i,:) =0.
1339 if (flag_qnwfa) rqnwfablten(i,:)=0.
1340 if (flag_qnifa) rqnifablten(i,:)=0.
1341 if (flag_qnbca) rqnbcablten(i,:)=0.
1343 dozone(i,kts:kte)=dozone1(kts:kte)
1344 if (icloud_bl > 0)
then
1345 qc_bl(i,kts:kte) =qc_bl1d(kts:kte)
1346 qi_bl(i,kts:kte) =qi_bl1d(kts:kte)
1347 cldfra_bl(i,kts:kte)=cldfra_bl1d(kts:kte)
1349 el_pbl(i,kts:kte)=el(kts:kte)
1350 qke(i,kts:kte) =qke1(kts:kte)
1351 tsq(i,kts:kte) =tsq1(kts:kte)
1352 qsq(i,kts:kte) =qsq1(kts:kte)
1353 cov(i,kts:kte) =cov1(kts:kte)
1354 sh3d(i,kts:kte) =sh(kts:kte)
1355 sm3d(i,kts:kte) =sm(kts:kte)
1357 if (tke_budget .eq. 1)
then
1361 qshear1(k) =4.*(ust(i)**3*phi_m/(karman*dz(i,k)))-qshear1(k+1)
1362 qbuoy1(k) =4.*(-ust(i)**3*zet/(karman*dz(i,k)))-qbuoy1(k+1)
1365 qshear(i,k)=0.5*(qshear1(k)+qshear1(k+1))
1366 qbuoy(i,k) =0.5*(qbuoy1(k)+qbuoy1(k+1))
1368 qdiss(i,k) =qdiss1(k)
1369 dqke(i,k) =(qke1(k)-dqke(i,k))*0.5/delt
1381 if (bl_mynn_output > 0)
then
1382 if (bl_mynn_edmf > 0)
then
1383 edmf_a(i,kts:kte) =edmf_a1(kts:kte)
1384 edmf_w(i,kts:kte) =edmf_w1(kts:kte)
1385 edmf_qt(i,kts:kte) =edmf_qt1(kts:kte)
1386 edmf_thl(i,kts:kte) =edmf_thl1(kts:kte)
1387 edmf_ent(i,kts:kte) =edmf_ent1(kts:kte)
1388 edmf_qc(i,kts:kte) =edmf_qc1(kts:kte)
1389 sub_thl3d(i,kts:kte)=sub_thl(kts:kte)
1390 sub_sqv3d(i,kts:kte)=sub_sqv(kts:kte)
1391 det_thl3d(i,kts:kte)=det_thl(kts:kte)
1392 det_sqv3d(i,kts:kte)=det_sqv(kts:kte)
1405 if ( debug_code .and. (i .eq. idbg))
THEN
1406 if ( abs(qfx(i))>.001)print*,&
1407 "SUSPICIOUS VALUES AT: i=",i,
" QFX=",qfx(i)
1408 if ( abs(hfx(i))>1100.)print*,&
1409 "SUSPICIOUS VALUES AT: i=",i,
" HFX=",hfx(i)
1411 IF ( sh(k) < 0. .OR. sh(k)> 200.)print*,&
1412 "SUSPICIOUS VALUES AT: i,k=",i,k,
" sh=",sh(k)
1413 IF ( abs(vt(k)) > 2.0 )print*,&
1414 "SUSPICIOUS VALUES AT: i,k=",i,k,
" vt=",vt(k)
1415 IF ( abs(vq(k)) > 7000.)print*,&
1416 "SUSPICIOUS VALUES AT: i,k=",i,k,
" vq=",vq(k)
1417 IF ( qke(i,k) < -1. .OR. qke(i,k)> 200.)print*,&
1418 "SUSPICIOUS VALUES AT: i,k=",i,k,
" qke=",qke(i,k)
1419 IF ( el_pbl(i,k) < 0. .OR. el_pbl(i,k)> 1500.)print*,&
1420 "SUSPICIOUS VALUES AT: i,k=",i,k,
" el_pbl=",el_pbl(i,k)
1421 IF ( exch_m(i,k) < 0. .OR. exch_m(i,k)> 2000.)print*,&
1422 "SUSPICIOUS VALUES AT: i,k=",i,k,
" exxch_m=",exch_m(i,k)
1423 IF (icloud_bl > 0)
then
1424 IF ( cldfra_bl(i,k) < 0.0 .OR. cldfra_bl(i,k)> 1.)
THEN
1425 print*,
"SUSPICIOUS VALUES: CLDFRA_BL=",cldfra_bl(i,k),
" qc_bl=",qc_bl(i,k)
1443 IF (bl_mynn_tkeadvect)
THEN
1448#ifdef HARDCODE_VERTICAL
1852 & rmo, flt, fltv, flq, &
1858 & Psig_bl, cldfra_bl1D, &
1859 & bl_mynn_mixlength, &
1864 integer,
intent(in) :: kts,kte
1866#ifdef HARDCODE_VERTICAL
1868# define kte HARDCODE_VERTICAL
1871 integer,
intent(in) :: bl_mynn_mixlength
1872 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz
1873 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
1874 real(kind_phys),
intent(in) :: rmo,flt,fltv,flq,Psig_bl,xland
1875 real(kind_phys),
intent(in) :: dx,zi
1876 real(kind_phys),
dimension(kts:kte),
intent(in) :: u1,v1, &
1877 &qke,vt,vq,cldfra_bl1D,edmf_w1,edmf_a1
1878 real(kind_phys),
dimension(kts:kte),
intent(out) :: qkw, el
1879 real(kind_phys),
dimension(kts:kte),
intent(in) :: dtv
1880 real(kind_phys):: elt,vsc
1881 real(kind_phys),
dimension(kts:kte),
intent(in) :: theta
1882 real(kind_phys),
dimension(kts:kte) :: qtke,elBLmin,elBLavg,thetaw
1883 real(kind_phys):: wt,wt2,zi2,h1,h2,hs,elBLmin0,elBLavg0,cldavg
1887 real(kind_phys):: cns, & !< for surface layer (els) in stable conditions
1888 alp1, & !< for turbulent length scale (elt)
1889 alp2, & !< for buoyancy length scale (elb)
1890 alp3, & !< for buoyancy enhancement factor of elb
1891 alp4, & !< for surface layer (els) in unstable conditions
1892 alp5, & !< for BouLac mixing length or above PBLH
1899 real(kind_phys),
parameter :: minzi = 300.
1900 real(kind_phys),
parameter :: maxdz = 750.
1903 real(kind_phys),
parameter :: mindz = 300.
1906 real(kind_phys),
parameter :: ZSLH = 100.
1907 real(kind_phys),
parameter :: CSL = 2.
1911 real(kind_phys):: afk,abk,zwk,zwk1,dzk,qdz,vflx,bv,tau_cloud, &
1912 & wstar,elb,els,elf,el_stab,el_mf,el_stab_mf,elb_mf, &
1913 & PBLH_PLUS_ENT,Uonset,Ugrid,wt_u,el_les
1914 real(kind_phys),
parameter :: ctau = 1000.
1919 SELECT CASE(bl_mynn_mixlength)
1931 zi2 = min(10000.,zw(kte-2))
1932 h1=max(0.3*zi2,mindz)
1936 qkw(kts) = sqrt(max(qke(kts),1.0e-10))
1938 afk = dz(k)/( dz(k)+dz(k-1) )
1940 qkw(k) = sqrt(max(qke(k)*abk+qke(k-1)*afk,1.0e-3))
1949 DO WHILE (zwk .LE. zi2+h1)
1950 dzk = 0.5*( dz(k)+dz(k-1) )
1951 qdz = max( qkw(k)-qmin, 0.03 )*dzk
1959 vflx = ( vt(kts)+1.0 )*flt +( vq(kts)+tv0 )*flq
1960 vsc = ( gtr*elt*max( vflx, 0.0 ) )**(1.0/3.0)
1970 IF ( dtv(k) .GT. 0.0 )
THEN
1971 bv = sqrt( gtr*dtv(k) )
1972 elb = alp2*qkw(k) / bv &
1973 & *( 1.0 + alp3/alp2*&
1974 &sqrt( vsc/( bv*elt ) ) )
1975 elf = alp2 * qkw(k)/bv
1983 IF ( rmo .GT. 0.0 )
THEN
1984 els = karman*zwk/(1.0+cns*min( zwk*rmo, zmax ))
1986 els = karman*zwk*( 1.0 - alp4* zwk*rmo )**0.2
1993 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
1995 el(k) = min(elb/( elb/elt+elb/els+1.0 ),elf)
2001 ugrid = sqrt(u1(kts)**2 + v1(kts)**2)
2003 wt_u = (1.0 - min(max(ugrid - uonset, 0.0)/30.0, 0.5))
2014 h1=max(0.3*zi2,300.)
2018 qtke(kts)=max(0.5*qke(kts), 0.01)
2019 thetaw(kts)=theta(kts)
2020 qkw(kts) = sqrt(max(qke(kts),1.0e-10))
2023 afk = dz(k)/( dz(k)+dz(k-1) )
2025 qkw(k) = sqrt(max(qke(k)*abk+qke(k-1)*afk,1.0e-3))
2026 qtke(k) = 0.5*(qkw(k)**2)
2027 thetaw(k)= theta(k)*abk + theta(k-1)*afk
2036 DO WHILE (zwk .LE. zi2+h1)
2037 dzk = 0.5*( dz(k)+dz(k-1) )
2038 qdz = min(max( qkw(k)-qmin, 0.03 ), 30.0)*dzk
2045 elt = min( max( alp1*elt/vsc, 10.), 400.)
2049 vsc = ( gtr*elt*max( vflx, 0.0 ) )**onethird
2056 CALL boulac_length(kts,kte,zw,dz,qtke,thetaw,elblmin,elblavg)
2062 IF ( dtv(k) .GT. 0.0 )
THEN
2063 bv = max( sqrt( gtr*dtv(k) ), 0.0001)
2064 elb = max(alp2*qkw(k), &
2065 & alp6*edmf_a1(k-1)*edmf_w1(k-1)) / bv &
2066 & *( 1.0 + alp3*sqrt( vsc/(bv*elt) ) )
2068 elf = 1.0 * qkw(k)/bv
2069 elblavg(k) = max(elblavg(k), alp6*edmf_a1(k-1)*edmf_w1(k-1)/bv)
2076 IF ( rmo .GT. 0.0 )
THEN
2077 els = karman*zwk/(1.0+cns*min( zwk*rmo, zmax ))
2079 els = karman*zwk*( 1.0 - alp4* zwk*rmo )**0.2
2083 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
2090 el(k) = sqrt( els**2/(1. + (els**2/elt**2)))
2091 el(k) = min(el(k), elb)
2092 el(k) = min(el(k), elf)
2093 el(k) = el(k)*(1.-wt) + alp5*elblavg(k)*wt
2096 el(k) = el(k)*psig_bl
2102 uonset = 3.5 + dz(kts)*0.1
2103 ugrid = sqrt(u1(kts)**2 + v1(kts)**2)
2117 h1=max(0.3*zi2,300.)
2121 qtke(kts)=max(0.5*qke(kts),0.01)
2122 qkw(kts) = sqrt(max(qke(kts),1.0e-4))
2125 afk = dz(k)/( dz(k)+dz(k-1) )
2127 qkw(k) = sqrt(max(qke(k)*abk+qke(k-1)*afk,1.0e-3))
2128 qtke(k) = 0.5*qkw(k)**2
2135 pblh_plus_ent = max(zi+h1, 100.)
2138 DO WHILE (zwk .LE. pblh_plus_ent)
2139 dzk = 0.5*( dz(k)+dz(k-1) )
2140 qdz = min(max( qkw(k)-qmin, 0.03 ), 30.0)*dzk
2147 elt = min( max(alp1*elt/vsc, 10.), 400.)
2151 vsc = ( gtr*elt*max( vflx, 0.0 ) )**onethird
2159 dzk = 0.5*( dz(k)+dz(k-1) )
2160 cldavg = 0.5*(cldfra_bl1d(k-1)+cldfra_bl1d(k))
2163 IF ( dtv(k) .GT. 0.0 )
THEN
2165 bv = max( sqrt( gtr*dtv(k) ), 0.001)
2167 elb_mf = max(alp2*qkw(k), &
2168 & alp6*edmf_a1(k-1)*edmf_w1(k-1)) / bv &
2169 & *( 1.0 + alp3*sqrt( vsc/( bv*elt ) ) )
2170 elb = min(max(alp5*qkw(k), alp6*edmf_a1(k)*edmf_w1(k))/bv, zwk)
2173 wstar = 1.25*(gtr*zi*max(vflx,1.0e-4))**onethird
2174 tau_cloud = min(max(ctau * wstar/grav, 30.), 150.)
2176 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
2177 tau_cloud = tau_cloud*(1.-wt) + 50.*wt
2178 elf = min(max(tau_cloud*sqrt(min(qtke(k),40.)), &
2179 & alp6*edmf_a1(k)*edmf_w1(k)/bv), zwk)
2200 wstar = 1.25*(gtr*zi*max(vflx,1.0e-4))**onethird
2201 tau_cloud = min(max(ctau * wstar/grav, 50.), 200.)
2203 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
2205 tau_cloud = tau_cloud*(1.-wt) + max(100.,dzk*0.25)*wt
2207 elb = min(tau_cloud*sqrt(min(qtke(k),40.)), zwk)
2212 elf = elf/(1. + (elf/800.))
2213 elb_mf = max(elb_mf, 0.01)
2216 IF ( rmo .GT. 0.0 )
THEN
2217 els = karman*zwk/(1.0+cns*min( zwk*rmo, zmax ))
2219 els = karman*zwk*( 1.0 - alp4* zwk*rmo )**0.2
2223 wt=.5*tanh((zwk - (zi2+h1))/h2) + .5
2226 el(k) = sqrt( els**2/(1. + (els**2/elt**2) +(els**2/elb_mf**2)))
2227 el(k) = el(k)*(1.-wt) + elf*wt
2230 el_les= min(els/(1. + (els/12.)), elb_mf)
2231 el(k) = el(k)*psig_bl + (1.-psig_bl)*el_les
2238#ifdef HARDCODE_VERTICAL
2622 & u, v, thl, thetav, ql, qw, &
2623 & qke, tsq, qsq, cov, &
2625 & rmo, flt, fltv, flq, &
2629 & Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, &
2630 & qWT1D,qSHEAR1D,qBUOY1D,qDISS1D, &
2632 & Psig_bl,Psig_shcu,cldfra_bl1D, &
2633 & bl_mynn_mixlength, &
2634 & edmf_w1,edmf_a1, &
2636 & spp_pbl,rstoch_col )
2640 integer,
intent(in) :: kts,kte
2642#ifdef HARDCODE_VERTICAL
2644# define kte HARDCODE_VERTICAL
2647 integer,
intent(in) :: bl_mynn_mixlength,tke_budget
2648 real(kind_phys),
intent(in) :: closure
2649 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz
2650 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
2651 real(kind_phys),
intent(in) :: rmo,flt,fltv,flq, &
2652 &Psig_bl,Psig_shcu,xland,dx,zi
2653 real(kind_phys),
dimension(kts:kte),
intent(in) :: u,v,thl,thetav,qw, &
2654 &ql,vt,vq,qke,tsq,qsq,cov,cldfra_bl1D,edmf_w1,edmf_a1, &
2657 real(kind_phys),
dimension(kts:kte),
intent(out) :: dfm,dfh,dfq, &
2658 &pdk,pdt,pdq,pdc,tcd,qcd,el
2660 real(kind_phys),
dimension(kts:kte),
intent(inout) :: &
2661 qWT1D,qSHEAR1D,qBUOY1D,qDISS1D
2662 real(kind_phys):: q3sq_old,dlsq1,qWTP_old,qWTP_new
2663 real(kind_phys):: dudz,dvdz,dTdz,upwp,vpwp,Tpwp
2665 real(kind_phys),
dimension(kts:kte) :: qkw,dtl,dqw,dtv,gm,gh,sm,sh
2669 real(kind_phys):: e6c,dzk,afk,abk,vtt,vqq, &
2670 &cw25,clow,cupp,gamt,gamq,smd,gamv,elq,elh
2672 real(kind_phys):: cldavg
2673 real(kind_phys),
dimension(kts:kte),
intent(in) :: theta
2675 real(kind_phys):: a2fac, duz, ri
2677 real:: auh,aum,adh,adm,aeh,aem,Req,Rsl,Rsl2, &
2678 gmelq,sm20,sh20,sm25max,sh25max,sm25min,sh25min, &
2679 sm_pbl,sh_pbl,zi2,wt,slht,wtpr
2681 DOUBLE PRECISION q2sq, t2sq, r2sq, c2sq, elsq, gmel, ghel
2682 DOUBLE PRECISION q3sq, t3sq, r3sq, c3sq, dlsq, qdiv
2683 DOUBLE PRECISION e1, e2, e3, e4, enum, eden, wden
2686 integer,
intent(in) :: spp_pbl
2687 real(kind_phys),
dimension(kts:kte) :: rstoch_col
2688 real(kind_phys):: Prnum, shb
2689 real(kind_phys),
parameter :: Prlimit = 5.0
2706 & u, v, thl, thetav, qw, &
2708 & dtl, dqw, dtv, gm, gh, sm, sh )
2713 & rmo, flt, fltv, flq, &
2719 & qkw,psig_bl,cldfra_bl1d, &
2720 & bl_mynn_mixlength, &
2725 dzk = 0.5 *( dz(k)+dz(k-1) )
2726 afk = dz(k)/( dz(k)+dz(k-1) )
2730 q2sq = b1*elsq*( sm(k)*gm(k)+sh(k)*gh(k) )
2732 sh20 = max(sh(k), 1e-5)
2733 sm20 = max(sm(k), 1e-5)
2734 sh(k)= max(sh(k), 1e-5)
2737 duz = ( u(k)-u(k-1) )**2 +( v(k)-v(k-1) )**2
2740 ri = -gh(k)/max( duz, 1.0e-10 )
2741 IF (ckmod .eq. 1)
THEN
2742 a2fac = 1./(1. + max(ri,0.0))
2753 prnum = min(0.76 + 4.0*max(ri,0.0), prlimit)
2762 IF ( debug_code )
THEN
2763 IF (sh(k)<0.0 .OR. sm(k)<0.0)
THEN
2764 print*,
"MYNN; mym_turbulence 2.0; sh=",sh(k),
" k=",k
2765 print*,
" gm=",gm(k),
" gh=",gh(k),
" sm=",sm(k)
2766 print*,
" q2sq=",q2sq,
" q3sq=",q3sq,
" q3/q2=",q3sq/q2sq
2767 print*,
" qke=",qke(k),
" el=",el(k),
" ri=",ri
2768 print*,
" PBLH=",zi,
" u=",u(k),
" v=",v(k)
2777 IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k)
2779 IF ( q3sq .LT. q2sq )
THEN
2781 qdiv = sqrt( q3sq/q2sq )
2798 sh(k) = sh(k) * qdiv
2799 sm(k) = sm(k) * qdiv
2816 e1 = q3sq - e1c*ghel*a2fac * qdiv**2
2817 e2 = q3sq - e2c*ghel*a2fac * qdiv**2
2818 e3 = e1 + e3c*ghel*a2fac**2 * qdiv**2
2819 e4 = e1 - e4c*ghel*a2fac * qdiv**2
2820 eden = e2*e4 + e3*e5c*gmel * qdiv**2
2821 eden = max( eden, 1.0d-20 )
2832 e1 = q3sq - e1c*ghel*a2fac
2833 e2 = q3sq - e2c*ghel*a2fac
2834 e3 = e1 + e3c*ghel*a2fac**2
2835 e4 = e1 - e4c*ghel*a2fac
2836 eden = e2*e4 + e3*e5c*gmel
2837 eden = max( eden, 1.0d-20 )
2841 sm(k) = q3sq*a1*( e3-3.0*c1*e4 )/eden
2845 sh(k) = q3sq*(a2*a2fac)*( e2+3.0*c1*e5c*gmel )/eden
2855 gmelq = max(gmel/q3sq, 1e-8)
2863 IF ( debug_code )
THEN
2864 IF ((sh(k)<sh25min .OR. sm(k)<sm25min .OR. &
2865 sh(k)>sh25max .OR. sm(k)>sm25max) )
THEN
2866 print*,
"In mym_turbulence 2.5: k=",k
2867 print*,
" sm=",sm(k),
" sh=",sh(k)
2868 print*,
" ri=",ri,
" Pr=",sm(k)/max(sh(k),1e-8)
2869 print*,
" gm=",gm(k),
" gh=",gh(k)
2870 print*,
" q2sq=",q2sq,
" q3sq=",q3sq, q3sq/q2sq
2871 print*,
" qke=",qke(k),
" el=",el(k)
2872 print*,
" PBLH=",zi,
" u=",u(k),
" v=",v(k)
2873 print*,
" SMnum=",q3sq*a1*( e3-3.0*c1*e4),
" SMdenom=",eden
2874 print*,
" SHnum=",q3sq*(a2*a2fac)*( e2+3.0*c1*e5c*gmel ),&
2880 IF ( sh(k) > sh25max ) sh(k) = sh25max
2881 IF ( sh(k) < sh25min ) sh(k) = sh25min
2893 shb = max(sh(k), 0.002)
2894 sm(k) = min(sm(k), prlimit*shb)
2897 IF ( closure .GE. 3.0 )
THEN
2898 t2sq = qdiv*b2*elsq*sh(k)*dtl(k)**2
2899 r2sq = qdiv*b2*elsq*sh(k)*dqw(k)**2
2900 c2sq = qdiv*b2*elsq*sh(k)*dtl(k)*dqw(k)
2901 t3sq = max( tsq(k)*abk+tsq(k-1)*afk, 0.0 )
2902 r3sq = max( qsq(k)*abk+qsq(k-1)*afk, 0.0 )
2903 c3sq = cov(k)*abk+cov(k-1)*afk
2906 c3sq = sign( min( abs(c3sq), sqrt(t3sq*r3sq) ), c3sq )
2908 vtt = 1.0 +vt(k)*abk +vt(k-1)*afk
2909 vqq = tv0 +vq(k)*abk +vq(k-1)*afk
2911 t2sq = vtt*t2sq +vqq*c2sq
2912 r2sq = vtt*c2sq +vqq*r2sq
2913 c2sq = max( vtt*t2sq+vqq*r2sq, 0.0d0 )
2914 t3sq = vtt*t3sq +vqq*c3sq
2915 r3sq = vtt*c3sq +vqq*r3sq
2916 c3sq = max( vtt*t3sq+vqq*r3sq, 0.0d0 )
2918 cw25 = e1*( e2 + 3.0*c1*e5c*gmel*qdiv**2 )/( 3.0*eden )
2922 IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k)
2927 auh = 27.*a1*((a2*a2fac)**2)*b2*(gtr)**2
2928 aum = 54.*(a1**2)*(a2*a2fac)*b2*c1*(gtr)
2929 adh = 9.*a1*((a2*a2fac)**2)*(12.*a1 + 3.*b2)*(gtr)**2
2930 adm = 18.*(a1**2)*(a2*a2fac)*(b2 - 3.*(a2*a2fac))*(gtr)
2932 aeh = (9.*a1*((a2*a2fac)**2)*b1 +9.*a1*((a2*a2fac)**2)* &
2933 (12.*a1 + 3.*b2))*(gtr)
2934 aem = 3.*a1*(a2*a2fac)*b1*(3.*(a2*a2fac) + 3.*b2*c1 + &
2935 (18.*a1*c1 - b2)) + &
2936 (18.)*(a1**2)*(a2*a2fac)*(b2 - 3.*(a2*a2fac))
2939 rsl = (auh + aum*req)/(3.*adh + 3.*adm*req)
2954 e2 = q3sq - e2c*ghel*a2fac * qdiv**2
2955 e3 = q3sq + e3c*ghel*a2fac**2 * qdiv**2
2956 e4 = q3sq - e4c*ghel*a2fac * qdiv**2
2957 eden = e2*e4 + e3 *e5c*gmel * qdiv**2
2962 wden = cc3*gtr**2 * dlsq**2/elsq * qdiv**2 &
2963 & *( e2*e4c*a2fac - e3c*e5c*gmel*a2fac**2 * qdiv**2 )
2965 IF ( wden .NE. 0.0 )
THEN
2967 clow = q3sq*( 0.12-cw25 )*eden/wden
2968 cupp = q3sq*( 0.76-cw25 )*eden/wden
2972 IF ( wden .GT. 0.0 )
THEN
2973 c3sq = min( max( c3sq, c2sq+clow ), c2sq+cupp )
2975 c3sq = max( min( c3sq, c2sq+clow ), c2sq+cupp )
2979 e1 = e2 + e5c*gmel * qdiv**2
2980 eden = max( eden, 1.0d-20 )
2985 e6c = 3.0*(a2*a2fac)*cc3*gtr * dlsq/elsq
2990 IF ( t2sq .GE. 0.0 )
THEN
2991 enum = max( qdiv*e6c*( t3sq-t2sq ), 0.0d0 )
2993 enum = min( qdiv*e6c*( t3sq-t2sq ), 0.0d0 )
2995 gamt =-e1 *
enum /eden
3000 IF ( r2sq .GE. 0.0 )
THEN
3001 enum = max( qdiv*e6c*( r3sq-r2sq ), 0.0d0 )
3003 enum = min( qdiv*e6c*( r3sq-r2sq ), 0.0d0 )
3005 gamq =-e1 *
enum /eden
3010 enum = max( qdiv*e6c*( c3sq-c2sq ), 0.0d0)
3014 smd = dlsq*enum*gtr/eden * qdiv**2 * (e3c*a2fac**2 + &
3015 & e4c*a2fac)*a1/(a2*a2fac)
3017 gamv = e1 *enum*gtr/eden
3025 IF ( debug_code )
THEN
3026 IF (sh(k)<-0.3 .OR. sm(k)<-0.3 .OR. &
3027 qke(k) < -0.1 .or. abs(smd) .gt. 2.0)
THEN
3028 print*,
" MYNN; mym_turbulence3.0; sh=",sh(k),
" k=",k
3029 print*,
" gm=",gm(k),
" gh=",gh(k),
" sm=",sm(k)
3030 print*,
" q2sq=",q2sq,
" q3sq=",q3sq,
" q3/q2=",q3sq/q2sq
3031 print*,
" qke=",qke(k),
" el=",el(k),
" ri=",ri
3032 print*,
" PBLH=",zi,
" u=",u(k),
" v=",v(k)
3047 cldavg = 0.5*(cldfra_bl1d(k-1) + cldfra_bl1d(k))
3048 IF (edmf_a1(k) > 0.001 .OR. cldavg > 0.02)
THEN
3050 sm(k) = max(sm(k), 0.03*min(10.*edmf_a1(k)*edmf_w1(k),1.0) )
3051 sh(k) = max(sh(k), 0.03*min(10.*edmf_a1(k)*edmf_w1(k),1.0) )
3053 sm(k) = max(sm(k), 0.05*min(cldavg,1.0) )
3054 sh(k) = max(sh(k), 0.05*min(cldavg,1.0) )
3062 pdk(k) = elq*( sm(k)*gm(k) &
3063 & +sh(k)*gh(k)+gamv ) + &
3065 pdt(k) = elh*( sh(k)*dtl(k)+gamt )*dtl(k)
3066 pdq(k) = elh*( sh(k)*dqw(k)+gamq )*dqw(k)
3067 pdc(k) = elh*( sh(k)*dtl(k)+gamt ) &
3069 & + elh*( sh(k)*dqw(k)+gamq )*dtl(k)*0.5
3076 dfm(k) = elq*sm(k) / dzk
3077 dfh(k) = elq*sh(k) / dzk
3084 IF (tke_budget .eq. 1)
THEN
3100 qshear1d(k) = elq*sm(k)*gm(k)
3108 qbuoy1d(k) = elq*(sh(k)*gh(k)+gamv)+0.5*tkeprodtd(k)
3131 tcd(k) = ( tcd(k+1)-tcd(k) )/( dzk )
3132 qcd(k) = ( qcd(k+1)-qcd(k) )/( dzk )
3135 if (spp_pbl==1)
then
3137 dfm(k)= dfm(k) + dfm(k)* rstoch_col(k) * 1.5 * max(exp(-max(zw(k)-8000.,0.0)/2000.),0.001)
3138 dfh(k)= dfh(k) + dfh(k)* rstoch_col(k) * 1.5 * max(exp(-max(zw(k)-8000.,0.0)/2000.),0.001)
3143#ifdef HARDCODE_VERTICAL
3603 & dx, dz, zw, xland, &
3604 & thl, qw, qv, qc, qi, qs, &
3607 & Sh, el, bl_mynn_cloudpdf, &
3608 & qc_bl1D, qi_bl1D, &
3611 & Vt, Vq, th, sgm, rmo, &
3612 & spp_pbl,rstoch_col )
3616 integer,
intent(in) :: kts,kte, bl_mynn_cloudpdf
3618#ifdef HARDCODE_VERTICAL
3620# define kte HARDCODE_VERTICAL
3623 real(kind_phys),
intent(in) :: HFX1,rmo,xland
3624 real(kind_phys),
intent(in) :: dx,pblh1
3625 real(kind_phys),
dimension(kts:kte),
intent(in) :: dz
3626 real(kind_phys),
dimension(kts:kte+1),
intent(in) :: zw
3627 real(kind_phys),
dimension(kts:kte),
intent(in) :: p,exner,thl,qw, &
3628 &qv,qc,qi,qs,tsq,qsq,cov,th
3630 real(kind_phys),
dimension(kts:kte),
intent(inout) :: vt,vq,sgm
3632 real(kind_phys),
dimension(kts:kte) :: alp,a,bet,b,ql,q1,RH
3633 real(kind_phys),
dimension(kts:kte),
intent(out) :: qc_bl1D,qi_bl1D, &
3635 DOUBLE PRECISION :: t3sq, r3sq, c3sq
3637 real(kind_phys):: qsl,esat,qsat,dqsl,cld0,q1k,qlk,eq1,qll, &
3638 &q2p,pt,rac,qt,t,xl,rsl,cpm,Fng,qww,alpha,beta,bb, &
3639 &ls,wt,wt2,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, &
3640 &qmq,qsat_tk,q1_rh,rh_hack,dzm1,zsl,maxqc
3641 real(kind_phys),
parameter :: qpct_sfc=0.025
3642 real(kind_phys),
parameter :: qpct_pbl=0.030
3643 real(kind_phys),
parameter :: qpct_trp=0.040
3644 real(kind_phys),
parameter :: rhcrit =0.83
3645 real(kind_phys),
parameter :: rhmax =1.02
3648 real(kind_phys):: erf
3651 real:: dth,dtl,dqw,dzk,els
3652 real(kind_phys),
dimension(kts:kte),
intent(in) :: Sh,el
3655 real(kind_phys) :: zagl,damp,PBLH2
3656 real(kind_phys) :: cfmax
3659 real(kind_phys) :: theta1, theta2, ht1, ht2
3663 integer,
intent(in) :: spp_pbl
3664 real(kind_phys),
dimension(kts:kte) :: rstoch_col
3665 real(kind_phys) :: qw_pert
3672 DO k = kte-3, kts, -1
3675 ht1 = 44307.692 * (1.0 - (p(k)/101325.)**0.190)
3676 ht2 = 44307.692 * (1.0 - (p(k+2)/101325.)**0.190)
3677 if ( (((theta2-theta1)/(ht2-ht1)) .lt. 10./1500. ) .AND. &
3678 & (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
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
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
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))
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)) + &
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
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)
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)
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
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
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
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
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
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)
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)+sqs(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
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.01)/0.03) + .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) +rhoz(k)*upa(k,i)*upw(k,i)
6509 edmf_qt(k) =edmf_qt(k) +rhoz(k)*upa(k,i)*upqt(k,i)
6510 edmf_thl(k)=edmf_thl(k)+rhoz(k)*upa(k,i)*upthl(k,i)
6511 edmf_ent(k)=edmf_ent(k)+rhoz(k)*upa(k,i)*ent(k,i)
6512 edmf_qc(k) =edmf_qc(k) +rhoz(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
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)
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...