module module_bl_shinhong contains subroutine shinhong(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, & znu,znw,p_top, & znt,ust,hpbl,psim,psih, & xland,hfx,qfx,wspd,br, & dt,kpbl2d, & exch_h, & u10,v10, & shinhong_tke_diag,tke_pbl,el_pbl,corf, & dx,dy, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & ctopo,ctopo2, & wstar,delta, & 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 ) :: shinhong_tke_diag real, intent(in ) :: dt,cp,g,rovcp,rovg,rd,xlv,rv real, intent(in ) :: ep1,ep2,karman real, intent(in ) :: dx,dy integer, dimension( ims:ime, jms:jme ) , & intent(out ) :: kpbl2d real, dimension( ims:ime, kms:kme, jms:jme ) , & intent(in ) :: u3d, & v3d real, dimension( ims:ime, kms:kme, jms:jme ) , & intent(in ) :: qv3d, & qc3d, & qi3d, & p3d, & pi3d, & th3d, & t3d, & dz8w 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 real, dimension( ims:ime, kms:kme, jms:jme ) , & intent(inout) :: tke_pbl, & el_pbl real, dimension( ims:ime, jms:jme ) , & intent(in ) :: xland, & hfx, & qfx, & corf, & br, & psfc real, dimension( ims:ime, jms:jme ) , & intent(in ) :: & psim, & psih real, dimension( ims:ime, jms:jme ) , & intent(inout) :: u10, & v10 real, dimension( ims:ime, jms:jme ) , & intent(inout) :: znt, & ust, & hpbl, & wspd logical, intent(in) :: flag_qi real, dimension( ims:ime, kms:kme, jms:jme ) , & intent(inout), optional :: rqiblten real, dimension( ims:ime, jms:jme ) , & intent(inout), optional :: wstar, & delta real, dimension( ims:ime, jms:jme ) , & intent(inout), optional :: regime real, dimension( ims:ime, jms:jme ) , & intent(in ), optional :: ctopo, & ctopo2 real, dimension( kms:kme ) , & intent(in ), optional :: znu, & znw real, optional, intent(in ) :: p_top 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 shinhong2d(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) & ,wstar=wstar(ims,j) & ,delta=delta(ims,j) & ,u10=u10(ims,j),v10=v10(ims,j) & ,ctopo=ctopo(ims,j),ctopo2=ctopo2(ims,j) & ,shinhong_tke_diag=shinhong_tke_diag & ,tke=tke_pbl(ims,kms,j),el_pbl=el_pbl(ims,kms,j) & ,corf=corf(ims,j) & ,dx=dx,dy=dy & ,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 shinhong subroutine shinhong2d(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, & wstar,delta, & shinhong_tke_diag, & tke,el_pbl,corf, & u10,v10, & ctopo,ctopo2, & dx,dy, & 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 :: 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 :: ckz = 0.001,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 real,parameter :: epsq2l = 0.01,c_1 = 1.0,gamcre = 0.224 real,parameter :: mltop = 1.0,sfcfracn1 = 0.075 real,parameter :: nlfrac = 0.7,enlfrac = -0.4 real,parameter :: a11 = 1.0,a12 = -1.15 real,parameter :: ezfac = 1.5 real,parameter :: cpent = -0.4,rigsmax = 100. real,parameter :: entfmin = 1.0, entfmax = 5.0 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 ) :: shinhong_tke_diag real, intent(in ) :: dt,rcl,cp,g,rovcp,rovg,rd,xlv,rv real, intent(in ) :: ep1,ep2,karman real, intent(in ) :: dx,dy integer, dimension( ims:ime ) , & intent(out ) :: kpbl1d real, dimension( ims:ime, kms:kme ) , & intent(in ) :: dz8w2d, & pi2d real, dimension( ims:ime, kms:kme ) , & intent(in ) :: ux, & vx real, dimension( ims:ime, kms:kme ) , & intent(in ) :: tx real, dimension( its:ite, kts:kte*ndiff ) , & intent(in ) :: qx real, dimension( its:ite, kts:kte+1 ) , & intent(in ) :: p2di real, dimension( its:ite, kts:kte ) , & intent(in ) :: p2d real, dimension( ims:ime, kms:kme ) , & intent(inout) :: utnp, & vtnp, & ttnp real, dimension( ims:ime, kms:kme ) , & intent(inout) :: exch_hx real, dimension( ims:ime, kms:kme ) , & intent(inout) :: tke, & el_pbl real, dimension( its:ite, kts:kte*ndiff ) , & intent(inout) :: qtnp real, dimension( ims:ime ) , & intent(in ) :: xland, & hfx, & qfx real, dimension( ims:ime ) , & intent(in ) :: br, & psim, & psih, & psfcpa real, dimension( ims:ime ) , & intent(in ) :: corf real, dimension( ims:ime ) , & intent(inout) :: ust, & hpbl, & znt real, dimension( ims:ime ) , & intent(inout) :: wspd real, dimension( ims:ime ) , & intent(inout) :: u10, & v10 real, dimension( ims:ime ) , & optional , & intent(in ) :: ctopo, & ctopo2 real, dimension( ims:ime ) , & optional , & intent(inout) :: regime real, dimension( its:ite ) , & intent(out ) :: wstar, & delta integer :: n,i,k,l,ic,is,nwmass integer :: klpbl, kqc, kqi integer :: lmh,lmxl real :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0 real :: ss,ri,qmean,tmean,alpha,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 :: delxy,pu1,pth1,pq1 real :: dex,hgame_c real :: zfacdx real :: amf1,amf2,bmf2,amf3,bmf3,amf4,bmf4,sflux0,snlflux0 real :: mlfrac,ezfrac,sfcfracn real :: uwst,uwstx,csfac real :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, & dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr, & prfac,prfac2,phim8z integer, dimension( its:ite ) :: kpbl real, dimension( its:ite ) :: hol real, dimension( its:ite ) :: deltaoh real, dimension( its:ite ) :: rigs, & enlfrac2, & cslen real, dimension( its:ite ) :: & rhox, & govrth, & zl1,thermal, & wscale, & hgamt,hgamq, & brdn,brup, & phim,phih, & dusfc,dvsfc, & dtsfc,dqsfc, & prpbl, & wspd1 real, dimension( its:ite ) :: & ust3, & wstar3, & hgamu,hgamv, & wm2, we, & bfxpbl, & hfxpbl,qfxpbl, & ufxpbl,vfxpbl, & dthvx real, dimension( its:ite ) :: & brcr, & sflux, & zol1, & brcr_sbro real, dimension( its:ite ) :: & efxpbl, & hpbl_cbl, & epshol, & ct real, dimension( its:ite, kts:kte ) :: xkzm,xkzh, & f1,f2, & r1,r2, & ad,au, & cu, & al, & xkzq, & zfac real, dimension( its:ite, kts:kte ) :: & thx,thvx, & del, & dza, & dzq, & xkzom, & xkzoh, & za real, dimension( its:ite, kts:kte ) :: & wscalek real, dimension( its:ite, kts:kte ) :: & xkzml,xkzhl, & zfacent,entfac real, dimension( its:ite, kts:kte ) :: & mf, & zfacmf, & entfacmf real, dimension( its:ite, kts:kte ) :: & q2x, & hgame2d, & tflux_e, & qflux_e, & tvflux_e real, dimension( its:ite, kts:kte+1 ) :: zq real, dimension( its:ite, kts:kte, ndiff ) :: r3,f3 real, dimension( kts:kte ) :: & uxk,vxk, & txk,thxk,thvxk, & q2xk, & hgame real, dimension( kts:kte ) :: & ps1d,pb1d,eps1d,pt1d, & xkze1d,eflx_l1d,eflx_nl1d, & ptke1 real, dimension( kts+1:kte ) :: & s2,gh,rig,el, & akmk,akhk, & mfk,ufxpblk,vfxpblk,qfxpblk real, dimension( kts:kte+1 ) :: zqk real, dimension( kts:kte*ndiff ) :: qxk logical, dimension( its:ite ) :: pblflg, & sfcflg, & stable logical, dimension( ndiff ) :: ifvmix klpbl = kte lmh = 1 lmxl = 1 cont=cp/g conq=xlv/g conw=1./g conwrc = conw*sqrt(rcl) conpr = bfac*karman*sfcfrac kqc = 1 + kte kqi = 1 + kte*2 nwmass = 3 ifvmix(:) = .true. do k = kts,kte do i = its,ite thx(i,k) = tx(i,k)/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) 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)*ux(i,1)+vx(i,1)*vx(i,1))+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 enddo do i = its,ite efxpbl(i) = 0.0 hpbl_cbl(i) = 0.0 epshol(i) = 0.0 ct(i) = 0.0 enddo do i = its,ite deltaoh(i) = 0.0 rigs(i) = 0.0 enlfrac2(i) = 0.0 cslen(i) = 0.0 enddo do k = kts,klpbl do i = its,ite wscalek(i,k) = 0.0 enddo enddo do k = kts,klpbl do i = its,ite zfac(i,k) = 0.0 enddo enddo do k = kts,kte do i = its,ite q2x(i,k) = 2.*tke(i,k) enddo enddo do k = kts,kte do i = its,ite el_pbl(i,k) = 0.0 hgame2d(i,k) = 0.0 tflux_e(i,k) = 0.0 qflux_e(i,k) = 0.0 tvflux_e(i,k) = 0.0 enddo enddo do k = kts,kte do i = its,ite mf(i,k) = 0.0 zfacmf(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) hpbl_cbl(i) = zq(i,1) zl1(i) = za(i,1) thermal(i)= thvx(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 epshol(i) = hol1 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) 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 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. if (wstar(i) .ne. 0) then uwst = abs(ust(i)/wstar(i)-0.5) uwstx = -80.*uwst+14. csfac = 0.5*(tanh(uwstx)+3.) else csfac = 1 endif cslen(i) = csfac*hpbl(i) endif enddo do i = its,ite hpbl_cbl(i) = hpbl(i) 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 delxy=sqrt(dx*dy) do i = its,ite pu1=pu(delxy,cslen(i)) pq1=pq(delxy,cslen(i)) if(pblflg(i)) then hgamu(i) = hgamu(i)*pu1 hgamv(i) = hgamv(i)*pu1 hgamq(i) = hgamq(i)*pq1 endif enddo delxy=sqrt(dx*dy) do i = its,ite if(pblflg(i)) then k = kpbl(i) - 1 prpbl(i) = 1.0 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) dthx = max(thx(i,k+1)-thx(i,k),tmin) dqx = min(qx(i,k+1)-qx(i,k),0.0) we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i))) hfxpbl(i) = we(i)*dthx pq1=pq(delxy,cslen(i)) qfxpbl(i) = we(i)*dqx*pq1 pu1=pu(delxy,cslen(i)) 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*pu1,-ust(i)*ust(i)) elseif(dux.lt.-tmin) then ufxpbl(i) = min(prpbl(i)*we(i)*dux*pu1,ust(i)*ust(i)) else ufxpbl(i) = 0.0 endif if(dvx.gt.tmin) then vfxpbl(i) = max(prpbl(i)*we(i)*dvx*pu1,-ust(i)*ust(i)) elseif(dvx.lt.-tmin) then vfxpbl(i) = min(prpbl(i)*we(i)*dvx*pu1,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.) delb = govrth(i)*dthvx(i) deltaoh(i) = d1*hpbl(i) + d2*wm2(i)/delb deltaoh(i) = max(ezfac*deltaoh(i),hpbl(i)-za(i,kpbl(i)-1)-1.) deltaoh(i) = min(deltaoh(i) ,hpbl(i)) if ((dux .ne. 0) .or. (dvx .ne. 0)) then rigs(i) = govrth(i)*dthvx(i)*deltaoh(i)/(dux**2.+dvx**2.) else rigs(i) = rigsmax endif rigs(i) = max(min(rigs(i), rigsmax),rimin) enlfrac2(i) = max(min(wm3/wstar3(i)/(1.+cpent/rigs(i)),entfmax), entfmin) enlfrac2(i) = enlfrac2(i)*enlfrac endif enddo do k = kts,klpbl do i = its,ite if(pblflg(i))then entfacmf(i,k) = sqrt(((zq(i,k+1)-hpbl(i))/deltaoh(i))**2.) endif if(pblflg(i).and.k.ge.kpbl(i))then entfac(i,k) = ((zq(i,k+1)-hpbl(i))/deltaoh(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 if(sfcflg(i)) then prfac = conpr prfac2 = 15.9*wstar3(i)/ust3(i)/(1.+4.*karman*wstar3(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 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.nwmass.ge.3)then if((qx(i,kqc+k-1)+qx(i,kqi+k-1)).gt.0.01e-3 & .and.(qx(i,kqc+k)+qx(i,kqi+k)).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)) alpha = xlv*qmean/rd/tmean chi = xlv*xlv*qmean/cp/rv/tmean/tmean ri = (1.+alpha)*(ri-g*g/ss/tmean/cp*((chi-alpha)/(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 i = its,ite deltaoh(i) = deltaoh(i)/hpbl(i) enddo delxy=sqrt(dx*dy) do i = its,ite mlfrac = mltop-deltaoh(i) ezfrac = mltop+deltaoh(i) zfacmf(i,1) = min(max((zq(i,2)/hpbl(i)),zfmin),1.) sfcfracn = max(sfcfracn1,zfacmf(i,1)) sflux0 = (a11+a12*sfcfracn)*sflux(i) snlflux0 = nlfrac*sflux0 amf1 = snlflux0/sfcfracn if (pblflg(i)) then amf2 = -snlflux0/(mlfrac-sfcfracn) bmf2 = -mlfrac*amf2 endif if ((deltaoh(i) .eq. 0) .and. (enlfrac2(i) .eq. 0)) then amf3 = 0. else amf3 = snlflux0*enlfrac2(i)/deltaoh(i) endif bmf3 = -amf3*mlfrac hfxpbl(i) = amf3+bmf3 pth1=pthnl(delxy,cslen(i)) hfxpbl(i) = hfxpbl(i)*pth1 do k = kts,klpbl zfacmf(i,k) = max((zq(i,k+1)/hpbl(i)),zfmin) if(pblflg(i).and.k.lt.kpbl(i)) then if(zfacmf(i,k).le.sfcfracn) then mf(i,k) = amf1*zfacmf(i,k) else if (zfacmf(i,k).le.mlfrac) then mf(i,k) = amf2*zfacmf(i,k)+bmf2 endif mf(i,k) = mf(i,k)+hfxpbl(i)*exp(-entfacmf(i,k)) mf(i,k) = mf(i,k)*pth1 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 delxy=sqrt(dx*dy) 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*(-mf(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 zfacdx=0.2*hpbl(i)/zq(i,k+1) delxy=sqrt(dx*dy)*max(zfacdx,1.0) pth1=pthl(delxy,hpbl(i)) if(pblflg(i).and.k.lt.kpbl(i)) then au(i,k) = au(i,k)*pth1 al(i,k) = al(i,k)*pth1 endif 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) if(k.eq.kte) then tflux_e(i,k) = ttend*dz8w2d(i,k) else tflux_e(i,k) = tflux_e(i,k+1) + ttend*dz8w2d(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. 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 zfacdx=0.2*hpbl(i)/zq(i,k+1) delxy=sqrt(dx*dy)*max(zfacdx,1.0) pq1=pq(delxy,hpbl(i)) if(pblflg(i).and.k.lt.kpbl(i)) then au(i,k) = au(i,k)*pq1 al(i,k) = al(i,k)*pq1 endif 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) if(k.eq.kte) then qflux_e(i,k) = qtend*dz8w2d(i,k) else qflux_e(i,k) = qflux_e(i,k+1) + qtend*dz8w2d(i,k) endif tvflux_e(i,k) = tflux_e(i,k) + qflux_e(i,k)*ep1*thx(i,k) enddo enddo do k = kts,kte do i = its,ite if(pblflg(i).and.k.lt.kpbl(i)) then hgame_c=c_1*0.2*2.5*(g/thvx(i,k))*wstar(i)/(0.25*(q2x(i,k+1)+q2x(i,k))) hgame_c=min(hgame_c,gamcre) if(k.eq.kte)then hgame2d(i,k)=hgame_c*0.5*tvflux_e(i,k)*hpbl(i) hgame2d(i,k)=max(hgame2d(i,k),0.0) else hgame2d(i,k)=hgame_c*0.5*(tvflux_e(i,k)+tvflux_e(i,k+1))*hpbl(i) hgame2d(i,k)=max(hgame2d(i,k),0.0) endif endif enddo enddo if(ndiff.ge.2) then do ic = 2,ndiff if(ifvmix(ic)) then 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 endif 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 ad(i,1) = 1.+ctopo(i)*ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 & *(wspd1(i)/wspd(i))**2 f1(i,1) = ux(i,1) f2(i,1) = vx(i,1) enddo delxy=sqrt(dx*dy) 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 zfacdx=0.2*hpbl(i)/zq(i,k+1) delxy=sqrt(dx*dy)*max(zfacdx,1.0) pu1=pu(delxy,hpbl(i)) if(pblflg(i).and.k.lt.kpbl(i)) then au(i,k) = au(i,k)*pu1 al(i,k) = al(i,k)*pu1 endif ad(i,k) = ad(i,k)-au(i,k) ad(i,k+1) = 1.-al(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 kpbl1d(i) = kpbl(i) enddo do i = its,ite u10(i) = ctopo2(i)*u10(i)+(1-ctopo2(i))*ux(i,1) v10(i) = ctopo2(i)*v10(i)+(1-ctopo2(i))*vx(i,1) enddo if (shinhong_tke_diag.eq.1) then tke_calculation: do i = its,ite do k = kts+1,kte s2(k) = 0.0 gh(k) = 0.0 rig(k) = 0.0 el(k) = 0.0 akmk(k) = 0.0 akhk(k) = 0.0 mfk(k) = 0.0 ufxpblk(k) = 0.0 vfxpblk(k) = 0.0 qfxpblk(k) = 0.0 enddo do k = kts,kte uxk(k) = 0.0 vxk(k) = 0.0 txk(k) = 0.0 thxk(k) = 0.0 thvxk(k) = 0.0 q2xk(k) = 0.0 hgame(k) = 0.0 ps1d(k) = 0.0 pb1d(k) = 0.0 eps1d(k) = 0.0 pt1d(k) = 0.0 xkze1d(k) = 0.0 eflx_l1d(k) = 0.0 eflx_nl1d(k) = 0.0 ptke1(k) = 1.0 enddo do k = kts,kte+1 zqk(k) = 0.0 enddo do k = kts,kte*ndiff qxk(k) = 0.0 enddo do k = kts,kte uxk(k) = ux(i,k) vxk(k) = vx(i,k) txk(k) = tx(i,k) thxk(k) = thx(i,k) thvxk(k) = thvx(i,k) q2xk(k) = q2x(i,k) hgame(k) = hgame2d(i,k) enddo do k = kts,kte-1 if(pblflg(i).and.k.le.kpbl(i)) then zfacdx = 0.2*hpbl(i)/za(i,k) delxy = sqrt(dx*dy)*max(zfacdx,1.0) ptke1(k+1) = ptke(delxy,hpbl(i)) endif enddo do k = kts,kte+1 zqk(k) = zq(i,k) enddo do k = kts,kte*ndiff qxk(k) = qx(i,k) enddo do k = kts+1,kte akmk(k) = xkzm(i,k-1) akhk(k) = xkzh(i,k-1) mfk(k) = mf(i,k-1)/xkzh(i,k-1) ufxpblk(k) = ufxpbl(i)*zfacent(i,k-1)/xkzm(i,k-1) vfxpblk(k) = vfxpbl(i)*zfacent(i,k-1)/xkzm(i,k-1) qfxpblk(k) = qfxpbl(i)*zfacent(i,k-1)/xkzq(i,k-1) enddo if(pblflg(i)) then k = kpbl(i) - 1 dex = 0.25*(q2xk(k+2)-q2xk(k)) efxpbl(i) = we(i)*dex endif delxy=sqrt(dx*dy) call mixlen(lmh,uxk,vxk,txk,thxk,qxk(kts),qxk(kte+1) & ,q2xk,zqk,ust(i),corf(i),epshol(i) & ,s2,gh,rig,el & ,hpbl(i),kpbl(i),lmxl,ct(i) & ,hgamu(i),hgamv(i),hgamq(i),pblflg(i) & ,mfk,ufxpblk,vfxpblk,qfxpblk & ,ep1,karman,cp & ,ids,ide,jds,jde,kds,kde & ,ims,ime,jms,jme,kms,kme & ,its,ite,jts,jte,kts,kte ) call prodq2(lmh,dt,ust(i),s2,rig,q2xk,el,zqk,akmk,akhk & ,uxk,vxk,thxk,thvxk & ,hgamu(i),hgamv(i),hgamq(i),delxy & ,hpbl(i),pblflg(i),kpbl(i) & ,mfk,ufxpblk,vfxpblk,qfxpblk & ,ep1 & ,ids,ide,jds,jde,kds,kde & ,ims,ime,jms,jme,kms,kme & ,its,ite,jts,jte,kts,kte ) call vdifq(lmh,dt,q2xk,el,zqk & ,akhk,ptke1 & ,hgame,hpbl(i),pblflg(i),kpbl(i) & ,efxpbl(i) & ,ids,ide,jds,jde,kds,kde & ,ims,ime,jms,jme,kms,kme & ,its,ite,jts,jte,kts,kte ) do k = kts,kte q2x(i,k) = amax1(q2xk(k),epsq2l) tke(i,k) = 0.5*q2x(i,k) if(k/=kts) el_pbl(i,k) = el(k) enddo enddo tke_calculation endif end subroutine shinhong2d 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 mixlen(lmh,u,v,t,the,q,cwm,q2,z,ustar,corf,epshol, & s2,gh,ri,el,hpbl,lpbl,lmxl,ct, & hgamu,hgamv,hgamq,pblflg, & mf,ufxpbl,vfxpbl,qfxpbl, & p608,vkarman,cp, & ids,ide,jds,jde,kds,kde, & ims,ime,jms,jme,kms,kme, & its,ite,jts,jte,kts,kte) implicit none real,parameter :: blckdr=0.0063,cn=0.75 real,parameter :: eps1=1.e-12,epsl=0.32,epsru=1.e-7,epsrs=1.e-7 real,parameter :: el0max=1000.,el0min=1.,elfc=0.23*0.5 real,parameter :: alph=0.30,beta=1./273.,g=9.81,btg=beta*g real,parameter :: a1=0.659888514560862645,a2x=0.6574209922667784586 real,parameter :: b1=11.87799326209552761,b2=7.226971804046074028 real,parameter :: c1=0.000830955950095854396 real,parameter :: adnh= 9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg real,parameter :: adnm=18.*a1*a1*a2x*(b2-3.*a2x)*btg real,parameter :: bdnh= 3.*a2x*(7.*a1+b2)*btg,bdnm= 6.*a1*a1 real,parameter :: aeqh=9.*a1*a2x*a2x*b1*btg*btg & +9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg real,parameter :: aeqm=3.*a1*a2x*b1*(3.*a2x+3.*b2*c1+18.*a1*c1-b2) & *btg+18.*a1*a1*a2x*(b2-3.*a2x)*btg real,parameter :: requ=-aeqh/aeqm real,parameter :: epsgh=1.e-9,epsgm=requ*epsgh real,parameter :: ubryl=(18.*requ*a1*a1*a2x*b2*c1*btg & +9.*a1*a2x*a2x*b2*btg*btg) & /(requ*adnm+adnh) real,parameter :: ubry=(1.+epsrs)*ubryl,ubry3=3.*ubry real,parameter :: aubh=27.*a1*a2x*a2x*b2*btg*btg-adnh*ubry3 real,parameter :: aubm=54.*a1*a1*a2x*b2*c1*btg -adnm*ubry3 real,parameter :: bubh=(9.*a1*a2x+3.*a2x*b2)*btg-bdnh*ubry3 real,parameter :: bubm=18.*a1*a1*c1 -bdnm*ubry3 real,parameter :: cubr=1.-ubry3,rcubr=1./cubr real,parameter :: elcbl=0.77 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte integer, intent(in ) :: lmh,lmxl,lpbl real, intent(in ) :: p608,vkarman,cp real, intent(in ) :: hpbl,corf,ustar,hgamu,hgamv,hgamq real, intent(inout) :: ct,epshol real, dimension( kts:kte ) , & intent(in ) :: cwm, & q, & q2, & t, & the, & u, & v real, dimension( kts+1:kte ) , & intent(in ) :: mf, & ufxpbl, & vfxpbl, & qfxpbl real, dimension( kts:kte+1 ) , & intent(in ) :: z real, dimension( kts+1:kte ) , & intent(out ) :: el, & ri, & gh, & s2 logical,intent(in) :: pblflg integer :: k,lpblm real :: suk,svk,elocp real :: a,aden,b,bden,aubr,bubr,blmx,el0,eloq2x,ghl,s2l, & qol2st,qol2un,qdzl,rdz,sq,srel,szq,tem,thm,vkrmz,rlambda, & rlb,rln,f real :: ckp real, dimension( kts:kte ) :: q1, & en2 real, dimension( kts+1:kte ) :: dth, & elm, & rel elocp=2.72e6/cp ct=0. do k = kts,kte q1(k) = 0. enddo do k = kts+1,kte dth(k) = the(k)-the(k-1) enddo do k = kts+2,kte if(dth(k)>0..and.dth(k-1)<=0.)then dth(k)=dth(k)+ct exit endif enddo do k = kte,kts+1,-1 rdz=2./(z(k+1)-z(k-1)) s2l=((u(k)-u(k-1))**2+(v(k)-v(k-1))**2)*rdz*rdz if(pblflg.and.k.le.lpbl)then suk=(u(k)-u(k-1))*rdz svk=(v(k)-v(k-1))*rdz s2l=(suk-hgamu/hpbl-ufxpbl(k))*suk+(svk-hgamv/hpbl-vfxpbl(k))*svk endif s2l=max(s2l,epsgm) s2(k)=s2l tem=(t(k)+t(k-1))*0.5 thm=(the(k)+the(k-1))*0.5 a=thm*p608 b=(elocp/tem-1.-p608)*thm ghl=(dth(k)*((q(k)+q(k-1)+cwm(k)+cwm(k-1))*(0.5*p608)+1.) & +(q(k)-q(k-1)+cwm(k)-cwm(k-1))*a & +(cwm(k)-cwm(k-1))*b)*rdz if(pblflg.and.k.le.lpbl)then ghl=ghl-mf(k)-(hgamq/hpbl+qfxpbl(k))*a endif if(abs(ghl)<=epsgh)ghl=epsgh en2(k)=ghl*g/thm gh(k)=ghl ri(k)=en2(k)/s2l enddo do k = kte,kts+1,-1 s2l=s2(k) ghl=gh(k) if(ghl>=epsgh)then if(s2l/ghl<=requ)then elm(k)=epsl else aubr=(aubm*s2l+aubh*ghl)*ghl bubr= bubm*s2l+bubh*ghl qol2st=(-0.5*bubr+sqrt(bubr*bubr*0.25-aubr*cubr))*rcubr eloq2x=1./qol2st elm(k)=max(sqrt(eloq2x*q2(k)),epsl) endif else aden=(adnm*s2l+adnh*ghl)*ghl bden= bdnm*s2l+bdnh*ghl qol2un=-0.5*bden+sqrt(bden*bden*0.25-aden) eloq2x=1./(qol2un+epsru) elm(k)=max(sqrt(eloq2x*q2(k)),epsl) endif enddo do k = lpbl,lmh,-1 q1(k)=sqrt(q2(k)) enddo szq=0. sq =0. do k = kte,kts+1,-1 qdzl=(q1(k)+q1(k-1))*(z(k)-z(k-1)) szq=(z(k)+z(k-1)-z(lmh)-z(lmh))*qdzl+szq sq=qdzl+sq enddo el0=min(alph*szq*0.5/sq,el0max) el0=max(el0 ,el0min) lpblm=min(lpbl+1,kte) do k = kte,lpblm,-1 el(k)=(z(k+1)-z(k-1))*elfc rel(k)=el(k)/elm(k) enddo epshol=min(epshol,0.0) ckp=elcbl*((1.0-8.0*epshol)**(1./3.)) if(lpbl>lmh)then do k = lpbl,lmh+1,-1 vkrmz=(z(k)-z(lmh))*vkarman if(pblflg) then vkrmz=ckp*(z(k)-z(lmh))*vkarman el(k)=vkrmz/(vkrmz/el0+1.) else el(k)=vkrmz/(vkrmz/el0+1.) endif rel(k)=el(k)/elm(k) enddo endif do k = lpbl-1,lmh+2,-1 srel=min(((rel(k-1)+rel(k+1))*0.5+rel(k))*0.5,rel(k)) el(k)=max(srel*elm(k),epsl) enddo f=max(corf,eps1) rlambda=f/(blckdr*ustar) do k = kte,kts+1,-1 if(en2(k)>=0.0)then vkrmz=(z(k)-z(lmh))*vkarman rlb=rlambda+1./vkrmz rln=sqrt(2.*en2(k)/q2(k))/cn el(k)=1./(rlb+rln) endif enddo end subroutine mixlen subroutine prodq2(lmh,dtturbl,ustar,s2,ri,q2,el,z,akm,akh, & uxk,vxk,thxk,thvxk, & hgamu,hgamv,hgamq,delxy, & hpbl,pblflg,kpbl, & mf,ufxpbl,vfxpbl,qfxpbl, & p608, & ids,ide,jds,jde,kds,kde, & ims,ime,jms,jme,kms,kme, & its,ite,jts,jte,kts,kte) implicit none real,parameter :: epsq2l = 0.01,c0 = 0.55,ceps = 16.6,g = 9.81 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte integer, intent(in ) :: lmh,kpbl real, intent(in ) :: p608,dtturbl,ustar real, intent(in ) :: hgamu,hgamv,hgamq,delxy,hpbl logical, intent(in ) :: pblflg real, dimension( kts:kte ) , & intent(in ) :: uxk, & vxk, & thxk, & thvxk real, dimension( kts+1:kte ) , & intent(in ) :: s2, & ri, & akm, & akh, & el, & mf, & ufxpbl, & vfxpbl, & qfxpbl real, dimension( kts:kte+1 ) , & intent(in ) :: z real, dimension( kts:kte ) , & intent(inout) :: q2 integer :: k real :: s2l,q2l,deltaz,akml,akhl,en2,pr,bpr,dis,rc02 real :: suk,svk,gthvk,govrthvk,pru,prv real :: thm,disel rc02=2.0/(c0*c0) main_integration: do k = kts+1,kte deltaz=0.5*(z(k+1)-z(k-1)) s2l=s2(k) q2l=q2(k) suk=(uxk(k)-uxk(k-1))/deltaz svk=(vxk(k)-vxk(k-1))/deltaz gthvk=(thvxk(k)-thvxk(k-1))/deltaz govrthvk=g/(0.5*(thvxk(k)+thvxk(k-1))) akml=akm(k) akhl=akh(k) en2=ri(k)*s2l thm=(thxk(k)+thxk(k-1))*0.5 if(pblflg.and.k.le.kpbl)then pru=(akml*(suk-hgamu/hpbl-ufxpbl(k)))*suk prv=(akml*(svk-hgamv/hpbl-vfxpbl(k)))*svk else pru=akml*suk*suk prv=akml*svk*svk endif pr=pru+prv if(pblflg.and.k.le.kpbl)then bpr=(akhl*(gthvk-mf(k)-(hgamq/hpbl+qfxpbl(k))*p608*thm))*govrthvk else bpr=akhl*gthvk*govrthvk endif disel=min(delxy,ceps*el(k)) dis=(q2l)**1.5/disel q2l=q2l+2.0*(pr-bpr-dis)*dtturbl q2(k)=amax1(q2l,epsq2l) enddo main_integration q2(kts)=amax1(rc02*ustar*ustar,epsq2l) end subroutine prodq2 subroutine vdifq(lmh,dtdif,q2,el,z, & akhk,ptke1, & hgame,hpbl,pblflg,kpbl, & efxpbl, & ids,ide,jds,jde,kds,kde, & ims,ime,jms,jme,kms,kme, & its,ite,jts,jte,kts,kte) implicit none real,parameter :: c_k=1.0,esq=5.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 ) :: lmh,kpbl real, intent(in ) :: dtdif,hpbl,efxpbl logical, intent(in ) :: pblflg real, dimension( kts:kte ) , & intent(in ) :: hgame, & ptke1 real, dimension( kts+1:kte ) , & intent(in ) :: el, & akhk real, dimension( kts:kte+1 ) , & intent(in ) :: z real, dimension( kts:kte ) , & intent(inout) :: q2 integer :: k real :: aden,akqs,bden,besh,besm,cden,cf,dtozs,ell,eloq2,eloq4 real :: elqdz,esh,esm,esqhf,ghl,gml,q1l,rden,rdz real :: zak real, dimension( kts+1:kte ) :: zfacentk real, dimension( kts+2:kte ) :: akq, & cm, & cr, & dtoz, & rsq2 esqhf=0.5*esq do k = kts+1,kte zak=0.5*(z(k)+z(k-1)) zfacentk(k)=(zak/hpbl)**3.0 enddo do k = kte,kts+2,-1 dtoz(k)=(dtdif+dtdif)/(z(k+1)-z(k-1)) akq(k)=c_k*(akhk(k)/(z(k+1)-z(k-1))+akhk(k-1)/(z(k)-z(k-2))) akq(k)=akq(k)*ptke1(k) cr(k)=-dtoz(k)*akq(k) enddo akqs=c_k*akhk(kts+1)/(z(kts+2)-z(kts)) akqs=akqs*ptke1(kts+1) cm(kte)=dtoz(kte)*akq(kte)+1. rsq2(kte)=q2(kte) do k = kte-1,kts+2,-1 cf=-dtoz(k)*akq(k+1)/cm(k+1) cm(k)=-cr(k+1)*cf+(akq(k+1)+akq(k))*dtoz(k)+1. rsq2(k)=-rsq2(k+1)*cf+q2(k) if(pblflg.and.k.lt.kpbl) then rsq2(k)=rsq2(k)-dtoz(k)*(2.0*hgame(k)/hpbl)*akq(k+1)*(z(k+1)-z(k)) & +dtoz(k)*(2.0*hgame(k-1)/hpbl)*akq(k)*(z(k)-z(k-1)) rsq2(k)=rsq2(k)-dtoz(k)*2.0*efxpbl*zfacentk(k+1) & +dtoz(k)*2.0*efxpbl*zfacentk(k) endif enddo dtozs=(dtdif+dtdif)/(z(kts+2)-z(kts)) cf=-dtozs*akq(lmh+2)/cm(lmh+2) if(pblflg.and.((lmh+1).lt.kpbl)) then q2(lmh+1)=(dtozs*akqs*q2(lmh)-rsq2(lmh+2)*cf+q2(lmh+1) & -dtozs*(2.0*hgame(lmh+1)/hpbl)*akq(lmh+2)*(z(lmh+2)-z(lmh+1)) & +dtozs*(2.0*hgame(lmh)/hpbl)*akqs*(z(lmh+1)-z(lmh))) q2(lmh+1)=q2(lmh+1)-dtozs*2.0*efxpbl*zfacentk(lmh+2) & +dtozs*2.0*efxpbl*zfacentk(lmh+1) q2(lmh+1)=q2(lmh+1)/((akq(lmh+2)+akqs)*dtozs-cr(lmh+2)*cf+1.) else q2(lmh+1)=(dtozs*akqs*q2(lmh)-rsq2(lmh+2)*cf+q2(lmh+1)) & /((akq(lmh+2)+akqs)*dtozs-cr(lmh+2)*cf+1.) endif do k = lmh+2,kte q2(k)=(-cr(k)*q2(k-1)+rsq2(k))/cm(k) enddo end subroutine vdifq subroutine shinhonginit(rublten,rvblten,rthblten,rqvblten, & rqcblten,rqiblten, & tke_pbl, & 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 real,parameter :: epsq2l = 0.01 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 real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: & tke_pbl 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. tke_pbl(i,k,j) = epsq2l/2. 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 shinhonginit function pu(d,h) implicit none real :: pu real,parameter :: pmin = 0.0,pmax = 1.0 real,parameter :: a1 = 1.0, a2 = 0.070, a3 = 1.0, a4 = 0.142, a5 = 0.071 real,parameter :: b1 = 2.0, b2 = 0.6666667 real :: d,h,doh,num,den if (h .ne. 0) then doh=d/h num=a1*(doh)**b1+a2*(doh)**b2 den=a3*(doh)**b1+a4*(doh)**b2+a5 pu=num/den else pu=1. endif pu=max(pu,pmin) pu=min(pu,pmax) return end function function pq(d,h) implicit none real :: pq real,parameter :: pmin = 0.0,pmax = 1.0 real,parameter :: a1 = 1.0, a2 = -0.098, a3 = 1.0, a4 = 0.106, a5 = 0.5 real,parameter :: b1 = 2.0 real :: d,h,doh,num,den if (h .ne. 0) then doh=d/h num=a1*(doh)**b1+a2 den=a3*(doh)**b1+a4 pq=a5*num/den+(1.-a5) else pq=1. endif pq=max(pq,pmin) pq=min(pq,pmax) return end function function pthnl(d,h) implicit none real :: pthnl real,parameter :: pmin = 0.0,pmax = 1.0 real,parameter :: a1 = 1.000, a2 = 0.936, a3 = -1.110, & a4 = 1.000, a5 = 0.312, a6 = 0.329, a7 = 0.243 real,parameter :: b1 = 2.0, b2 = 0.875 real :: d,h,doh,num,den if (h .ne. 0) then doh=d/h num=a1*(doh)**b1+a2*(doh)**b2+a3 den=a4*(doh)**b1+a5*(doh)**b2+a6 pthnl=a7*num/den+(1.-a7) else pthnl=1. endif pthnl=max(pthnl,pmin) pthnl=min(pthnl,pmax) return end function function pthl(d,h) implicit none real :: pthl real,parameter :: pmin = 0.0,pmax = 1.0 real,parameter :: a1 = 1.000, a2 = 0.870, a3 = -0.913, & a4 = 1.000, a5 = 0.153, a6 = 0.278, a7 = 0.280 real,parameter :: b1 = 2.0, b2 = 0.5 real :: d,h,doh,num,den if (h .ne. 0) then doh=d/h num=a1*(doh)**b1+a2*(doh)**b2+a3 den=a4*(doh)**b1+a5*(doh)**b2+a6 pthl=a7*num/den+(1.-a7) else pthl=1. endif pthl=max(pthl,pmin) pthl=min(pthl,pmax) return end function function ptke(d,h) implicit none real :: ptke real,parameter :: pmin = 0.0,pmax = 1.0 real,parameter :: a1 = 1.000, a2 = 0.070, & a3 = 1.000, a4 = 0.142, a5 = 0.071 real,parameter :: b1 = 2.0, b2 = 0.6666667 real :: d,h,doh,num,den if (h .ne. 0) then doh=d/h num=a1*(doh)**b1+a2*(doh)**b2 den=a3*(doh)**b1+a4*(doh)**b2+a5 ptke=num/den else ptke=1. endif ptke=max(ptke,pmin) ptke=min(ptke,pmax) return end function end module module_bl_shinhong