MODULE module_bl_mynn USE module_model_constants, only: & &karman, g, p1000mb, & &cp, r_d, r_v, rcp, xlv, xlf, xls, & &svp1, svp2, svp3, svpt0, ep_1, ep_2, rvovrd, & &cpv, cliq, cice USE module_state_description, only: param_first_scalar, & &p_qc, p_qr, p_qi, p_qs, p_qg, p_qnc, p_qni USE module_wrf_error IMPLICIT NONE REAL, PARAMETER :: cphm_st=5.0, cphm_unst=16.0, & cphh_st=5.0, cphh_unst=16.0 REAL, PARAMETER :: xlvcp=xlv/cp, xlscp=(xlv+xlf)/cp, ev=xlv, rd=r_d, & &rk=cp/rd, svp11=svp1*1.e3, p608=ep_1, ep_3=1.-ep_2 REAL, PARAMETER :: tref=300.0 REAL, PARAMETER :: TKmin=253.0 REAL, PARAMETER :: tv0=p608*tref, tv1=(1.+p608)*tref, gtr=g/tref REAL, PARAMETER :: & &vk = karman, & &pr = 0.74, & &g1 = 0.235, & &b1 = 24.0, & &b2 = 15.0, & &c2 = 0.729, & &c3 = 0.340, & &c4 = 0.0, & &c5 = 0.2, & &a1 = b1*( 1.0-3.0*g1 )/6.0, & &c1 = g1 -1.0/( 3.0*a1*2.88449914061481660), & &a2 = a1*( g1-c1 )/( g1*pr ), & &g2 = b2/b1*( 1.0-c3 ) +2.0*a1/b1*( 3.0-2.0*c2 ) REAL, PARAMETER :: & &cc2 = 1.0-c2, & &cc3 = 1.0-c3, & &e1c = 3.0*a2*b2*cc3, & &e2c = 9.0*a1*a2*cc2, & &e3c = 9.0*a2*a2*cc2*( 1.0-c5 ), & &e4c = 12.0*a1*a2*cc2, & &e5c = 6.0*a1*a1 REAL, PARAMETER :: qmin=0.0, zmax=1.0, Sqfac=3.0 REAL, PARAMETER :: gno=1.0 REAL, PARAMETER :: gpw=5./3., qcgmin=1.e-8, qkemin=1.e-12 REAL, PARAMETER :: rr2=0.7071068, rrp=0.3989423 REAL, PARAMETER :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0 REAL, PARAMETER :: CKmod=1. REAL, PARAMETER :: scaleaware=1. INTEGER, PARAMETER :: bl_mynn_mixchem = 0 REAL, PARAMETER:: J0= .611583699E03 REAL, PARAMETER:: J1= .444606896E02 REAL, PARAMETER:: J2= .143177157E01 REAL, PARAMETER:: J3= .264224321E-1 REAL, PARAMETER:: J4= .299291081E-3 REAL, PARAMETER:: J5= .203154182E-5 REAL, PARAMETER:: J6= .702620698E-8 REAL, PARAMETER:: J7= .379534310E-11 REAL, PARAMETER:: J8=-.321582393E-13 REAL, PARAMETER:: K0= .609868993E03 REAL, PARAMETER:: K1= .499320233E02 REAL, PARAMETER:: K2= .184672631E01 REAL, PARAMETER:: K3= .402737184E-1 REAL, PARAMETER:: K4= .565392987E-3 REAL, PARAMETER:: K5= .521693933E-5 REAL, PARAMETER:: K6= .307839583E-7 REAL, PARAMETER:: K7= .105785160E-9 REAL, PARAMETER:: K8= .161444444E-12 INTEGER :: mynn_level CHARACTER*128 :: mynn_message INTEGER, PARAMETER :: kdebug=27 CONTAINS SUBROUTINE mym_initialize ( & & kts,kte, & & dz, zw, & & u, v, thl, qw, & & zi, theta, sh, & & ust, rmo, el, & & Qke, Tsq, Qsq, Cov, Psig_bl, cldfra_bl1D, & & bl_mynn_mixlength, & & edmf_a1,edmf_qc1,bl_mynn_edmf, & & spp_pbl,rstoch_col) INTEGER, INTENT(IN) :: kts,kte INTEGER, INTENT(IN) :: bl_mynn_mixlength,bl_mynn_edmf REAL, INTENT(IN) :: ust, rmo, Psig_bl REAL, DIMENSION(kts:kte), INTENT(in) :: dz REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,cldfra_bl1D,& edmf_a1,edmf_qc1 REAL, DIMENSION(kts:kte), INTENT(out) :: tsq,qsq,cov REAL, DIMENSION(kts:kte), INTENT(inout) :: el,qke REAL, DIMENSION(kts:kte) :: & &ql,pdk,pdt,pdq,pdc,dtl,dqw,dtv,& &gm,gh,sm,sh,qkw,vt,vq INTEGER :: k,l,lmax REAL :: phm,vkz,elq,elv,b1l,b2l,pmz=1.,phh=1.,flt=0.,flq=0.,tmpq REAL :: zi REAL, DIMENSION(kts:kte) :: theta REAL, DIMENSION(kts:kte) :: rstoch_col INTEGER ::spp_pbl DO k = kts,kte ql(k) = 0.0 vt(k) = 0.0 vq(k) = 0.0 END DO CALL mym_level2 ( kts,kte,& & dz, & & u, v, thl, qw, & & ql, vt, vq, & & dtl, dqw, dtv, gm, gh, sm, sh ) el (kts) = 0.0 qke(kts) = ust**2 * ( b1*pmz )**(2.0/3.0) phm = phh*b2 / ( b1*pmz )**(1.0/3.0) tsq(kts) = phm*( flt/ust )**2 qsq(kts) = phm*( flq/ust )**2 cov(kts) = phm*( flt/ust )*( flq/ust ) DO k = kts+1,kte vkz = vk*zw(k) el (k) = vkz/( 1.0 + vkz/100.0 ) qke(k) = 0.0 tsq(k) = 0.0 qsq(k) = 0.0 cov(k) = 0.0 END DO lmax = 5 DO l = 1,lmax CALL mym_length ( & & kts,kte, & & dz, zw, & & rmo, flt, flq, & & vt, vq, & & qke, & & dtv, & & el, & & zi,theta, & & qkw,Psig_bl,cldfra_bl1D,bl_mynn_mixlength,& & edmf_a1,edmf_qc1,bl_mynn_edmf,& & spp_pbl,rstoch_col) DO k = kts+1,kte elq = el(k)*qkw(k) pdk(k) = elq*( sm(k)*gm (k)+& &sh(k)*gh (k) ) pdt(k) = elq* sh(k)*dtl(k)**2 pdq(k) = elq* sh(k)*dqw(k)**2 pdc(k) = elq* sh(k)*dtl(k)*dqw(k) END DO vkz = vk*0.5*dz(kts) elv = 0.5*( el(kts+1)+el(kts) ) / vkz qke(kts) = ust**2 * ( b1*pmz*elv )**(2.0/3.0) phm = phh*b2 / ( b1*pmz/elv**2 )**(1.0/3.0) tsq(kts) = phm*( flt/ust )**2 qsq(kts) = phm*( flq/ust )**2 cov(kts) = phm*( flt/ust )*( flq/ust ) DO k = kts+1,kte-1 b1l = b1*0.25*( el(k+1)+el(k) ) tmpq=MAX(b1l*( pdk(k+1)+pdk(k) ),qkemin) qke(k) = tmpq**(2.0/3.0) IF ( qke(k) .LE. 0.0 ) THEN b2l = 0.0 ELSE b2l = b2*( b1l/b1 ) / SQRT( qke(k) ) END IF tsq(k) = b2l*( pdt(k+1)+pdt(k) ) qsq(k) = b2l*( pdq(k+1)+pdq(k) ) cov(k) = b2l*( pdc(k+1)+pdc(k) ) END DO END DO qke(kte)=qke(kte-1) tsq(kte)=tsq(kte-1) qsq(kte)=qsq(kte-1) cov(kte)=cov(kte-1) END SUBROUTINE mym_initialize SUBROUTINE mym_level2 (kts,kte,& & dz, & & u, v, thl, qw, & & ql, vt, vq, & & dtl, dqw, dtv, gm, gh, sm, sh ) INTEGER, INTENT(IN) :: kts,kte REAL, DIMENSION(kts:kte), INTENT(in) :: dz REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,ql,vt,vq REAL, DIMENSION(kts:kte), INTENT(out) :: & &dtl,dqw,dtv,gm,gh,sm,sh INTEGER :: k REAL :: rfc,f1,f2,rf1,rf2,smc,shc,& &ri1,ri2,ri3,ri4,duz,dtz,dqz,vtt,vqq,dtq,dzk,afk,abk,ri,rf REAL :: a2den rfc = g1/( g1+g2 ) f1 = b1*( g1-c1 ) +3.0*a2*( 1.0 -c2 )*( 1.0-c5 ) & & +2.0*a1*( 3.0-2.0*c2 ) f2 = b1*( g1+g2 ) -3.0*a1*( 1.0 -c2 ) rf1 = b1*( g1-c1 )/f1 rf2 = b1* g1 /f2 smc = a1 /a2* f1/f2 shc = 3.0*a2*( g1+g2 ) ri1 = 0.5/smc ri2 = rf1*smc ri3 = 4.0*rf2*smc -2.0*ri2 ri4 = ri2**2 DO k = kts+1,kte dzk = 0.5 *( dz(k)+dz(k-1) ) afk = dz(k)/( dz(k)+dz(k-1) ) abk = 1.0 -afk duz = ( u(k)-u(k-1) )**2 +( v(k)-v(k-1) )**2 duz = duz /dzk**2 dtz = ( thl(k)-thl(k-1) )/( dzk ) dqz = ( qw(k)-qw(k-1) )/( dzk ) vtt = 1.0 +vt(k)*abk +vt(k-1)*afk vqq = tv0 +vq(k)*abk +vq(k-1)*afk dtq = vtt*dtz +vqq*dqz dtl(k) = dtz dqw(k) = dqz dtv(k) = dtq gm (k) = duz gh (k) = -dtq*gtr ri = -gh(k)/MAX( duz, 1.0e-10 ) IF (CKmod .eq. 1) THEN a2den = 1. + MAX(ri,0.0) ELSE a2den = 1. + 0.0 ENDIF rfc = g1/( g1+g2 ) f1 = b1*( g1-c1 ) +3.0*(a2/a2den)*( 1.0 -c2 )*( 1.0-c5 ) & & +2.0*a1*( 3.0-2.0*c2 ) f2 = b1*( g1+g2 ) -3.0*a1*( 1.0 -c2 ) rf1 = b1*( g1-c1 )/f1 rf2 = b1* g1 /f2 smc = a1 /(a2/a2den)* f1/f2 shc = 3.0*(a2/a2den)*( g1+g2 ) ri1 = 0.5/smc ri2 = rf1*smc ri3 = 4.0*rf2*smc -2.0*ri2 ri4 = ri2**2 rf = MIN( ri1*( ri+ri2-SQRT(ri**2-ri3*ri+ri4) ), rfc ) sh (k) = shc*( rfc-rf )/( 1.0-rf ) sm (k) = smc*( rf1-rf )/( rf2-rf ) * sh(k) END DO RETURN END SUBROUTINE mym_level2 SUBROUTINE mym_length ( & & kts,kte, & & dz, zw, & & rmo, flt, flq, & & vt, vq, & & qke, & & dtv, & & el, & & zi,theta, & & qkw,Psig_bl,cldfra_bl1D,bl_mynn_mixlength,& & edmf_a1,edmf_qc1,bl_mynn_edmf,& & spp_pbl,rstoch_col) INTEGER, INTENT(IN) :: kts,kte INTEGER, INTENT(IN) :: bl_mynn_mixlength,bl_mynn_edmf REAL, DIMENSION(kts:kte), INTENT(in) :: dz REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw REAL, INTENT(in) :: rmo,flt,flq,Psig_bl REAL, DIMENSION(kts:kte), INTENT(IN) :: qke,vt,vq,cldfra_bl1D,& edmf_a1,edmf_qc1 REAL, DIMENSION(kts:kte), INTENT(out) :: qkw, el REAL, DIMENSION(kts:kte), INTENT(in) :: dtv REAL :: elt,vsc REAL, DIMENSION(kts:kte), INTENT(IN) :: theta REAL, DIMENSION(kts:kte) :: qtke,elBLmin,elBLavg,thetaw REAL :: wt,wt2,zi,zi2,h1,h2,hs,elBLmin0,elBLavg0 REAL :: cns, & alp1, & alp2, & alp3, & alp4, & alp5 REAL, PARAMETER :: minzi = 300. REAL, PARAMETER :: maxdz = 750. REAL, PARAMETER :: mindz = 300. REAL, PARAMETER :: ZSLH = 100. REAL, PARAMETER :: CSL = 2. REAL :: z_m INTEGER :: i,j,k REAL :: afk,abk,zwk,zwk1,dzk,qdz,vflx,bv,tau_cloud,elb,els,els1,elf, & & el_stab,el_unstab,el_mf,el_stab_mf,elb_mf,PBLH_PLUS_ENT INTEGER, INTENT(IN) :: spp_pbl REAL, DIMENSION(kts:kte), INTENT(in) :: rstoch_col IF ( bl_mynn_mixlength .EQ. 0 ) THEN cns = 2.7 alp1 = 0.23 alp2 = 1.0 alp3 = 5.0 alp4 = 100. alp5 = 0.4 ELSEIF ( bl_mynn_mixlength .EQ. 1 ) THEN cns = 2.3 alp1 = 0.23 alp2 = 0.6 alp3 = 3.0 alp4 = 20. alp5 = 0.4 ELSEIF ( bl_mynn_mixlength .GE. 2 ) THEN cns = 3.5 alp1 = 0.23 alp2 = 0.25 alp3 = 3.0 alp4 = 10. alp5 = 0.3 ENDIF IF ( bl_mynn_mixlength .EQ. 0 ) THEN zi2=10000. ELSE zi2=MAX(zi,minzi) ENDIF h1=MAX(0.3*zi2,mindz) h1=MIN(h1,maxdz) h2=h1/2.0 qtke(kts)=MAX(qke(kts)/2.,0.01) thetaw(kts)=theta(kts) qkw(kts) = SQRT(MAX(qke(kts),1.0e-10)) DO k = kts+1,kte afk = dz(k)/( dz(k)+dz(k-1) ) abk = 1.0 -afk qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk,1.0e-3)) qtke(k) = (qkw(k)**2.)/2. thetaw(k)= theta(k)*abk + theta(k-1)*afk END DO elt = 1.0e-5 vsc = 1.0e-5 IF ( bl_mynn_mixlength .EQ. 2 ) THEN PBLH_PLUS_ENT = MAX(zi, 100.) ELSE PBLH_PLUS_ENT = zi2+h1 ENDIF k = kts+1 zwk = zw(k) DO WHILE (zwk .LE. PBLH_PLUS_ENT) dzk = 0.5*( dz(k)+dz(k-1) ) qdz = MAX( qkw(k)-qmin, 0.03 )*dzk elt = elt +qdz*zwk vsc = vsc +qdz k = k+1 zwk = zw(k) END DO elt = MAX(alp1*elt/vsc, 10.) vflx = ( vt(kts)+1.0 )*flt +( vq(kts)+tv0 )*flq vsc = ( gtr*elt*MAX( vflx, 0.0 ) )**(1.0/3.0) el(kts) = 0.0 IF ( bl_mynn_mixlength .EQ. 1 ) THEN CALL boulac_length(kts,kte,zw,dz,qtke,thetaw,elBLmin,elBLavg) ENDIF DO k = kts+1,kte zwk = zw(k) IF (k .EQ. kts+1) zwk1=zwk IF ( dtv(k) .GT. 0.0 ) THEN bv = SQRT( gtr*dtv(k) ) IF ( bl_mynn_mixlength .EQ. 0 ) THEN elb = alp2*qkw(k) / bv & & *( 1.0 + alp3/alp2*& &SQRT( vsc/( bv*elt ) ) ) elf = alp2 * qkw(k)/bv ELSE IF ( bl_mynn_mixlength .EQ. 1 ) THEN elb = alp2*qkw(k) / bv & & *( 1.0 + alp3/alp2*& &SQRT( vsc/( bv*elt ) ) ) elb = MIN(elb, zwk) elf = alp2 * qkw(k)/bv ELSE IF ( bl_mynn_mixlength .GE. 2 ) THEN elb_mf = alp2*qkw(k) / bv & & *( 1.0 + alp3/alp2*& &SQRT( vsc/( bv*elt ) ) ) elb = MIN(alp2*qkw(k)/bv, zwk) elf = elb IF (zwk > zi .AND. elf > 500.) THEN CALL boulac_length0(k,kts,kte,zw,dz,qtke,thetaw,elBLmin0,elBLavg0) elf = alp5*elBLavg0 ENDIF END IF ELSE IF ( bl_mynn_mixlength .LE. 1 ) THEN elb = 1.0e10 elf = elb ELSE IF ( bl_mynn_mixlength .GE. 2 ) THEN tau_cloud = MIN(MAX(0.5*zi/((gtr*zi*MAX(flt,1.0e-4))**(1.0/3.0)),10.),100.) wt=.5*TANH((zwk - (zi2+h1))/h2) + .5 tau_cloud = tau_cloud*(1.-wt) + 50.*wt elb = MIN(tau_cloud*SQRT(MIN(qtke(k),20.)), zwk) elf = elb elb_mf = elb IF (zwk > zi .AND. elf > 500.) THEN CALL boulac_length0(k,kts,kte,zw,dz,qtke,thetaw,elBLmin0,elBLavg0) elf = alp5*elBLavg0*(1.-cldfra_bl1D(k)) + elf*cldfra_bl1D(k) END IF elf = elf*(1.-cldfra_bl1D(k)) + elb*cldfra_bl1D(k) END IF END IF z_m = MAX(0.,zwk - 4.) IF ( rmo .GT. 0.0 ) THEN els = vk*zwk/(1.0+cns*MIN( zwk*rmo, zmax )) els1 = vk*z_m/(1.0+cns*MIN( zwk*rmo, zmax )) ELSE els = vk*zwk*( 1.0 - alp4* zwk*rmo )**0.2 els1 = vk*z_m*( 1.0 - alp4* zwk*rmo )**0.2 END IF wt=.5*TANH((zwk - (zi2+h1))/h2) + .5 IF ( bl_mynn_mixlength .EQ. 0 ) THEN el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf) ELSE IF ( bl_mynn_mixlength .EQ. 1 ) THEN el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf) el(k) = el(k)*(1.-wt) + alp5*elBLmin(k)*wt ELSE IF ( bl_mynn_mixlength .GE. 2 ) THEN el_unstab = els/(1. + (els1/elt)) el_stab = MIN(el_unstab, elb) el_stab_mf = MIN(el_unstab, elb_mf) IF (bl_mynn_edmf > 0 .AND. edmf_a1(kts)>0.0) THEN el(k) = el_stab_mf ELSE el(k) = el_stab ENDIF el(k) = el(k)*(1.-wt) + elf*wt END IF IF ( bl_mynn_mixlength .GE. 1 ) THEN el(k) = el(k)*Psig_bl END IF if (spp_pbl==1) then if (k.lt.25) then el(k)= el(k) + el(k)* rstoch_col(k) * 1.5 * MAX(exp(-MAX(zwk-3000.,0.0)/2000.),0.01) endif endif IF ( wrf_at_debug_level(3000) ) THEN IF (el(k) > 1000.) THEN WRITE ( mynn_message , FMT='(A,F7.0,I5,4F7.0)' ) & ' MYNN; mym_length; LARGE el,k,elb,elt,elf,tau:'& , el(k),k,elb,elt,elf,tau_cloud CALL wrf_debug ( 0 , mynn_message ) ENDIF ENDIF END DO RETURN END SUBROUTINE mym_length SUBROUTINE boulac_length0(k,kts,kte,zw,dz,qtke,theta,lb1,lb2) INTEGER, INTENT(IN) :: k,kts,kte REAL, DIMENSION(kts:kte), INTENT(IN) :: qtke,dz,theta REAL, INTENT(OUT) :: lb1,lb2 REAL, DIMENSION(kts:kte+1), INTENT(IN) :: zw INTEGER :: izz, found REAL :: dlu,dld REAL :: dzt, zup, beta, zup_inf, bbb, tl, zdo, zdo_sup, zzz zup=0. dlu=zw(kte+1)-zw(k)-dz(k)/2. zzz=0. zup_inf=0. beta=g/theta(k) if (k .lt. kte) then found = 0 izz=k DO WHILE (found .EQ. 0) if (izz .lt. kte) then dzt=dz(izz) zup=zup-beta*theta(k)*dzt zup=zup+beta*(theta(izz+1)+theta(izz))*dzt/2. zzz=zzz+dzt if (qtke(k).lt.zup .and. qtke(k).ge.zup_inf) then bbb=(theta(izz+1)-theta(izz))/dzt if (bbb .ne. 0.) then tl=(-beta*(theta(izz)-theta(k)) + & & sqrt( max(0.,(beta*(theta(izz)-theta(k)))**2. + & & 2.*bbb*beta*(qtke(k)-zup_inf))))/bbb/beta else if (theta(izz) .ne. theta(k))then tl=(qtke(k)-zup_inf)/(beta*(theta(izz)-theta(k))) else tl=0. endif endif dlu=zzz-dzt+tl found =1 endif zup_inf=zup izz=izz+1 ELSE found = 1 ENDIF ENDDO endif zdo=0. zdo_sup=0. dld=zw(k) zzz=0. if (k .gt. kts) then found = 0 izz=k DO WHILE (found .EQ. 0) if (izz .gt. kts) then dzt=dz(izz-1) zdo=zdo+beta*theta(k)*dzt zdo=zdo-beta*(theta(izz-1)+theta(izz))*dzt/2. zzz=zzz+dzt if (qtke(k).lt.zdo .and. qtke(k).ge.zdo_sup) then bbb=(theta(izz)-theta(izz-1))/dzt if (bbb .ne. 0.) then tl=(beta*(theta(izz)-theta(k))+ & & sqrt( max(0.,(beta*(theta(izz)-theta(k)))**2. + & & 2.*bbb*beta*(qtke(k)-zdo_sup))))/bbb/beta else if (theta(izz) .ne. theta(k)) then tl=(qtke(k)-zdo_sup)/(beta*(theta(izz)-theta(k))) else tl=0. endif endif dld=zzz-dzt+tl found = 1 endif zdo_sup=zdo izz=izz-1 ELSE found = 1 ENDIF ENDDO endif dld = min(dld,zw(k+1)) lb1 = min(dlu,dld) dlu=MAX(0.1,MIN(dlu,1000.)) dld=MAX(0.1,MIN(dld,1000.)) lb2 = sqrt(dlu*dld) if (k .eq. kte) then lb1 = 0. lb2 = 0. endif END SUBROUTINE boulac_length0 SUBROUTINE boulac_length(kts,kte,zw,dz,qtke,theta,lb1,lb2) INTEGER, INTENT(IN) :: kts,kte REAL, DIMENSION(kts:kte), INTENT(IN) :: qtke,dz,theta REAL, DIMENSION(kts:kte), INTENT(OUT) :: lb1,lb2 REAL, DIMENSION(kts:kte+1), INTENT(IN) :: zw INTEGER :: iz, izz, found REAL, DIMENSION(kts:kte) :: dlu,dld REAL, PARAMETER :: Lmax=2000. REAL :: dzt, zup, beta, zup_inf, bbb, tl, zdo, zdo_sup, zzz do iz=kts,kte zup=0. dlu(iz)=zw(kte+1)-zw(iz)-dz(iz)/2. zzz=0. zup_inf=0. beta=g/theta(iz) if (iz .lt. kte) then found = 0 izz=iz DO WHILE (found .EQ. 0) if (izz .lt. kte) then dzt=dz(izz) zup=zup-beta*theta(iz)*dzt zup=zup+beta*(theta(izz+1)+theta(izz))*dzt/2. zzz=zzz+dzt if (qtke(iz).lt.zup .and. qtke(iz).ge.zup_inf) then bbb=(theta(izz+1)-theta(izz))/dzt if (bbb .ne. 0.) then tl=(-beta*(theta(izz)-theta(iz)) + & & sqrt( max(0.,(beta*(theta(izz)-theta(iz)))**2. + & & 2.*bbb*beta*(qtke(iz)-zup_inf))))/bbb/beta else if (theta(izz) .ne. theta(iz))then tl=(qtke(iz)-zup_inf)/(beta*(theta(izz)-theta(iz))) else tl=0. endif endif dlu(iz)=zzz-dzt+tl found =1 endif zup_inf=zup izz=izz+1 ELSE found = 1 ENDIF ENDDO endif zdo=0. zdo_sup=0. dld(iz)=zw(iz) zzz=0. if (iz .gt. kts) then found = 0 izz=iz DO WHILE (found .EQ. 0) if (izz .gt. kts) then dzt=dz(izz-1) zdo=zdo+beta*theta(iz)*dzt zdo=zdo-beta*(theta(izz-1)+theta(izz))*dzt/2. zzz=zzz+dzt if (qtke(iz).lt.zdo .and. qtke(iz).ge.zdo_sup) then bbb=(theta(izz)-theta(izz-1))/dzt if (bbb .ne. 0.) then tl=(beta*(theta(izz)-theta(iz))+ & & sqrt( max(0.,(beta*(theta(izz)-theta(iz)))**2. + & & 2.*bbb*beta*(qtke(iz)-zdo_sup))))/bbb/beta else if (theta(izz) .ne. theta(iz)) then tl=(qtke(iz)-zdo_sup)/(beta*(theta(izz)-theta(iz))) else tl=0. endif endif dld(iz)=zzz-dzt+tl found = 1 endif zdo_sup=zdo izz=izz-1 ELSE found = 1 ENDIF ENDDO endif dld(iz) = min(dld(iz),zw(iz+1)) lb1(iz) = min(dlu(iz),dld(iz)) dlu(iz)=MAX(0.1,MIN(dlu(iz),1000.)) dld(iz)=MAX(0.1,MIN(dld(iz),1000.)) lb2(iz) = sqrt(dlu(iz)*dld(iz)) lb1(iz) = lb1(iz)/(1. + (lb1(iz)/Lmax)) lb2(iz) = lb2(iz)/(1. + (lb2(iz)/Lmax)) if (iz .eq. kte) then lb1(kte) = lb1(kte-1) lb2(kte) = lb2(kte-1) endif ENDDO END SUBROUTINE boulac_length SUBROUTINE mym_turbulence ( & & kts,kte, & & levflag, & & dz, zw, & & u, v, thl, ql, qw, & & qke, tsq, qsq, cov, & & vt, vq, & & rmo, flt, flq, & & zi,theta, & & sh, & & El, & & Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, & & qWT1D,qSHEAR1D,qBUOY1D,qDISS1D, & & bl_mynn_tkebudget, & & Psig_bl,Psig_shcu,cldfra_bl1D,bl_mynn_mixlength,& & edmf_a1,edmf_qc1,bl_mynn_edmf,& & spp_pbl,rstoch_col) INTEGER, INTENT(IN) :: kts,kte INTEGER, INTENT(IN) :: levflag,bl_mynn_mixlength,bl_mynn_edmf REAL, DIMENSION(kts:kte), INTENT(in) :: dz REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw REAL, INTENT(in) :: rmo,flt,flq,Psig_bl,Psig_shcu REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,& &ql,vt,vq,qke,tsq,qsq,cov,cldfra_bl1D,edmf_a1,edmf_qc1 REAL, DIMENSION(kts:kte), INTENT(out) :: dfm,dfh,dfq,& &pdk,pdt,pdq,pdc,tcd,qcd,el REAL, DIMENSION(kts:kte), INTENT(inout) :: & qWT1D,qSHEAR1D,qBUOY1D,qDISS1D REAL :: q3sq_old,dlsq1,qWTP_old,qWTP_new REAL :: dudz,dvdz,dTdz,& upwp,vpwp,Tpwp INTEGER, INTENT(in) :: bl_mynn_tkebudget REAL, DIMENSION(kts:kte) :: qkw,dtl,dqw,dtv,gm,gh,sm,sh INTEGER :: k REAL :: e6c,dzk,afk,abk,vtt,vqq,& &cw25,clow,cupp,gamt,gamq,smd,gamv,elq,elh REAL :: zi REAL, DIMENSION(kts:kte), INTENT(in) :: theta REAL :: a2den, duz, ri, HLmod REAL:: auh,aum,adh,adm,aeh,aem,Req,Rsl,Rsl2 DOUBLE PRECISION q2sq, t2sq, r2sq, c2sq, elsq, gmel, ghel DOUBLE PRECISION q3sq, t3sq, r3sq, c3sq, dlsq, qdiv DOUBLE PRECISION e1, e2, e3, e4, enum, eden, wden INTEGER, INTENT(IN) :: spp_pbl REAL, DIMENSION(KTS:KTE) :: rstoch_col REAL :: prlimit CALL mym_level2 (kts,kte,& & dz, & & u, v, thl, qw, & & ql, vt, vq, & & dtl, dqw, dtv, gm, gh, sm, sh ) CALL mym_length ( & & kts,kte, & & dz, zw, & & rmo, flt, flq, & & vt, vq, & & qke, & & dtv, & & el, & & zi,theta, & & qkw,Psig_bl,cldfra_bl1D,bl_mynn_mixlength, & & edmf_a1,edmf_qc1,bl_mynn_edmf, & & spp_pbl,rstoch_col) DO k = kts+1,kte dzk = 0.5 *( dz(k)+dz(k-1) ) afk = dz(k)/( dz(k)+dz(k-1) ) abk = 1.0 -afk elsq = el (k)**2 q2sq = b1*elsq*( sm(k)*gm(k)+sh(k)*gh(k) ) q3sq = qkw(k)**2 duz = ( u(k)-u(k-1) )**2 +( v(k)-v(k-1) )**2 duz = duz /dzk**2 ri = -gh(k)/MAX( duz, 1.0e-10 ) IF (CKmod .eq. 1) THEN a2den = 1. + MAX(ri,0.0) ELSE a2den = 1. + 0.0 ENDIF gmel = gm (k)*elsq ghel = gh (k)*elsq IF ( wrf_at_debug_level(3000) ) THEN IF (sh(k)<0.0 .OR. sm(k)<0.0) THEN WRITE ( mynn_message , FMT='(A,F8.1,A,I6)' ) & " MYNN; mym_turbulence2.0; sh=",sh(k)," k=",k CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " gm=",gm(k)," gh=",gh(k)," sm=",sm(k) CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " q2sq=",q2sq," q3sq=",q3sq," q3/q2=",q3sq/q2sq CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " qke=",qke(k)," el=",el(k)," ri=",ri CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " PBLH=",zi," u=",u(k)," v=",v(k) CALL wrf_debug ( 0 , mynn_message ) ENDIF ENDIF IF (CKmod .eq. 1) THEN HLmod = q2sq -1. ELSE HLmod = q3sq ENDIF dlsq = elsq IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k) IF ( q3sq .LT. q2sq ) THEN qdiv = SQRT( q3sq/q2sq ) sm(k) = sm(k) * qdiv sh(k) = sh(k) * qdiv e1 = q3sq - e1c*ghel/a2den * qdiv**2 e2 = q3sq - e2c*ghel/a2den * qdiv**2 e3 = e1 + e3c*ghel/(a2den**2) * qdiv**2 e4 = e1 - e4c*ghel/a2den * qdiv**2 eden = e2*e4 + e3*e5c*gmel * qdiv**2 eden = MAX( eden, 1.0d-20 ) ELSE e1 = q3sq - e1c*ghel/a2den e2 = q3sq - e2c*ghel/a2den e3 = e1 + e3c*ghel/(a2den**2) e4 = e1 - e4c*ghel/a2den eden = e2*e4 + e3*e5c*gmel eden = MAX( eden, 1.0d-20 ) qdiv = 1.0 sm(k) = q3sq*a1*( e3-3.0*c1*e4 )/eden sh(k) = q3sq*(a2/a2den)*( e2+3.0*c1*e5c*gmel )/eden END IF IF ( wrf_at_debug_level(3000) ) THEN IF (sh(k)<0.0 .OR. sm(k)<0.0 .OR. & sh(k) > 0.76*b2 .or. (sm(k)**2*gm(k) .gt. .44**2)) THEN WRITE ( mynn_message , FMT='(A,F8.1,A,I6)' ) & " MYNN; mym_turbulence2.5; sh=",sh(k)," k=",k CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " gm=",gm(k)," gh=",gh(k)," sm=",sm(k) CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " q2sq=",q2sq," q3sq=",q3sq," q3/q2=",q3sq/q2sq CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " qke=",qke(k)," el=",el(k)," ri=",ri CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " PBLH=",zi," u=",u(k)," v=",v(k) CALL wrf_debug ( 0 , mynn_message ) ENDIF ENDIF IF ( levflag .EQ. 3 ) THEN t2sq = qdiv*b2*elsq*sh(k)*dtl(k)**2 r2sq = qdiv*b2*elsq*sh(k)*dqw(k)**2 c2sq = qdiv*b2*elsq*sh(k)*dtl(k)*dqw(k) t3sq = MAX( tsq(k)*abk+tsq(k-1)*afk, 0.0 ) r3sq = MAX( qsq(k)*abk+qsq(k-1)*afk, 0.0 ) c3sq = cov(k)*abk+cov(k-1)*afk c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq ) vtt = 1.0 +vt(k)*abk +vt(k-1)*afk vqq = tv0 +vq(k)*abk +vq(k-1)*afk t2sq = vtt*t2sq +vqq*c2sq r2sq = vtt*c2sq +vqq*r2sq c2sq = MAX( vtt*t2sq+vqq*r2sq, 0.0d0 ) t3sq = vtt*t3sq +vqq*c3sq r3sq = vtt*c3sq +vqq*r3sq c3sq = MAX( vtt*t3sq+vqq*r3sq, 0.0d0 ) cw25 = e1*( e2 + 3.0*c1*e5c*gmel*qdiv**2 )/( 3.0*eden ) dlsq = elsq IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k) auh = 27.*a1*((a2/a2den)**2)*b2*(g/tref)**2 aum = 54.*(a1**2)*(a2/a2den)*b2*c1*(g/tref) adh = 9.*a1*((a2/a2den)**2)*(12.*a1 + 3.*b2)*(g/tref)**2 adm = 18.*(a1**2)*(a2/a2den)*(b2 - 3.*(a2/a2den))*(g/tref) aeh = (9.*a1*((a2/a2den)**2)*b1 +9.*a1*((a2/a2den)**2)* & (12.*a1 + 3.*b2))*(g/tref) aem = 3.*a1*(a2/a2den)*b1*(3.*(a2/a2den) + 3.*b2*c1 + & (18.*a1*c1 - b2)) + & (18.)*(a1**2)*(a2/a2den)*(b2 - 3.*(a2/a2den)) Req = -aeh/aem Rsl = (auh + aum*Req)/(3.*adh + 3.*adm*Req) Rsl = .12 Rsl2= 1.0 - 2.*Rsl e2 = q3sq - e2c*ghel/a2den * qdiv**2 e3 = q3sq + e3c*ghel/(a2den**2) * qdiv**2 e4 = q3sq - e4c*ghel/a2den * qdiv**2 eden = e2*e4 + e3 *e5c*gmel * qdiv**2 wden = cc3*gtr**2 * dlsq**2/elsq * qdiv**2 & & *( e2*e4c/a2den - e3c*e5c*gmel/(a2den**2) * qdiv**2 ) IF ( wden .NE. 0.0 ) THEN clow = q3sq*( Rsl -cw25 )*eden/wden cupp = q3sq*( Rsl2-cw25 )*eden/wden IF ( wden .GT. 0.0 ) THEN c3sq = MIN( MAX( c3sq, c2sq+clow ), c2sq+cupp ) ELSE c3sq = MAX( MIN( c3sq, c2sq+clow ), c2sq+cupp ) END IF END IF e1 = e2 + e5c*gmel * qdiv**2 eden = MAX( eden, 1.0d-20 ) e6c = 3.0*(a2/a2den)*cc3*gtr * dlsq/elsq IF ( t2sq .GE. 0.0 ) THEN enum = MAX( qdiv*e6c*( t3sq-t2sq ), 0.0d0 ) ELSE enum = MIN( qdiv*e6c*( t3sq-t2sq ), 0.0d0 ) ENDIF gamt =-e1 *enum /eden IF ( r2sq .GE. 0.0 ) THEN enum = MAX( qdiv*e6c*( r3sq-r2sq ), 0.0d0 ) ELSE enum = MIN( qdiv*e6c*( r3sq-r2sq ), 0.0d0 ) ENDIF gamq =-e1 *enum /eden enum = MAX( qdiv*e6c*( c3sq-c2sq ), 0.0d0) smd = dlsq*enum*gtr/eden * qdiv**2 * (e3c/(a2den**2) + & & e4c/a2den)*a1/(a2/a2den) gamv = e1 *enum*gtr/eden sm(k) = sm(k) +smd qdiv = 1.0 IF ( wrf_at_debug_level(3000) ) THEN IF (sh(k)<-0.3 .OR. sm(k)<-0.3 .OR. & qke(k) < -0.1 .or. ABS(smd) .gt. 2.0) THEN WRITE ( mynn_message , FMT='(A,F8.1,A,I6)' ) & " MYNN; mym_turbulence3.0; sh=",sh(k)," k=",k CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " gm=",gm(k)," gh=",gh(k)," sm=",sm(k) CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " q2sq=",q2sq," q3sq=",q3sq," q3/q2=",q3sq/q2sq CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " qke=",qke(k)," el=",el(k)," ri=",ri CALL wrf_debug ( 0 , mynn_message ) WRITE ( mynn_message , FMT='(A,F6.1,A,F6.1,A,F6.1)' ) & " PBLH=",zi," u=",u(k)," v=",v(k) CALL wrf_debug ( 0 , mynn_message ) ENDIF ENDIF ELSE gamt = 0.0 gamq = 0.0 gamv = 0.0 END IF if (spp_pbl==1) then prlimit = MIN(MAX(1.,2.5 + 5.0*rstoch_col(k)), 10.) IF(sm(k) > sh(k)*Prlimit) THEN sm(k) = sh(k)*Prlimit ENDIF ENDIF elq = el(k)*qkw(k) elh = elq*qdiv pdk(k) = elq*( sm(k)*gm(k) & & +sh(k)*gh(k)+gamv ) pdt(k) = elh*( sh(k)*dtl(k)+gamt )*dtl(k) pdq(k) = elh*( sh(k)*dqw(k)+gamq )*dqw(k) pdc(k) = elh*( sh(k)*dtl(k)+gamt )& &*dqw(k)*0.5 & &+elh*( sh(k)*dqw(k)+gamq )*dtl(k)*0.5 tcd(k) = elq*gamt qcd(k) = elq*gamq dfm(k) = elq*sm(k) / dzk dfh(k) = elq*sh(k) / dzk dfq(k) = dfm(k) IF ( bl_mynn_tkebudget == 1) THEN dudz = ( u(k)-u(k-1) )/dzk dvdz = ( v(k)-v(k-1) )/dzk dTdz = ( thl(k)-thl(k-1) )/dzk upwp = -elq*sm(k)*dudz vpwp = -elq*sm(k)*dvdz Tpwp = -elq*sh(k)*dTdz Tpwp = SIGN(MAX(ABS(Tpwp),1.E-6),Tpwp) IF ( k .EQ. kts+1 ) THEN qWT1D(kts)=0. q3sq_old =0. qWTP_old =0. dlsq1 = MAX(el(kts)**2,1.0) IF ( q3sq_old/dlsq1 .LT. -gh(k) ) q3sq_old = -dlsq1*gh(k) ENDIF qWTP_new = elq*Sqfac*sm(k)*(q3sq - q3sq_old)/dzk qWT1D(k) = 0.5*(qWTP_new - qWTP_old)/dzk qWTP_old = elq*Sqfac*sm(k)*(q3sq - q3sq_old)/dzk q3sq_old = q3sq qSHEAR1D(k) = elq*sm(k)*gm(k) qBUOY1D(k) = elq*(sh(k)*(-dTdz*g/thl(k)) + gamv) qDISS1D(k) = (q3sq**(3./2.))/(b1*MAX(el(k),1.)) ENDIF END DO dfm(kts) = 0.0 dfh(kts) = 0.0 dfq(kts) = 0.0 tcd(kts) = 0.0 qcd(kts) = 0.0 tcd(kte) = 0.0 qcd(kte) = 0.0 DO k = kts,kte-1 dzk = dz(k) tcd(k) = ( tcd(k+1)-tcd(k) )/( dzk ) qcd(k) = ( qcd(k+1)-qcd(k) )/( dzk ) END DO IF ( bl_mynn_tkebudget == 1) THEN qWT1D(kts)=0. qSHEAR1D(kts)=qSHEAR1D(kts+1) qBUOY1D(kts)=qBUOY1D(kts+1) qDISS1D(kts)=qDISS1D(kts+1) ENDIF RETURN END SUBROUTINE mym_turbulence SUBROUTINE mym_predict (kts,kte,& & levflag, & & delt,& & dz, & & ust, flt, flq, pmz, phh, & & el, dfq, & & pdk, pdt, pdq, pdc,& & qke, tsq, qsq, cov, & & s_aw,s_awqke,bl_mynn_edmf_tke & &) INTEGER, INTENT(IN) :: kts,kte INTEGER, INTENT(IN) :: levflag INTEGER, INTENT(IN) :: bl_mynn_edmf_tke REAL, INTENT(IN) :: delt REAL, DIMENSION(kts:kte), INTENT(IN) :: dz, dfq,el REAL, DIMENSION(kts:kte), INTENT(INOUT) :: pdk, pdt, pdq, pdc REAL, INTENT(IN) :: flt, flq, ust, pmz, phh REAL, DIMENSION(kts:kte), INTENT(INOUT) :: qke,tsq, qsq, cov REAL, DIMENSION(kts:kte+1), INTENT(INOUT) :: s_awqke,s_aw INTEGER :: k,nz REAL, DIMENSION(kts:kte) :: qkw, bp, rp, df3q REAL :: vkz,pdk1,phm,pdt1,pdq1,pdc1,b1l,b2l,onoff REAL, DIMENSION(kts:kte) :: dtz REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d nz=kte-kts+1 IF (bl_mynn_edmf_tke == 0) THEN onoff=0.0 ELSE onoff=1.0 ENDIF vkz = vk*0.5*dz(kts) DO k = kts,kte qkw(k) = SQRT( MAX( qke(k), 0.0 ) ) df3q(k)=Sqfac*dfq(k) dtz(k)=delt/dz(k) END DO pdk1 = 2.0*ust**3*pmz/( vkz ) phm = 2.0/ust *phh/( vkz ) pdt1 = phm*flt**2 pdq1 = phm*flq**2 pdc1 = phm*flt*flq pdk(kts) = pdk1 -pdk(kts+1) pdt(kts) = pdt(kts+1) pdq(kts) = pdq(kts+1) pdc(kts) = pdc(kts+1) DO k = kts,kte-1 b1l = b1*0.5*( el(k+1)+el(k) ) bp(k) = 2.*qkw(k) / b1l rp(k) = pdk(k+1) + pdk(k) END DO DO k=kts,kte-1 a(k-kts+1)=-dtz(k)*df3q(k) + 0.5*dtz(k)*s_aw(k)*onoff b(k-kts+1)=1. + dtz(k)*(df3q(k)+df3q(k+1)) & + 0.5*dtz(k)*(s_aw(k)-s_aw(k+1))*onoff + bp(k)*delt c(k-kts+1)=-dtz(k)*df3q(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff d(k-kts+1)=rp(k)*delt + qke(k) + dtz(k)*(s_awqke(k)-s_awqke(k+1))*onoff ENDDO a(nz)=-1. b(nz)=1. c(nz)=0. d(nz)=0. CALL tridiag(nz,a,b,c,d) DO k=kts,kte qke(k)=d(k-kts+1) ENDDO IF ( levflag .EQ. 3 ) THEN DO k = kts,kte-1 b2l = b2*0.5*( el(k+1)+el(k) ) bp(k) = 2.*qkw(k) / b2l rp(k) = pdt(k+1) + pdt(k) END DO DO k=kts,kte-1 a(k-kts+1)=-dtz(k)*dfq(k) b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt c(k-kts+1)=-dtz(k)*dfq(k+1) d(k-kts+1)=rp(k)*delt + tsq(k) ENDDO a(nz)=-1. b(nz)=1. c(nz)=0. d(nz)=0. CALL tridiag(nz,a,b,c,d) DO k=kts,kte tsq(k)=d(k-kts+1) ENDDO DO k = kts,kte-1 b2l = b2*0.5*( el(k+1)+el(k) ) bp(k) = 2.*qkw(k) / b2l rp(k) = pdq(k+1) +pdq(k) END DO DO k=kts,kte-1 a(k-kts+1)=-dtz(k)*dfq(k) b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt c(k-kts+1)=-dtz(k)*dfq(k+1) d(k-kts+1)=rp(k)*delt + qsq(k) ENDDO a(nz)=-1. b(nz)=1. c(nz)=0. d(nz)=0. CALL tridiag(nz,a,b,c,d) DO k=kts,kte qsq(k)=d(k-kts+1) ENDDO DO k = kts,kte-1 b2l = b2*0.5*( el(k+1)+el(k) ) bp(k) = 2.*qkw(k) / b2l rp(k) = pdc(k+1) + pdc(k) END DO DO k=kts,kte-1 a(k-kts+1)=-dtz(k)*dfq(k) b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt c(k-kts+1)=-dtz(k)*dfq(k+1) d(k-kts+1)=rp(k)*delt + cov(k) ENDDO a(nz)=-1. b(nz)=1. c(nz)=0. d(nz)=0. CALL tridiag(nz,a,b,c,d) DO k=kts,kte cov(k)=d(k-kts+1) ENDDO ELSE DO k = kts,kte-1 IF ( qkw(k) .LE. 0.0 ) THEN b2l = 0.0 ELSE b2l = b2*0.25*( el(k+1)+el(k) )/qkw(k) END IF tsq(k) = b2l*( pdt(k+1)+pdt(k) ) qsq(k) = b2l*( pdq(k+1)+pdq(k) ) cov(k) = b2l*( pdc(k+1)+pdc(k) ) END DO tsq(kte)=tsq(kte-1) qsq(kte)=qsq(kte-1) cov(kte)=cov(kte-1) END IF END SUBROUTINE mym_predict SUBROUTINE mym_condensation (kts,kte, & & dx, dz, & & thl, qw, & & p,exner, & & tsq, qsq, cov, & & Sh, el, bl_mynn_cloudpdf,& & qc_bl1D, cldfra_bl1D, & & PBLH1,HFX1, & & Vt, Vq, th, sgm) INTEGER, INTENT(IN) :: kts,kte, bl_mynn_cloudpdf REAL, INTENT(IN) :: dx,PBLH1,HFX1 REAL, DIMENSION(kts:kte), INTENT(IN) :: dz REAL, DIMENSION(kts:kte), INTENT(IN) :: p,exner, thl, qw, & &tsq, qsq, cov, th REAL, DIMENSION(kts:kte), INTENT(INOUT) :: vt,vq,sgm REAL, DIMENSION(kts:kte) :: qmq,alp,a,bet,b,ql,q1,cld,RH REAL, DIMENSION(kts:kte), INTENT(OUT) :: qc_bl1D,cldfra_bl1D DOUBLE PRECISION :: t3sq, r3sq, c3sq REAL :: p2a,qsl,esat,qsat,tlk,qsat_tl,dqsl,cld0,q1k,eq1,qll,& &q2p,pt,rac,qt,t,xl,rsl,cpm,cdhdz,Fng,qww,alpha,beta,bb,ls_min,ls,wt INTEGER :: i,j,k REAL :: erf REAL::dth,dtl,dqw,dzk REAL, DIMENSION(kts:kte), INTENT(IN) :: Sh,el REAL::zagl,cld9,damp,edown,RHcrit,RHmean,RHsum,RHnum,Hshcu,PBLH2,ql_limit REAL, PARAMETER :: Hfac = 3.0 REAL, PARAMETER :: HFXmin = 50.0 REAL :: RH_00L, RH_00O, phi_dz, lfac REAL, PARAMETER :: cdz = 2.0 REAL, PARAMETER :: mdz = 1.5 REAL :: theta1, theta2, ht1, ht2 INTEGER :: k_tropo DO k = kte-3, kts, -1 theta1 = th(k) theta2 = th(k+2) ht1 = 44307.692 * (1.0 - (p(k)/101325.)**0.190) ht2 = 44307.692 * (1.0 - (p(k+2)/101325.)**0.190) if ( (((theta2-theta1)/(ht2-ht1)) .lt. 10./1500. ) .AND. & & (ht1.lt.19000.) .and. (ht1.gt.4000.) ) then goto 86 endif ENDDO 86 continue k_tropo = MAX(kts+2, k+2) zagl = 0. DO k = kts,kte-1 t = th(k)*exner(k) esat = esat_blend(t) qsl=ep_2*esat/(p(k)-ep_3*esat) dqsl = qsl*ep_2*ev/( rd*t**2 ) RH(k)=MAX(MIN(1.0,qw(k)/MAX(1.E-8,qsl)),0.001) alp(k) = 1.0/( 1.0+dqsl*xlvcp ) bet(k) = dqsl*p2a IF (bl_mynn_cloudpdf == 0) THEN t3sq = MAX( tsq(k), 0.0 ) r3sq = MAX( qsq(k), 0.0 ) c3sq = cov(k) c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq ) r3sq = r3sq +bet(k)**2*t3sq -2.0*bet(k)*c3sq qmq(k) = qw(k) -qsl sgm(k) = SQRT( MAX( r3sq, 1.0d-10 )) q1(k) = qmq(k) / sgm(k) cld(k) = 0.5*( 1.0+erf( q1(k)*rr2 ) ) ELSE IF (bl_mynn_cloudpdf == 1 .OR. bl_mynn_cloudpdf == -1) THEN if (k .eq. kts) then dzk = 0.5*dz(k) else dzk = 0.5*( dz(k) + dz(k-1) ) end if dth = 0.5*(thl(k+1)+thl(k)) - 0.5*(thl(k)+thl(MAX(k-1,kts))) dqw = 0.5*(qw(k+1) + qw(k)) - 0.5*(qw(k) + qw(MAX(k-1,kts))) sgm(k) = SQRT( MAX( (alp(k)**2 * MAX(el(k)**2,0.1) * & b2 * MAX(Sh(k),0.03))/4. * & (dqw/dzk - bet(k)*(dth/dzk ))**2 , 1.0e-10) ) qmq(k) = qw(k) -qsl q1(k) = qmq(k) / sgm(k) cld(k) = 0.5*( 1.0+erf( q1(k)*rr2 ) ) ELSE IF (bl_mynn_cloudpdf == 2 .OR. bl_mynn_cloudpdf == -2) THEN xl = xl_blend(t) tlk = thl(k)*(p(k)/p1000mb)**rcp qsat_tl = qsat_blend(tlk,p(k)) rsl = xl*qsat_tl / (r_v*tlk**2) cpm = cp + qw(k)*cpv a(k) = 1./(1. + xl*rsl/cpm) qmq(k) = a(k) * (qw(k) - qsat_tl) b(k) = a(k)*rsl dtl = 0.5*(thl(k+1)*(p(k+1)/p1000mb)**rcp + tlk) & & - 0.5*(tlk + thl(MAX(k-1,kts))*(p(MAX(k-1,kts))/p1000mb)**rcp) dqw = 0.5*(qw(k+1) + qw(k)) - 0.5*(qw(k) + qw(MAX(k-1,kts))) if (k .eq. kts) then dzk = 0.5*dz(k) else dzk = 0.5*( dz(k) + dz(k-1) ) end if cdhdz = dtl/dzk + (g/cpm)*(1.+qw(k)) zagl = zagl + dz(k) ls_min = MIN(MAX(zagl,25.),300.) lfac=MIN(4.25+dx/4000.,6.) ls = MAX(MIN(lfac*el(k),900.),ls_min) sgm(k) = MAX(1.e-10, 0.225*ls*SQRT(MAX(0., & & (a(k)*dqw/dzk)**2 & & -2*a(k)*b(k)*cdhdz*dqw/dzk & & +b(k)**2 * cdhdz**2))) q1(k) = qmq(k) / sgm(k) cld(k) = MAX(0., MIN(1., 0.5+0.36*ATAN(1.55*q1(k)))) ENDIF END DO zagl = 0. RHsum=0. RHnum=0. RHmean=0.1 damp =0 PBLH2=MAX(10.,PBLH1) DO k = kts,kte-1 t = th(k)*exner(k) q1k = q1(k) zagl = zagl + dz(k) IF ( bl_mynn_cloudpdf <= 1 .AND. bl_mynn_cloudpdf >= -1) THEN IF (zagl < PBLH2 .AND. PBLH2 > 400.) THEN RHsum=RHsum+RH(k) RHnum=RHnum+1.0 RHmean=RHsum/RHnum ENDIF RHcrit = 1. - 0.35*(1.0 - (MAX(250.- MAX(HFX1,HFXmin),0.0)/200.)**2) if(HFX1 > HFXmin)then cld9=MIN(MAX(0., (rh(k)-RHcrit)/(1.1-RHcrit)), 1.)**2 else cld9=0.0 endif edown=PBLH2*.1 Hshcu=200. + (RHmean+0.5)**1.5*MAX(HFX1,0.)*Hfac if(zagl < PBLH2-edown)then damp=MIN(1.0,exp(-ABS(((PBLH2-edown)-zagl)/edown))) elseif(zagl >= PBLH2-edown .AND. zagl < PBLH2+Hshcu)then damp=1. elseif (zagl >= PBLH2+Hshcu)then damp=MIN(1.0,exp(-ABS((zagl-(PBLH2+Hshcu))/500.))) endif cldfra_bl1D(k)=cld9*damp eq1 = rrp*EXP( -0.5*q1k*q1k ) qll = MAX( cldfra_bl1D(k)*q1k + eq1, 0.0 ) ql (k) = alp(k)*sgm(k)*qll if(cldfra_bl1D(k)>0.01 .and. ql(k)<1.E-6)ql(k)=1.E-6 qc_bl1D(k)=ql(k)*damp eq1 = rrp*EXP( -0.5*q1k*q1k ) qll = MAX( cld(k)*q1k + eq1, 0.0 ) ql (k) = alp(k)*sgm(k)*qll q2p = xlvcp/exner(k) pt = thl(k) +q2p*ql(k) qt = 1.0 +p608*qw(k) -(1.+p608)*ql(k) rac = alp(k)*( cld(k)-qll*eq1 )*( q2p*qt-(1.+p608)*pt ) vt(k) = qt-1.0 -rac*bet(k) vq(k) = p608*pt-tv0 +rac ELSE IF ( bl_mynn_cloudpdf == 2 .OR. bl_mynn_cloudpdf == -2) THEN IF (q1k < 0.) THEN ql (k) = sgm(k)*EXP(1.2*q1k-1) ELSE IF (q1k > 2.) THEN ql (k) = sgm(k)*q1k ELSE ql (k) = sgm(k)*(EXP(-1.) + 0.66*q1k + 0.086*q1k**2) ENDIF if ((cld(k) .gt. 0.) .or. (ql(k) .gt. 0.)) then if (k .le. k_tropo) then ql_limit = 0.005 * qsat_blend( th(k)*exner(k), p(k) ) ql(k) = MIN( ql(k), ql_limit ) else cld(k) = 0. ql(k) = 0. endif endif Fng = 1. xl = xl_blend(t) bb = b(k)*t/th(k) qww = 1.+0.61*qw(k) alpha = 0.61*th(k) beta = (th(k)/t)*(xl/cp) - 1.61*th(k) vt(k) = qww - cld(k)*beta*bb*Fng - 1. vq(k) = alpha + cld(k)*beta*a(k)*Fng - tv0 if (zagl .lt. PBLH2+1000.) cld(k) = MIN( 1., 1.8*cld(k) ) cldfra_bl1D(k) = cld(k) qc_bl1D(k) = ql(k) ENDIF IF (bl_mynn_cloudpdf .LT. 0) THEN cldfra_bl1D(k) = 0.0 qc_bl1D(k) = 0.0 ENDIF IF (QC_BL1D(k) < 1E-6 .AND. ABS(CLDFRA_BL1D(k)) > 0.01) QC_BL1D(k)= 1E-6 IF (CLDFRA_BL1D(k) < 1E-2)THEN CLDFRA_BL1D(k)=0. QC_BL1D(k)=0. ENDIF END DO cld(kte) = cld(kte-1) ql(kte) = ql(kte-1) vt(kte) = vt(kte-1) vq(kte) = vq(kte-1) qc_bl1D(kte)=0. cldfra_bl1D(kte)=0. RETURN END SUBROUTINE mym_condensation SUBROUTINE mynn_tendencies(kts,kte, & &levflag,grav_settling, & &delt,dz, & &u,v,th,tk,qv,qc,qi,qni,qnc, & &p,exner, & &thl,sqv,sqc,sqi,sqw, & &ust,flt,flq,flqv,flqc,wspd,qcg, & &uoce,voce, & &tsq,qsq,cov, & &tcd,qcd, & &dfm,dfh,dfq, & &Du,Dv,Dth,Dqv,Dqc,Dqi,Dqni, & &vdfg1, & &s_aw,s_awthl,s_awqt,s_awqv,s_awqc, & &s_awu,s_awv, & &FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC, & &cldfra_bl1d, & &bl_mynn_cloudmix, & &bl_mynn_mixqt, & &bl_mynn_edmf, & &bl_mynn_edmf_mom ) INTEGER, INTENT(in) :: kts,kte INTEGER, INTENT(in) :: grav_settling,levflag INTEGER, INTENT(in) :: bl_mynn_cloudmix,bl_mynn_mixqt,& bl_mynn_edmf,bl_mynn_edmf_mom LOGICAL, INTENT(IN) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC REAL,DIMENSION(kts:kte+1), INTENT(in) :: s_aw,s_awthl,s_awqt,& s_awqv,s_awqc,s_awu,s_awv REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,th,tk,qv,qc,qi,qni,qnc,& &p,exner,dfq,dz,tsq,qsq,cov,tcd,qcd,cldfra_bl1d REAL, DIMENSION(kts:kte), INTENT(inout) :: thl,sqw,sqv,sqc,sqi,& &dfm,dfh REAL, DIMENSION(kts:kte), INTENT(inout) :: du,dv,dth,dqv,dqc,dqi,& &dqni REAL, INTENT(IN) :: delt,ust,flt,flq,flqv,flqc,wspd,uoce,voce,qcg REAL, DIMENSION(kts:kte) :: dtz,vt,vq,dfhc,dfmc REAL, DIMENSION(kts:kte) :: sqv2,sqc2,sqi2,sqw2,qni2 REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d REAL :: rhs,gfluxm,gfluxp,dztop,maxdfh,mindfh REAL :: grav_settling2,vdfg1 REAL :: t,esat,qsl,onoff INTEGER :: k,kk,nz nz=kte-kts+1 dztop=.5*(dz(kte)+dz(kte-1)) IF (bl_mynn_edmf_mom == 0) THEN onoff=0.0 ELSE onoff=1.0 ENDIF maxdfh=maxval(dfh(1:15)) mindfh=maxdfh*0.01 DO k=kts,kte dtz(k)=delt/dz(k) IF (dfm(k) > dfh(k)) THEN IF (qc(k) > 1.e-6 .OR. & qi(k) > 1.e-6 .OR. & cldfra_bl1D(k) > 0.05 ) THEN dfh(k)= MAX(dfh(k),dfm(k)*0.5) ENDIF ENDIF IF(bl_mynn_edmf==2 .AND. k > 1 .AND. s_aw(k)>0.0) THEN dfh(k)=MAX(mindfh,dfh(k)) dfm(k)=MAX(mindfh,dfm(k)) ENDIF ENDDO k=kts a(1)=0. b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff c(1)=-dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff d(1)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff DO k=kts+1,kte-1 kk=k-kts+1 a(kk)=-dtz(k)*dfm(k) + 0.5*dtz(k)*s_aw(k)*onoff b(kk)=1. + dtz(k)*(dfm(k)+dfm(k+1)) + 0.5*dtz(k)*(s_aw(k)-s_aw(k+1))*onoff c(kk)=-dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff d(kk)=u(k) + dtz(k)*(s_awu(k)-s_awu(k+1))*onoff ENDDO a(nz)=0 b(nz)=1. c(nz)=0. d(nz)=u(kte) CALL tridiag(nz,a,b,c,d) DO k=kts,kte du(k)=(d(k-kts+1)-u(k))/delt ENDDO k=kts a(1)=0. b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff c(1)=-dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff d(1)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff DO k=kts+1,kte-1 kk=k-kts+1 a(kk)=-dtz(k)*dfm(k) + 0.5*dtz(k)*s_aw(k)*onoff b(kk)=1. + dtz(k)*(dfm(k)+dfm(k+1)) + 0.5*dtz(k)*(s_aw(k)-s_aw(k+1))*onoff c(kk)=-dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff d(kk)=v(k) + dtz(k)*(s_awv(k)-s_awv(k+1))*onoff ENDDO a(nz)=0 b(nz)=1. c(nz)=0. d(nz)=v(kte) CALL tridiag(nz,a,b,c,d) DO k=kts,kte dv(k)=(d(k-kts+1)-v(k))/delt ENDDO k=kts a(1)=0. b(1)=1.+dtz(k)*dfh(k+1) - 0.5*dtz(k)*s_aw(k+1) c(1)=-dtz(k)*(dfh(k+1) + 0.5*s_aw(k+1)) rhs= tcd(k) d(1)=thl(k) + dtz(k)*flt + rhs*delt -dtz(1)*s_awthl(kts+1) DO k=kts+1,kte-1 kk=k-kts+1 a(kk)=-dtz(k)*(dfh(k) - 0.5*s_aw(k)) b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1) + 0.5*(s_aw(k)-s_aw(k+1))) c(kk)=-dtz(k)*(dfh(k+1) + 0.5*s_aw(k+1)) rhs= tcd(k) d(kk)=thl(k) + rhs*delt + dtz(k)*(s_awthl(k)-s_awthl(k+1)) ENDDO a(nz)=0. b(nz)=1. c(nz)=0. d(nz)=thl(kte) CALL tridiag(nz,a,b,c,d) DO k=kts,kte thl(k)=d(k-kts+1) ENDDO IF (bl_mynn_mixqt > 0) THEN k=kts a(1)=0. b(1)=1.+dtz(k)*dfh(k+1) - 0.5*dtz(k)*s_aw(k+1) c(1)= -dtz(k)*(dfh(k+1) + 0.5*s_aw(k+1)) rhs= qcd(k) d(1)=sqw(k) + dtz(k)*flq + rhs*delt - dtz(k)*s_awqt(k+1) DO k=kts+1,kte-1 kk=k-kts+1 a(kk)=-dtz(k)*(dfh(k) - 0.5*s_aw(k)) b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1) + 0.5*(s_aw(k)-s_aw(k+1))) c(kk)=-dtz(k)*(dfh(k+1) + 0.5*s_aw(k+1)) rhs= qcd(k) d(kk)=sqw(k) + rhs*delt + dtz(k)*(s_awqt(k)-s_awqt(k+1)) ENDDO a(nz)=0. b(nz)=1. c(nz)=0. d(nz)=sqw(kte) CALL tridiag(nz,a,b,c,d) DO k=kts,kte sqw2(k)=d(k-kts+1) ENDDO ELSE sqw2=sqw ENDIF IF (bl_mynn_mixqt == 0) THEN IF (bl_mynn_cloudmix > 0 .AND. FLAG_QC) THEN k=kts a(1)=0. b(1)=1.+dtz(k)*dfh(k+1) - 0.5*dtz(k)*s_aw(k+1) c(1)=-dtz(k)*dfh(k+1) - 0.5*dtz(k)*s_aw(k+1) rhs = qcd(k) d(1)=sqc(k) + dtz(k)*flqc + rhs*delt -dtz(k)*s_awqc(k+1) DO k=kts+1,kte-1 kk=k-kts+1 a(kk)=-dtz(k)*dfh(k) + 0.5*dtz(k)*s_aw(k) b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1)) + 0.5*dtz(k)*(s_aw(k)-s_aw(k+1)) c(kk)=-dtz(k)*dfh(k+1) - 0.5*dtz(k)*s_aw(k+1) rhs = qcd(k) d(kk)=sqc(k) + rhs*delt + dtz(k)*(s_awqc(k)-s_awqc(k+1)) ENDDO a(nz)=0. b(nz)=1. c(nz)=0. d(nz)=sqc(kte) CALL tridiag(nz,a,b,c,d) DO k=kts,kte sqc2(k)=d(k-kts+1) ENDDO ELSE sqc2=sqc ENDIF ENDIF IF (bl_mynn_mixqt == 0) THEN k=kts a(1)=0. b(1)=1.+dtz(k)*dfh(k+1) - 0.5*dtz(k)*s_aw(k+1) c(1)=-dtz(k)*dfh(k+1) - 0.5*dtz(k)*s_aw(k+1) d(1)=sqv(k) + dtz(k)*flqv + qcd(k)*delt - dtz(k)*s_awqv(k+1) DO k=kts+1,kte-1 kk=k-kts+1 a(kk)=-dtz(k)*dfh(k) + 0.5*dtz(k)*s_aw(k) b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1)) + 0.5*dtz(k)*(s_aw(k)-s_aw(k+1)) c(kk)=-dtz(k)*dfh(k+1) - 0.5*dtz(k)*s_aw(k+1) d(kk)=sqv(k) + qcd(k)*delt + dtz(k)*(s_awqv(k)-s_awqv(k+1)) ENDDO a(nz)=0. b(nz)=1. c(nz)=0. d(nz)=sqv(kte) CALL tridiag(nz,a,b,c,d) DO k=kts,kte sqv2(k)=d(k-kts+1) ENDDO ELSE sqv2=sqv ENDIF IF (bl_mynn_cloudmix > 0 .AND. FLAG_QI) THEN k=kts a(1)=0. b(1)=1.+dtz(k)*dfh(k+1) c(1)=-dtz(k)*dfh(k+1) d(1)=sqi(k) DO k=kts+1,kte-1 kk=k-kts+1 a(kk)=-dtz(k)*dfh(k) b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1)) c(kk)=-dtz(k)*dfh(k+1) d(kk)=sqi(k) ENDDO a(nz)=0. b(nz)=1. c(nz)=0. d(nz)=sqi(kte) CALL tridiag(nz,a,b,c,d) DO k=kts,kte sqi2(k)=d(k-kts+1) ENDDO ELSE sqi2=sqi ENDIF qni2=qni DO k=kts,kte IF (bl_mynn_mixqt > 0) THEN t = th(k)*exner(k) esat=esat_blend(t) qsl=ep_2*esat/(p(k)-ep_3*esat) IF (FLAG_QI) THEN sqi2(k) = MAX(0., sqi2(k)) sqc2(k) = MAX(0., sqw2(k) - sqi2(k) - qsl) sqv2(k) = MAX(0., sqw2(k) - sqc2(k) - sqi2(k)) ELSE sqi2(k) = 0.0 sqc2(k) = MAX(0., sqw2(k) - qsl) sqv2(k) = MAX(0., sqw2(k) - sqc2(k)) ENDIF ENDIF Dqv(k)=(sqv2(k)/(1.-sqv2(k)) - qv(k))/delt IF (bl_mynn_cloudmix > 0 .AND. FLAG_QC) THEN Dqc(k)=(sqc2(k)/(1.-sqc2(k)) - qc(k))/delt IF(Dqc(k)*delt + qc(k) < 0.) THEN Dqc(k)=-qc(k)/delt ENDIF ELSE Dqc(k)=0. ENDIF IF (bl_mynn_cloudmix > 0 .AND. FLAG_QI) THEN Dqi(k)=(sqi2(k)/(1.-sqi2(k)) - qi(k))/delt IF(Dqi(k)*delt + qi(k) < 0.) THEN Dqi(k)=-qi(k)/delt ENDIF IF (FLAG_QNI) THEN Dqni(k)=(qni2(k)-qni(k))/delt IF(Dqni(k)*delt + qni(k) < 0.)Dqni(k)=-qni(k)/delt ELSE Dqni(k)=0. ENDIF ELSE Dqi(k)=0. Dqni(k)=0. ENDIF IF (FLAG_QI) THEN Dth(k)=(thl(k) + xlvcp/exner(k)*sqc(k) & & + xlscp/exner(k)*sqi(k) & & - th(k))/delt ELSE Dth(k)=(thl(k)+xlvcp/exner(k)*sqc2(k) - th(k))/delt ENDIF ENDDO END SUBROUTINE mynn_tendencies SUBROUTINE retrieve_exchange_coeffs(kts,kte,& &dfm,dfh,dfq,dz,& &K_m,K_h,K_q) INTEGER , INTENT(in) :: kts,kte REAL, DIMENSION(KtS:KtE), INTENT(in) :: dz,dfm,dfh,dfq REAL, DIMENSION(KtS:KtE), INTENT(out) :: & &K_m, K_h, K_q INTEGER :: k REAL :: dzk K_m(kts)=0. K_h(kts)=0. K_q(kts)=0. DO k=kts+1,kte dzk = 0.5 *( dz(k)+dz(k-1) ) K_m(k)=dfm(k)*dzk K_h(k)=dfh(k)*dzk K_q(k)=Sqfac*dfq(k)*dzk ENDDO END SUBROUTINE retrieve_exchange_coeffs SUBROUTINE tridiag(n,a,b,c,d) INTEGER, INTENT(in):: n REAL, DIMENSION(n), INTENT(in) :: a,b REAL, DIMENSION(n), INTENT(inout) :: c,d INTEGER :: i REAL :: p REAL, DIMENSION(n) :: q c(n)=0. q(1)=-c(1)/b(1) d(1)=d(1)/b(1) DO i=2,n p=1./(b(i)+a(i)*q(i-1)) q(i)=-c(i)*p d(i)=(d(i)-a(i)*d(i-1))*p ENDDO DO i=n-1,1,-1 d(i)=d(i)+q(i)*d(i+1) ENDDO END SUBROUTINE tridiag SUBROUTINE mynn_bl_driver( & &initflag,grav_settling, & &delt,dz,dx,znt, & &u,v,w,th,qv,qc,qi,qni,qnc, & &p,exner,rho,T3D, & &xland,ts,qsfc,qcg,ps, & &ust,ch,hfx,qfx,rmol,wspd, & &uoce,voce, & &vdfg, & &Qke,tke_pbl, & &qke_adv,bl_mynn_tkeadvect, & &Tsq,Qsq,Cov, & &RUBLTEN,RVBLTEN,RTHBLTEN, & &RQVBLTEN,RQCBLTEN,RQIBLTEN, & &RQNIBLTEN, & &exch_h,exch_m, & &Pblh,kpbl, & &el_pbl, & &dqke,qWT,qSHEAR,qBUOY,qDISS, & &wstar,delta, & &bl_mynn_tkebudget, & &bl_mynn_cloudpdf,Sh3D, & &bl_mynn_mixlength, & &icloud_bl,qc_bl,cldfra_bl, & &bl_mynn_edmf, & &bl_mynn_edmf_mom,bl_mynn_edmf_tke, & &bl_mynn_edmf_part, & &bl_mynn_cloudmix,bl_mynn_mixqt, & &edmf_a,edmf_w,edmf_qt, & &edmf_thl,edmf_ent,edmf_qc, & &spp_pbl,pattern_spp_pbl, & &FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC & &,IDS,IDE,JDS,JDE,KDS,KDE & &,IMS,IME,JMS,JME,KMS,KME & &,ITS,ITE,JTS,JTE,KTS,KTE) INTEGER, INTENT(in) :: initflag INTEGER, INTENT(in) :: grav_settling INTEGER, INTENT(in) :: bl_mynn_tkebudget INTEGER, INTENT(in) :: bl_mynn_cloudpdf INTEGER, INTENT(in) :: bl_mynn_mixlength INTEGER, INTENT(in) :: bl_mynn_edmf LOGICAL, INTENT(IN) :: bl_mynn_tkeadvect INTEGER, INTENT(in) :: bl_mynn_edmf_mom INTEGER, INTENT(in) :: bl_mynn_edmf_tke INTEGER, INTENT(in) :: bl_mynn_edmf_part INTEGER, INTENT(in) :: bl_mynn_cloudmix INTEGER, INTENT(in) :: bl_mynn_mixqt INTEGER, INTENT(in) :: icloud_bl LOGICAL, INTENT(IN) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC INTEGER,INTENT(IN) :: & & IDS,IDE,JDS,JDE,KDS,KDE & &,IMS,IME,JMS,JME,KMS,KME & &,ITS,ITE,JTS,JTE,KTS,KTE REAL, INTENT(in) :: delt,dx REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(in) :: dz,& &u,v,w,th,qv,p,exner,rho,T3D REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), OPTIONAL, INTENT(in)::& &qc,qi,qni,qnc REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(in) :: xland,ust,& &ch,rmol,ts,qsfc,qcg,ps,hfx,qfx, wspd,uoce,voce, vdfg,znt REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: & &Qke,Tsq,Qsq,Cov, & &tke_pbl, & &qke_adv REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: & &RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN,& &RQIBLTEN,RQNIBLTEN REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: & &exch_h,exch_m REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), OPTIONAL, INTENT(inout) :: & & edmf_a,edmf_w,edmf_qt,edmf_thl,edmf_ent,edmf_qc REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(inout) :: & &Pblh,wstar,delta REAL, DIMENSION(IMS:IME,JMS:JME) :: & &Psig_bl,Psig_shcu INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: & &KPBL REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: & &el_pbl REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: & &qWT,qSHEAR,qBUOY,qDISS,dqke REAL, DIMENSION(KTS:KTE) :: qWT1,qSHEAR1,qBUOY1,qDISS1,dqke1 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME) :: K_q,Sh3D REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: & &qc_bl,cldfra_bl REAL, DIMENSION(KTS:KTE) :: qc_bl1D,cldfra_bl1D,& qc_bl1D_old,cldfra_bl1D_old INTEGER :: ITF,JTF,KTF, IMD,JMD INTEGER :: i,j,k REAL, DIMENSION(KTS:KTE) :: thl,tl,sqv,sqc,sqi,sqw,& &El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, & &Vt, Vq, sgm REAL, DIMENSION(KTS:KTE) :: thetav,sh,u1,v1,w1,p1,ex1,dz1,th1,tk1,rho1,& & qke1,tsq1,qsq1,cov1,qv1,qi1,qc1,du1,dv1,dth1,dqv1,dqc1,dqi1, & & k_m1,k_h1,k_q1,qni1,dqni1,qnc1 REAL, DIMENSION(KTS:KTE) :: dth1mf,dqv1mf,dqc1mf,du1mf,dv1mf REAL, DIMENSION(KTS:KTE) :: edmf_a1,edmf_w1,edmf_qt1,edmf_thl1,& edmf_ent1,edmf_qc1 REAL,DIMENSION(KTS:KTE+1) :: s_aw1,s_awthl1,s_awqt1,& s_awqv1,s_awqc1,s_awu1,s_awv1,s_awqke1 REAL, DIMENSION(KTS:KTE+1) :: zw REAL :: cpm,sqcg,flt,flq,flqv,flqc,pmz,phh,exnerg,zet,& &afk,abk,ts_decay,th_sfc real,parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001 real,parameter :: h1 = 0.33333335, h2 = 0.6666667 REAL :: govrth, sflux, bfx0, wstar3, wm2, wm3, delb INTEGER, SAVE :: levflag INTEGER, INTENT(IN) ::spp_pbl REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN),OPTIONAL ::pattern_spp_pbl REAL, DIMENSION(KTS:KTE) :: rstoch_col IF ( wrf_at_debug_level(3000) ) THEN WRITE ( mynn_message , FMT='(A)' ) & 'in MYNN driver; at beginning' CALL wrf_debug ( 0 , mynn_message ) ENDIF IMD=(IMS+IME)/2 JMD=(JMS+JME)/2 JTF=MIN0(JTE,JDE-1) ITF=MIN0(ITE,IDE-1) KTF=MIN0(KTE,KDE-1) levflag=mynn_level IF (bl_mynn_edmf > 0) THEN edmf_a(its:ite,kts:kte,jts:jte)=0. edmf_w(its:ite,kts:kte,jts:jte)=0. edmf_qt(its:ite,kts:kte,jts:jte)=0. edmf_thl(its:ite,kts:kte,jts:jte)=0. edmf_ent(its:ite,kts:kte,jts:jte)=0. edmf_qc(its:ite,kts:kte,jts:jte)=0. ENDIF IF (initflag > 0) THEN Sh3D(its:ite,kts:kte,jts:jte)=0. el_pbl(its:ite,kts:kte,jts:jte)=0. tsq(its:ite,kts:kte,jts:jte)=0. qsq(its:ite,kts:kte,jts:jte)=0. cov(its:ite,kts:kte,jts:jte)=0. dqc1(kts:kte)=0.0 dqi1(kts:kte)=0.0 dqni1(kts:kte)=0.0 qc_bl1D(kts:kte)=0.0 cldfra_bl1D(kts:kte)=0.0 qc_bl1D_old(kts:kte)=0.0 cldfra_bl1D_old(kts:kte)=0.0 edmf_a1(kts:kte)=0.0 edmf_qc1(kts:kte)=0.0 sgm(kts:kte)=0.0 vt(kts:kte)=0.0 vq(kts:kte)=0.0 DO j=JTS,JTF DO i=ITS,ITF DO k=KTS,KTF dz1(k)=dz(i,k,j) u1(k) = u(i,k,j) v1(k) = v(i,k,j) w1(k) = w(i,k,j) th1(k)=th(i,k,j) tk1(k)=T3D(i,k,j) rho1(k)=rho(i,k,j) sqc(k)=qc(i,k,j)/(1.+qc(i,k,j)) sqv(k)=qv(i,k,j)/(1.+qv(i,k,j)) thetav(k)=th(i,k,j)*(1.+0.61*sqv(k)) IF (PRESENT(qi) .AND. FLAG_QI ) THEN sqi(k)=qi(i,k,j)/(1.+qi(i,k,j)) sqw(k)=sqv(k)+sqc(k)+sqi(k) thl(k)=th(i,k,j)- xlvcp/exner(i,k,j)*sqc(k) & & - xlscp/exner(i,k,j)*sqi(k) ELSE sqi(k)=0.0 sqw(k)=sqv(k)+sqc(k) thl(k)=th(i,k,j)-xlvcp/exner(i,k,j)*sqc(k) ENDIF IF (k==kts) THEN zw(k)=0. ELSE zw(k)=zw(k-1)+dz(i,k-1,j) ENDIF exch_m(i,k,j)=0. exch_h(i,k,j)=0. K_q(i,k,j)=0. qke(i,k,j)=0.1-MIN(zw(k)*0.001, 0.0) qke1(k)=qke(i,k,j) el(k)=el_pbl(i,k,j) sh(k)=Sh3D(i,k,j) tsq1(k)=tsq(i,k,j) qsq1(k)=qsq(i,k,j) cov1(k)=cov(i,k,j) if (spp_pbl==1) then rstoch_col(k)=pattern_spp_pbl(i,k,j) else rstoch_col(k)=0.0 endif IF ( bl_mynn_tkebudget == 1) THEN qWT(i,k,j)=0. qSHEAR(i,k,j)=0. qBUOY(i,k,j)=0. qDISS(i,k,j)=0. dqke(i,k,j)=0. ENDIF ENDDO zw(kte+1)=zw(kte)+dz(i,kte,j) CALL GET_PBLH(KTS,KTE,PBLH(i,j),thetav,& & Qke1,zw,dz1,xland(i,j),KPBL(i,j)) IF (scaleaware > 0.) THEN CALL SCALE_AWARE(dx,PBLH(i,j),Psig_bl(i,j),Psig_shcu(i,j)) ELSE Psig_bl(i,j)=1.0 Psig_shcu(i,j)=1.0 ENDIF CALL mym_initialize ( & &kts,kte, & &dz1, zw, u1, v1, thl, sqv, & &PBLH(i,j), th1, sh, & &ust(i,j), rmol(i,j), & &el, Qke1, Tsq1, Qsq1, Cov1, & &Psig_bl(i,j), cldfra_bl1D, & &bl_mynn_mixlength, & &edmf_a1,edmf_qc1,bl_mynn_edmf,& &spp_pbl,rstoch_col ) DO k=KTS,KTE el_pbl(i,k,j)=el(k) sh3d(i,k,j)=sh(k) qke(i,k,j)=qke1(k) tsq(i,k,j)=tsq1(k) qsq(i,k,j)=qsq1(k) cov(i,k,j)=cov1(k) IF (bl_mynn_tkeadvect) THEN qke_adv(i,k,j)=qke1(k) ENDIF ENDDO ENDDO ENDDO ENDIF IF (bl_mynn_tkeadvect) THEN qke=qke_adv ENDIF DO j=JTS,JTF DO i=ITS,ITF DO k=KTS,KTF IF ( bl_mynn_tkebudget == 1) THEN dqke(i,k,j)=qke(i,k,j) END IF dz1(k)= dz(i,k,j) u1(k) = u(i,k,j) v1(k) = v(i,k,j) w1(k) = w(i,k,j) th1(k)= th(i,k,j) tk1(k)=T3D(i,k,j) rho1(k)=rho(i,k,j) qv1(k)= qv(i,k,j) qc1(k)= qc(i,k,j) sqv(k)= qv(i,k,j)/(1.+qv(i,k,j)) sqc(k)= qc(i,k,j)/(1.+qc(i,k,j)) IF(icloud_bl > 0)cldfra_bl1D_old(k)=cldfra_bl(i,k,j) IF(icloud_bl > 0)qc_bl1D_old(k)=qc_bl(i,k,j) dqc1(k)=0.0 dqi1(k)=0.0 dqni1(k)=0.0 IF(PRESENT(qi) .AND. FLAG_QI)THEN qi1(k)= qi(i,k,j) sqi(k)= qi(i,k,j)/(1.+qi(i,k,j)) sqw(k)= sqv(k)+sqc(k)+sqi(k) thl(k)= th(i,k,j) - xlvcp/exner(i,k,j)*sqc(k) & & - xlscp/exner(i,k,j)*sqi(k) ELSE qi1(k)=0.0 sqi(k)=0.0 sqw(k)= sqv(k)+sqc(k) thl(k)= th(i,k,j)-xlvcp/exner(i,k,j)*sqc(k) ENDIF IF (PRESENT(qni) .AND. FLAG_QNI ) THEN qni1(k)=qni(i,k,j) ELSE qni1(k)=0.0 ENDIF IF (PRESENT(qnc) .AND. FLAG_QNC ) THEN qnc1(k)=qnc(i,k,j) ELSE qnc1(k)=0.0 ENDIF thetav(k)=th(i,k,j)*(1.+0.608*sqv(k)) p1(k) = p(i,k,j) ex1(k)= exner(i,k,j) el(k) = el_pbl(i,k,j) qke1(k)=qke(i,k,j) sh(k) = sh3d(i,k,j) tsq1(k)=tsq(i,k,j) qsq1(k)=qsq(i,k,j) cov1(k)=cov(i,k,j) if (spp_pbl==1) then rstoch_col(k)=pattern_spp_pbl(i,k,j) else rstoch_col(k)=0.0 endif edmf_a1(k)=0.0 edmf_qc1(k)=0.0 s_aw1(k)=0. s_awthl1(k)=0. s_awqt1(k)=0. s_awqv1(k)=0. s_awqc1(k)=0. s_awu1(k)=0. s_awv1(k)=0. s_awqke1(k)=0. IF (k==kts) THEN zw(k)=0. ELSE zw(k)=zw(k-1)+dz(i,k-1,j) ENDIF ENDDO zw(kte+1)=zw(kte)+dz(i,kte,j) s_aw1(kte+1)=0. s_awthl1(kte+1)=0. s_awqt1(kte+1)=0. s_awqv1(kte+1)=0. s_awqc1(kte+1)=0. s_awu1(kte+1)=0. s_awv1(kte+1)=0. s_awqke1(kte+1)=0. CALL GET_PBLH(KTS,KTE,PBLH(i,j),thetav,& & Qke1,zw,dz1,xland(i,j),KPBL(i,j)) IF (scaleaware > 0.) THEN CALL SCALE_AWARE(dx,PBLH(i,j),Psig_bl(i,j),Psig_shcu(i,j)) ELSE Psig_bl(i,j)=1.0 Psig_shcu(i,j)=1.0 ENDIF sqcg= 0.0 cpm=cp*(1.+0.84*qv(i,kts,j)) exnerg=(ps(i,j)/p1000mb)**rcp flt = hfx(i,j)/( rho(i,kts,j)*cpm ) & & +xlvcp*vdfg(i,j)*(sqc(kts)/exner(i,kts,j)- sqcg/exnerg) flq = qfx(i,j)/ rho(i,kts,j) & & -vdfg(i,j)*(sqc(kts) - sqcg ) flqv = qfx(i,j)/rho(i,kts,j) flqc = -vdfg(i,j)*(sqc(kts) - sqcg ) th_sfc = ts(i,j)/ex1(kts) zet = 0.5*dz(i,kts,j)*rmol(i,j) if ( zet >= 0.0 ) then pmz = 1.0 + (cphm_st-1.0) * zet phh = 1.0 + cphh_st * zet else pmz = 1.0/ (1.0-cphm_unst*zet)**0.25 - zet phh = 1.0/SQRT(1.0-cphh_unst*zet) end if govrth = g/th1(kts) sflux = hfx(i,j)/rho(i,kts,j)/cpm + & qfx(i,j)/rho(i,kts,j)*ep_1*th1(kts) bfx0 = max(sflux,0.) wstar3 = (govrth*bfx0*pblh(i,j)) wstar(i,j) = wstar3**h1 wm3 = wstar3 + 5.*ust(i,j)**3. wm2 = wm3**h2 delb = govrth*d3*pblh(i,j) delta(i,j) = min(d1*pblh(i,j) + d2*wm2/delb, 100.) CALL mym_condensation ( kts,kte, & &dx,dz1,thl,sqw,p1,ex1, & &tsq1, qsq1, cov1, & &Sh,el,bl_mynn_cloudpdf, & &qc_bl1D,cldfra_bl1D, & &PBLH(i,j),HFX(i,j), & &Vt, Vq, th1, sgm ) IF (bl_mynn_edmf == 1) THEN CALL StEM_mf( & &kts,kte,delt,zw,dz1,p1, & &bl_mynn_edmf_mom, & &bl_mynn_edmf_tke, & &u1,v1,w1,th1,thl,thetav,tk1, & &sqw,sqv,sqc,qke1, & &ex1,Vt,Vq,sgm, & &ust(i,j),flt,flq,flqv,flqc, & &PBLH(i,j),KPBL(i,j),DX, & &xland(i,j),th_sfc, & & edmf_a1,edmf_w1,edmf_qt1, & & edmf_thl1,edmf_ent1,edmf_qc1, & & s_aw1,s_awthl1,s_awqt1, & & s_awqv1,s_awqc1,s_awu1,s_awv1, & & s_awqke1, & & qc_bl1D,cldfra_bl1D, & & FLAG_QI,FLAG_QC, & & Psig_shcu(i,j), & & spp_pbl,rstoch_col & ) ELSEIF (bl_mynn_edmf == 2) THEN CALL temf_mf( & &kts,kte,delt,zw,p1,ex1, & &u1,v1,w1,th1,thl,thetav, & &sqw,sqv,sqc,qke1, & &ust(i,j),flt,flq,flqv,flqc, & &hfx(i,j),qfx(i,j),ts(i,j), & &pblh(i,j),rho1,dfh,dx,znt(i,j),ep_2, & & edmf_a1,edmf_w1,edmf_qt1, & & edmf_thl1,edmf_ent1,edmf_qc1, & & s_aw1,s_awthl1,s_awqt1, & & s_awqv1,s_awqc1, & & s_awu1,s_awv1,s_awqke1, & & qc_bl1D,cldfra_bl1D & &,FLAG_QI,FLAG_QC & &,Psig_shcu(i,j) & &,spp_pbl,rstoch_col & &,i,j,ids,ide,jds,jde & ) ENDIF CALL mym_turbulence ( & &kts,kte,levflag, & &dz1, zw, u1, v1, thl, sqc, sqw, & &qke1, tsq1, qsq1, cov1, & &vt, vq, & &rmol(i,j), flt, flq, & &PBLH(i,j),th1, & &Sh,el, & &Dfm,Dfh,Dfq, & &Tcd,Qcd,Pdk, & &Pdt,Pdq,Pdc, & &qWT1,qSHEAR1,qBUOY1,qDISS1, & &bl_mynn_tkebudget, & &Psig_bl(i,j),Psig_shcu(i,j), & &cldfra_bl1D,bl_mynn_mixlength, & &edmf_a1,edmf_qc1,bl_mynn_edmf, & &spp_pbl,rstoch_col) IF (bl_mynn_edmf_part > 0 .AND. bl_mynn_edmf > 0) THEN DO k=kts,kte dfm(k)=dfm(k) * (1. - 0.5*edmf_a1(k)) dfh(k)=dfh(k) * (1. - 0.5*edmf_a1(k)) dfq(k)=dfq(k) * (1. - 0.5*edmf_a1(k)) ENDDO ENDIF CALL mym_predict (kts,kte,levflag, & &delt, dz1, & &ust(i,j), flt, flq, pmz, phh, & &el, dfq, pdk, pdt, pdq, pdc, & &Qke1, Tsq1, Qsq1, Cov1, & &s_aw1, s_awqke1, bl_mynn_edmf_tke) CALL mynn_tendencies(kts,kte, & &levflag,grav_settling, & &delt, dz1, & &u1, v1, th1, tk1, qv1, qc1, qi1, & &qni1,qnc1, & &p1, ex1, thl, sqv, sqc, sqi, sqw,& &ust(i,j),flt,flq,flqv,flqc, & &wspd(i,j),qcg(i,j), & &uoce(i,j),voce(i,j), & &tsq1, qsq1, cov1, & &tcd, qcd, & &dfm, dfh, dfq, & &Du1, Dv1, Dth1, Dqv1, & &Dqc1, Dqi1, Dqni1, & &vdfg(i,j), & &s_aw1,s_awthl1,s_awqt1, & &s_awqv1,s_awqc1,s_awu1,s_awv1, & &FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC,& &cldfra_bl1d, & &bl_mynn_cloudmix, & &bl_mynn_mixqt, & &bl_mynn_edmf, & &bl_mynn_edmf_mom) CALL retrieve_exchange_coeffs(kts,kte,& &dfm, dfh, dfq, dz1,& &K_m1, K_h1, K_q1) DO k=KTS,KTF exch_m(i,k,j)=K_m1(k) exch_h(i,k,j)=K_h1(k) K_q(i,k,j)=K_q1(k) RUBLTEN(i,k,j)=du1(k) RVBLTEN(i,k,j)=dv1(k) RTHBLTEN(i,k,j)=dth1(k) RQVBLTEN(i,k,j)=dqv1(k) IF(bl_mynn_cloudmix > 0)THEN IF (PRESENT(qc) .AND. FLAG_QC) RQCBLTEN(i,k,j)=dqc1(k) IF (PRESENT(qi) .AND. FLAG_QI) RQIBLTEN(i,k,j)=dqi1(k) IF (PRESENT(qni) .AND. FLAG_QNI) RQNIBLTEN(i,k,j)=dqni1(k) ELSE IF (PRESENT(qc) .AND. FLAG_QC) RQCBLTEN(i,k,j)=0. IF (PRESENT(qi) .AND. FLAG_QI) RQIBLTEN(i,k,j)=0. IF (PRESENT(qni) .AND. FLAG_QNI) RQNIBLTEN(i,k,j)=0. ENDIF IF(icloud_bl > 0)THEN qc_bl(i,k,j)=qc_bl1D(k) cldfra_bl(i,k,j)=cldfra_bl1D(k) IF (CLDFRA_BL(i,k,j) > cldfra_bl1D_old(k)) THEN ELSE ts_decay = MIN( 1800., 3.*dx/MAX(SQRT(u1(k)**2 + v1(k)**2), 1.0) ) cldfra_bl(i,k,j)= MAX(cldfra_bl1D(k), cldfra_bl1D_old(k)-(0.25*delt/ts_decay)) IF (cldfra_bl(i,k,j) > 0.1) THEN IF (QC_BL(i,k,j) < 1E-5)QC_BL(i,k,j)= MAX(qc_bl1D_old(k), 1E-5) ELSE CLDFRA_BL(i,k,j)= 0. QC_BL(i,k,j) = 0. ENDIF ENDIF if (spp_pbl==1) then cldfra_bl(i,k,j)= cldfra_bl(i,k,j)*(1.0-rstoch_col(k)) IF ((cldfra_bl(i,k,j) > 1.0) .or. (cldfra_bl(i,k,j) < 0.0)) then cldfra_bl(i,k,j)=MAX(MIN(cldfra_bl(i,k,j),1.0),0.0) ENDIF ENDIF IF (icloud_bl > 0) THEN IF (QC_BL(i,k,j) < 1E-6 .AND. ABS(CLDFRA_BL(i,k,j)) > 0.1)QC_BL(i,k,j)= 1E-6 IF (CLDFRA_BL(i,k,j) < 1E-2)CLDFRA_BL(i,k,j)= 0. ENDIF ENDIF el_pbl(i,k,j)=el(k) qke(i,k,j)=qke1(k) tsq(i,k,j)=tsq1(k) qsq(i,k,j)=qsq1(k) cov(i,k,j)=cov1(k) sh3d(i,k,j)=sh(k) IF ( bl_mynn_tkebudget == 1) THEN dqke(i,k,j) = (qke1(k)-dqke(i,k,j))*0.5 qWT(i,k,j) = qWT1(k)*delt qSHEAR(i,k,j)= qSHEAR1(k)*delt qBUOY(i,k,j) = qBUOY1(k)*delt qDISS(i,k,j) = qDISS1(k)*delt ENDIF IF (bl_mynn_edmf > 0) THEN edmf_a(i,k,j)=edmf_a1(k) edmf_w(i,k,j)=edmf_w1(k) edmf_qt(i,k,j)=edmf_qt1(k) edmf_thl(i,k,j)=edmf_thl1(k) edmf_ent(i,k,j)=edmf_ent1(k) edmf_qc(i,k,j)=edmf_qc1(k) ENDIF IF ( wrf_at_debug_level(3000) ) THEN IF ( sh(k) < 0. .OR. sh(k)> 200.)print*,"**SUSPICIOUS VALUES AT: i,j,k=",i,j,k," sh=",sh(k) IF ( qke(i,k,j) < -1. .OR. qke(i,k,j)> 200.)print*,"**SUSPICIOUS VALUES AT: i,j,k=",i,j,k," qke=",qke(i,k,j) IF ( el_pbl(i,k,j) < 0. .OR. el_pbl(i,k,j)> 2000.)print*,"**SUSPICIOUS VALUES AT: i,j,k=",i,j,k," el_pbl=",el_pbl(i,k,j) IF ( ABS(vt(k)) > 0.8 )print*,"**SUSPICIOUS VALUES AT: i,j,k=",i,j,k," vt=",vt(k) IF ( ABS(vq(k)) > 6000.)print*,"**SUSPICIOUS VALUES AT: i,j,k=",i,j,k," vq=",vq(k) IF ( exch_m(i,k,j) < 0. .OR. exch_m(i,k,j)> 2000.)print*,"**SUSPICIOUS VALUES AT: i,j,k=",i,j,k," exxch_m=",exch_m(i,k,j) IF ( vdfg(i,j) < 0. .OR. vdfg(i,j)>5. )print*,"**SUSPICIOUS VALUES AT: i,j,k=",i,j,k," vdfg=",vdfg(i,j) IF ( ABS(QFX(i,j))>.001)print*, "**SUSPICIOUS VALUES AT: i,j=",i,j," QFX=",QFX(i,j) IF ( ABS(HFX(i,j))>1000.)print*, "**SUSPICIOUS VALUES AT: i,j=",i,j," HFX=",HFX(i,j) IF (icloud_bl > 0) then IF( cldfra_bl(i,k,j) < 0.0 .OR. cldfra_bl(i,k,j)> 1.)THEN PRINT*,"**SUSPICIOUS VALUES: CLDFRA_BL=",cldfra_bl(i,k,j)," qc_bl=",QC_BL(i,k,j) ENDIF ENDIF ENDIF ENDDO tke_pbl(i,kts,j) = 0.5*MAX(qke(i,kts,j),1.0e-10) DO k = kts+1,kte afk = dz1(k)/( dz1(k)+dz1(k-1) ) abk = 1.0 -afk tke_pbl(i,k,j) = 0.5*MAX(qke(i,k,j)*abk+qke(i,k-1,j)*afk,1.0e-3) ENDDO ENDDO ENDDO IF (bl_mynn_tkeadvect) THEN qke_adv=qke ENDIF END SUBROUTINE mynn_bl_driver SUBROUTINE mynn_bl_init_driver( & &RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, & &RQCBLTEN,RQIBLTEN & &,QKE,TKE_PBL,EXCH_H & &,RESTART,ALLOWED_TO_READ,LEVEL & &,IDS,IDE,JDS,JDE,KDS,KDE & &,IMS,IME,JMS,JME,KMS,KME & &,ITS,ITE,JTS,JTE,KTS,KTE) LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART INTEGER,INTENT(IN) :: LEVEL INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, & & IMS,IME,JMS,JME,KMS,KME, & & ITS,ITE,JTS,JTE,KTS,KTE REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: & &RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, & &RQCBLTEN,RQIBLTEN,& &QKE,TKE_PBL,EXCH_H INTEGER :: I,J,K,ITF,JTF,KTF JTF=MIN0(JTE,JDE-1) KTF=MIN0(KTE,KDE-1) ITF=MIN0(ITE,IDE-1) IF(.NOT.RESTART)THEN DO J=JTS,JTF DO K=KTS,KTF DO I=ITS,ITF RUBLTEN(i,k,j)=0. RVBLTEN(i,k,j)=0. RTHBLTEN(i,k,j)=0. RQVBLTEN(i,k,j)=0. if( p_qc >= param_first_scalar ) RQCBLTEN(i,k,j)=0. if( p_qi >= param_first_scalar ) RQIBLTEN(i,k,j)=0. TKE_PBL(i,k,j)=0. EXCH_H(i,k,j)=0. ENDDO ENDDO ENDDO ENDIF mynn_level=level END SUBROUTINE mynn_bl_init_driver SUBROUTINE GET_PBLH(KTS,KTE,zi,thetav1D,qke1D,zw1D,dz1D,landsea,kzi) INTEGER,INTENT(IN) :: KTS,KTE REAL, INTENT(OUT) :: zi REAL, INTENT(IN) :: landsea REAL, DIMENSION(KTS:KTE), INTENT(IN) :: thetav1D, qke1D, dz1D REAL, DIMENSION(KTS:KTE+1), INTENT(IN) :: zw1D REAL :: PBLH_TKE,qtke,qtkem1,wt,maxqke,TKEeps,minthv REAL :: delt_thv REAL, PARAMETER :: sbl_lim = 250. REAL, PARAMETER :: sbl_damp = 500. INTEGER :: I,J,K,kthv,ktke,kzi,kzi2 kzi = 2 kzi2= 2 k = kts+1 kthv = 1 ktke = 1 maxqke = 0. minthv = 9.E9 DO WHILE (zw1D(k) .LE. sbl_lim) qtke =MAX(Qke1D(k),0.) IF (maxqke < qtke) then maxqke = qtke ktke = k ENDIF IF (minthv > thetav1D(k)) then minthv = thetav1D(k) kthv = k ENDIF k = k+1 ENDDO TKEeps = maxqke/40. TKEeps = MAX(TKEeps,0.02) zi=0. k = kthv+1 IF((landsea-1.5).GE.0)THEN delt_thv = 0.75 ELSE delt_thv = 1.25 ENDIF zi=0. k = kthv+1 DO WHILE (zi .EQ. 0.) IF (thetav1D(k) .GE. (minthv + delt_thv))THEN kzi = MAX(k-1,1) zi = zw1D(k) - dz1D(k-1)* & & MIN((thetav1D(k)-(minthv + delt_thv))/ & & MAX(thetav1D(k)-thetav1D(k-1),1E-6),1.0) ENDIF k = k+1 IF (k .EQ. kte-1) zi = zw1D(kts+1) ENDDO PBLH_TKE=0. k = ktke+1 DO WHILE (PBLH_TKE .EQ. 0.) qtke =MAX(Qke1D(k)/2.,0.) qtkem1=MAX(Qke1D(k-1)/2.,0.) IF (qtke .LE. TKEeps) THEN kzi2 = MAX(k-1,1) PBLH_TKE = zw1D(k) - dz1D(k-1)* & & MIN((TKEeps-qtke)/MAX(qtkem1-qtke, 1E-6), 1.0) PBLH_TKE = MAX(PBLH_TKE,zw1D(kts+1)) ENDIF k = k+1 IF (k .EQ. kte-1) PBLH_TKE = zw1D(kts+1) ENDDO PBLH_TKE = MIN(PBLH_TKE,zi+500.) PBLH_TKE = MAX(PBLH_TKE,MAX(zi-500.,10.)) wt=.5*TANH((zi - sbl_lim)/sbl_damp) + .5 IF (maxqke <= 0.05) THEN ELSE zi=PBLH_TKE*(1.-wt) + zi*wt ENDIF kzi = MAX(INT(kzi2*(1.-wt) + kzi*wt),1) END SUBROUTINE GET_PBLH SUBROUTINE StEM_mf( & & kts,kte,dt,zw,dz,p, & & momentum_opt, & & tke_opt, & & u,v,w,th,thl,thv,tk, & & qt,qv,qc,qke, & & exner,vt,vq,sgm, & & ust,flt,flq,flqv,flqc, & & pblh,kpbl,DX,landsea,ts, & & edmf_a,edmf_w, & & edmf_qt,edmf_thl, & & edmf_ent,edmf_qc, & & s_aw,s_awthl,s_awqt, & & s_awqv,s_awqc, & & s_awu,s_awv,s_awqke, & & qc_bl1d,cldfra_bl1d, & &F_QC,F_QI, & &Psig_shcu, & &spp_pbl,rstoch_col) INTEGER, INTENT(IN) :: spp_pbl REAL, DIMENSION(KTS:KTE) :: rstoch_col INTEGER, INTENT(IN) :: KTS,KTE,momentum_opt,tke_opt,kpbl REAL,DIMENSION(KTS:KTE), INTENT(IN) :: U,V,W,TH,THL,TK,QT,QV,QC,& THV,P,qke,exner,dz REAL,DIMENSION(KTS:KTE+1), INTENT(IN) :: ZW REAL, INTENT(IN) :: DT,UST,FLT,FLQ,FLQV,FLQC,PBLH,& DX,Psig_shcu,landsea,ts LOGICAL, OPTIONAL :: F_QC,F_QI REAL,DIMENSION(KTS:KTE), INTENT(OUT) :: edmf_a,edmf_w, & & edmf_qt,edmf_thl, edmf_ent,edmf_qc INTEGER :: nup2,ktop REAL :: maxmf REAL,DIMENSION(KTS:KTE+1) :: s_aw, & s_awthl, & s_awqt, & s_awqv, & s_awqc, & s_awu, & s_awv, & s_awqke REAL,DIMENSION(KTS:KTE), INTENT(INOUT) :: qc_bl1d,cldfra_bl1d INTEGER, PARAMETER :: NUP=10 REAL,DIMENSION(KTS:KTE+1,1:NUP) :: UPW,UPTHL,UPQT,UPQC,UPQV, & UPA,UPU,UPV,UPTHV,UPQKE REAL,DIMENSION(KTS:KTE,1:NUP) :: ENT,ENTf INTEGER,DIMENSION(KTS:KTE,1:NUP) :: ENTi INTEGER :: K,I,k50 REAL :: fltv,wstar,qstar,thstar,sigmaW,sigmaQT,sigmaTH,z0, & pwmin,pwmax,wmin,wmax,wlv,wtv,Psig_w,maxw,maxqc,wpbl REAL :: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,Wn2,EntEXP,EntW REAL,PARAMETER :: & &Wa=2./3., & &Wb=0.002,& &Wc=1.5 REAL,PARAMETER :: & & L0=100.,& & ENT0=0.1 REAL, PARAMETER :: Atot = 0.1 REAL, PARAMETER :: lmax = 1000. REAL, PARAMETER :: dl = 100. REAL, PARAMETER :: dcut = 1.0 REAL :: d REAL :: cn,c,l,n,an2,hux,maxwidth,wspd_pbl,cloud_base,width_flx REAL :: ERF INTEGER, PARAMETER :: cldfra_opt = 0 REAL :: xcldfra, UMF_new, dcf, ktop_dcf, kbcon_dcf REAL, PARAMETER :: coef_p = 0.25, coef_gamm = 0.49, coef_alph = 100. REAL :: satvp,rhgrid,h2oliq LOGICAL :: superadiabatic REAL,DIMENSION(KTS:KTE), INTENT(INOUT) :: vt, vq, sgm REAL :: sigq,xl,tlk,qsat_tl,rsl,cpm,a,qmq,mf_cf,diffqt,& Fng,qww,alpha,beta,bb,f,pt,t,q2p,b9 REAL :: csigma,acfac INTEGER :: overshoot REAL :: bvf, Frz UPW=0. UPTHL=0. UPTHV=0. UPQT=0. UPA=0. UPU=0. UPV=0. UPQC=0. UPQV=0. UPQKE=0. ENT=0.001 edmf_a =0. edmf_w =0. edmf_qt =0. edmf_thl=0. edmf_ent=0. edmf_qc =0. s_aw=0. s_awthl=0. s_awqt=0. s_awqv=0. s_awqc=0. s_awu=0. s_awv=0. s_awqke=0. k = 1 maxw = 0.0 cloud_base = 9000.0 DO WHILE (ZW(k) < pblh + 500.) wpbl = w(k) IF(w(k) < 0.)wpbl = 2.*w(k) maxw = MAX(maxw,ABS(wpbl)) IF(ZW(k)<=50.)k50=k IF(qc(k)>1E-5 .AND. cloud_base == 9000.0)THEN cloud_base = 0.5*(ZW(k)+ZW(k+1)) ENDIF k = k + 1 ENDDO maxw = MAX(0.,maxw - 0.5) Psig_w = MAX(0.0, 1.0 - maxw/0.5) Psig_w = MIN(Psig_w, Psig_shcu) fltv = flt + svp1*flq IF(Psig_w == 0.0 .and. fltv > 0.0) fltv = -1.*fltv superadiabatic = .false. IF((landsea-1.5).GE.0)THEN hux = -0.002 ELSE hux = -0.005 ENDIF DO k=1,MAX(1,k50-1) IF (k == 1) then IF ((th(k)-ts)/(0.5*dz(k)) < hux) THEN superadiabatic = .true. ELSE superadiabatic = .false. exit ENDIF ELSE IF ((th(k)-th(k-1))/(0.5*(dz(k)+dz(k-1))) < hux) THEN superadiabatic = .true. ELSE superadiabatic = .false. exit ENDIF ENDIF ENDDO NUP2 = max(1,min(NUP,INT(dx*dcut/dl))) wspd_pbl=SQRT(MAX(u(kpbl)**2 + v(kpbl)**2, 0.01)) maxwidth = 0.9*PBLH - MIN(15.*MAX(wspd_pbl - 7.5, 0.), 0.3*PBLH) maxwidth = MIN(maxwidth,0.5*cloud_base) IF((landsea-1.5).LT.0)THEN IF (cloud_base .LT. 2000.) THEN width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.120)/0.03) + .5),1000.), 0.) ELSE width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.085)/0.04) + .5),1000.), 0.) ENDIF maxwidth = MIN(maxwidth,width_flx) ENDIF NUP2 = MIN(MAX(INT((maxwidth - MOD(maxwidth,100.))/100), 0), NUP2) ktop = 0 maxmf= 0.0 IF ( fltv > 0.002 .AND. NUP2 .GE. 1 .AND. superadiabatic) then cn = 0. d=-1.9 do I=1,NUP2 l = dl*I cn = cn + l**d * (l*l)/(dx*dx) * dl enddo C = Atot/cn An2 = 0. do I=1,NUP2 l = dl*I N = C*l**d UPA(1,I) = N*l*l/(dx*dx) * dl acfac = .5*tanh((fltv - 0.07)/0.09) + .5 UPA(1,I)=UPA(1,I)*acfac An2 = An2 + UPA(1,I) end do z0=50. pwmin=0.1 pwmax=0.5 wstar=max(1.E-2,(g/thv(1)*fltv*pblh)**(1./3.)) qstar=max(flq,1.0E-5)/wstar thstar=flt/wstar IF((landsea-1.5).GE.0)THEN csigma = 1.34 ELSE csigma = 1.34 ENDIF sigmaW =1.34*wstar*(z0/pblh)**(1./3.)*(1 - 0.8*z0/pblh) sigmaQT=csigma*qstar*(z0/pblh)**(-1./3.) sigmaTH=csigma*thstar*(z0/pblh)**(-1./3.) wmin=sigmaW*pwmin wmax=sigmaW*pwmax acfac = .5*tanh((fltv - 0.08)/0.07) + .5 DO I=1,NUP2 wlv=wmin+(wmax-wmin)/NUP2*(i-1) wtv=wmin+(wmax-wmin)/NUP2*i UPW(1,I)=wmin + REAL(i)/REAL(NUP)*(wmax-wmin) UPU(1,I)=U(1) UPV(1,I)=V(1) UPQC(1,I)=0 UPQT(1,I) = 0. UPTHV(1,I)= 0. UPTHL(1,I)= 0. k50=1 DO k=1,k50 UPQT(1,I) = UPQT(1,I) +QT(k) +0.58*UPW(1,I)*sigmaQT/sigmaW *acfac UPTHV(1,I)= UPTHV(1,I)+THV(k)+0.58*UPW(1,I)*sigmaTH/sigmaW *acfac UPTHL(1,I)= UPTHL(1,I)+THL(k)+0.58*UPW(1,I)*sigmaTH/sigmaW *acfac ENDDO UPQT(1,I) = UPQT(1,I)/REAL(k50) UPTHV(1,I)= UPTHV(1,I)/REAL(k50) UPTHL(1,I)= UPTHL(1,I)/REAL(k50) UPQKE(1,I)= QKE(1) ENDDO DO I=1,NUP2 QCn = 0. overshoot = 0 l = dl*I DO k=KTS+1,KTE ENT(k,i) = 0.75/(MIN(MAX(UPW(K-1,I),0.75),1.25)*l) ENT(k,i) = max(ENT(k,i),0.0006) ENT(k,i) = max(ENT(k,i),0.05/ZW(k)) IF(ZW(k) >= pblh+1000.) ENT(k,i) =ENT(k,i) + (ZW(k)-(pblh+1000.)) * 5.0E-6 IF(UPW(K-1,I) > 1.5) ENT(k,i) = ENT(k,i) + 0.004*(UPW(K-1,I) - 1.5) ENT(k,i) = min(ENT(k,i),0.9/(ZW(k)-ZW(k-1))) EntExp= ENT(K,I)*(ZW(k)-ZW(k-1)) QTn =UPQT(k-1,I) *(1.-EntExp) + QT(k-1)*EntExp THLn=UPTHL(k-1,I)*(1.-EntExp) + THL(k-1)*EntExp Un =UPU(k-1,I) *(1.-EntExp) + U(k-1)*EntExp Vn =UPV(k-1,I) *(1.-EntExp) + V(k-1)*EntExp QKEn=UPQKE(k-1,I)*(1.-EntExp) + QKE(k-1)*EntExp call condensation_edmf(QTn,THLn,(P(K)+P(K-1))/2.,ZW(k),THVn,QCn) B=g*(0.5*(THVn+UPTHV(k-1,I))/THV(k-1) - 1.0) EntW=exp(-2.*(Wb+Wc*ENT(K,I))*(ZW(k)-ZW(k-1))) Wn2=UPW(K-1,I)**2*EntW + (1.-EntW)*0.5*Wa*B/(Wb+Wc*ENT(K,I)) Wn2=MAX(Wn2,0.0) IF (fltv > 0.05 .AND. Wn2 <= 0 .AND. overshoot == 0) THEN overshoot = 1 IF ( THV(k)-THV(k-1) .GT. 0.0 ) THEN bvf = SQRT( gtr*(THV(k)-THV(k-1))/(0.5*(dz(k)+dz(k-1))) ) Frz = UPW(K-1,I)/(bvf*0.5*(dz(k)+dz(k-1))) IF ( Frz >= 0.5 ) Wn2 = MIN(Frz,1.0)*(UPW(K-1,I)**2) ENDIF ELSEIF (fltv > 0.05 .AND. overshoot == 1) THEN Wn2 = 0.0 ENDIF Wn2=Wn2*EXP(-MAX(ZW(k)-(pblh+1500.),0.0)/1000.) IF(ZW(k) >= pblh+3000.)Wn2=0. IF (fltv < 0.06) THEN IF(ZW(k) >= pblh-200. .AND. qc(k) > 1e-5 .AND. I > 2) Wn2=0. ENDIF IF (Wn2 > 0.) THEN UPW(K,I)=sqrt(Wn2) UPTHV(K,I)=THVn UPTHL(K,I)=THLn UPQT(K,I)=QTn UPQC(K,I)=QCn UPU(K,I)=Un UPV(K,I)=Vn UPQKE(K,I)=QKEn UPA(K,I)=UPA(K-1,I) ktop = MAX(ktop,k) ELSE exit END IF ENDDO ENDDO ELSE NUP2=0. END IF ktop=MIN(ktop,KTE-1) IF(nup2 > 0) THEN DO k=KTS,ktop DO I=1,NUP2 edmf_a(K)=edmf_a(K)+UPA(K+1,I) edmf_w(K)=edmf_w(K)+UPA(K+1,I)*UPW(K+1,I) edmf_qt(K)=edmf_qt(K)+UPA(K+1,I)*UPQT(K+1,I) edmf_thl(K)=edmf_thl(K)+UPA(K+1,I)*UPTHL(K+1,I) edmf_ent(K)=edmf_ent(K)+UPA(K+1,I)*ENT(K+1,I) edmf_qc(K)=edmf_qc(K)+UPA(K+1,I)*UPQC(K+1,I) ENDDO IF (edmf_a(k)>0.) THEN edmf_w(k)=edmf_w(k)/edmf_a(k) edmf_qt(k)=edmf_qt(k)/edmf_a(k) edmf_thl(k)=edmf_thl(k)/edmf_a(k) edmf_ent(k)=edmf_ent(k)/edmf_a(k) edmf_qc(k)=edmf_qc(k)/edmf_a(k) edmf_a(k)=edmf_a(k)*Psig_w IF(edmf_a(k)*edmf_w(k) > maxmf) maxmf = edmf_a(k)*edmf_w(k) ENDIF ENDDO DO k=KTS,ktop DO I=1,NUP2 s_aw(k) = s_aw(K) + UPA(K,I)*UPW(K,I)*Psig_w * (1.0+rstoch_col(k)) s_awthl(k)= s_awthl(K) + UPA(K,i)*UPW(K,I)*UPTHL(K,I)*Psig_w * (1.0+rstoch_col(k)) s_awqt(k) = s_awqt(K) + UPA(K,i)*UPW(K,I)*UPQT(K,I)*Psig_w * (1.0+rstoch_col(k)) s_awqc(k) = s_awqc(K) + UPA(K,i)*UPW(K,I)*UPQC(K,I)*Psig_w * (1.0+rstoch_col(k)) IF (momentum_opt > 0) THEN s_awu(k) = s_awu(K) + UPA(K,i)*UPW(K,I)*UPU(K,I)*Psig_w * (1.0+rstoch_col(k)) s_awv(k) = s_awv(K) + UPA(K,i)*UPW(K,I)*UPV(K,I)*Psig_w * (1.0+rstoch_col(k)) ENDIF IF (tke_opt > 0) THEN s_awqke(k)= s_awqke(K) + UPA(K,i)*UPW(K,I)*UPQKE(K,I)*Psig_w * (1.0+rstoch_col(k)) ENDIF ENDDO s_awqv(k) = s_awqt(k) - s_awqc(k) ENDDO DO K=KTS,ktop IF (cldfra_opt == 0) THEN IF(edmf_qc(k)>0.0)THEN satvp = 3.80*exp(17.27*(th(k)-273.)/ & (th(k)-36.))/(.01*p(k)) rhgrid = max(.01,MIN( 1., qv(k) /satvp)) xl = xl_blend(tk(k)) tlk = thl(k)*(p(k)/p1000mb)**rcp qsat_tl = qsat_blend(tlk,p(k)) rsl = xl*qsat_tl / (r_v*tlk**2) cpm = cp + qt(k)*cpv a = 1./(1. + xl*rsl/cpm) b9 = a*rsl q2p = xlvcp/exner(k) pt = thl(k) +q2p*edmf_qc(k) bb = b9*tk(k)/pt qww = 1.+0.61*qt(k) alpha = 0.61*pt t = th(k)*exner(k) beta = pt*xl/(t*cp) - 1.61*pt if (a > 0.0) then f = MIN(1.0/a, 4.0) else f = 1.0 endif sigq = 6.E-3 * edmf_a(k) * edmf_w(k) * f sigq = SQRT(sigq**2 + sgm(k)**2) qmq = a * (qt(k) - qsat_tl) mf_cf = min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.01),1.0) IF (rhgrid >= .93) THEN cldfra_bl1d(k) = MAX(mf_cf, cldfra_bl1d(k)) IF (cldfra_bl1d(k) > edmf_a(k)) THEN qc_bl1d(k) = edmf_qc(k)*edmf_a(k)/cldfra_bl1d(k) ELSE cldfra_bl1d(k)=edmf_a(k) qc_bl1d(k) = edmf_qc(k) ENDIF ELSE IF (mf_cf > edmf_a(k)) THEN cldfra_bl1d(k) = mf_cf qc_bl1d(k) = edmf_qc(k)*edmf_a(k)/mf_cf ELSE cldfra_bl1d(k)=edmf_a(k) qc_bl1d(k) = edmf_qc(k) ENDIF ENDIF Fng = 2.05 vt(k) = qww - MIN(0.3,cldfra_bl1D(k))*beta*bb*Fng - 1. vq(k) = alpha + MIN(0.3,cldfra_bl1D(k))*beta*a*Fng - tv0 ENDIF ELSEIF(cldfra_opt == 1) THEN qc_bl1d(k) = MAX(qc_bl1d(k), edmf_qc(k)) if(F_qc .and. .not. F_qi)then satvp = 3.80*exp(17.27*(th(k)-273.)/ & (th(k)-36.))/(.01*p(k)) rhgrid = max(.1,MIN( .95, qv(k) /satvp)) h2oliq=1000.*qc_bl1d(k) satvp=1000.*satvp cldfra_bl1d(k)=(1.-exp(-coef_alph*h2oliq/& ((1.-rhgrid)*satvp)**coef_gamm))*(rhgrid**coef_p) cldfra_bl1d(k)=max(0.0,MIN(1.,cldfra_bl1d(k))) elseif(F_qc .and. F_qi)then satvp = 3.80*exp(17.27*(th(k)-273.)/ & (th(k)-36.))/(.01*p(k)) rhgrid = max(.1,MIN( .95, qv(k) /satvp)) h2oliq=1000.*qc_bl1d(k) satvp=1000.*satvp cldfra_bl1d(k)=(1.-exp(-coef_alph*h2oliq/& ((1.-rhgrid)*satvp)**coef_gamm))*(rhgrid**coef_p) cldfra_bl1d(k)=max(0.0,MIN(1.,cldfra_bl1d(k))) endif ENDIF ENDDO ENDIF IF (ktop > 0) THEN maxqc = maxval(edmf_qc(1:ktop)) IF ( maxqc < 1.E-8) maxmf = -1.0*maxmf ENDIF IF (edmf_w(1) > 4.0) THEN print *,'flq:',flq,' fltv:',fltv print *,'pblh:',pblh,' wstar:',wstar print *,'sigmaW=',sigmaW,' sigmaTH=',sigmaTH,' sigmaQT=',sigmaQT print *,' edmf_a',edmf_a(1:14) print *,' edmf_w',edmf_w(1:14) print *,' edmf_qt:',edmf_qt(1:14) print *,' edmf_thl:',edmf_thl(1:14) ENDIF END SUBROUTINE StEM_MF subroutine Poisson(istart,iend,jstart,jend,mu,POI) integer, intent(in) :: istart,iend,jstart,jend real,dimension(istart:iend,jstart:jend),intent(in) :: MU integer, dimension(istart:iend,jstart:jend), intent(out) :: POI integer :: i,j do i=istart,iend do j=jstart,jend call random_Poisson(mu(i,j),.true.,POI(i,j)) enddo enddo end subroutine Poisson subroutine init_random_seed() implicit none integer, allocatable :: seed(:) integer :: i, n, un, istat, dt(8), pid integer :: t call random_seed(size = n) allocate(seed(n)) un=191 open(unit=un, file="/dev/urandom", access="stream", & form="unformatted", action="read", status="old", iostat=istat) if (istat == 0) then read(un) seed close(un) else call system_clock(t) if (t == 0) then call date_and_time(values=dt) t = dt(6) * 60 & + dt(7) end if pid = 666 + MOD(t,10) t = ieor(t, int(pid, kind(t))) do i = 1, n seed(i) = lcg(t) end do end if call random_seed(put=seed) contains function lcg(s) integer :: lcg integer :: s if (s == 0) then s = 1047 else s = mod(s, 71) end if s = mod(s * 23, 17) lcg = int(mod(s, int(s/3.5))) end function lcg end subroutine init_random_seed subroutine random_Poisson(mu,first,ival) REAL, INTENT(IN) :: mu LOGICAL, INTENT(IN) :: first INTEGER :: ival REAL :: b1_s, b2_s, c, c0, c1_s, c2_s, c3_s, del, difmuk, e, fk, fx, fy, g_s, & omega, px, py, t, u, v, x, xx REAL, SAVE :: s, d, p, q, p0 INTEGER :: j, k, kflag LOGICAL, SAVE :: full_init INTEGER, SAVE :: l, m REAL, SAVE :: pp(35) REAL, PARAMETER :: a0 = -.5, a1_s = .3333333, a2_s = -.2500068, a3 = .2000118, & a4 = -.1661269, a5 = .1421878, a6 = -0.1384794, & a7 = .1250060 REAL, PARAMETER :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040., & 40320., 362880. /) difmuk = 0. fk = 1.0 u = 0. IF (mu > 10.0) THEN IF (first) THEN s = SQRT(mu) d = 6.0*mu*mu l = mu - 1.1484 full_init = .false. END IF g_s = mu + s*random_normal() IF (g_s > 0.0) THEN ival = g_s IF (ival>=l) RETURN fk = ival difmuk = mu - fk CALL RANDOM_NUMBER(u) IF (d*u >= difmuk*difmuk*difmuk) RETURN END IF IF (.NOT. full_init) THEN omega = .3989423/s b1_s = .4166667E-1/mu b2_s = .3*b1_s*b1_s c3_s = .1428571*b1_s*b2_s c2_s = b2_s - 15.*c3_s c1_s = b1_s - 6.*b2_s + 45.*c3_s c0 = 1. - b1_s + 3.*b2_s - 15.*c3_s c = .1069/mu full_init = .true. END IF IF (g_s < 0.0) GO TO 50 kflag = 0 GO TO 70 40 IF (fy-u*fy <= py*EXP(px-fx)) RETURN 50 e = random_exponential() CALL RANDOM_NUMBER(u) u = u + u - one t = 1.8 + SIGN(e, u) IF (t <= (-.6744)) GO TO 50 ival = mu + s*t fk = ival difmuk = mu - fk kflag = 1 GO TO 70 60 IF (c*ABS(u) > py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50 RETURN 70 IF (ival>=10) GO TO 80 px = -mu py = mu**ival/fact(MAX(ival+1,1)) GO TO 110 80 del = .8333333E-1/fk del = del - 4.8*del*del*del v = difmuk/fk IF (ABS(v)>0.25) THEN px = fk*LOG(one + v) - difmuk - del ELSE px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2_s)*v+a1_s)*v+a0) - del END IF py = .3989423/SQRT(fk) 110 x = (half - difmuk)/s xx = x*x fx = -half*xx fy = omega* (((c3_s*xx + c2_s)*xx + c1_s)*xx + c0) IF (kflag <= 0) GO TO 40 GO TO 60 ELSE IF (first) THEN m = MAX(1, INT(mu)) l = 0 p = EXP(-mu) q = p p0 = p END IF DO CALL RANDOM_NUMBER(u) ival = 0 IF (u <= p0) RETURN IF (l == 0) GO TO 150 j = 1 IF (u > 0.458) j = MIN(l, m) DO k = j, l IF (u <= pp(k)) GO TO 180 END DO IF (l == 35) CYCLE 150 l = l + 1 DO k = l, 35 p = p*mu / k q = q + p pp(k) = q IF (u <= q) GO TO 170 END DO l = 35 END DO 170 l = k 180 ival = k RETURN END IF RETURN END subroutine random_Poisson FUNCTION random_normal() RESULT(fn_val) REAL :: fn_val REAL :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472, & r1 = 0.27597, r2 = 0.27846, u, v, x, y, q DO CALL RANDOM_NUMBER(u) CALL RANDOM_NUMBER(v) v = 1.7156 * (v - half) x = u - s y = ABS(v) - t q = x**2 + y*(a*y - b*x) IF (q < r1) EXIT IF (q > r2) CYCLE IF (v**2 < -4.0*LOG(u)*u**2) EXIT END DO fn_val = v/u RETURN END FUNCTION random_normal FUNCTION random_exponential() RESULT(fn_val) REAL :: fn_val REAL :: r DO CALL RANDOM_NUMBER(r) IF (r > zero) EXIT END DO fn_val = -LOG(r) RETURN END FUNCTION random_exponential subroutine condensation_edmf(QT,THL,P,zagl,THV,QC) real,intent(in) :: QT,THL,P,zagl real,intent(out) :: THV real,intent(inout):: QC integer :: niter,i real :: diff,exn,t,th,qs,qcold niter=50 diff=2.e-5 EXN=(P/p1000mb)**rcp do i=1,NITER T=EXN*THL + xlv/cp*QC QS=qsat_blend(T,P) QCOLD=QC QC=0.5*QC + 0.5*MAX((QT-QS),0.) if (abs(QC-QCOLD) 1.0) Psig_bl=1.0 If(Psig_bl < 0.0) Psig_bl=0.0 If(Psig_shcu > 1.0) Psig_shcu=1.0 If(Psig_shcu < 0.0) Psig_shcu=0.0 END SUBROUTINE SCALE_AWARE FUNCTION esat_blend(t) IMPLICIT NONE REAL, INTENT(IN):: t REAL :: esat_blend,XC,ESL,ESI,chi XC=MAX(-80.,t-273.16) IF (t .GE. 273.16) THEN esat_blend = J0+XC*(J1+XC*(J2+XC*(J3+XC*(J4+XC*(J5+XC*(J6+XC*(J7+XC*J8))))))) ELSE IF (t .LE. 253.) THEN esat_blend = K0+XC*(K1+XC*(K2+XC*(K3+XC*(K4+XC*(K5+XC*(K6+XC*(K7+XC*K8))))))) ELSE ESL = J0+XC*(J1+XC*(J2+XC*(J3+XC*(J4+XC*(J5+XC*(J6+XC*(J7+XC*J8))))))) ESI = K0+XC*(K1+XC*(K2+XC*(K3+XC*(K4+XC*(K5+XC*(K6+XC*(K7+XC*K8))))))) chi = (273.16-t)/20.16 esat_blend = (1.-chi)*ESL + chi*ESI END IF END FUNCTION esat_blend FUNCTION qsat_blend(t, P) IMPLICIT NONE REAL, INTENT(IN):: t, P REAL :: qsat_blend,XC,ESL,ESI,RSLF,RSIF,chi XC=MAX(-80.,t-273.16) IF (t .GE. 273.16) THEN ESL = J0+XC*(J1+XC*(J2+XC*(J3+XC*(J4+XC*(J5+XC*(J6+XC*(J7+XC*J8))))))) qsat_blend = 0.622*ESL/(P-ESL) ELSE IF (t .LE. 253.) THEN ESI = K0+XC*(K1+XC*(K2+XC*(K3+XC*(K4+XC*(K5+XC*(K6+XC*(K7+XC*K8))))))) qsat_blend = 0.622*ESI/(P-ESI) ELSE ESL = J0+XC*(J1+XC*(J2+XC*(J3+XC*(J4+XC*(J5+XC*(J6+XC*(J7+XC*J8))))))) ESI = K0+XC*(K1+XC*(K2+XC*(K3+XC*(K4+XC*(K5+XC*(K6+XC*(K7+XC*K8))))))) RSLF = 0.622*ESL/(P-ESL) RSIF = 0.622*ESI/(P-ESI) chi = (273.16-t)/20.16 qsat_blend = (1.-chi)*RSLF + chi*RSIF END IF END FUNCTION qsat_blend FUNCTION xl_blend(t) IMPLICIT NONE REAL, INTENT(IN):: t REAL :: xl_blend,xlvt,xlst,chi IF (t .GE. 273.16) THEN xl_blend = xlv + (cpv-cliq)*(t-273.16) ELSE IF (t .LE. 253.) THEN xl_blend = xls + (cpv-cice)*(t-273.16) ELSE xlvt = xlv + (cpv-cliq)*(t-273.16) xlst = xls + (cpv-cice)*(t-273.16) chi = (273.16-t)/20.16 xl_blend = (1.-chi)*xlvt + chi*xlst END IF END FUNCTION xl_blend SUBROUTINE temf_mf( & & kts,kte,dt,zw,p,pi1d, & & u,v,w,th,thl,thv,qt,qv,qc,& & qke,ust,flt,flq,flqv,flqc,& & hfx,qfx,tsk, & & pblh,rho,dfh,dx,znt,ep_2, & & edmf_a,edmf_w,edmf_qt, & & edmf_thl,edmf_ent,edmf_qc,& & s_aw,s_awthl,s_awqt, & & s_awqv,s_awqc, & & s_awu,s_awv,s_awqke, & & qc_bl1d,cldfra_bl1d, & &F_QC,F_QI,psig, & &spp_pbl,rstoch_col, & &ii,jj,ids,ide,jds,jde) INTEGER, INTENT(IN) :: spp_pbl REAL, DIMENSION(kts:kte), INTENT(in) :: rstoch_col INTEGER, INTENT(IN) :: kts,kte,ii,jj,ids,ide,jds,jde REAL,DIMENSION(kts:kte), INTENT(IN) :: u,v,w,th,thl,qt,qv,qc,thv,p,pi1d REAL,DIMENSION(kts:kte), INTENT(IN) :: qke REAL,DIMENSION(kts:kte+1), INTENT(IN) :: zw REAL,DIMENSION(kts:kte), INTENT(IN) :: rho REAL,DIMENSION(kts:kte), INTENT(IN) :: dfh REAL, INTENT(IN) :: dt,ust,flt,flq,flqv,flqc,hfx,qfx,tsk,pblh,dx,znt,ep_2,psig LOGICAL, OPTIONAL :: f_qc,f_qi REAL,DIMENSION(kts:kte), INTENT(OUT) :: & & edmf_a,edmf_w,edmf_qt, & & edmf_thl,edmf_ent,edmf_qc REAL,DIMENSION(kts:kte+1) :: s_aw, & s_awthl, & s_awqt, & s_awqv, & s_awqc, & s_awu, & s_awv, & s_awqke REAL,DIMENSION(kts:kte), INTENT(INOUT) :: qc_bl1d,cldfra_bl1d real, parameter :: CM = 0.03 real, parameter :: Cdelt = 0.006 real, parameter :: Cw = 0.5 real, parameter :: Cc = 3.0 real, parameter :: lasymp = 200.0 real, parameter :: hmax = 4000.0 integer, parameter :: Nupd = 8 integer :: i, k, kt, nu integer:: h0idx real:: h0 real:: wstr, ang, wm real, dimension( Nupd) :: hd,lcl,hct,ht real:: convection_TKE_surface_src, sfcFTE real:: sfcTHVF real:: z0t integer, dimension( Nupd) :: hdidx,lclidx,hctidx,htidx integer:: hmax_idx integer:: tval real, dimension( kts:kte) :: zm, zt, dzm, dzt real, dimension( kts:kte) :: thetal, qtot real, dimension( kts:kte) :: u_temf, v_temf real, dimension( kts:kte) :: rv, rl, rt real, dimension( kts:kte) :: chi_poisson, gam real, dimension( kts:kte) :: dthdz real, dimension( kts:kte) :: lepsmin real, dimension( kts:kte) :: thetav real, dimension( kts:kte) :: dmoist_qtdz real, dimension( kts:kte) :: B, Bmoist real, dimension( kts:kte, Nupd) :: epsmf, deltmf, dMdz real, dimension( kts:kte, Nupd) :: UUPD, VUPD real, dimension( kts:kte, Nupd) :: thetavUPD, qlUPD, TEUPD real, dimension( kts:kte, Nupd) :: thetavUPDmoist, wUPD_dry real, dimension( kts:kte, Nupd) :: dthUPDdz, dwUPDdz real, dimension( kts:kte, Nupd) :: dwUPDmoistdz real, dimension( kts:kte, Nupd) :: dUUPDdz, dVUPDdz, dTEUPDdz real, dimension( kts:kte, Nupd) :: TUPD, rstUPD, rUPD, rlUPD, qstUPD real, dimension( kts:kte, Nupd) :: MUPD, wUPD, qtUPD, thlUPD, qcUPD real, dimension( kts:kte, Nupd) :: aUPD, cldfraUPD, aUPDt real, dimension( kts:kte) :: N2, S, Ri, beta, ftau, fth, ratio real, dimension( kts:kte) :: TKE, TE2 real, dimension( kts:kte) :: ustrtilde, linv, leps real, dimension( kts:kte) :: km, kh real, dimension( kts:kte) :: Fz, QFK, uwk, vwk real, dimension( kts:kte) :: km_conv, kh_conv, lconv real, dimension( kts:kte) :: alpha2, beta2 real, dimension( kts:kte) :: THVF, buoy_src, srcs real, dimension( kts:kte) :: beta1 real, dimension( kts:kte) :: MFCth real Cepsmf real red_fact real, dimension( kts:kte) :: edmf_u, edmf_v, edmf_qke integer:: bdy_dist,taper_dist real:: taper real, dimension( kts:kte, Nupd) :: & shf_temfx, qf_temfx, uw_temfx, vw_temfx , & mf_temfx real, dimension( Nupd) :: hd_temfx, lcl_temfx, hct_temfx, cfm_temfx logical is_convective real, dimension( kts:kte) :: sigq, qst, satdef real :: sigq2, rst, cldfra_sum, psig_w, maxw IF ( wrf_at_debug_level(3000) ) THEN WRITE ( mynn_message , FMT='(A)' ) & ' MYNN; in TEMF_MF, beginning' CALL wrf_debug ( 0 , mynn_message ) ENDIF s_aw = 0. s_awthl= 0. s_awqt = 0. s_awqv = 0. s_awqc = 0. s_awu = 0. s_awv = 0. s_awqke= 0. edmf_a = 0. edmf_w = 0. edmf_qt= 0. edmf_thl=0. edmf_ent=0. edmf_qc= 0. edmf_u=0. edmf_v=0. edmf_qke=0. z0t = znt do k = kts,kte rv(k) = qv(k) / (1.-qv(k)) rl(k) = qc(k) / (1.-qc(k)) rt(k) = qt(k) thetal(k) = thl(k) qtot(k) = qt(k) thetav(k) = thv(k) end do do k = kts,kte u_temf(k) = u(k) v_temf(k) = v(k) end do k = 1 maxw = 0.0 DO WHILE (ZW(k) < pblh + 500.) maxw = MAX(maxw,ABS(W(k))) k = k+1 ENDDO maxw = MAX(0.,maxw - 0.5) Psig_w = MAX(0.0, 1.0 - maxw/0.5) Psig_w = MIN(Psig_w, Psig) zm(1) = znt dzt(1) = zw(2) - zm(1) zt(1) = (zw(2) - znt) / 2. do kt = kts+1,kte zm(kt) = zw(kt) zt(kt) = (zm(kt) + zw(kt+1)) / 2. dzm(kt) = zt(kt) - zt(kt-1) dzt(kt) = zw(kt+1) - zw(kt) end do dzm(1) = dzm(2) dthdz(1) = (thetal(2)-thetal(1)) / (zt(1) * log10(zm(2)/z0t)) sfcTHVF = hfx/(rho(1)*cp) * (1.+0.608*(qv(1)+qc(1))) + 0.608*thetav(1)*qfx h0idx = 1 h0 = zm(1) lepsmin(kts) = 0. hmax_idx = kte-1 do k = kts+1,kte-1 lepsmin(k) = 0. dthdz(k) = (thetal(k+1) - thetal(k)) / dzt(k) if (thetav(k) > thetav(1) .AND. h0idx .EQ. 1) then if (zm(k) < hmax) then h0idx = k h0 = zm(k) else h0idx = k h0 = hmax end if end if if (zm(k) > hmax) then hmax_idx = min(hmax_idx,k) end if end do dthdz(kte) = dthdz(kte-1) if ( hfx > 0.) then wstr = (g * h0 / thetav(2) * hfx/(rho(1)*cp) ) ** (1./3.) bdy_dist = min( min((ii-ids),(ide-ii)) , min((jj-jds),(jde-jj)) ) taper_dist = 5 if (bdy_dist .LE. taper_dist) then taper = max(0., min( 1., real(bdy_dist) / real(taper_dist) ) ) wstr = wstr * taper end if else wstr = 0. end if IF ( wrf_at_debug_level(3000) ) THEN WRITE ( mynn_message , FMT='(A,F5.1,F6.1,F5.1,F5.1,F5.1)' ) & ' MYNN; in TEMF_MF: wstr,hfx,dtdz1,dtdz2,h0:', wstr,hfx,dthdz(1),dthdz(2),h0 CALL wrf_debug ( 0 , mynn_message ) ENDIF is_convective = wstr > 0. .AND. dthdz(1)<0. .AND. dthdz(2)<0. if ( is_convective) then IF ( wrf_at_debug_level(3000) ) THEN WRITE ( mynn_message , FMT='(A)' ) & ' MYNN; in TEMF_MF: inconvective branch' CALL wrf_debug ( 0 , mynn_message ) ENDIF Cepsmf = 1.0 / max(200.,h0) Cepsmf = max(Cepsmf,0.0010) do nu = 1,Nupd do k = kts,kte epsmf(k,nu) = Cepsmf * (1+0.25*(nu-1)) end do thlUPD(1,nu) = thetal(1) + Cw*wstr qtUPD(1,nu) = qtot(1) + 0.0*qfx/wstr rUPD(1,nu) = qtUPD(1,nu) / (1. - qtUPD(1,nu)) wUPD(1,nu) = Cw * wstr wUPD_dry(1,nu) = Cw * wstr UUPD(1,nu) = u_temf(1) VUPD(1,nu) = v_temf(1) thetavUPD(1,nu) = thlUPD(1,nu) * (1. + 0.608*qtUPD(1,nu)) thetavUPDmoist(1,nu) = thetavUPD(1,nu) TEUPD(1,nu) = qke(1) + g / thetav(1) * sfcTHVF qlUPD(1,nu) = qc(1) TUPD(1,nu) = thlUPD(1,nu) * pi1d(1) rstUPD(1,nu) = qsat_blend(TUPD(1,nu),p(1)) rlUPD(1,nu) = 0. do k = 2,kte if ( k < hmax_idx) then dthUPDdz(k-1,nu) = -epsmf(k,nu) * (thlUPD(k-1,nu) - thetal(k-1)) thlUPD(k,nu) = thlUPD(k-1,nu) + dthUPDdz(k-1,nu) * dzm(k-1) dmoist_qtdz(k-1) = -epsmf(k,nu) * (qtUPD(k-1,nu) - qtot(k-1)) qtUPD(k,nu) = qtUPD(k-1,nu) + dmoist_qtdz(k-1) * dzm(k-1) thetavUPD(k,nu) = thlUPD(k,nu) * (1. + 0.608*qtUPD(k,nu)) B(k-1) = g * (thetavUPD(k,nu) - thetav(k)) / thetav(k) if ( wUPD_dry(k-1,nu) < 1e-15 ) then wUPD_dry(k,nu) = 0. else dwUPDdz(k-1,nu) = -2. *epsmf(k,nu)*wUPD_dry(k-1,nu) + 0.33*B(k-1)/wUPD_dry(k-1,nu) wUPD_dry(k,nu) = wUPD_dry(k-1,nu) + dwUPDdz(k-1,nu) * dzm(k-1) end if dUUPDdz(k-1,nu) = -epsmf(k,nu) * (UUPD(k-1,nu) - u_temf(k-1)) UUPD(k,nu) = UUPD(k-1,nu) + dUUPDdz(k-1,nu) * dzm(k-1) dVUPDdz(k-1,nu) = -epsmf(k,nu) * (VUPD(k-1,nu) - v_temf(k-1)) VUPD(k,nu) = VUPD(k-1,nu) + dVUPDdz(k-1,nu) * dzm(k-1) dTEUPDdz(k-1,nu) = -epsmf(k,nu) * (TEUPD(k-1,nu) - qke(k-1)) TEUPD(k,nu) = TEUPD(k-1,nu) + dTEUPDdz(k-1,nu) * dzm(k-1) rUPD(k,nu) = qtUPD(k,nu) / (1. - qtUPD(k,nu)) TUPD(k,nu) = thlUPD(k,nu) * pi1d(k) rstUPD(k,nu) = qsat_blend(TUPD(k,nu),p(k-1)) beta1(k) = 0.622 * (xlv/(r_d*TUPD(k,nu))) * (xlv/(cp*TUPD(k,nu))) rstUPD(k,nu) = rstUPD(k,nu) * (1.0+beta1(k)*rUPD(k,nu)) / (1.0+beta1(k)*rstUPD(k,nu)) qstUPD(k,nu) = rstUPD(k,nu) / (1. + rstUPD(k,nu)) if (rUPD(k,nu) > rstUPD(k,nu)) then rlUPD(k,nu) = rUPD(k,nu) - rstUPD(k,nu) qlUPD(k,nu) = rlUPD(k,nu) / (1. + rlUPD(k,nu)) thetavUPDmoist(k,nu) = (thlUPD(k,nu) + ((xlv/cp)*qlUPD(k,nu)/pi1d(k))) * & (1. + 0.608*qstUPD(k,nu) - qlUPD(k,nu)) else rlUPD(k,nu) = 0. qlUPD(k,nu) = qc(k-1) thetavUPDmoist(k,nu) = thlUPD(k,nu) * (1. + 0.608*qtUPD(k,nu)) end if Bmoist(k-1) = g * (thetavUPDmoist(k,nu) - thetav(k)) / thetav(k) if ( wUPD(k-1,nu) < 1e-15 ) then wUPD(k,nu) = 0. else dwUPDmoistdz(k-1,nu) = -2. *epsmf(k,nu)*wUPD(k-1,nu) + 0.33*Bmoist(k-1)/wUPD(k-1,nu) wUPD(k,nu) = wUPD(k-1,nu) + dwUPDmoistdz(k-1,nu) * dzm(k-1) end if else thlUPD(k,nu) = thetal(k) qtUPD(k,nu) = qtot(k) wUPD_dry(k,nu) = 0. UUPD(k,nu) = u_temf(k) VUPD(k,nu) = v_temf(k) TEUPD(k,nu) = qke(k) qlUPD(k,nu) = qc(k-1) wUPD(k,nu) = 0. end if IF ( wrf_at_debug_level(3000) ) THEN IF ( ABS(wUPD(k,nu))>10. ) THEN WRITE ( mynn_message , FMT='(A,2I3)' ) & ' MYNN, in TEMF_MF, huge w at (nu,k):', nu,k CALL wrf_debug ( 0 , mynn_message ) print *," thlUPD(1:k,nu) = ", thlUPD(1:k,nu) print *," wUPD(1:k,nu) = ", wUPD(1:k,nu) print *," Bmoist(1:k-1) = ", Bmoist(1:k-1) print *," epsmf(1:k,nu) = ", epsmf(1:k,nu) ENDIF ENDIF ENDDO if (wUPD_dry(1,nu) == 0.) then hdidx(nu) = 1 else hdidx(nu) = kte do k = 2,kte if (wUPD_dry(k,nu) <= 0. .OR. zm(k) > hmax) then hdidx(nu) = k exit end if end do end if 100 hd(nu) = zm(hdidx(nu)) lclidx(nu) = kte do k = kts,kte if ( k < hmax_idx .AND. rUPD(k,nu) > rstUPD(k,nu)) then lclidx(nu) = k exit end if end do 200 lcl(nu) = zm(lclidx(nu)) if (hd(nu) > lcl(nu)) then if (wUPD(1,nu) == 0.) then hctidx(nu) = 1 else hctidx(nu) = kte do k = 2,kte if (wUPD(k,nu) <= 0. .OR. zm(k) > hmax) then hctidx(nu) = k exit end if end do end if 300 hct(nu) = zm(hctidx(nu)) if (hctidx(nu) <= hdidx(nu)+1) then hct(nu) = hd(nu) hctidx(nu) = hdidx(nu) else end if else hct(nu) = hd(nu) hctidx(nu) = hdidx(nu) end if ht(nu) = max(hd(nu),hct(nu)) htidx(nu) = max(hdidx(nu),hctidx(nu)) do k = 1,kte if (zm(k) < 0.9*ht(nu)) then tval = 1 else if (zm(k) >= 0.9*ht(nu) .AND. zm(k) <= 1.0*ht(nu)) then tval = 1. - ((zm(k) - 0.9*ht(nu)) / (1.0*ht(nu) - 0.9*ht(nu))) else tval = 0. end if thlUPD(k,nu) = tval * thlUPD(k,nu) + (1-tval)*thetal(k) thetavUPD(k,nu) = tval * thetavUPD(k,nu) + (1-tval)*thetav(k) qtUPD(k,nu) = tval * qtUPD(k,nu) + (1-tval) * qtot(k) if (k > 1) then qlUPD(k,nu) = tval * qlUPD(k,nu) + (1-tval) * qc(k-1) end if UUPD(k,nu) = tval * UUPD(k,nu) + (1-tval) * u_temf(k) VUPD(k,nu) = tval * VUPD(k,nu) + (1-tval) * v_temf(k) TEUPD(k,nu) = tval * TEUPD(k,nu) + (1-tval) * qke(k) if (zm(k) > ht(nu)) then wUPD(k,nu) = 0. dwUPDmoistdz(k,nu) = 0. wUPD_dry(k,nu) = 0. dwUPDdz(k,nu) = 0. end if end do deltmf(:,nu) = epsmf(:,nu) mf_temfx(1,nu) = CM * wstr / Nupd mf_temfx(1,nu) = min(mf_temfx(1,nu),0.2/Nupd) do kt = 2,kte-1 dMdz(kt,nu) = (epsmf(kt,nu) - deltmf(kt,nu)) * mf_temfx(kt-1,nu) * dzt(kt) mf_temfx(kt,nu) = mf_temfx(kt-1,nu) + dMdz(kt,nu) mf_temfx(kt,nu) = max(mf_temfx(kt,nu),0.0) IF ( wrf_at_debug_level(3000) ) THEN IF ( mf_temfx(kt,nu)>=0.2/NUPD ) THEN WRITE ( mynn_message , FMT='(A,2I3)' ) & ' MYNN, in TEMF_MF, huge MF at (nu,k):', nu,kt CALL wrf_debug ( 0 , mynn_message ) print *," mf_temfx(1:kt,nu) = ", mf_temfx(1:kt,nu) ENDIF ENDIF end do mf_temfx(kte,nu) = 0. aUPD(1,nu) = 0.06 / Nupd do k = 2,kte if (wUPD(k-1,nu) >= 1.0e-5 .AND. wUPD(k,nu) >= 1.0e-5) then aUPD(k,nu) = ((mf_temfx(k-1,nu)+mf_temfx(k,nu))/2.0) / & ((wUPD(k-1,nu)+wUPD(k,nu))/2.0) else aUPD(k,nu) = 0.0 end if sigq2 = aUPD(k,nu) * (qtUPD(k,nu)-qtot(k)) if (sigq2 > 0.0) then sigq(k) = sqrt(sigq2) else sigq(k) = 0.0 end if rst = qsat_blend(th(k-1)*pi1d(k-1),p(k-1)) qst(k) = rst / (1. + rst) satdef(k) = qtot(k) - qst(k) if (satdef(k) <= 0.0) then if (sigq(k) > 1.0e-15) then cldfraUPD(k,nu) = max(0.5 + 0.36 * atan(1.55*(satdef(k)/sigq(k))),0.0) / Nupd else cldfraUPD(k,nu) = 0.0 end if else cldfraUPD(k,nu) = 1.0 / Nupd end if if (zm(k) < lcl(nu)) then cldfraUPD(k,nu) = 0.0 end if end do end do cfm_temfx = 0.0 do k = 2,kte cldfra_sum = 0.0 edmf_a(k) = 0.0 edmf_w(k) = 0.0 edmf_thl(k) = 0.0 edmf_qt(k) = 0.0 edmf_qc(k) = 0.0 edmf_u(k) = 0.0 edmf_v(k) = 0.0 edmf_qke(k) = 0.0 edmf_ent(k) = 0.0 do nu = 1,Nupd aUPDt(k,nu) = mf_temfx(k,nu) / wUPD(k,nu) if (aUPDt(k,nu) >= 1.0e-3 .AND. wUPD(k,nu) >= 1.0e-5) then edmf_a(k) = edmf_a(k) + aUPDt(k,nu) edmf_w(k) = edmf_w(k) + aUPDt(k,nu)*wUPD(k,nu) edmf_thl(k) = edmf_thl(k) + aUPDt(k,nu)*thlUPD(k,nu) edmf_qt(k) = edmf_qt(k) + aUPDt(k,nu)*qtUPD(k,nu) edmf_qc(k) = edmf_qc(k) + aUPDt(k,nu)*qlUPD(k,nu) edmf_u(k) = edmf_u(k) + aUPDt(k,nu)*UUPD(k,nu) edmf_v(k) = edmf_v(k) + aUPDt(k,nu)*VUPD(k,nu) edmf_qke(k) = edmf_qke(k) + aUPDt(k,nu)*TEUPD(k,nu) edmf_ent(k) = edmf_ent(k) + aUPDt(k,nu)*epsmf(k,nu) cldfra_sum = cldfra_sum + cldfraUPD(k,nu) end if end do IF ( wrf_at_debug_level(3000) ) THEN ENDIF if (edmf_a(k)>1.e-3) then edmf_w(k)=edmf_w(k)/edmf_a(k) edmf_qt(k)=edmf_qt(k)/edmf_a(k) edmf_thl(k)=edmf_thl(k)/edmf_a(k) edmf_ent(k)=edmf_ent(k)/edmf_a(k) edmf_qc(k)=edmf_qc(k)/edmf_a(k) edmf_u(k)=edmf_u(k)/edmf_a(k) edmf_v(k)=edmf_v(k)/edmf_a(k) edmf_qke(k)=edmf_qke(k)/edmf_a(k) if (edmf_qc(k) > 0.0) then IF (cldfra_sum > edmf_a(k)) THEN cldfra_bl1d(k) = cldfra_sum qc_bl1d(k) = edmf_qc(k)*edmf_a(k)/cldfra_sum ELSE cldfra_bl1d(k)=edmf_a(k) qc_bl1d(k) = edmf_qc(k) ENDIF endif endif if (zt(k) <= hmax) then cfm_temfx = max(cldfra_bl1d(k),cfm_temfx) end if end do do k=kts,kte if (edmf_a(k)>1.0e-3) then s_aw(k) = edmf_a(k)*edmf_w(k)*psig_w * (1.0+rstoch_col(k)) s_awthl(k)= edmf_a(k)*edmf_w(k)*edmf_thl(k)*psig_w * (1.0+rstoch_col(k)) s_awqt(k) = edmf_a(k)*edmf_w(k)*edmf_qt(k)*psig_w * (1.0+rstoch_col(k)) s_awqc(k) = edmf_a(k)*edmf_w(k)*edmf_qc(k)*psig_w * (1.0+rstoch_col(k)) s_awqv(k) = s_awqt(k) - s_awqc(k) s_awu(k) = edmf_a(k)*edmf_w(k)*edmf_u(k)*psig_w * (1.0+rstoch_col(k)) s_awv(k) = edmf_a(k)*edmf_w(k)*edmf_v(k)*psig_w * (1.0+rstoch_col(k)) s_awqke(k) = edmf_a(k)*edmf_w(k)*edmf_qke(k)*psig_w * (1.0+rstoch_col(k)) endif edmf_a(k)=edmf_a(k)*psig_w enddo else edmf_a = 0. edmf_w = 0. edmf_qt = 0. edmf_thl = 0. edmf_ent = 0. edmf_u = 0. edmf_v = 0. edmf_qke = 0. s_aw = 0. s_awthl= 0. s_awqt = 0. s_awqv = 0. s_awqc = 0. s_awu = 0. s_awv = 0. s_awqke= 0. edmf_qc(1) = qc(1) do k = kts+1,kte-1 edmf_qc(k) = qc(k-1) end do end if END SUBROUTINE temf_mf real function rsat_temf(p,T,ep2) implicit none real p, T, ep2 real temp, x real, parameter :: c0 = 0.6105851e+3 real, parameter :: c1 = 0.4440316e+2 real, parameter :: c2 = 0.1430341e+1 real, parameter :: c3 = 0.2641412e-1 real, parameter :: c4 = 0.2995057e-3 real, parameter :: c5 = 0.2031998e-5 real, parameter :: c6 = 0.6936113e-8 real, parameter :: c7 = 0.2564861e-11 real, parameter :: c8 = -0.3704404e-13 temp = T - 273.15 x =c0+temp*(c1+temp*(c2+temp*(c3+temp*(c4+temp*(c5+temp*(c6+temp*(c7+temp*c8))))))) rsat_temf = ep2*x/(p-x) return end function rsat_temf END MODULE module_bl_mynn