module bl_ysu use ccpp_kind_types,only: kind_phys implicit none private public:: bl_ysu_run, & bl_ysu_init, & bl_ysu_finalize contains subroutine bl_ysu_init(errmsg,errflg) character(len=*),intent(out):: errmsg integer,intent(out):: errflg errmsg = 'bl_ysu_init OK' errflg = 0 end subroutine bl_ysu_init subroutine bl_ysu_finalize(errmsg,errflg) character(len=*),intent(out):: errmsg integer,intent(out):: errflg errmsg = 'bl_ysu_finalize OK' errflg = 0 end subroutine bl_ysu_finalize subroutine bl_ysu_run(ux,vx,tx,qvx,qcx,qix,nmix,qmix,p2d,p2di,pi2d, & f_qc,f_qi, & utnp,vtnp,ttnp,qvtnp,qctnp,qitnp,qmixtnp, & cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, & dz8w2d,psfcpa, & znt,ust,hpbl,dusfc,dvsfc,dtsfc,dqsfc,psim,psih, & xland,hfx,qfx,wspd,br, & dt,kpbl1d, & exch_hx,exch_mx, & wstar,delta, & u10,v10, & uox,vox, & rthraten, & ysu_topdown_pblmix, & ctopo,ctopo2, & a_u,a_v,a_t,a_q,a_e, & b_u,b_v,b_t,b_q,b_e, & sfk,vlk,dlu,dlg,frcurb, & flag_bep, & its,ite,kte,kme, & errmsg,errflg & ) implicit none real(kind=kind_phys),parameter :: xkzminm = 0.1,xkzminh = 0.01 real(kind=kind_phys),parameter :: xkzmin = 0.01,xkzmax = 1000.,rimin = -100. real(kind=kind_phys),parameter :: rlam = 30.,prmin = 0.25,prmax = 4. real(kind=kind_phys),parameter :: brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4 real(kind=kind_phys),parameter :: afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0 real(kind=kind_phys),parameter :: phifac = 8.,sfcfrac = 0.1 real(kind=kind_phys),parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001 real(kind=kind_phys),parameter :: h1 = 0.33333335, h2 = 0.6666667 real(kind=kind_phys),parameter :: zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16. real(kind=kind_phys),parameter :: tmin=1.e-2 real(kind=kind_phys),parameter :: gamcrt = 3.,gamcrq = 2.e-3 real(kind=kind_phys),parameter :: xka = 2.4e-5 integer,parameter :: imvdif = 1 real(kind=kind_phys),parameter :: rcl = 1.0 integer,parameter :: kts=1, kms=1 integer, intent(in ) :: its,ite,kte,kme logical, intent(in) :: ysu_topdown_pblmix integer, intent(in) :: nmix real(kind=kind_phys), intent(in ) :: dt,cp,g,rovcp,rovg,rd,xlv,rv real(kind=kind_phys), intent(in ) :: ep1,ep2,karman logical, intent(in ) :: f_qc, f_qi real(kind=kind_phys), dimension( its:,: ) , & intent(in) :: dz8w2d, & pi2d real(kind=kind_phys), dimension( its:,: ) , & intent(in ) :: tx, & qvx, & qcx, & qix real(kind=kind_phys), dimension( its:,:,: ) , & intent(in ) :: qmix real(kind=kind_phys), dimension( its:,: ) , & intent(out ) :: utnp, & vtnp, & ttnp, & qvtnp, & qctnp, & qitnp real(kind=kind_phys), dimension( its:,:,: ) , & intent(out ) :: qmixtnp real(kind=kind_phys), dimension( its:,: ) , & intent(in ) :: p2di real(kind=kind_phys), dimension( its:,: ) , & intent(in ) :: p2d real(kind=kind_phys), dimension( its: ) , & intent(out ) :: hpbl real(kind=kind_phys), dimension( its: ) , & intent(out ), optional :: dusfc, & dvsfc, & dtsfc, & dqsfc real(kind=kind_phys), dimension( its: ) , & intent(in ) :: ust, & znt real(kind=kind_phys), dimension( its: ) , & intent(in ) :: xland, & hfx, & qfx real(kind=kind_phys), dimension( its: ), intent(in ) :: wspd real(kind=kind_phys), dimension( its: ), intent(in ) :: br real(kind=kind_phys), dimension( its: ), intent(in ) :: psim, & psih real(kind=kind_phys), dimension( its: ), intent(in ) :: psfcpa integer, dimension( its: ), intent(out ) :: kpbl1d real(kind=kind_phys), dimension( its:,: ) , & intent(in ) :: ux, & vx, & rthraten real(kind=kind_phys), dimension( its: ) , & optional , & intent(in ) :: ctopo, & ctopo2 logical, intent(in ) :: flag_bep real(kind=kind_phys), dimension( its:,: ) , & optional , & intent(in ) :: a_u, & a_v, & a_t, & a_q, & a_e, & b_u, & b_v, & b_t, & b_q, & b_e, & sfk, & vlk, & dlu, & dlg real(kind=kind_phys), dimension( its: ) , & optional , & intent(in ) :: frcurb character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg real(kind=kind_phys), dimension( its:ite ) :: hol real(kind=kind_phys), dimension( its:ite, kms:kme ) :: zq real(kind=kind_phys), dimension( its:ite, kts:kte ) :: & thx,thvx,thlix, & del, & dza, & dzq, & xkzom, & xkzoh, & za real(kind=kind_phys), dimension( its:ite ) :: & rhox, & govrth, & zl1,thermal, & wscale, & hgamt,hgamq, & brdn,brup, & phim,phih, & prpbl, & wspd1,thermalli real(kind=kind_phys), dimension( its:ite, kts:kte ) :: xkzh,xkzm,xkzq, & f1,f2, & r1,r2, & ad,au, & cu, & al, & zfac, & rhox2, & hgamt2, & ad1,adm,adv real(kind=kind_phys), dimension( its:ite, kts:kte ) , & intent(out ) :: exch_hx, & exch_mx real(kind=kind_phys), dimension( its:ite ) , & intent(inout) :: u10, & v10 real(kind=kind_phys), dimension( its:ite ), optional , & intent(in ) :: uox, & vox real(kind=kind_phys), dimension( its:ite ) :: uoxl, & voxl real(kind=kind_phys), dimension( its:ite ) :: & brcr, & sflux, & zol1, & brcr_sbro real(kind=kind_phys), dimension( its:ite, kts:kte) :: r3,f3 integer, dimension( its:ite ) :: kpbl,kpblold logical, dimension( its:ite ) :: pblflg, & sfcflg, & stable, & cloudflg logical :: definebrup integer :: n,i,k,l,ic,is,kk integer :: klpbl real(kind=kind_phys) :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0 real(kind=kind_phys) :: ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri real(kind=kind_phys) :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz real(kind=kind_phys) :: utend,vtend,ttend,qtend real(kind=kind_phys) :: dtstep,govrthv real(kind=kind_phys) :: cont, conq, conw, conwrc real(kind=kind_phys), dimension( its:ite, kts:kte ) :: wscalek,wscalek2 real(kind=kind_phys), dimension( its:ite ), intent(out) :: wstar, & delta real(kind=kind_phys), dimension( its:ite, kts:kte ) :: xkzml,xkzhl, & zfacent,entfac real(kind=kind_phys), dimension( its:ite, kts:kte ) :: qcxl, & qixl real(kind=kind_phys), dimension( its:ite ) :: ust3, & wstar3, & wstar3_2, & hgamu,hgamv, & wm2, we, & bfxpbl, & hfxpbl,qfxpbl, & ufxpbl,vfxpbl, & dthvx real(kind=kind_phys) :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, & dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr, & prfac,prfac2,phim8z,radsum,tmp1,templ,rvls,temps,ent_eff, & rcldb,bruptmp,radflux,vconvlim,vconvnew,fluxc,vconvc,vconv real(kind=kind_phys), dimension( its:ite, kts:kte ) :: fric, & tke_ysu,& el_ysu,& shear_ysu,& buoy_ysu real(kind=kind_phys), dimension( its:ite) :: pblh_ysu,& vconvfx real(kind=kind_phys) :: bepswitch real(kind=kind_phys), dimension( its:ite, kts:kte ) :: & a_u2d,a_v2d,a_t2d,a_q2d,a_e2d,b_u2d,b_v2d,b_t2d,b_q2d,b_e2d, & sfk2d,vlk2d,dlu2d,dlg2d real(kind=kind_phys), dimension( its:ite ) :: & frc_urb1d real(kind=kind_phys), dimension( kts:kte ) :: thvx_1d,tke_1d,dzq_1d real(kind=kind_phys), dimension( kts:kte+1) :: zq_1d klpbl = kte cont=cp/g conq=xlv/g conw=1./g conwrc = conw*sqrt(rcl) conpr = bfac*karman*sfcfrac if(f_qc) then do k = kts,kte do i = its,ite qcxl(i,k) = qcx(i,k) enddo enddo else do k = kts,kte do i = its,ite qcxl(i,k) = 0. enddo enddo endif if(f_qi) then do k = kts,kte do i = its,ite qixl(i,k) = qix(i,k) enddo enddo else do k = kts,kte do i = its,ite qixl(i,k) = 0. enddo enddo endif do k = kts,kte do i = its,ite thx(i,k) = tx(i,k)/pi2d(i,k) thlix(i,k) = (tx(i,k)-xlv*qcxl(i,k)/cp-2.834E6*qixl(i,k)/cp)/pi2d(i,k) enddo enddo do k = kts,kte do i = its,ite tvcon = (1.+ep1*qvx(i,k)) thvx(i,k) = thx(i,k)*tvcon enddo enddo if ( present(uox) .and. present(vox) ) then do i =its,ite uoxl(i) = uox(i) voxl(i) = vox(i) enddo else do i =its,ite uoxl(i) = 0 voxl(i) = 0 enddo endif do i = its,ite tvcon = (1.+ep1*qvx(i,1)) rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon) govrth(i) = g/thx(i,1) enddo if(present(a_u) .and. present(a_v) .and. present(a_t) .and. & present(a_q) .and. present(a_t) .and. present(a_e) .and. & present(b_u) .and. present(b_v) .and. present(b_t) .and. & present(b_q) .and. present(b_e) .and. present(dlg) .and. & present(dlu) .and. present(sfk) .and. present(vlk) .and. & present(frcurb) .and. flag_bep) then bepswitch=1.0 do k = kts, kte do i = its,ite a_u2d(i,k) = a_u(i,k) a_v2d(i,k) = a_v(i,k) a_t2d(i,k) = a_t(i,k) a_q2d(i,k) = a_q(i,k) a_e2d(i,k) = a_e(i,k) b_u2d(i,k) = b_u(i,k) b_v2d(i,k) = b_v(i,k) b_t2d(i,k) = b_t(i,k) b_q2d(i,k) = b_q(i,k) b_e2d(i,k) = b_e(i,k) dlg2d(i,k) = dlg(i,k) dlu2d(i,k) = dlu(i,k) vlk2d(i,k) = vlk(i,k) sfk2d(i,k) = sfk(i,k) enddo enddo do i = its, ite frc_urb1d(i) = frcurb(i) enddo else bepswitch=0.0 do k = kts, kte do i = its,ite a_u2d(i,k) = 0.0 a_v2d(i,k) = 0.0 a_t2d(i,k) = 0.0 a_q2d(i,k) = 0.0 a_e2d(i,k) = 0.0 b_u2d(i,k) = 0.0 b_v2d(i,k) = 0.0 b_t2d(i,k) = 0.0 b_q2d(i,k) = 0.0 b_e2d(i,k) = 0.0 dlg2d(i,k) = 0.0 dlu2d(i,k) = 0.0 vlk2d(i,k) = 1.0 sfk2d(i,k) = 1.0 enddo enddo do i = its, ite frc_urb1d(i) = 0.0 enddo endif do i = its,ite zq(i,1) = 0. enddo do k = kts,kte do i = its,ite zq(i,k+1) = dz8w2d(i,k)+zq(i,k) tvcon = (1.+ep1*qvx(i,k)) rhox2(i,k) = p2d(i,k)/(rd*tx(i,k)*tvcon) enddo enddo do k = kts,kte do i = its,ite za(i,k) = 0.5*(zq(i,k)+zq(i,k+1)) dzq(i,k) = zq(i,k+1)-zq(i,k) del(i,k) = p2di(i,k)-p2di(i,k+1) enddo enddo do i = its,ite dza(i,1) = za(i,1) enddo do k = kts+1,kte do i = its,ite dza(i,k) = za(i,k)-za(i,k-1) enddo enddo do k = kts,kte do i = its,ite exch_hx(i,k) = 0. exch_mx(i,k) = 0. xkzh(i,k) = 0. xkzhl(i,k) = 0. xkzm(i,k) = 0. xkzml(i,k) = 0. xkzq(i,k) = 0. enddo enddo do i = its,ite wspd1(i) = sqrt( (ux(i,1)-uoxl(i))*(ux(i,1)-uoxl(i)) + (vx(i,1)-voxl(i))*(vx(i,1)-voxl(i)) )+1.e-9 enddo dtstep = dt dt2 = 2.*dtstep rdt = 1./dt2 do i = its,ite bfxpbl(i) = 0.0 hfxpbl(i) = 0.0 qfxpbl(i) = 0.0 ufxpbl(i) = 0.0 vfxpbl(i) = 0.0 hgamu(i) = 0.0 hgamv(i) = 0.0 delta(i) = 0.0 wstar3_2(i) = 0.0 enddo do k = kts,klpbl do i = its,ite wscalek(i,k) = 0.0 wscalek2(i,k) = 0.0 enddo enddo do k = kts,klpbl do i = its,ite zfac(i,k) = 0.0 enddo enddo do k = kts,klpbl-1 do i = its,ite xkzom(i,k) = xkzminm xkzoh(i,k) = xkzminh enddo enddo do i = its,ite if(present(dusfc)) dusfc(i) = 0. if(present(dvsfc)) dvsfc(i) = 0. if(present(dtsfc)) dtsfc(i) = 0. if(present(dqsfc)) dqsfc(i) = 0. enddo do i = its,ite hgamt(i) = 0. hgamq(i) = 0. wscale(i) = 0. we(i) = 0. kpbl(i) = 1 hpbl(i) = zq(i,1) zl1(i) = za(i,1) thermal(i)= thvx(i,1) thermalli(i) = thlix(i,1) pblflg(i) = .true. sfcflg(i) = .true. sflux(i) = hfx(i)/rhox(i)/cp + qfx(i)/rhox(i)*ep1*thx(i,1) if(br(i).gt.0.0) sfcflg(i) = .false. enddo do i = its,ite stable(i) = .false. brup(i) = br(i) brcr(i) = brcr_ub enddo do k = 2,klpbl do i = its,ite if(.not.stable(i))then brdn(i) = brup(i) spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.) brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2 kpbl(i) = k stable(i) = brup(i).gt.brcr(i) endif enddo enddo do i = its,ite k = kpbl(i) if(brdn(i).ge.brcr(i))then brint = 0. elseif(brup(i).le.brcr(i))then brint = 1. else brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i)) endif hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1)) if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1 if(kpbl(i).le.1) pblflg(i) = .false. enddo do i = its,ite fm = psim(i) fh = psih(i) zol1(i) = max(br(i)*fm*fm/fh,rimin) if(sfcflg(i))then zol1(i) = min(zol1(i),-zfmin) else zol1(i) = max(zol1(i),zfmin) endif hol1 = zol1(i)*hpbl(i)/zl1(i)*sfcfrac if(sfcflg(i))then phim(i) = (1.-aphi16*hol1)**(-1./4.) phih(i) = (1.-aphi16*hol1)**(-1./2.) bfx0 = max(sflux(i),0.) hfx0 = max(hfx(i)/rhox(i)/cp,0.) qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.) wstar3(i) = (govrth(i)*bfx0*hpbl(i)) wstar(i) = (wstar3(i))**h1 else phim(i) = (1.+aphi5*hol1) phih(i) = phim(i) wstar(i) = 0. wstar3(i) = 0. endif ust3(i) = ust(i)**3. wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1 wscale(i) = min(wscale(i),ust(i)*aphi16) wscale(i) = max(wscale(i),ust(i)/aphi5) enddo do i = its,ite if(sfcflg(i).and.sflux(i).gt.0.0)then gamfac = bfac/rhox(i)/wscale(i) hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt) hgamq(i) = min(gamfac*qfx(i),gamcrq) vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0) thermalli(i)= thermalli(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0) hgamt(i) = max(hgamt(i),0.0) hgamq(i) = max(hgamq(i),0.0) brint = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.) hgamu(i) = brint*ux(i,1) hgamv(i) = brint*vx(i,1) else pblflg(i) = .false. endif enddo do i = its,ite if(pblflg(i))then kpbl(i) = 1 hpbl(i) = zq(i,1) endif enddo do i = its,ite if(pblflg(i))then stable(i) = .false. brup(i) = br(i) brcr(i) = brcr_ub endif enddo do k = 2,klpbl do i = its,ite if(.not.stable(i).and.pblflg(i))then brdn(i) = brup(i) spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.) brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2 kpbl(i) = k stable(i) = brup(i).gt.brcr(i) endif enddo enddo if (ysu_topdown_pblmix)then do i = its,ite kpblold(i) = kpbl(i) definebrup=.false. do k = kpblold(i), kte-1 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.) bruptmp = (thlix(i,k)-thermalli(i))*(g*za(i,k)/thlix(i,1))/spdk2 stable(i) = bruptmp.ge.brcr(i) if (definebrup) then kpbl(i) = k brup(i) = bruptmp definebrup=.false. endif if (.not.stable(i)) then brdn(i)=bruptmp definebrup=.true. pblflg(i)=.true. endif enddo enddo endif do i = its,ite if(pblflg(i)) then k = kpbl(i) if(brdn(i).ge.brcr(i))then brint = 0. elseif(brup(i).le.brcr(i))then brint = 1. else brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i)) endif hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1)) if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1 if(kpbl(i).le.1) pblflg(i) = .false. endif enddo do i = its,ite if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then brup(i) = br(i) stable(i) = .false. else stable(i) = .true. endif enddo do i = its,ite if((.not.stable(i)).and.((xland(i)-1.5).ge.0))then wspd10 = u10(i)*u10(i) + v10(i)*v10(i) wspd10 = sqrt(wspd10) ross = wspd10 / (cori*znt(i)) brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3) endif enddo do i = its,ite if(.not.stable(i))then if((xland(i)-1.5).ge.0)then brcr(i) = brcr_sbro(i) else brcr(i) = brcr_sb endif endif enddo do k = 2,klpbl do i = its,ite if(.not.stable(i))then brdn(i) = brup(i) spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.) brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2 kpbl(i) = k stable(i) = brup(i).gt.brcr(i) endif enddo enddo do i = its,ite if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then k = kpbl(i) if(brdn(i).ge.brcr(i))then brint = 0. elseif(brup(i).le.brcr(i))then brint = 1. else brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i)) endif hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1)) if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1 if(kpbl(i).le.1) pblflg(i) = .false. endif enddo do i = its,ite cloudflg(i)=.false. if(pblflg(i)) then k = kpbl(i) - 1 wm3 = wstar3(i) + 5. * ust3(i) wm2(i) = wm3**h2 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i) dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin) we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i))) if((qcxl(i,k)+qixl(i,k)).gt.0.01e-3.and.ysu_topdown_pblmix)then if ( kpbl(i) .ge. 2) then cloudflg(i)=.true. templ=thlix(i,k)*(p2di(i,k+1)/100000)**rovcp rvls=100.*6.112*EXP(17.67*(templ-273.16)/(templ-29.65))*(ep2/p2di(i,k+1)) temps=templ + ((qvx(i,k)+qcxl(i,k))-rvls)/(cp/xlv + & ep2*xlv*rvls/(rd*templ**2)) rvls=100.*6.112*EXP(17.67*(temps-273.15)/(temps-29.65))*(ep2/p2di(i,k+1)) rcldb=max((qvx(i,k)+qcxl(i,k))-rvls,0.) dthvx(i) = (thlix(i,k+2)+thx(i,k+2)*ep1*(qvx(i,k+2)+qcxl(i,k+2))) & - (thlix(i,k) + thx(i,k) *ep1*(qvx(i,k) +qcxl(i,k))) dthvx(i) = max(dthvx(i),0.1) tmp1 = xlv/cp * rcldb/(pi2d(i,k)*dthvx(i)) ent_eff = 0.2 * 8. * tmp1 +0.2 radsum=0. do kk = 1,kpbl(i)-1 radflux=rthraten(i,kk)*pi2d(i,kk) radflux=radflux*cp/g*(p2di(i,kk)-p2di(i,kk+1)) if (radflux < 0.0 ) radsum=abs(radflux)+radsum enddo radsum=max(radsum,0.0) bfx0 = max(max(sflux(i),0.0)-radsum/rhox2(i,k)/cp,0.) bfx0 = max(sflux(i),0.0) wm3 = (govrth(i)*bfx0*hpbl(i))+5. * ust3(i) wm2(i) = wm3**h2 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i) dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin) we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i))) bfx0 = max(radsum/rhox2(i,k)/cp-max(sflux(i),0.0),0.) bfx0 = max(radsum/rhox2(i,k)/cp,0.) wm3 = (g/thvx(i,k)*bfx0*hpbl(i)) wm2(i) = wm2(i)+wm3**h2 bfxpbl(i) = - ent_eff * bfx0 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),0.1) we(i) = we(i) + max(bfxpbl(i)/dthvx(i),-sqrt(wm3**h2)) bfx0 = max(radsum/rhox2(i,k)/cp,0.) wstar3_2(i) = (g/thvx(i,k)*bfx0*hpbl(i)) wscale(i) = (ust3(i)+phifac*karman*(wstar3(i)+wstar3_2(i))*0.5)**h1 wscale(i) = min(wscale(i),ust(i)*aphi16) wscale(i) = max(wscale(i),ust(i)/aphi5) gamfac = bfac/rhox(i)/wscale(i) hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt) hgamq(i) = min(gamfac*qfx(i),gamcrq) gamfac = bfac/rhox2(i,k)/wscale(i) hgamt2(i,k) = min(gamfac*radsum/cp,gamcrt) hgamt(i) = max(hgamt(i),0.0) + max(hgamt2(i,k),0.0) brint = -15.9*ust(i)*ust(i)/wspd(i)*(wstar3(i)+wstar3_2(i))/(wscale(i)**4.) hgamu(i) = brint*ux(i,1) hgamv(i) = brint*vx(i,1) endif endif prpbl(i) = 1.0 dthx = max(thx(i,k+1)-thx(i,k),tmin) dqx = min(qvx(i,k+1)-qvx(i,k),0.0) hfxpbl(i) = we(i)*dthx qfxpbl(i) = we(i)*dqx dux = ux(i,k+1)-ux(i,k) dvx = vx(i,k+1)-vx(i,k) if(dux.gt.tmin) then ufxpbl(i) = max(prpbl(i)*we(i)*dux,-ust(i)*ust(i)) elseif(dux.lt.-tmin) then ufxpbl(i) = min(prpbl(i)*we(i)*dux,ust(i)*ust(i)) else ufxpbl(i) = 0.0 endif if(dvx.gt.tmin) then vfxpbl(i) = max(prpbl(i)*we(i)*dvx,-ust(i)*ust(i)) elseif(dvx.lt.-tmin) then vfxpbl(i) = min(prpbl(i)*we(i)*dvx,ust(i)*ust(i)) else vfxpbl(i) = 0.0 endif delb = govrth(i)*d3*hpbl(i) delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.) endif enddo do k = kts,klpbl do i = its,ite if(pblflg(i).and.k.ge.kpbl(i))then entfac(i,k) = ((zq(i,k+1)-hpbl(i))/delta(i))**2. else entfac(i,k) = 1.e30 endif enddo enddo do k = kts,klpbl do i = its,ite if(k.lt.kpbl(i)) then zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.) zfacent(i,k) = (1.-zfac(i,k))**3. wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1 wscalek2(i,k) = (phifac*karman*wstar3_2(i)*(zfac(i,k)))**h1 if(sfcflg(i)) then prfac = conpr prfac2 = 15.9*(wstar3(i)+wstar3_2(i))/ust3(i)/(1.+4.*karman*(wstar3(i)+wstar3_2(i))/ust3(i)) prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2. else prfac = 0. prfac2 = 0. prnumfac = 0. phim8z = 1.+aphi5*zol1(i)*zq(i,k+1)/zl1(i) wscalek(i,k) = ust(i)/phim8z wscalek(i,k) = max(wscalek(i,k),0.001) endif prnum0 = (phih(i)/phim(i)+prfac) prnum0 = max(min(prnum0,prmax),prmin) xkzm(i,k) = wscalek(i,k) *karman* zq(i,k+1) * zfac(i,k)**pfac+ & wscalek2(i,k)*karman*(hpbl(i)-zq(i,k+1))*(1-zfac(i,k))**pfac if (k.eq.kpbl(i)-1.and.cloudflg(i).and.we(i).lt.0.0) then xkzm(i,k) = 0.0 endif prnum = 1. + (prnum0-1.)*exp(prnumfac) xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac) prnum0 = prnum0/(1.+prfac2*karman*sfcfrac) prnum = 1. + (prnum0-1.)*exp(prnumfac) xkzh(i,k) = xkzm(i,k)/prnum xkzm(i,k) = xkzm(i,k)+xkzom(i,k) xkzh(i,k) = xkzh(i,k)+xkzoh(i,k) xkzq(i,k) = xkzq(i,k)+xkzoh(i,k) xkzm(i,k) = min(xkzm(i,k),xkzmax) xkzh(i,k) = min(xkzh(i,k),xkzmax) xkzq(i,k) = min(xkzq(i,k),xkzmax) endif enddo enddo do k = kts,kte-1 do i = its,ite if(k.ge.kpbl(i)) then ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k)) & +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k))) & /(dza(i,k+1)*dza(i,k+1))+1.e-9 govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k))) ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1)) if(imvdif.eq.1)then if((qcxl(i,k)+qixl(i,k)).gt.0.01e-3.and. & (qcxl(i,k+1)+qixl(i,k+1)).gt.0.01e-3)then qmean = 0.5*(qvx(i,k)+qvx(i,k+1)) tmean = 0.5*(tx(i,k)+tx(i,k+1)) alph = xlv*qmean/rd/tmean chi = xlv*xlv*qmean/cp/rv/tmean/tmean ri = (1.+alph)*(ri-g*g/ss/tmean/cp*((chi-alph)/(1.+chi))) endif endif zk = karman*zq(i,k+1) rlamdz = min(max(0.1*dza(i,k+1),rlam),300.) rlamdz = min(dza(i,k+1),rlamdz) rl2 = (zk*rlamdz/(rlamdz+zk))**2 dk = rl2*sqrt(ss) if(ri.lt.0.)then ri = max(ri, rimin) sri = sqrt(-ri) xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri)) xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri)) else xkzh(i,k) = dk/(1+5.*ri)**2 prnum = 1.0+2.1*ri prnum = min(prnum,prmax) xkzm(i,k) = xkzh(i,k)*prnum endif xkzm(i,k) = xkzm(i,k)+xkzom(i,k) xkzh(i,k) = xkzh(i,k)+xkzoh(i,k) xkzm(i,k) = min(xkzm(i,k),xkzmax) xkzh(i,k) = min(xkzh(i,k),xkzmax) xkzml(i,k) = xkzm(i,k) xkzhl(i,k) = xkzh(i,k) endif enddo enddo do k = kts,kte do i = its,ite au(i,k) = 0. al(i,k) = 0. ad(i,k) = 0. f1(i,k) = 0. enddo enddo do i = its,ite ad(i,1) = 1. f1(i,1) = thx(i,1)-300.+(1.0-bepswitch)*hfx(i)/cont/del(i,1)*dt2 enddo do k = kts,kte-1 do i = its,ite dtodsd = sfk2d(i,k)*dt2/del(i,k) dtodsu = sfk2d(i,k)*dt2/del(i,k+1) dsig = p2d(i,k)-p2d(i,k+1) rdz = 1./dza(i,k+1) tem1 = dsig*xkzh(i,k)*rdz if(pblflg(i).and.k.lt.kpbl(i)) then dsdzt = tem1*(-hgamt(i)/hpbl(i)-hfxpbl(i)*zfacent(i,k)/xkzh(i,k)) f1(i,k) = f1(i,k)+dtodsd*dsdzt f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k)) xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k)) xkzh(i,k) = max(xkzh(i,k),xkzoh(i,k)) xkzh(i,k) = min(xkzh(i,k),xkzmax) f1(i,k+1) = thx(i,k+1)-300. else f1(i,k+1) = thx(i,k+1)-300. endif tem1 = dsig*xkzh(i,k)*rdz dsdz2 = tem1*rdz au(i,k) = -dtodsd*dsdz2/vlk2d(i,k) al(i,k) = -dtodsu*dsdz2/vlk2d(i,k) ad(i,k) = ad(i,k)-au(i,k) ad(i,k+1) = 1.-al(i,k) exch_hx(i,k+1) = xkzh(i,k) enddo enddo do k = kts,kte do i = its,ite ad(i,k) = ad(i,k) - a_t2d(i,k)*dt2 f1(i,k) = f1(i,k) + b_t2d(i,k)*dt2 enddo enddo do k = kts,kte do i = its,ite cu(i,k) = au(i,k) r1(i,k) = f1(i,k) enddo enddo call tridin_ysu(al,ad,cu,r1,au,f1,its,ite,kts,kte,1) do k = kte,kts,-1 do i = its,ite ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k) ttnp(i,k) = ttend if(present(dtsfc)) dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)/pi2d(i,k) enddo enddo do i = its,ite do k = kts,kte-1 if(k .ge. kpbl(i)) xkzq(i,k) = xkzh(i,k) enddo enddo do i = its,ite do k = kts,kte au(i,k) = 0. al(i,k) = 0. ad(i,k) = 0. f1(i,k) = 0. r1(i,k) = 0. enddo k = 1 ad(i,1) = 1. f1(i,1) = qvx(i,1)+(1.0-bepswitch)*qfx(i)*g/del(i,1)*dt2 do k = kts,kte-1 dtodsd = sfk2d(i,k)*dt2/del(i,k) dtodsu = sfk2d(i,k)*dt2/del(i,k+1) dsig = p2d(i,k)-p2d(i,k+1) rdz = 1./dza(i,k+1) tem1 = dsig*xkzq(i,k)*rdz if(pblflg(i).and.k.lt.kpbl(i)) then dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzq(i,k)) f1(i,k) = f1(i,k)+dtodsd*dsdzq f1(i,k+1) = qvx(i,k+1)-dtodsu*dsdzq elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k)) xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k)) xkzq(i,k) = max(xkzq(i,k),xkzoh(i,k)) xkzq(i,k) = min(xkzq(i,k),xkzmax) f1(i,k+1) = qvx(i,k+1) else f1(i,k+1) = qvx(i,k+1) endif tem1 = dsig*xkzq(i,k)*rdz dsdz2 = tem1*rdz au(i,k) = -dtodsd*dsdz2/vlk2d(i,k) al(i,k) = -dtodsu*dsdz2/vlk2d(i,k) ad(i,k) = ad(i,k)-au(i,k) ad(i,k+1) = 1.-al(i,k) enddo do k = kts,kte adv(i,k) = ad(i,k) - a_q2d(i,k)*dt2 f1(i,k) = f1(i,k) + b_q2d(i,k)*dt2 enddo do k = kts,kte cu(i,k) = au(i,k) r1(i,k) = f1(i,k) enddo enddo call tridin_ysu(al,adv,cu,r1,au,f1,its,ite,kts,kte,1) do i = its,ite do k = kte,kts,-1 qtend = (f1(i,k)-qvx(i,k))*rdt qvtnp(i,k) = qtend if(present(dqsfc)) dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k) enddo enddo if(f_qc) then do i = its,ite do k = kts,kte f1(i,k) = qcxl(i,k) r1(i,k) = f1(i,k) enddo enddo call tridin_ysu(al,ad,cu,r1,au,f1,its,ite,kts,kte,1) do i = its,ite do k = kte,kts,-1 qtend = (f1(i,k)-qcxl(i,k))*rdt qctnp(i,k) = qtend enddo enddo endif if(f_qi) then do i = its,ite do k = kts,kte f1(i,k) = qixl(i,k) r1(i,k) = f1(i,k) enddo enddo call tridin_ysu(al,ad,cu,r1,au,f1,its,ite,kts,kte,1) do i = its,ite do k = kte,kts,-1 qtend = (f1(i,k)-qixl(i,k))*rdt qitnp(i,k) = qtend enddo enddo endif do n = 1, nmix do i = its,ite do k = kts,kte f1(i,k) = qmix(i,k,n) r1(i,k) = f1(i,k) enddo enddo call tridin_ysu(al,ad,cu,r1,au,f1,its,ite,kts,kte,1) do i = its,ite do k = kte,kts,-1 qtend = (f1(i,k)-qmix(i,k,n))*rdt qmixtnp(i,k,n) = qtend enddo enddo enddo do i = its,ite do k = kts,kte au(i,k) = 0. al(i,k) = 0. ad(i,k) = 0. f1(i,k) = 0. f2(i,k) = 0. enddo enddo do i = its,ite do k= kts, kte-1 shear_ysu(i,k)=xkzm(i,k)*((-hgamu(i)/hpbl(i)+(ux(i,k+1)-ux(i,k))/dza(i,k+1))*(ux(i,k+1)-ux(i,k))/dza(i,k+1) & + (-hgamv(i)/hpbl(i)+(vx(i,k+1)-vx(i,k))/dza(i,k+1))*(vx(i,k+1)-vx(i,k))/dza(i,k+1)) buoy_ysu(i,k)=xkzh(i,k)*g*(1.0/thx(i,k))*(-hgamt(i)/hpbl(i)+(thx(i,k+1)-thx(i,k))/dza(i,k+1)) zk = karman*zq(i,k+1) if (k.ge.kpbl(i)) then rlamdz = min(max(0.1*dza(i,k+1),rlam),300.) rlamdz = min(dza(i,k+1),rlamdz) else rlamdz = 150.0 endif el_ysu(i,k) = zk*rlamdz/(rlamdz+zk) tke_ysu(i,k)=16.6*el_ysu(i,k)*(shear_ysu(i,k)-buoy_ysu(i,k)) if(tke_ysu(i,k).le.0) then tke_ysu(i,k)=0.0 else tke_ysu(i,k)=(tke_ysu(i,k))**0.66 endif enddo do k = kts,kte thvx_1d(k) = thvx(i,k) tke_1d(k) = tke_ysu(i,k) zq_1d(k) = zq(i,k) dzq_1d(k) = dzq(i,k) enddo zq_1d(kte+1) = zq(i,kte+1) call get_pblh(kts,kte,pblh_ysu(i),thvx_1d,tke_1d,zq_1d,dzq_1d,xland(i)) if (xland(i).lt.1.5) then fluxc = max(sflux(i),0.0) vconvc=1. VCONV = vconvc*(g/thvx(i,1)*pblh_ysu(i)*fluxc)**.33 else VCONV = 0. endif vconvfx(i) = vconv fric(i,1)=ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 & *(wspd1(i)/wspd(i))**2 if(present(ctopo)) then vconvnew=0.9*vconvfx(i)+1.5*(max((pblh_ysu(i)-500)/1000.0,0.0)) vconvlim = min(vconvnew,1.0) ad(i,1) = 1.+fric(i,1)*vconvlim+ctopo(i)*fric(i,1)*(1-vconvlim) ad(i,1) = ad(i,1) - bepswitch*frc_urb1d(i)* & (fric(i,1)*vconvlim+ctopo(i)*fric(i,1)*(1-vconvlim)) else ad(i,1) = 1.+fric(i,1) endif f1(i,1) = ux(i,1)+uoxl(i)*ust(i)**2*rhox(i)*g/del(i,1)*dt2/wspd1(i)*(wspd1(i)/wspd(i))**2 f2(i,1) = vx(i,1)+voxl(i)*ust(i)**2*rhox(i)*g/del(i,1)*dt2/wspd1(i)*(wspd1(i)/wspd(i))**2 enddo do k = kts,kte-1 do i = its,ite dtodsd = sfk2d(i,k)*dt2/del(i,k) dtodsu = sfk2d(i,k)*dt2/del(i,k+1) dsig = p2d(i,k)-p2d(i,k+1) rdz = 1./dza(i,k+1) tem1 = dsig*xkzm(i,k)*rdz if(pblflg(i).and.k.lt.kpbl(i))then dsdzu = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k)) dsdzv = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k)) f1(i,k) = f1(i,k)+dtodsd*dsdzu f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu f2(i,k) = f2(i,k)+dtodsd*dsdzv f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then xkzm(i,k) = prpbl(i)*xkzh(i,k) xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k)) xkzm(i,k) = max(xkzm(i,k),xkzom(i,k)) xkzm(i,k) = min(xkzm(i,k),xkzmax) f1(i,k+1) = ux(i,k+1) f2(i,k+1) = vx(i,k+1) else f1(i,k+1) = ux(i,k+1) f2(i,k+1) = vx(i,k+1) endif tem1 = dsig*xkzm(i,k)*rdz dsdz2 = tem1*rdz au(i,k) = -dtodsd*dsdz2/vlk2d(i,k) al(i,k) = -dtodsu*dsdz2/vlk2d(i,k) ad(i,k) = ad(i,k)-au(i,k) ad(i,k+1) = 1.-al(i,k) exch_mx(i,k+1) = xkzm(i,k) enddo enddo do k = kts,kte do i = its,ite ad1(i,k) = ad(i,k) end do end do do k = kts,kte do i = its,ite ad(i,k) = ad(i,k) - a_u2d(i,k)*dt2 ad1(i,k) = ad1(i,k) - a_v2d(i,k)*dt2 f1(i,k) = f1(i,k) + b_u2d(i,k)*dt2 f2(i,k) = f2(i,k) + b_v2d(i,k)*dt2 enddo enddo do k = kts,kte do i = its,ite cu(i,k) = au(i,k) r1(i,k) = f1(i,k) r2(i,k) = f2(i,k) enddo enddo call tridi2n(al,ad,ad1,cu,r1,r2,au,f1,f2,its,ite,kts,kte,1) do k = kte,kts,-1 do i = its,ite utend = (f1(i,k)-ux(i,k))*rdt vtend = (f2(i,k)-vx(i,k))*rdt utnp(i,k) = utend vtnp(i,k) = vtend if(present(dusfc)) dusfc(i) = dusfc(i) + utend*conwrc*del(i,k) if(present(dvsfc)) dvsfc(i) = dvsfc(i) + vtend*conwrc*del(i,k) enddo enddo do i = its,ite if(present(ctopo).and.present(ctopo2)) then u10(i) = ctopo2(i)*u10(i)+(1-ctopo2(i))*ux(i,1) v10(i) = ctopo2(i)*v10(i)+(1-ctopo2(i))*vx(i,1) endif enddo do i = its,ite kpbl1d(i) = kpbl(i) enddo errmsg = 'bl_ysu_run OK' errflg = 0 end subroutine bl_ysu_run subroutine tridi2n(cl,cm,cm1,cu,r1,r2,au,f1,f2,its,ite,kts,kte,nt) implicit none integer, intent(in ) :: its,ite, kts,kte, nt real(kind=kind_phys), dimension( its:ite, kts+1:kte+1 ) , & intent(in ) :: cl real(kind=kind_phys), dimension( its:ite, kts:kte ) , & intent(in ) :: cm, & cm1, & r1 real(kind=kind_phys), dimension( its:ite, kts:kte,nt ) , & intent(in ) :: r2 real(kind=kind_phys), dimension( its:ite, kts:kte ) , & intent(inout) :: au, & cu, & f1 real(kind=kind_phys), dimension( its:ite, kts:kte,nt ) , & intent(inout) :: f2 real(kind=kind_phys) :: fk integer :: i,k,l,n,it l = ite n = kte do i = its,l fk = 1./cm(i,1) au(i,1) = fk*cu(i,1) f1(i,1) = fk*r1(i,1) enddo do it = 1,nt do i = its,l fk = 1./cm1(i,1) f2(i,1,it) = fk*r2(i,1,it) enddo enddo do k = kts+1,n-1 do i = its,l fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) au(i,k) = fk*cu(i,k) f1(i,k) = fk*(r1(i,k)-cl(i,k)*f1(i,k-1)) enddo enddo do it = 1,nt do k = kts+1,n-1 do i = its,l fk = 1./(cm1(i,k)-cl(i,k)*au(i,k-1)) f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it)) enddo enddo enddo do i = its,l fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) f1(i,n) = fk*(r1(i,n)-cl(i,n)*f1(i,n-1)) enddo do it = 1,nt do i = its,l fk = 1./(cm1(i,n)-cl(i,n)*au(i,n-1)) f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it)) enddo enddo do k = n-1,kts,-1 do i = its,l f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1) enddo enddo do it = 1,nt do k = n-1,kts,-1 do i = its,l f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it) enddo enddo enddo end subroutine tridi2n subroutine tridin_ysu(cl,cm,cu,r2,au,f2,its,ite,kts,kte,nt) implicit none integer, intent(in ) :: its,ite, kts,kte, nt real(kind=kind_phys), dimension( its:ite, kts+1:kte+1 ) , & intent(in ) :: cl real(kind=kind_phys), dimension( its:ite, kts:kte ) , & intent(in ) :: au, & cm, & cu real(kind=kind_phys), dimension( its:ite, kts:kte,nt ) , & intent(in ) :: r2 real(kind=kind_phys), dimension( its:ite, kts:kte,nt ) , & intent(inout) :: f2 real(kind=kind_phys) :: fk real(kind=kind_phys), dimension( its:ite, kts:kte ) :: aul integer :: i,k,l,n,it l = ite n = kte do i = its,ite do k = kts,kte aul(i,k) = 0. enddo enddo do it = 1,nt do i = its,l fk = 1./cm(i,1) aul(i,1) = fk*cu(i,1) f2(i,1,it) = fk*r2(i,1,it) enddo enddo do it = 1,nt do k = kts+1,n-1 do i = its,l fk = 1./(cm(i,k)-cl(i,k)*aul(i,k-1)) aul(i,k) = fk*cu(i,k) f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it)) enddo enddo enddo do it = 1,nt do i = its,l fk = 1./(cm(i,n)-cl(i,n)*aul(i,n-1)) f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it)) enddo enddo do it = 1,nt do k = n-1,kts,-1 do i = its,l f2(i,k,it) = f2(i,k,it)-aul(i,k)*f2(i,k+1,it) enddo enddo enddo end subroutine tridin_ysu subroutine get_pblh(kts,kte,zi,thetav1d,qke1d,zw1d,dz1d,landsea) integer,intent(in) :: kts,kte real(kind=kind_phys), intent(out) :: zi real(kind=kind_phys), intent(in) :: landsea real(kind=kind_phys), dimension(kts:kte), intent(in) :: thetav1d, qke1d, dz1d real(kind=kind_phys), dimension(kts:kte+1), intent(in) :: zw1d real(kind=kind_phys) :: pblh_tke,qtke,qtkem1,wt,maxqke,tkeeps,minthv real(kind=kind_phys) :: delt_thv real(kind=kind_phys), parameter :: sbl_lim = 200. real(kind=kind_phys), parameter :: sbl_damp = 400. integer :: i,j,k,kthv,ktke k = kts+1 kthv = 1 ktke = 1 maxqke = 0. minthv = 9.e9 do while (zw1d(k) .le. 500.) 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.025) tkeeps = min(tkeeps,0.25) zi=0. k = kthv+1 if((landsea-1.5).ge.0)then delt_thv = 0.75 else delt_thv = 1.5 endif zi=0. k = kthv+1 do while (zi .eq. 0.) if (thetav1d(k) .ge. (minthv + delt_thv))then 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 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 wt=.5*tanh((zi - sbl_lim)/sbl_damp) + .5 zi=pblh_tke*(1.-wt) + zi*wt end subroutine get_pblh end module bl_ysu