module module_bl_ysu contains subroutine ysu(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, & rublten,rvblten,rthblten, & rqvblten,rqcblten,rqiblten,flag_qi, & cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, & dz8w,psfc, & znt,ust,hpbl,psim,psih, & xland,hfx,qfx,wspd,br, & dt,kpbl2d, & exch_h,exch_m, & wstar,delta, & u10,v10, & uoce,voce, & rthraten,ysu_topdown_pblmix, & ctopo,ctopo2, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & regime & ) implicit none integer,parameter :: ndiff = 3 real,parameter :: rcl = 1.0 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte integer, intent(in) :: ysu_topdown_pblmix real, intent(in ) :: dt,cp,g,rovcp,rovg,rd,xlv,rv real, intent(in ) :: ep1,ep2,karman real, dimension( ims:ime, kms:kme, jms:jme ) , & intent(in ) :: qv3d, & qc3d, & qi3d, & p3d, & pi3d, & th3d, & t3d, & dz8w, & rthraten real, dimension( ims:ime, kms:kme, jms:jme ) , & intent(in ) :: p3di real, dimension( ims:ime, kms:kme, jms:jme ) , & intent(inout) :: rublten, & rvblten, & rthblten, & rqvblten, & rqcblten real, dimension( ims:ime, kms:kme, jms:jme ) , & intent(inout) :: exch_h, & exch_m real, dimension( ims:ime, jms:jme ) , & intent(inout) :: wstar real, dimension( ims:ime, jms:jme ) , & intent(inout) :: delta real, dimension( ims:ime, jms:jme ) , & intent(inout) :: u10, & v10 real, dimension( ims:ime, jms:jme ) , & intent(in ) :: uoce, & voce real, dimension( ims:ime, jms:jme ) , & intent(in ) :: xland, & hfx, & qfx, & br, & psfc real, dimension( ims:ime, jms:jme ) , & intent(in ) :: & psim, & psih real, dimension( ims:ime, jms:jme ) , & intent(inout) :: znt, & ust, & hpbl, & wspd real, dimension( ims:ime, kms:kme, jms:jme ) , & intent(in ) :: u3d, & v3d integer, dimension( ims:ime, jms:jme ) , & intent(out ) :: kpbl2d logical, intent(in) :: flag_qi real, dimension( ims:ime, jms:jme ) , & optional , & intent(inout) :: regime real, dimension( ims:ime, kms:kme, jms:jme ) , & optional , & intent(inout) :: rqiblten real, dimension( ims:ime, jms:jme ) , & optional , & intent(in ) :: ctopo, & ctopo2 integer :: i,j,k real, dimension( its:ite, kts:kte*ndiff ) :: rqvbl2dt, & qv2d real, dimension( its:ite, kts:kte ) :: pdh real, dimension( its:ite, kts:kte+1 ) :: pdhi real, dimension( its:ite ) :: & dusfc, & dvsfc, & dtsfc, & dqsfc qv2d(its:ite,:) = 0.0 do j = jts,jte do k = kts,kte+1 do i = its,ite if(k.le.kte)pdh(i,k) = p3d(i,k,j) pdhi(i,k) = p3di(i,k,j) enddo enddo do k = kts,kte do i = its,ite qv2d(i,k) = qv3d(i,k,j) qv2d(i,k+kte) = qc3d(i,k,j) if(present(rqiblten)) qv2d(i,k+kte+kte) = qi3d(i,k,j) enddo enddo call ysu2d(J=j,ux=u3d(ims,kms,j),vx=v3d(ims,kms,j) & ,tx=t3d(ims,kms,j) & ,qx=qv2d(its,kts) & ,p2d=pdh(its,kts),p2di=pdhi(its,kts) & ,pi2d=pi3d(ims,kms,j) & ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) & ,ttnp=rthblten(ims,kms,j),qtnp=rqvbl2dt(its,kts),ndiff=ndiff & ,cp=cp,g=g,rovcp=rovcp,rd=rd,rovg=rovg & ,xlv=xlv,rv=rv & ,ep1=ep1,ep2=ep2,karman=karman & ,dz8w2d=dz8w(ims,kms,j) & ,psfcpa=psfc(ims,j),znt=znt(ims,j),ust=ust(ims,j) & ,hpbl=hpbl(ims,j) & ,regime=regime(ims,j),psim=psim(ims,j) & ,psih=psih(ims,j),xland=xland(ims,j) & ,hfx=hfx(ims,j),qfx=qfx(ims,j) & ,wspd=wspd(ims,j),br=br(ims,j) & ,dusfc=dusfc,dvsfc=dvsfc,dtsfc=dtsfc,dqsfc=dqsfc & ,dt=dt,rcl=1.0,kpbl1d=kpbl2d(ims,j) & ,exch_hx=exch_h(ims,kms,j) & ,exch_mx=exch_m(ims,kms,j) & ,wstar=wstar(ims,j) & ,delta=delta(ims,j) & ,u10=u10(ims,j),v10=v10(ims,j) & ,uox=uoce(ims,j),vox=voce(ims,j) & ,rthraten=rthraten(ims,kms,j),p2diORG=p3di(ims,kms,j) & ,ysu_topdown_pblmix=ysu_topdown_pblmix & ,ctopo=ctopo(ims,j),ctopo2=ctopo2(ims,j) & ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde & ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme & ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte ) do k = kts,kte do i = its,ite rthblten(i,k,j) = rthblten(i,k,j)/pi3d(i,k,j) rqvblten(i,k,j) = rqvbl2dt(i,k) rqcblten(i,k,j) = rqvbl2dt(i,k+kte) if(present(rqiblten)) rqiblten(i,k,j) = rqvbl2dt(i,k+kte+kte) enddo enddo enddo end subroutine ysu subroutine ysu2d(j,ux,vx,tx,qx,p2d,p2di,pi2d, & utnp,vtnp,ttnp,qtnp,ndiff, & cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, & dz8w2d,psfcpa, & znt,ust,hpbl,psim,psih, & xland,hfx,qfx,wspd,br, & dusfc,dvsfc,dtsfc,dqsfc, & dt,rcl,kpbl1d, & exch_hx,exch_mx, & wstar,delta, & u10,v10, & uox,vox, & rthraten,p2diORG, & ysu_topdown_pblmix, & ctopo,ctopo2, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & regime & ) implicit none real,parameter :: xkzminm = 0.1,xkzminh = 0.01 real,parameter :: xkzmin = 0.01,xkzmax = 1000.,rimin = -100. real,parameter :: rlam = 30.,prmin = 0.25,prmax = 4. real,parameter :: brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4 real,parameter :: afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0 real,parameter :: phifac = 8.,sfcfrac = 0.1 real,parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001 real,parameter :: h1 = 0.33333333, h2 = 0.6666667 real,parameter :: zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16. real,parameter :: tmin=1.e-2 real,parameter :: gamcrt = 3.,gamcrq = 2.e-3 real,parameter :: xka = 2.4e-5 integer,parameter :: imvdif = 1 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & j,ndiff integer, intent(in) :: ysu_topdown_pblmix real, intent(in ) :: dt,rcl,cp,g,rovcp,rovg,rd,xlv,rv real, intent(in ) :: ep1,ep2,karman real, dimension( ims:ime, kms:kme ), & intent(in) :: dz8w2d, & pi2d, & p2diorg real, dimension( ims:ime, kms:kme ) , & intent(in ) :: tx real, dimension( its:ite, kts:kte*ndiff ) , & intent(in ) :: qx real, dimension( ims:ime, kms:kme ) , & intent(inout) :: utnp, & vtnp, & ttnp real, dimension( its:ite, kts:kte*ndiff ) , & intent(inout) :: qtnp real, dimension( its:ite, kts:kte+1 ) , & intent(in ) :: p2di real, dimension( its:ite, kts:kte ) , & intent(in ) :: p2d real, dimension( ims:ime ) , & intent(inout) :: ust, & hpbl, & znt real, dimension( ims:ime ) , & intent(in ) :: xland, & hfx, & qfx real, dimension( ims:ime ), intent(inout) :: wspd real, dimension( ims:ime ), intent(in ) :: br real, dimension( ims:ime ), intent(in ) :: psim, & psih real, dimension( ims:ime ), intent(in ) :: psfcpa integer, dimension( ims:ime ), intent(out ) :: kpbl1d real, dimension( ims:ime, kms:kme ) , & intent(in ) :: ux, & vx, & rthraten real, dimension( ims:ime ) , & optional , & intent(in ) :: ctopo, & ctopo2 real, dimension( ims:ime ) , & optional , & intent(inout) :: regime real, dimension( its:ite ) :: hol real, dimension( its:ite, kts:kte+1 ) :: zq real, dimension( its:ite, kts:kte ) :: & thx,thvx,thlix, & del, & dza, & dzq, & xkzom, & xkzoh, & za real, dimension( its:ite ) :: & rhox, & govrth, & zl1,thermal, & wscale, & hgamt,hgamq, & brdn,brup, & phim,phih, & dusfc,dvsfc, & dtsfc,dqsfc, & prpbl, & wspd1,thermalli real, dimension( its:ite, kts:kte ) :: xkzm,xkzh, & f1,f2, & r1,r2, & ad,au, & cu, & al, & xkzq, & zfac, & rhox2, & hgamt2 real, dimension( ims:ime, kms:kme ) , & intent(inout) :: exch_hx, & exch_mx real, dimension( ims:ime ) , & intent(inout) :: u10, & v10 real, dimension( ims:ime ) , & intent(in ) :: uox, & vox real, dimension( its:ite ) :: & brcr, & sflux, & zol1, & brcr_sbro real, dimension( its:ite, kts:kte, ndiff) :: 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, ktrace1, ktrace2, ktrace3 real :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0 real :: ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri real :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz real :: utend,vtend,ttend,qtend real :: dtstep,govrthv real :: cont, conq, conw, conwrc real, dimension( its:ite, kts:kte ) :: wscalek,wscalek2 real, dimension( ims:ime ) :: wstar real, dimension( ims:ime ) :: delta real, dimension( its:ite, kts:kte ) :: xkzml,xkzhl, & zfacent,entfac real, dimension( its:ite ) :: ust3, & wstar3, & wstar3_2, & hgamu,hgamv, & wm2, we, & bfxpbl, & hfxpbl,qfxpbl, & ufxpbl,vfxpbl, & dthvx real :: 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, dimension( ims:ime, kms:kme ) :: fric, & tke_ysu,& el_ysu,& shear_ysu,& buoy_ysu real, dimension( ims:ime ) :: pblh_ysu,& vconvfx klpbl = kte cont=cp/g conq=xlv/g conw=1./g conwrc = conw*sqrt(rcl) conpr = bfac*karman*sfcfrac ktrace1 = 0 ktrace2 = 0 + kte ktrace3 = 0 + kte*2 do k = kts,kte do i = its,ite thx(i,k) = tx(i,k)/pi2d(i,k) thlix(i,k) = (tx(i,k)-xlv*qx(i,ktrace2+k)/cp-2.834E6*qx(i,ktrace3+k)/cp)/pi2d(i,k) enddo enddo do k = kts,kte do i = its,ite tvcon = (1.+ep1*qx(i,k)) thvx(i,k) = thx(i,k)*tvcon enddo enddo do i = its,ite tvcon = (1.+ep1*qx(i,1)) rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon) govrth(i) = g/thx(i,1) enddo 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*qx(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 utnp(its:ite,:) = 0. vtnp(its:ite,:) = 0. ttnp(its:ite,:) = 0. qtnp(its:ite,:) = 0. do i = its,ite wspd1(i) = sqrt( (ux(i,1)-uox(i))*(ux(i,1)-uox(i)) + (vx(i,1)-vox(i))*(vx(i,1)-vox(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 dusfc(i) = 0. dvsfc(i) = 0. dtsfc(i) = 0. dqsfc(i) = 0. enddo do i = its,ite hgamt(i) = 0. hgamq(i) = 0. wscale(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.eq.1)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((qx(i,ktrace2+k)+qx(i,ktrace3+k)).gt.0.01e-3.and.ysu_topdown_pblmix.eq.1)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 + ((qx(i,k)+qx(i,ktrace2+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((qx(i,k)+qx(i,ktrace2+k))-rvls,0.) dthvx(i) = (thlix(i,k+2)+thx(i,k+2)*ep1*(qx(i,k+2)+qx(i,ktrace2+k+2))) & - (thlix(i,k) + thx(i,k) *ep1*(qx(i,k) +qx(i,ktrace2+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*(p2diORG(i,kk)-p2diORG(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(qx(i,k+1)-qx(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.and.ndiff.ge.3)then if((qx(i,ktrace2+k)+qx(i,ktrace3+k)).gt.0.01e-3.and.(qx(i & ,ktrace2+k+1)+qx(i,ktrace3+k+1)).gt.0.01e-3)then qmean = 0.5*(qx(i,k)+qx(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.+hfx(i)/cont/del(i,1)*dt2 enddo do k = kts,kte-1 do i = its,ite dtodsd = dt2/del(i,k) dtodsu = 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 al(i,k) = -dtodsu*dsdz2 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 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) = ttnp(i,k)+ttend dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)/pi2d(i,k) enddo enddo do k = kts,kte do i = its,ite au(i,k) = 0. al(i,k) = 0. ad(i,k) = 0. enddo enddo do ic = 1,ndiff do i = its,ite do k = kts,kte f3(i,k,ic) = 0. enddo enddo enddo do i = its,ite ad(i,1) = 1. f3(i,1,1) = qx(i,1)+qfx(i)*g/del(i,1)*dt2 enddo if(ndiff.ge.2) then do ic = 2,ndiff is = (ic-1) * kte do i = its,ite f3(i,1,ic) = qx(i,1+is) enddo enddo endif do k = kts,kte-1 do i = its,ite if(k.ge.kpbl(i)) then xkzq(i,k) = xkzh(i,k) endif enddo enddo do k = kts,kte-1 do i = its,ite dtodsd = dt2/del(i,k) dtodsu = 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)) f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq f3(i,k+1,1) = qx(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) f3(i,k+1,1) = qx(i,k+1) else f3(i,k+1,1) = qx(i,k+1) endif tem1 = dsig*xkzq(i,k)*rdz dsdz2 = tem1*rdz au(i,k) = -dtodsd*dsdz2 al(i,k) = -dtodsu*dsdz2 ad(i,k) = ad(i,k)-au(i,k) ad(i,k+1) = 1.-al(i,k) enddo enddo if(ndiff.ge.2) then do ic = 2,ndiff is = (ic-1) * kte do k = kts,kte-1 do i = its,ite f3(i,k+1,ic) = qx(i,k+1+is) enddo enddo enddo endif do k = kts,kte do i = its,ite cu(i,k) = au(i,k) enddo enddo do ic = 1,ndiff do k = kts,kte do i = its,ite r3(i,k,ic) = f3(i,k,ic) enddo enddo enddo call tridin_ysu(al,ad,cu,r3,au,f3,its,ite,kts,kte,ndiff) do k = kte,kts,-1 do i = its,ite qtend = (f3(i,k,1)-qx(i,k))*rdt qtnp(i,k) = qtnp(i,k)+qtend dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k) enddo enddo if(ndiff.ge.2) then do ic = 2,ndiff is = (ic-1) * kte do k = kte,kts,-1 do i = its,ite qtend = (f3(i,k,ic)-qx(i,k+is))*rdt qtnp(i,k+is) = qtnp(i,k+is)+qtend enddo enddo enddo endif 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 CALL GET_PBLH(KTS,KTE,pblh_ysu(i),thvx(i,kts:kte),& & tke_ysu(i,kts:kte),zq(i,kts:kte+1),dzq(i,kts:kte),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) else ad(i,1) = 1.+fric(i,1) endif f1(i,1) = ux(i,1)+uox(i)*ust(i)**2*rhox(i)*g/del(i,1)*dt2/wspd1(i)*(wspd1(i)/wspd(i))**2 f2(i,1) = vx(i,1)+vox(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 = dt2/del(i,k) dtodsu = 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 al(i,k) = -dtodsu*dsdz2 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 cu(i,k) = au(i,k) r1(i,k) = f1(i,k) r2(i,k) = f2(i,k) enddo enddo call tridi1n(al,ad,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) = utnp(i,k)+utend vtnp(i,k) = vtnp(i,k)+vtend dusfc(i) = dusfc(i) + utend*conwrc*del(i,k) 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 end subroutine ysu2d subroutine tridi1n(cl,cm,cu,r1,r2,au,f1,f2,its,ite,kts,kte,nt) implicit none integer, intent(in ) :: its,ite, kts,kte, nt real, dimension( its:ite, kts+1:kte+1 ) , & intent(in ) :: cl real, dimension( its:ite, kts:kte ) , & intent(in ) :: cm, & r1 real, dimension( its:ite, kts:kte,nt ) , & intent(in ) :: r2 real, dimension( its:ite, kts:kte ) , & intent(inout) :: au, & cu, & f1 real, dimension( its:ite, kts:kte,nt ) , & intent(inout) :: f2 real :: 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./cm(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./(cm(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./(cm(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 tridi1n 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, dimension( its:ite, kts+1:kte+1 ) , & intent(in ) :: cl real, dimension( its:ite, kts:kte ) , & intent(in ) :: cm real, dimension( its:ite, kts:kte,nt ) , & intent(in ) :: r2 real, dimension( its:ite, kts:kte ) , & intent(inout) :: au, & cu real, dimension( its:ite, kts:kte,nt ) , & intent(inout) :: f2 real :: fk integer :: i,k,l,n,it l = ite n = kte do it = 1,nt do i = its,l fk = 1./cm(i,1) au(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)*au(i,k-1)) au(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)*au(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)-au(i,k)*f2(i,k+1,it) enddo enddo enddo end subroutine tridin_ysu subroutine ysuinit(rublten,rvblten,rthblten,rqvblten, & rqcblten,rqiblten,p_qi,p_first_scalar, & restart, allowed_to_read, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) implicit none logical , intent(in) :: restart, allowed_to_read integer , intent(in) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte integer , intent(in) :: p_qi,p_first_scalar real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: & rublten, & rvblten, & rthblten, & rqvblten, & rqcblten, & rqiblten 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. rqcblten(i,k,j) = 0. enddo enddo enddo endif if (p_qi .ge. p_first_scalar .and. .not.restart) then do j = jts,jtf do k = kts,ktf do i = its,itf rqiblten(i,k,j) = 0. enddo enddo enddo endif end subroutine ysuinit SUBROUTINE GET_PBLH(KTS,KTE,zi,thetav1D,qke1D,zw1D,dz1D,landsea) 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 = 200. REAL, 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 module_bl_ysu