!> \file shinhongvdif.F90 !! This file contains the CCPP-compliant Shinhong (saYSU) scheme which computes !! subgrid vertical turbulence mixing using traditional K-profile method !! Please refer to (Shin and Hong, 2013,2015). !! !! Subroutine 'shinhongvdif_run' computes subgrid vertical turbulence mixing !! using scale-aware YSU K-profile method !! !---------------------------------------------------------------------- module shinhongvdif contains subroutine shinhongvdif_init (shinhong,errmsg,errflg) logical, intent(in) :: shinhong character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! Initialize CCPP error handling variables errmsg = '' errflg = 0 ! Consistency checks if (.not. shinhong) then write(errmsg,fmt='(*(a))') 'Logic error: shinhong = .false.' errflg = 1 return end if end subroutine shinhongvdif_init !> \defgroup SHINHONG FV3GFS shinhongvdif_run Main !! \brief This subroutine contains all of the logic for the !! scale-aware Shinhong scheme. !! !> \section arg_table_shinhongvdif_run Argument Table !! \htmlinclude shinhongvdif_run.html !! !------------------------------------------------------------------------------- subroutine shinhongvdif_run(im,km,ux,vx,tx,qx,p2d,p2di,pi2d,karman, & utnp,vtnp,ttnp,qtnp,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, & dx,lssav,ldiag3d, & flag_for_pbl_generic_tend,ntoz,ntqv,dtend,dtidx, & index_of_process_pbl,index_of_temperature,index_of_x_wind, & index_of_y_wind,errmsg,errflg ) use machine , only : kind_phys ! !------------------------------------------------------------------------------- implicit none !------------------------------------------------------------------------------- ! ! the shinhongpbl (shin and hong 2015) is based on the les study of shin ! and hong (2013). the major ingredients of the shinhongpbl are ! 1) the prescribed nonlocal heat transport profile fit to the les and ! 2) inclusion of explicit scale dependency functions for vertical ! transport in convective pbl. ! so, the shinhongpbl works at the gray zone resolution of convective pbl. ! note that honnert et al. (2011) first suggested explicit scale dependency ! function, and shin and hong (2013) further classified the function by ! stability (u*/w*) in convective pbl and calculated the function for ! nonlocal and local transport separately. ! vertical mixing in the stable boundary layer and free atmosphere follows ! hong (2010) and hong et al. (2006), same as the ysupbl scheme. ! ! shinhongpbl: ! coded and implemented by hyeyum hailey shin (ncar) ! summer 2014 ! ! ysupbl: ! coded by song-you hong (yonsei university) and implemented by ! song-you hong (yonsei university) and jimy dudhia (ncar) ! summer 2002 ! ! references: ! shin and hong (2015) mon. wea. rev. ! shin and hong (2013) j. atmos. sci. ! honnert, masson, and couvreux (2011) j. atmos. sci. ! hong (2010) quart. j. roy. met. soc ! hong, noh, and dudhia (2006), mon. wea. rev. ! !------------------------------------------------------------------------------- ! real(kind=kind_phys),parameter :: xkzminm = 0.1,xkzminh = 0.01 real(kind=kind_phys),parameter :: 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 :: ckz = 0.001,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),intent(in) :: karman real(kind=kind_phys),parameter :: corf=0.000073 real(kind=kind_phys),parameter :: rcl = 1.0 integer,parameter :: imvdif = 1 integer,parameter :: shinhong_tke_diag = 0 ! ! tunable parameters for tke ! real(kind=kind_phys),parameter :: epsq2l = 0.01,c_1 = 1.0,gamcre = 0.224 ! ! tunable parameters for prescribed nonlocal transport profile ! real(kind=kind_phys),parameter :: mltop = 1.0,sfcfracn1 = 0.075 real(kind=kind_phys),parameter :: nlfrac = 0.7,enlfrac = -0.4 real(kind=kind_phys),parameter :: a11 = 1.0,a12 = -1.15 real(kind=kind_phys),parameter :: ezfac = 1.5 real(kind=kind_phys),parameter :: cpent = -0.4,rigsmax = 100. real(kind=kind_phys),parameter :: entfmin = 1.0, entfmax = 5.0 ! 1D in integer, intent(in ) :: im,km,ntrac,ndiff,ntcw,ntiw,ntoz real(kind=kind_phys), intent(in ) :: g,cp,rd,rv,ep1,ep2,xlv,dt logical, intent(in ) :: lssav, ldiag3d, flag_for_pbl_generic_tend ! 3D in real(kind=kind_phys), dimension(:,:) , & intent(in ) :: phil, & pi2d, & p2d, & ux, & vx, & tx real(kind=kind_phys), dimension(:,:,:) , & intent(in ) :: qx real(kind=kind_phys), dimension(:,:) , & intent(in ) :: p2di, & phii ! 3D in&out real(kind=kind_phys), dimension(:,:) , & intent(inout) :: utnp, & vtnp, & ttnp real(kind=kind_phys), dimension(:,:,:) , & intent(inout) :: qtnp ! 2D in integer, dimension(:) , & intent(in ) :: landmask real(kind=kind_phys), dimension(:) , & intent(in ) :: heat, & evap, & br, & psim, & psih, & psfcpa, & stress, & zorl, & wspd, & u10, & v10, & dx ! 2D: out integer, dimension(:) , & intent(out ) :: kpbl1d real(kind=kind_phys), dimension(:) , & intent(out ) :: hpbl, & dusfc, & dvsfc, & dtsfc, & dqsfc real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:) integer, intent(in) :: dtidx(:,:), index_of_process_pbl, ntqv, & index_of_x_wind, index_of_y_wind, index_of_temperature ! Index within dtend third dimension for tendency of interest: integer :: idtend ! error messages character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! ! local vars ! integer :: n,i,k,l,ic integer :: klpbl integer :: lmh,lmxl,kts,kte,its,ite ! real(kind=kind_phys) :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0 real(kind=kind_phys) :: ss,ri,qmean,tmean,alpha,chi,zk,rl2,dk,sri real(kind=kind_phys) :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz real(kind=kind_phys) :: utend,vtend,ttend,qtend real(kind=kind_phys) :: dtstep,govrthv real(kind=kind_phys) :: cont, conq, conw, conwrc real(kind=kind_phys) :: delxy,pu1,pth1,pq1 real(kind=kind_phys) :: dex,hgame_c real(kind=kind_phys) :: zfacdx real(kind=kind_phys) :: amf1,amf2,bmf2,amf3,bmf3,amf4,bmf4,sflux0,snlflux0 real(kind=kind_phys) :: mlfrac,ezfrac,sfcfracn real(kind=kind_phys) :: uwst,uwstx,csfac 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 ! integer, dimension(im) :: kpbl real(kind=kind_phys), dimension(im) :: hol real(kind=kind_phys), dimension(im) :: deltaoh real(kind=kind_phys), dimension(im) :: rigs, & enlfrac2, & cslen real(kind=kind_phys), dimension(im) :: & rhox, & govrth, & zl1,thermal, & wscale, & hgamt,hgamq, & brdn,brup, & phim,phih, & prpbl, & wspd1, & ust,hfx,qfx,znt, & xland real(kind=kind_phys), dimension(im) :: & ust3, & wstar3, & wstar,delta, & hgamu,hgamv, & wm2, we, & bfxpbl, & hfxpbl,qfxpbl, & ufxpbl,vfxpbl, & dthvx real(kind=kind_phys), dimension(im) :: & brcr, & sflux, & zol1, & brcr_sbro real(kind=kind_phys), dimension(im) :: & efxpbl, & hpbl_cbl, & epshol, & ct ! real(kind=kind_phys), dimension(im,km) :: & xkzm,xkzh, & f1,f2, & r1,r2, & ad,au, & cu, & al, & xkzq, & zfac real(kind=kind_phys), dimension(im,km) :: & thx,thvx, & del, & dza, & dzq, & xkzom, & xkzoh, & za real(kind=kind_phys), dimension(im,km) :: & wscalek real(kind=kind_phys), dimension(im,km) :: & xkzml,xkzhl, & zfacent,entfac real(kind=kind_phys), dimension(im,km) :: & mf, & zfacmf, & entfacmf real(kind=kind_phys), dimension(im,km) :: & q2x, & hgame2d, & tflux_e, & qflux_e, & tvflux_e real(kind=kind_phys), dimension( im, km+1 ) :: zq real(kind=kind_phys), dimension( im, km, ndiff ) :: r3,f3 ! real(kind=kind_phys), dimension( km ) :: & uxk,vxk, & txk,thxk,thvxk, & q2xk, & hgame real(kind=kind_phys), dimension( km ) :: & ps1d,pb1d,eps1d,pt1d, & xkze1d,eflx_l1d,eflx_nl1d, & ptke1 real(kind=kind_phys), dimension( 2:km ) :: & s2,gh,rig,el, & akmk,akhk, & mfk,ufxpblk,vfxpblk,qfxpblk real(kind=kind_phys), dimension( km+1 ) :: zqk real(kind=kind_phys), dimension(im,km) :: dz8w2d ! logical, dimension(im) :: pblflg, & sfcflg, & stable logical, dimension( ndiff ) :: ifvmix ! !------------------------------------------------------------------------------- ! Initialize CCPP error handling variables errmsg = '' errflg = 0 its = 1 ite = im kts = 1 kte = km klpbl = kte lmh = 1 lmxl = 1 ! cont=cp/g conq=xlv/g conw=1./g conwrc = conw*sqrt(rcl) conpr = bfac*karman*sfcfrac ! change xland values do i=its,ite if(landmask(i).eq.0) then !ocean xland(i) = 2 else xland(i) = 1 !land end if end do ! ! k-start index for cloud and rain ! 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,1)) thvx(i,k) = thx(i,k)*tvcon enddo enddo ! do i = its,ite 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 enddo ! !-----compute the height of full- and half-sigma levels above ground ! level, and the layer thicknesses. ! do i = its,ite zq(i,1) = 0. enddo ! do k = kts,kte do i = its,ite zq(i,k+1) = phii(i,k+1)*conw za(i,k) = phil(i,k)*conw enddo enddo ! do k = kts,kte do i = its,ite dzq(i,k) = zq(i,k+1)-zq(i,k) del(i,k) = p2di(i,k)-p2di(i,k+1) dz8w2d(i,k)=dzq(i,k) enddo enddo ! do i = its,ite dza(i,1) = za(i,1) enddo ! do k = kts+1,kte do i = its,ite dza(i,k) = za(i,k)-za(i,k-1) enddo enddo ! do i = its,ite wspd1(i) = sqrt(ux(i,1)*ux(i,1)+vx(i,1)*vx(i,1))+1.e-9 enddo ! write(0,*)"===CALLING shinhong; 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*,"dz8w2d:",dz8w2d(1,1),dz8w2d(1,2),dz8w2d(1,km) ! print *,"del:",del(1,1),del(1,2),del(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*,"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*,"dx,u10,v10:",dx(1),u10(1),v10(1) ! print*,"psfcpa,cp:",psfcpa(1),cp ! print*,"ntrac,ndiff,ntcw,ntiw:",ntrac,ndiff,ntcw,ntiw ! !---- compute vertical diffusion ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! compute preliminary variables ! 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) = 1.e-4 enddo enddo ! do k = kts,kte do i = its,ite 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 ! ! compute the first guess of pbl height ! 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 ! ! compute the surface variables for pbl height estimation ! under unstable conditions ! 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 ! ! enhance the pbl height by considering the thermal ! 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. uwst = abs(ust(i)/wstar(i)-0.5) uwstx = -80.*uwst+14. csfac = 0.5*(tanh(uwstx)+3.) cslen(i) = csfac*hpbl(i) endif enddo ! ! stable boundary layer ! 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 ! ! scale dependency for nonlocal momentum and moisture transport ! do i = its,ite pu1=pu(dx(i),cslen(i)) pq1=pq(dx(i),cslen(i)) if(pblflg(i)) then hgamu(i) = hgamu(i)*pu1 hgamv(i) = hgamv(i)*pu1 hgamq(i) = hgamq(i)*pq1 endif enddo ! ! estimate the entrainment parameters ! 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,1)-qx(i,k,1),0.0) we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i))) hfxpbl(i) = we(i)*dthx pq1=pq(dx(i),cslen(i)) qfxpbl(i) = we(i)*dqx*pq1 ! pu1=pu(dx(i),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)) rigs(i) = govrth(i)*dthvx(i)*deltaoh(i)/(dux**2.+dvx**2.) 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 ! ! compute diffusion coefficients below pbl ! 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 ! ! compute diffusion coefficients over pbl (free atmosphere) ! 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)) ! in cloud 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 qmean = 0.5*(qx(i,k,1)+qx(i,k+1,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 ! 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 ! ! prescribe nonlocal heat transport below pbl ! do i = its,ite deltaoh(i) = deltaoh(i)/hpbl(i) enddo ! 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 amf2 = -snlflux0/(mlfrac-sfcfracn) bmf2 = -mlfrac*amf2 amf3 = snlflux0*enlfrac2(i)/deltaoh(i) bmf3 = -amf3*mlfrac hfxpbl(i) = amf3+bmf3 pth1=pthnl(dx(i),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 ! ! compute tridiagonal matrix elements for heat ! 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*(-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 ! ! scale dependency for local heat transport ! zfacdx=0.2*hpbl(i)/zq(i,k+1) delxy=dx(i)*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) enddo enddo ! ! copies here to avoid duplicate input args for tridin ! 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) ! ! recover tendencies of heat ! 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) 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 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 = 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,1)+qfx(i)*g/del(i,1)*dt2 enddo ! if(ndiff.ge.2) then do ic = 2,ndiff do i = its,ite f3(i,1,ic) = qx(i,1,ic) 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,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 ! ! scale dependency for local moisture transport ! zfacdx=0.2*hpbl(i)/zq(i,k+1) delxy=dx(i)*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 do k = kts,kte-1 do i = its,ite 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 = 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 ! ! solve tridiagonal problem for moisture, clouds, and gases ! call tridin_ysu(al,ad,cu,r3,au,f3,its,ite,kts,kte,ndiff) ! ! recover tendencies of heat and moisture ! do k = kte,kts,-1 do i = its,ite 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) 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 if(lssav .and. ldiag3d .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*rdt*(f3(:,:,1)-qx(:,:,1)) endif endif ! print*,"qtnp:",maxval(qtnp(:,:,1)),minval(qtnp(:,:,1)) ! 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 do k = kte,kts,-1 do i = its,ite qtend = (f3(i,k,ic)-qx(i,k,ic))*rdt qtnp(i,k,ic) = qtnp(i,k,ic)+qtend enddo enddo endif enddo if(lssav .and. ldiag3d .and. ntoz>0 .and. & & .not. flag_for_pbl_generic_tend) then idtend=dtidx(ntoz+100,index_of_process_pbl) if(idtend>=1) then dtend(:,:,idtend) = dtend(:,:,idtend) + qtend*(f3(:,:,ntoz)-qx(:,:,ntoz)) endif endif endif ! ! compute tridiagonal matrix elements for momentum ! 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.+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 ! 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 ! ! scale dependency for local momentum transport ! zfacdx=0.2*hpbl(i)/zq(i,k+1) delxy=dx(i)*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 ! ! copies here to avoid duplicate input args for tridin ! 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 ! ! solve tridiagonal problem for momentum ! call tridi1n(al,ad,cu,r1,r2,au,f1,f2,its,ite,kts,kte,1) ! ! recover tendencies of momentum ! 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 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*rdt*(f1-ux) endif idtend=dtidx(index_of_y_wind,index_of_process_pbl) if(idtend>=1) then dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*rdt*(f2-vx) endif endif ! do i = its,ite kpbl1d(i) = kpbl(i) enddo ! !---- calculate sgs tke which is consistent with shinhongpbl algorithm ! 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 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 = dx(i)*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+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 ! !---- find the mixing length ! call mixlen(lmh,uxk,vxk,txk,thxk,qx(i,kts:kte,1),qx(i,kts:kte,ntcw) & ,q2xk,zqk,ust(i),corf,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 & ,kts,kte ) ! !---- solve for the production/dissipation of the turbulent kinetic energy ! call prodq2(lmh,dt,ust(i),s2,rig,q2xk,el,zqk,akmk,akhk & ,uxk,vxk,thxk,thvxk & ,hgamu(i),hgamv(i),hgamq(i),dx(i) & ,hpbl(i),pblflg(i),kpbl(i) & ,mfk,ufxpblk,vfxpblk,qfxpblk & ,ep1 & ,kts,kte ) ! ! !---- carry out the vertical diffusion of turbulent kinetic energy ! call vdifq(lmh,dt,q2xk,el,zqk & ,akhk,ptke1 & ,hgame,hpbl(i),pblflg(i),kpbl(i) & ,efxpbl(i) & ,kts,kte ) ! !---- save the new tke and mixing length. ! do k = kts,kte q2x(i,k) = amax1(q2xk(k),epsq2l) enddo ! enddo tke_calculation endif ! !---- end of tke calculation ! ! !---- end of vertical diffusion ! end subroutine shinhongvdif_run !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- subroutine tridi1n(cl,cm,cu,r1,r2,au,f1,f2,its,ite,kts,kte,nt) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- ! integer, intent(in ) :: its,ite, kts,kte, nt ! real(kind=kind_phys), dimension( its:ite, kts+1:kte+1 ) , & intent(in ) :: cl ! real(kind=kind_phys), dimension( its:ite, kts:kte ) , & intent(in ) :: cm, & r1 real(kind=kind_phys), dimension( its:ite, kts:kte,nt ) , & intent(in ) :: r2 ! real(kind=kind_phys), dimension( its:ite, kts:kte ) , & intent(inout) :: au, & cu, & f1 real(kind=kind_phys), dimension( its:ite, kts:kte,nt ) , & intent(inout) :: f2 ! real(kind=kind_phys) :: fk integer :: i,k,l,n,it ! !------------------------------------------------------------------------------- ! l = ite n = kte ! do i = its,l fk = 1./cm(i,1) au(i,1) = fk*cu(i,1) f1(i,1) = fk*r1(i,1) enddo ! do it = 1,nt do i = its,l fk = 1./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) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- ! integer, intent(in ) :: its,ite, kts,kte, nt ! real(kind=kind_phys), dimension( its:ite, kts+1:kte+1 ) , & intent(in ) :: cl ! real(kind=kind_phys), dimension( its:ite, kts:kte ) , & intent(in ) :: cm real(kind=kind_phys), dimension( its:ite, kts:kte,nt) , & intent(in ) :: r2 ! real(kind=kind_phys), dimension( its:ite, kts:kte ) , & intent(inout) :: au, & cu real(kind=kind_phys), dimension( its:ite, kts:kte,nt) , & intent(inout) :: f2 ! real(kind=kind_phys) :: fk integer :: i,k,l,n,it ! !------------------------------------------------------------------------------- ! l = ite n = kte ! do 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, & kts,kte) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- ! qnse model constants !------------------------------------------------------------------------------- real(kind=kind_phys),parameter :: blckdr=0.0063,cn=0.75 real(kind=kind_phys),parameter :: eps1=1.e-12,epsl=0.32,epsru=1.e-7,epsrs=1.e-7 real(kind=kind_phys),parameter :: el0max=1000.,el0min=1.,elfc=0.23*0.5 real(kind=kind_phys),parameter :: alph=0.30,beta=1./273.,g=9.81,btg=beta*g real(kind=kind_phys),parameter :: a1=0.659888514560862645,a2x=0.6574209922667784586 real(kind=kind_phys),parameter :: b1=11.87799326209552761,b2=7.226971804046074028 real(kind=kind_phys),parameter :: c1=0.000830955950095854396 real(kind=kind_phys),parameter :: adnh= 9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg real(kind=kind_phys),parameter :: adnm=18.*a1*a1*a2x*(b2-3.*a2x)*btg real(kind=kind_phys),parameter :: bdnh= 3.*a2x*(7.*a1+b2)*btg,bdnm= 6.*a1*a1 !------------------------------------------------------------------------------- ! free term in the equilibrium equation for (l/q)**2 !------------------------------------------------------------------------------- real(kind=kind_phys),parameter :: aeqh=9.*a1*a2x*a2x*b1*btg*btg & +9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg real(kind=kind_phys),parameter :: aeqm=3.*a1*a2x*b1*(3.*a2x+3.*b2*c1+18.*a1*c1-b2) & *btg+18.*a1*a1*a2x*(b2-3.*a2x)*btg !------------------------------------------------------------------------------- ! forbidden turbulence area !------------------------------------------------------------------------------- real(kind=kind_phys),parameter :: requ=-aeqh/aeqm real(kind=kind_phys),parameter :: epsgh=1.e-9,epsgm=requ*epsgh !------------------------------------------------------------------------------- ! near isotropy for shear turbulence, ww/q2 lower limit !------------------------------------------------------------------------------- real(kind=kind_phys),parameter :: ubryl=(18.*requ*a1*a1*a2x*b2*c1*btg & +9.*a1*a2x*a2x*b2*btg*btg) & /(requ*adnm+adnh) real(kind=kind_phys),parameter :: ubry=(1.+epsrs)*ubryl,ubry3=3.*ubry real(kind=kind_phys),parameter :: aubh=27.*a1*a2x*a2x*b2*btg*btg-adnh*ubry3 real(kind=kind_phys),parameter :: aubm=54.*a1*a1*a2x*b2*c1*btg -adnm*ubry3 real(kind=kind_phys),parameter :: bubh=(9.*a1*a2x+3.*a2x*b2)*btg-bdnh*ubry3 real(kind=kind_phys),parameter :: bubm=18.*a1*a1*c1 -bdnm*ubry3 real(kind=kind_phys),parameter :: cubr=1.-ubry3,rcubr=1./cubr !------------------------------------------------------------------------------- ! k profile constants !------------------------------------------------------------------------------- real(kind=kind_phys),parameter :: elcbl=0.77 !------------------------------------------------------------------------------- ! integer, intent(in ) :: kts,kte integer, intent(in ) :: lmh,lmxl,lpbl ! real(kind=kind_phys), intent(in ) :: p608,vkarman,cp real(kind=kind_phys), intent(in ) :: hpbl,corf,ustar,hgamu,hgamv,hgamq real(kind=kind_phys), intent(inout) :: ct,epshol ! real(kind=kind_phys), dimension( kts:kte ) , & intent(in ) :: cwm, & q, & q2, & t, & the, & u, & v ! real(kind=kind_phys), dimension( kts+1:kte ) , & intent(in ) :: mf, & ufxpbl, & vfxpbl, & qfxpbl ! real(kind=kind_phys), dimension( kts:kte+1 ) , & intent(in ) :: z ! real(kind=kind_phys), dimension( kts+1:kte ) , & intent(out ) :: el, & ri, & gh, & s2 ! logical,intent(in) :: pblflg ! ! local vars ! integer :: k,lpblm real(kind=kind_phys) :: suk,svk,elocp real(kind=kind_phys) :: 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(kind=kind_phys) :: ckp real(kind=kind_phys), dimension( kts:kte ) :: q1, & en2 real(kind=kind_phys), 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 ! ! compute local gradient richardson number ! 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 ! s**2 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 ! dtheta/dz 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 ! n**2 gh(k)=ghl ri(k)=en2(k)/s2l enddo ! ! find maximum mixing lengths and the level of the pbl top ! 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) ! repsr1/qol2un 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 ! ! computation of asymptotic l in blackadar formula ! el0=min(alph*szq*0.5/sq,el0max) el0=max(el0 ,el0min) ! ! above the pbl top ! 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 ! ! inside the pbl ! 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 ! ! mixing length for the qnse model in stable case ! f=max(corf,eps1) rlambda=f/(blckdr*ustar) do k = kte,kts+1,-1 if(en2(k)>=0.0)then ! stable case 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, & kts,kte) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- ! real(kind=kind_phys),parameter :: epsq2l = 0.01,c0 = 0.55,ceps = 16.6,g = 9.81 ! integer, intent(in ) :: kts,kte integer, intent(in ) :: lmh,kpbl ! real(kind=kind_phys), intent(in ) :: p608,dtturbl,ustar real(kind=kind_phys), intent(in ) :: hgamu,hgamv,hgamq,delxy,hpbl ! logical, intent(in ) :: pblflg ! real(kind=kind_phys), dimension( : ) , & intent(in ) :: uxk, & vxk, & thxk, & thvxk real(kind=kind_phys), dimension( : ) , & intent(in ) :: s2, & ri, & akm, & akh, & el, & mf, & ufxpbl, & vfxpbl, & qfxpbl ! real(kind=kind_phys), dimension( : ) , & intent(in ) :: z ! real(kind=kind_phys), dimension( : ) , & intent(inout) :: q2 ! ! local vars ! integer :: k ! real(kind=kind_phys) :: s2l,q2l,deltaz,akml,akhl,en2,pr,bpr,dis,rc02 real(kind=kind_phys) :: suk,svk,gthvk,govrthvk,pru,prv real(kind=kind_phys) :: thm,disel ! !------------------------------------------------------------------------------- ! rc02=2.0/(c0*c0) ! ! start of production/dissipation loop ! 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 !n**2 thm=(thxk(k)+thxk(k-1))*0.5 ! ! turbulence production term ! 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 ! ! buoyancy production ! 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 ! ! dissipation ! disel=min(delxy,ceps*el(k)) dis=(q2l)**1.5/disel ! q2l=q2l+2.0*(pr-bpr-dis)*dtturbl q2(k)=amax1(q2l,epsq2l) ! ! end of production/dissipation loop ! enddo main_integration ! ! lower boundary condition for q2 ! q2(kts)=amax1(rc02*ustar*ustar,epsq2l) ! end subroutine prodq2 !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- subroutine vdifq(lmh,dtdif,q2,el,z, & akhk,ptke1, & hgame,hpbl,pblflg,kpbl, & efxpbl, & kts,kte) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- ! real(kind=kind_phys),parameter :: c_k=1.0,esq=5.0 ! integer, intent(in ) :: kts,kte integer, intent(in ) :: lmh,kpbl ! real(kind=kind_phys), intent(in ) :: dtdif,hpbl,efxpbl ! logical, intent(in ) :: pblflg ! real(kind=kind_phys), dimension( : ) , & intent(in ) :: hgame, & ptke1 real(kind=kind_phys), dimension( : ) , & intent(in ) :: el, & akhk real(kind=kind_phys), dimension( : ) , & intent(in ) :: z ! real(kind=kind_phys), dimension( : ) , & intent(inout) :: q2 ! ! local vars ! integer :: k ! real(kind=kind_phys) :: aden,akqs,bden,besh,besm,cden,cf,dtozs,ell,eloq2,eloq4 real(kind=kind_phys) :: elqdz,esh,esm,esqhf,ghl,gml,q1l,rden,rdz real(kind=kind_phys) :: zak ! real(kind=kind_phys), dimension( kts+1:kte ) :: zfacentk real(kind=kind_phys), dimension( kts+2:kte ) :: akq, & cm, & cr, & dtoz, & rsq2 ! !------------------------------------------------------------------------------- ! ! vertical turbulent diffusion ! esqhf=0.5*esq do k = kts+1,kte zak=0.5*(z(k)+z(k-1)) !zak of vdifq = za(k-1) of shinhong2d 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 !------------------------------------------------------------------------------- function pu(d,h) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- real(kind=kind_phys) :: pu real(kind=kind_phys),parameter :: pmin = 0.0,pmax = 1.0 real(kind=kind_phys),parameter :: a1 = 1.0, a2 = 0.070, a3 = 1.0, a4 = 0.142, a5 = 0.071 real(kind=kind_phys),parameter :: b1 = 2.0, b2 = 0.6666667 real(kind=kind_phys) :: d,h,doh,num,den ! doh=d/h num=a1*(doh)**b1+a2*(doh)**b2 den=a3*(doh)**b1+a4*(doh)**b2+a5 pu=num/den pu=max(pu,pmin) pu=min(pu,pmax) ! return end function !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- function pq(d,h) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- real(kind=kind_phys) :: pq real(kind=kind_phys),parameter :: pmin = 0.0,pmax = 1.0 real(kind=kind_phys),parameter :: a1 = 1.0, a2 = -0.098, a3 = 1.0, a4 = 0.106, a5 = 0.5 real(kind=kind_phys),parameter :: b1 = 2.0 real(kind=kind_phys) :: d,h,doh,num,den ! doh=d/h num=a1*(doh)**b1+a2 den=a3*(doh)**b1+a4 pq=a5*num/den+(1.-a5) pq=max(pq,pmin) pq=min(pq,pmax) ! return end function !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- function pthnl(d,h) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- real(kind=kind_phys) :: pthnl real(kind=kind_phys),parameter :: pmin = 0.0,pmax = 1.0 real(kind=kind_phys),parameter :: a1 = 1.000, a2 = 0.936, a3 = -1.110, & a4 = 1.000, a5 = 0.312, a6 = 0.329, a7 = 0.243 real(kind=kind_phys),parameter :: b1 = 2.0, b2 = 0.875 real(kind=kind_phys) :: d,h,doh,num,den ! 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) pthnl=max(pthnl,pmin) pthnl=min(pthnl,pmax) ! return end function !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- function pthl(d,h) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- real(kind=kind_phys) :: pthl real(kind=kind_phys),parameter :: pmin = 0.0,pmax = 1.0 real(kind=kind_phys),parameter :: a1 = 1.000, a2 = 0.870, a3 = -0.913, & a4 = 1.000, a5 = 0.153, a6 = 0.278, a7 = 0.280 real(kind=kind_phys),parameter :: b1 = 2.0, b2 = 0.5 real(kind=kind_phys) :: d,h,doh,num,den ! 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) pthl=max(pthl,pmin) pthl=min(pthl,pmax) ! return end function !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- function ptke(d,h) !------------------------------------------------------------------------------- use machine , only : kind_phys implicit none !------------------------------------------------------------------------------- real(kind=kind_phys) :: ptke real(kind=kind_phys),parameter :: pmin = 0.0,pmax = 1.0 real(kind=kind_phys),parameter :: a1 = 1.000, a2 = 0.070, & a3 = 1.000, a4 = 0.142, a5 = 0.071 real(kind=kind_phys),parameter :: b1 = 2.0, b2 = 0.6666667 real(kind=kind_phys) :: d,h,doh,num,den ! doh=d/h num=a1*(doh)**b1+a2*(doh)**b2 den=a3*(doh)**b1+a4*(doh)**b2+a5 ptke=num/den ptke=max(ptke,pmin) ptke=min(ptke,pmax) ! return end function !------------------------------------------------------------------------------- end module shinhongvdif