!> \file ysuvdif.F90 !! This file contains the CCPP-compliant YSU scheme which computes !! subgrid vertical turbulence mixing using traditional K-profile method !! Please refer to (Hong, Noh and Dudhia, 2006, MWR). !! !! Subroutine 'ysuvdif_run' computes subgrid vertical turbulence mixing !! using YSU K-profile method !! !---------------------------------------------------------------------- module ysuvdif contains subroutine ysuvdif_init (do_ysu,errmsg,errflg) logical, intent(in) :: do_ysu character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! Initialize CCPP error handling variables errmsg = '' errflg = 0 ! Consistency checks if (.not. do_ysu) then write(errmsg,fmt='(*(a))') 'Logic error: do_ysu = .false.' errflg = 1 return end if end subroutine ysuvdif_init !> \defgroup YSU FV3GFS ysuvdif_run Main !! \brief This subroutine contains all of the logic for the !! YSU scheme. !! !> \section arg_table_ysuvdif_run Argument Table !! \htmlinclude ysuvdif_run.html !! !------------------------------------------------------------------------------- subroutine ysuvdif_run(im,km,ux,vx,tx,qx,p2d,p2di,pi2d,karman, & utnp,vtnp,ttnp,qtnp, & swh,hlw,xmu,ntrac,ndiff,ntcw,ntiw, & phii,phil,psfcpa, & zorl,stress,hpbl,psim,psih, & landmask,heat,evap,wspd,br, & g,rd,cp,rv,ep1,ep2,xlv, & dusfc,dvsfc,dtsfc,dqsfc, & dt,kpbl1d,u10,v10,lssav,ldiag3d,qdiag3d, & flag_for_pbl_generic_tend,ntoz,ntqv,dtend,dtidx, & index_of_temperature,index_of_x_wind,index_of_y_wind, & index_of_process_pbl,errmsg,errflg ) use machine , only : kind_phys ! !------------------------------------------------------------------------------- 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 real(kind=kind_phys),parameter :: rcl = 1.0 real(kind=kind_phys),intent(in) :: karman integer,parameter :: imvdif = 1 integer,parameter :: ysu_topdown_pblmix = 1 ! !------------------------------------------------------------------------------------- ! input variables integer, intent(in ) :: im,km,ntrac,ndiff,ntcw,ntiw,ntoz real(kind=kind_phys), intent(in ) :: g,cp,rd,rv,ep1,ep2,xlv,dt real(kind=kind_phys), dimension( :,: ), & intent(in) :: pi2d,p2d,phil,ux,vx,swh,hlw,tx real(kind=kind_phys), dimension( :,:,: ) , & intent(in ) :: qx real(kind=kind_phys), dimension( :,: ) , & intent(in ) :: p2di,phii real(kind=kind_phys), dimension( : ) , & intent(in) :: stress,zorl,heat,evap,wspd,br,psim,psih,psfcpa, & u10,v10,xmu integer, dimension(:) ,& intent(in ) :: landmask logical, intent(in ) :: lssav, ldiag3d, qdiag3d, & flag_for_pbl_generic_tend ! !---------------------------------------------------------------------------------- ! input/output variables ! real(kind=kind_phys), dimension( :,: ) , & intent(inout) :: utnp,vtnp,ttnp real(kind=kind_phys), dimension( :,:,: ) , & intent(inout) :: qtnp real(kind=kind_phys), optional, intent(inout) :: dtend(:,:,:) integer, intent(in) :: dtidx(:,:), ntqv, index_of_temperature, & index_of_x_wind, index_of_y_wind, index_of_process_pbl ! !--------------------------------------------------------------------------------- ! output variables integer, dimension( : ), intent(out ) :: kpbl1d real(kind=kind_phys), dimension( : ), & intent(out) :: hpbl real(kind=kind_phys), dimension( : ), & intent(out) :: dusfc,dvsfc, dtsfc,dqsfc ! error messages character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! !-------------------------------------------------------------------------------- ! ! local vars ! real(kind=kind_phys), dimension( im ) :: hol real(kind=kind_phys), dimension( im, km+1 ) :: zq ! real(kind=kind_phys), dimension( im, km ) :: & thx,thvx,thlix, & del, & dza, & dzq, & xkzom, & xkzoh, & za ! real(kind=kind_phys), dimension( im ) :: & rhox, & govrth, & zl1,thermal, & wscale, & hgamt,hgamq, & brdn,brup, & phim,phih, & prpbl, & wspd1,thermalli ! real(kind=kind_phys), dimension( im, km ) :: xkzm,xkzh, & f1,f2, & r1,r2, & ad,au, & cu, & al, & xkzq, & zfac, & rhox2, & hgamt2 ! real(kind=kind_phys), dimension( im ) :: & brcr, & sflux, & zol1, & brcr_sbro ! real(kind=kind_phys), dimension( im ) :: xland real(kind=kind_phys), dimension( im ) :: ust real(kind=kind_phys), dimension( im ) :: hfx real(kind=kind_phys), dimension( im ) :: qfx real(kind=kind_phys), dimension( im ) :: znt real(kind=kind_phys), dimension( im ) :: uox real(kind=kind_phys), dimension( im ) :: vox ! real(kind=kind_phys), dimension( im, km, ndiff) :: r3,f3 integer, dimension( im ) :: kpbl,kpblold ! logical, dimension( im ) :: pblflg, & sfcflg, & stable, & cloudflg logical :: definebrup ! integer :: n,i,k,l,ic,is,kk integer :: klpbl, ktrace1, ktrace2, ktrace3 ! ! 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, rovcp ! real(kind=kind_phys), dimension( im, km ) :: wscalek,wscalek2 real(kind=kind_phys), dimension( im ) :: wstar real(kind=kind_phys), dimension( im ) :: delta real(kind=kind_phys), dimension( im, km ) :: xkzml,xkzhl, & zfacent,entfac real(kind=kind_phys), dimension( im ) :: 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 integer :: idtend ! !------------------------------------------------------------------------------- ! ! Initialize CCPP error handling variables errmsg = '' errflg = 0 klpbl = km ! rovcp=rd/cp cont=cp/g conq=xlv/g conw=1./g conwrc = conw*sqrt(rcl) conpr = bfac*karman*sfcfrac ! ! change xland values do i=1,im if(landmask(i).eq.0) then !ocean xland(i) = 2 else xland(i) = 1 !land end if end do ! do k = 1,km do i = 1,im thx(i,k) = tx(i,k)/pi2d(i,k) thlix(i,k) = (tx(i,k)-xlv*qx(i,k,ntcw)/cp-2.834E6*qx(i,k,ntiw)/cp)/pi2d(i,k) enddo enddo ! do k = 1,km do i = 1,im tvcon = (1.+ep1*qx(i,k,1)) thvx(i,k) = thx(i,k)*tvcon enddo enddo ! do i = 1,im tvcon = (1.+ep1*qx(i,1,1)) rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon) govrth(i) = g/thx(i,1) hfx(i) = heat(i)*rhox(i)*cp ! reset to the variable in WRF qfx(i) = evap(i)*rhox(i) ! reset to the variable in WRF ust(i) = sqrt(stress(i)) ! reset to the variable in WRF znt(i) = 0.01*zorl(i) ! reset to the variable in WRF uox(i) = 0.0 vox(i) = 0.0 enddo ! !-----compute the height of full- and half-sigma levels above ground ! level, and the layer thicknesses. ! do i = 1,im zq(i,1) = 0. enddo ! do k = 1,km do i = 1,im zq(i,k+1) = phii(i,k+1)*conw tvcon = (1.+ep1*qx(i,k,1)) rhox2(i,k) = p2d(i,k)/(rd*tx(i,k)*tvcon) enddo enddo ! do k = 1,km do i = 1,im za(i,k) = phil(i,k)*conw dzq(i,k) = zq(i,k+1)-zq(i,k) del(i,k) = p2di(i,k)-p2di(i,k+1) enddo enddo ! do i = 1,im dza(i,1) = za(i,1) enddo ! do k = 2,km do i = 1,im dza(i,k) = za(i,k)-za(i,k-1) enddo enddo ! write(0,*)"===CALLING ysu; input:" ! print*,"t:",tx(1,1),tx(1,2),tx(1,km) ! print*,"u:",ux(1,1),ux(1,2),ux(1,km) ! print*,"v:",vx(1,1),vx(1,2),vx(1,km) ! print*,"q:",qx(1,1,1),qx(1,2,1),qx(1,km,1) ! print*,"exner:",pi2d(1,1),pi2d(1,2),pi2d(1,km) ! print*,"phii:",zq(1,1),zq(1,2),zq(1,km+1) ! print*,"phil:",za(1,1),za(1,2),za(1,km) ! print*,"p2d:",p2d(1,1),p2d(1,2),p2d(1,km) ! print*,"p2di:",p2di(1,1),p2di(1,2),p2di(1,km+1) ! print *,"del:",del(1,1),del(1,2),del(1,km) ! print*,"znt,ust,wspd:",znt(1),ust(1),wspd(1) ! print*,"hfx,qfx,xland:",hfx(1),qfx(1),xland(1) ! print*,"rd,rv,g:",rd,rv,g ! print*,"ep1,ep2,xlv:",ep1,ep2,xlv ! print*,"br,psim,psih:",br(1),psim(1),psih(1) ! print*,"u10,v10:",u10(1),v10(1) ! print*,"psfcpa,cp:",psfcpa(1),cp ! print*,"ntrac,ndiff,ntcw,ntiw:",ntrac,ndiff,ntcw,ntiw ! ! !-----initialize vertical tendencies and ! ! utnp(:,:) = 0. ! vtnp(:,:) = 0. ! ttnp(:,:) = 0. ! qtnp(:,:,:) = 0. ! do i = 1,im 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 ! !---- compute vertical diffusion ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! compute preliminary variables ! dtstep = dt dt2 = 2.*dtstep rdt = 1./dt2 ! do i = 1,im 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 = 1,klpbl do i = 1,im wscalek(i,k) = 0.0 wscalek2(i,k) = 0.0 enddo enddo ! do k = 1,klpbl do i = 1,im zfac(i,k) = 0.0 enddo enddo do k = 1,klpbl-1 do i = 1,im xkzom(i,k) = xkzminm xkzoh(i,k) = xkzminh enddo enddo ! do i = 1,im dusfc(i) = 0. dvsfc(i) = 0. dtsfc(i) = 0. dqsfc(i) = 0. enddo ! do i = 1,im 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 ! ! compute the first guess of pbl height ! do i = 1,im stable(i) = .false. brup(i) = br(i) brcr(i) = brcr_ub enddo ! do k = 2,klpbl do i = 1,im 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 = 1,im 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 = 1,im 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 ! ! compute the surface variables for pbl height estimation ! under unstable conditions ! do i = 1,im 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 ! ! enhance the pbl height by considering the thermal ! do i = 1,im if(pblflg(i))then kpbl(i) = 1 hpbl(i) = zq(i,1) endif enddo ! do i = 1,im if(pblflg(i))then stable(i) = .false. brup(i) = br(i) brcr(i) = brcr_ub endif enddo ! do k = 2,klpbl do i = 1,im 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 ! ! enhance pbl by theta-li ! if (ysu_topdown_pblmix.eq.1)then do i = 1,im kpblold(i) = kpbl(i) definebrup=.false. do k = kpblold(i), km-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 !overwrite brup brdn values brdn(i)=bruptmp definebrup=.true. pblflg(i)=.true. endif enddo enddo endif do i = 1,im 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 ! ! stable boundary layer ! do i = 1,im 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 = 1,im 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 = 1,im 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 = 1,im 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 = 1,im 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 ! ! estimate the entrainment parameters ! do i = 1,im 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,k,ntcw)+qx(i,k,ntiw)).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 is ws at full level rvls=100.*6.112*EXP(17.67*(templ-273.16)/(templ-29.65))*(ep2/p2di(i,k+1)) temps=templ + ((qx(i,k,1)+qx(i,k,ntcw))-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,1)+qx(i,k,ntcw))-rvls,0.) !entrainment efficiency dthvx(i) = (thlix(i,k+2)+thx(i,k+2)*ep1*(qx(i,k+2,1)+qx(i,k+2,ntcw))) & - (thlix(i,k) + thx(i,k) *ep1*(qx(i,k,1) +qx(i,k,ntcw))) 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=swh(i,kk)*xmu(i)+hlw(i,kk) !radiative heating rate temp/s radflux=radflux*cp/g*(p2di(i,kk)-p2di(i,kk+1)) ! converts temp/s to W/m^2 if (radflux < 0.0 ) radsum=abs(radflux)+radsum enddo radsum=max(radsum,0.0) !recompute entrainment from sfc thermals 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))) !entrainment from PBL top thermals 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)) ! this is wstar3(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)) !wstar3_2 bfx0 = max(radsum/rhox2(i,k)/cp,0.) wstar3_2(i) = (g/thvx(i,k)*bfx0*hpbl(i)) !recompute hgamt 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,1)-qx(i,k,1),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 = 1,klpbl do i = 1,im 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 ! ! compute diffusion coefficients below pbl ! do k = 1,klpbl do i = 1,im 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 !Do not include xkzm at kpbl-1 since it changes entrainment 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 ! ! compute diffusion coefficients over pbl (free atmosphere) ! do k = 1,km-1 do i = 1,im 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.ntcw.ge.2.and.ntiw.ge.2)then if((qx(i,k,ntcw)+qx(i,k,ntiw)).gt.0.01e-3.and.(qx(i & ,k+1,ntcw)+qx(i,k+1,ntiw)).gt.0.01e-3)then ! in cloud qmean = 0.5*(qx(i,k,1)+qx(i,k+1,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 ! unstable regime 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 ! stable regime 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 ! ! compute tridiagonal matrix elements for heat ! do k = 1,km do i = 1,im au(i,k) = 0. al(i,k) = 0. ad(i,k) = 0. f1(i,k) = 0. enddo enddo ! do i = 1,im ad(i,1) = 1. f1(i,1) = thx(i,1)-300.+hfx(i)/cont/del(i,1)*dt2 enddo ! do k = 1,km-1 do i = 1,im 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) enddo enddo ! ! copies here to avoid duplicate input args for tridin ! do k = 1,km do i = 1,im cu(i,k) = au(i,k) r1(i,k) = f1(i,k) enddo enddo ! call tridin_ysu(al,ad,cu,r1,au,f1,im,km,1) ! ! recover tendencies of heat ! do k = km,1,-1 do i = 1,im 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) enddo enddo if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend) then idtend = dtidx(index_of_temperature,index_of_process_pbl) if(idtend>=1) then dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f1-thx+300.)*rdt*pi2d endif endif ! ! compute tridiagonal matrix elements for moisture, clouds, and gases ! do k = 1,km do i = 1,im au(i,k) = 0. al(i,k) = 0. ad(i,k) = 0. enddo enddo ! do ic = 1,ndiff do i = 1,im do k = 1,km f3(i,k,ic) = 0. enddo enddo enddo ! do i = 1,im ad(i,1) = 1. f3(i,1,1) = qx(i,1,1)+qfx(i)*g/del(i,1)*dt2 enddo ! if(ndiff.ge.2) then do ic = 2,ndiff do i = 1,im f3(i,1,ic) = qx(i,1,ic) enddo enddo endif ! do k = 1,km-1 do i = 1,im if(k.ge.kpbl(i)) then xkzq(i,k) = xkzh(i,k) endif enddo enddo ! do k = 1,km-1 do i = 1,im 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,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,1) else f3(i,k+1,1) = qx(i,k+1,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 do k = 1,km-1 do i = 1,im f3(i,k+1,ic) = qx(i,k+1,ic) enddo enddo enddo endif ! ! copies here to avoid duplicate input args for tridin ! do k = 1,km do i = 1,im cu(i,k) = au(i,k) enddo enddo ! do ic = 1,ndiff do k = 1,km do i = 1,im r3(i,k,ic) = f3(i,k,ic) enddo enddo enddo ! ! solve tridiagonal problem for moisture, clouds, and gases ! call tridin_ysu(al,ad,cu,r3,au,f3,im,km,ndiff) ! ! recover tendencies of heat and moisture ! do k = km,1,-1 do i = 1,im qtend = (f3(i,k,1)-qx(i,k,1))*rdt qtnp(i,k,1) = qtnp(i,k,1)+qtend dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k) enddo enddo if(lssav .and. ldiag3d .and. qdiag3d .and. .not. flag_for_pbl_generic_tend) then idtend = dtidx(ntqv+100,index_of_process_pbl) if(idtend>=1) then dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f3(:,:,1)-qx(:,:,1))*rdt endif endif ! if(ndiff.ge.2) then do ic = 2,ndiff do k = km,1,-1 do i = 1,im qtend = (f3(i,k,ic)-qx(i,k,ic))*rdt qtnp(i,k,ic) = qtnp(i,k,ic)+qtend enddo enddo enddo if(lssav .and. ldiag3d .and. ntoz>0 .and. qdiag3d .and. & & .not. flag_for_pbl_generic_tend) then idtend = dtidx(100+ntoz,index_of_process_pbl) if(idtend>=1) then dtend(:,:,idtend) = dtend(:,:,idtend) + f3(:,:,ntoz)-qx(:,:,ntoz) endif endif endif ! ! compute tridiagonal matrix elements for momentum ! do i = 1,im do k = 1,km au(i,k) = 0. al(i,k) = 0. ad(i,k) = 0. f1(i,k) = 0. f2(i,k) = 0. enddo enddo ! do i = 1,im ad(i,1) = 1.+ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 & *(wspd1(i)/wspd(i))**2 f1(i,1) = ux(i,1)+uox(i)*ust(i)**2*g/del(i,1)*dt2/wspd1(i) f2(i,1) = vx(i,1)+vox(i)*ust(i)**2*g/del(i,1)*dt2/wspd1(i) enddo ! do k = 1,km-1 do i = 1,im 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) enddo enddo ! ! copies here to avoid duplicate input args for tridin ! do k = 1,km do i = 1,im cu(i,k) = au(i,k) r1(i,k) = f1(i,k) r2(i,k) = f2(i,k) enddo enddo ! ! solve tridiagonal problem for momentum ! call tridi1n(al,ad,cu,r1,r2,au,f1,f2,im,km,1) ! ! recover tendencies of momentum ! do k = km,1,-1 do i = 1,im 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 if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend) then idtend = dtidx(index_of_x_wind,index_of_process_pbl) if(idtend>=1) then dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f1-ux)*rdt endif idtend = dtidx(index_of_y_wind,index_of_process_pbl) if(idtend>=1) then dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f2-vx)*rdt endif endif ! !---- end of vertical diffusion ! do i = 1,im kpbl1d(i) = kpbl(i) enddo ! ! end subroutine ysuvdif_run !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- subroutine tridi1n(cl,cm,cu,r1,r2,au,f1,f2,im,km,nt) use machine , only : kind_phys !------------------------------------------------------------------------------- implicit none !------------------------------------------------------------------------------- ! integer, intent(in ) :: im, km, nt ! real(kind=kind_phys), dimension( im, 2:km+1 ) , & intent(in ) :: cl ! real(kind=kind_phys), dimension( im, km ) , & intent(in ) :: cm, & r1 real(kind=kind_phys), dimension( im, km,nt ) , & intent(in ) :: r2 ! real(kind=kind_phys), dimension( im, km ) , & intent(inout) :: au, & cu, & f1 real(kind=kind_phys), dimension( im, km,nt ) , & intent(inout) :: f2 ! real(kind=kind_phys) :: fk integer :: i,k,l,n,it ! !------------------------------------------------------------------------------- ! l = im n = km ! do i = 1,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 = 1,l fk = 1./cm(i,1) f2(i,1,it) = fk*r2(i,1,it) enddo enddo ! do k = 2,n-1 do i = 1,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 = 2,n-1 do i = 1,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 = 1,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 = 1,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,1,-1 do i = 1,l f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1) enddo enddo ! do it = 1,nt do k = n-1,1,-1 do i = 1,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,im,km,nt) use machine , only : kind_phys !------------------------------------------------------------------------------- implicit none !------------------------------------------------------------------------------- ! integer, intent(in ) :: im, km, nt ! real(kind=kind_phys), dimension( im, 2:km+1 ) , & intent(in ) :: cl ! real(kind=kind_phys), dimension( im, km ) , & intent(in ) :: cm real(kind=kind_phys), dimension( im, km,nt ) , & intent(in ) :: r2 ! real(kind=kind_phys), dimension( im, km ) , & intent(inout) :: au, & cu real(kind=kind_phys), dimension( im, km,nt ) , & intent(inout) :: f2 ! real(kind=kind_phys) :: fk integer :: i,k,l,n,it ! !------------------------------------------------------------------------------- ! l = im n = km ! do it = 1,nt do i = 1,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 = 2,n-1 do i = 1,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 = 1,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,1,-1 do i = 1,l f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it) enddo enddo enddo ! end subroutine tridin_ysu !------------------------------------------------------------------------------- end module ysuvdif !-------------------------------------------------------------------------------