Skip to content

Syncing MYNN-EDMF + bug fix #1996

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
weiwangncar merged 10 commits into from
Feb 1, 2024
Merged

Syncing MYNN-EDMF + bug fix #1996

Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
e74b61c
Updates to MYNN-EDMF scheme
joeolson42 Nov 17, 2023
2d06ef6
patch for missing rsqblten array
joeolson42 Nov 18, 2023
a5b14a7
Remove an empty line in module_bl_mynn_common.F to trigger the reg test.
weiwangncar Dec 5, 2023
6b37be0
Merge branch 'wrf-model:develop' into develop
joeolson42 Dec 19, 2023
cfcc342
bug fix from Rob Fovell
joeolson42 Dec 19, 2023
ece28d1
test kind_phys as defined in module_gfs_machine.F -- expect all singl…
joeolson42 Jan 9, 2024
9f67c6d
Adding new routine to define kind_phys for CCPP-influenced modules.
joeolson42 Jan 10, 2024
647511f
Merge branch 'wrf-model:develop' into develop
joeolson42 Jan 23, 2024
60a824e
bug fixes and tweaks from the CCPP version and Laura Fowler
joeolson42 Jan 25, 2024
6dba638
Merge branch 'wrf-model:develop' into develop
joeolson42 Jan 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
166 changes: 77 additions & 89 deletions phys/module_bl_mynn.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -232,7 +232,7 @@
! bl_mynn_cloudpdf = 2 (Chab-Becht).
! Removed WRF_CHEM dependencies.
! Many miscellaneous tweaks.
! v4.5.2 / CCPP
! v4.6 / CCPP
! Some code optimization. Removed many conditions from loops. Redesigned the mass-
! flux scheme to use 8 plumes instead of a variable n plumes. This results in
! the removal of the output variable "nudprafts" and adds maxwidth and ztop_plume.
Expand All @@ -242,6 +242,7 @@
! Now outputs all SGS cloud mixing ratios as grid-mean values, not in-cloud. This
! results in a change in the pre-radiation code to no longer multiply mixing ratios
! by cloud fractions.
! Bug fix for the momentum transport.
! Lots of code cleanup: removal of test code, comments, changing text case, etc.
! Many misc tuning/tweaks.
!
Expand Down Expand Up @@ -1114,7 +1115,7 @@ SUBROUTINE mynn_bl_driver( &
!! Calculate the buoyancy production of TKE from cloud-top cooling when
!! \p bl_mynn_topdown =1.
if (bl_mynn_topdown.eq.1) then
call topdown_cloudrad(kts,kte,dz1,zw, &
call topdown_cloudrad(kts,kte,dz1,zw,fltv, &
&xland(i),kpbl(i),PBLH(i), &
&sqc,sqi,sqw,thl,th1,ex1,p1,rho1,thetav, &
&cldfra_bl1D,rthraten(i,:), &
Expand Down Expand Up @@ -2001,9 +2002,9 @@ SUBROUTINE mym_length ( &
uonset= 15.
wt_u = (1.0 - min(max(ugrid - uonset, 0.0)/30.0, 0.5))
cns = 2.7 !was 3.5
alp1 = 0.22
alp1 = 0.23
alp2 = 0.3
alp3 = 2.0 * wt_u !taper off bouyancy enhancement in shear-driven pbls
alp3 = 2.5 * wt_u !taper off bouyancy enhancement in shear-driven pbls
alp4 = 5.0
alp5 = 0.3
alp6 = 50.
Expand Down Expand Up @@ -2059,12 +2060,12 @@ SUBROUTINE mym_length ( &

! ** Length scale limited by the buoyancy effect **
IF ( dtv(k) .GT. 0.0 ) THEN
bv = max( sqrt( gtr*dtv(k) ), 0.001)
bv = max( sqrt( gtr*dtv(k) ), 0.0001)
elb = MAX(alp2*qkw(k), &
& alp6*edmf_a1(k-1)*edmf_w1(k-1)) / bv &
& *( 1.0 + alp3*SQRT( vsc/(bv*elt) ) )
elb = MIN(elb, zwk)
elf = 0.80 * qkw(k)/bv
elf = 1.0 * qkw(k)/bv
elBLavg(k) = MAX(elBLavg(k), alp6*edmf_a1(k-1)*edmf_w1(k-1)/bv)
ELSE
elb = 1.0e10
Expand All @@ -2084,8 +2085,10 @@ SUBROUTINE mym_length ( &
!add blending to use BouLac mixing length in free atmos;
!defined relative to the PBLH (zi) + transition layer (h1)
!el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf)
!try squared-blending
el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2)))
!try squared-blending - but take out elb (makes it underdiffusive)
!el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2)))
el(k) = sqrt( els**2/(1. + (els**2/elt**2)))
el(k) = min(el(k), elb)
el(k) = MIN (el(k), elf)
el(k) = el(k)*(1.-wt) + alp5*elBLavg(k)*wt

Expand Down Expand Up @@ -3633,13 +3636,13 @@ SUBROUTINE mym_condensation (kts,kte, &

real(kind_phys):: qsl,esat,qsat,dqsl,cld0,q1k,qlk,eq1,qll, &
&q2p,pt,rac,qt,t,xl,rsl,cpm,Fng,qww,alpha,beta,bb, &
&ls,wt,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, &
&ls,wt,wt2,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, &
&qmq,qsat_tk,q1_rh,rh_hack,dzm1,zsl,maxqc
real(kind_phys), parameter :: qpct_sfc=0.025
real(kind_phys), parameter :: qpct_pbl=0.030
real(kind_phys), parameter :: qpct_trp=0.040
real(kind_phys), parameter :: rhcrit =0.83 !for cloudpdf = 2
real(kind_phys), parameter :: rhmax =1.01 !for cloudpdf = 2
real(kind_phys), parameter :: rhmax =1.02 !for cloudpdf = 2
integer :: i,j,k

real(kind_phys):: erf
Expand Down Expand Up @@ -3864,25 +3867,18 @@ SUBROUTINE mym_condensation (kts,kte, &
!Add condition for falling/settling into low-RH layers, so at least
!some cloud fraction is applied for all qc, qs, and qi.
rh_hack= rh(k)
wt2 = min(max( zagl - pblh2, 0.0 )/300., 1.0)
!ensure adequate RH & q1 when qi is at least 1e-9 (above the PBLH)
if (qi(k)>1.e-9 .and. zagl .gt. pblh2) then
rh_hack =min(rhmax, rhcrit + 0.07*(9.0 + log10(qi(k))))
if ((qi(k)+qs(k))>1.e-9 .and. (zagl .gt. pblh2)) then
rh_hack =min(rhmax, rhcrit + wt2*0.045*(9.0 + log10(qi(k)+qs(k))))
rh(k) =max(rh(k), rh_hack)
!add rh-based q1
q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
q1(k) =max(q1_rh, q1(k) )
endif
!ensure adequate RH & q1 when qc is at least 1e-6
if (qc(k)>1.e-6) then
rh_hack =min(rhmax, rhcrit + 0.09*(6.0 + log10(qc(k))))
rh(k) =max(rh(k), rh_hack)
!add rh-based q1
q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
q1(k) =max(q1_rh, q1(k) )
endif
!ensure adequate RH & q1 when qs is at least 1e-8 (above the PBLH)
if (qs(k)>1.e-8 .and. zagl .gt. pblh2) then
rh_hack =min(rhmax, rhcrit + 0.07*(8.0 + log10(qs(k))))
!ensure adequate rh & q1 when qc is at least 1e-6 (above the PBLH)
if (qc(k)>1.e-6 .and. (zagl .gt. pblh2)) then
rh_hack =min(rhmax, rhcrit + wt2*0.08*(6.0 + log10(qc(k))))
rh(k) =max(rh(k), rh_hack)
!add rh-based q1
q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit)
Expand Down Expand Up @@ -3994,7 +3990,7 @@ SUBROUTINE mym_condensation (kts,kte, &
fac_damp = min(zagl * 0.0025, 1.0)
!cld_factor = 1.0 + fac_damp*MAX(0.0, ( RH(k) - 0.75 ) / 0.26 )**1.9 !HRRRv4
!cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.25 )**2, 0.3)
cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.35)
cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.37)
cldfra_bl1D(K) = min( 1., cld_factor*cldfra_bl1D(K) )
enddo

Expand Down Expand Up @@ -4182,38 +4178,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &

k=kts

!original approach (drag in b-vector):
! a(1)=0.
! b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff
! c(1)=-dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff
! d(1)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff + &
! sub_u(k)*delt + det_u(k)*delt

!rho-weighted (drag in b-vector):
a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) &
& - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
& - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff - &
& dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt

!rho-weighted with drag term moved out of b-array
! a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
! b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
! c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
! d(k)=u(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*uoce*ust**2/wspd - &
! !!!d(k)=u(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*uoce*ust**2/wspd - &
! & dtz(k)*rhoinv(k)*s_awu(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt

DO k=kts+1,kte-1
a(k)= -dtz(k)*kmdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + &
& 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff + &
& sub_u(k)*delt + det_u(k)*delt
ENDDO
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=u(k) + dtz(k)*uoce*ust**2/wspd &
& - dtz(k)*rhoinv(k)*s_awu(k+1)*onoff &
& + dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff &
& + sub_u(k)*delt + det_u(k)*delt

do k=kts+1,kte-1
a(k)= -dtz(k)*kmdz(k)*rhoinv(k) &
& + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff &
& + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
b(k)=1.+ dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) &
& + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff &
& + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
c(k)= - dtz(k)*kmdz(k+1)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff &
& - dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff &
& + sub_u(k)*delt + det_u(k)*delt
enddo

!! no flux at the top
! a(kte)=-1.
Expand Down Expand Up @@ -4248,37 +4239,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i, &

k=kts

!original approach (drag in b-vector):
! a(1)=0.
! b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff
! c(1)= - dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff
! d(1)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff + &
! sub_v(k)*delt + det_v(k)*delt

!rho-weighted (drag in b-vector):
a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) &
& - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + &
& sub_v(k)*delt + det_v(k)*delt

!rho-weighted with drag term moved out of b-array
! a(k)= -dtz(k)*kmdz(k)*rhoinv(k)
! b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
! c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
! d(k)=v(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*voce*ust**2/wspd - &
! !!!d(k)=v(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*voce*ust**2/wspd - &
! & dtz(k)*rhoinv(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + sub_v(k)*delt + det_v(k)*delt

DO k=kts+1,kte-1
a(k)= -dtz(k)*kmdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + &
& 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff + &
& sub_v(k)*delt + det_v(k)*delt
ENDDO
b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=v(k) + dtz(k)*voce*ust**2/wspd &
& - dtz(k)*rhoinv(k)*s_awv(k+1)*onoff &
& + dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff &
& + sub_v(k)*delt + det_v(k)*delt

do k=kts+1,kte-1
a(k)= -dtz(k)*kmdz(k)*rhoinv(k) &
& + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff &
& + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff
b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) &
& + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff &
& + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff
c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) &
& - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff &
& - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff
d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff &
& - dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff &
& + sub_v(k)*delt + det_v(k)*delt
enddo

!! no flux at the top
! a(kte)=-1.
Expand Down Expand Up @@ -7637,7 +7624,8 @@ FUNCTION phih(zet)

END FUNCTION phih
! ==================================================================
SUBROUTINE topdown_cloudrad(kts,kte,dz1,zw,xland,kpbl,PBLH, &
SUBROUTINE topdown_cloudrad(kts,kte, &
&dz1,zw,fltv,xland,kpbl,PBLH, &
&sqc,sqi,sqw,thl,th1,ex1,p1,rho1,thetav, &
&cldfra_bl1D,rthraten, &
&maxKHtopdown,KHtopdown,TKEprodTD )
Expand All @@ -7648,15 +7636,15 @@ SUBROUTINE topdown_cloudrad(kts,kte,dz1,zw,xland,kpbl,PBLH, &
thl,th1,ex1,p1,rho1,thetav,cldfra_bl1D
real(kind_phys), dimension(kts:kte), intent(in) :: rthraten
real(kind_phys), dimension(kts:kte+1), intent(in) :: zw
real(kind_phys), intent(in) :: pblh
real(kind_phys), intent(in) :: pblh,fltv
real(kind_phys), intent(in) :: xland
integer , intent(in) :: kpbl
!output
real(kind_phys), intent(out) :: maxKHtopdown
real(kind_phys), dimension(kts:kte), intent(out) :: KHtopdown,TKEprodTD
!local
real(kind_phys), dimension(kts:kte) :: zfac,wscalek2,zfacent
real(kind_phys) :: bfx0,sflux,wm2,wm3,h1,h2,bfxpbl,dthvx,tmp1
real(kind_phys) :: bfx0,wm2,wm3,bfxpbl,dthvx,tmp1
real(kind_phys) :: temps,templ,zl1,wstar3_2
real(kind_phys) :: ent_eff,radsum,radflux,we,rcldb,rvls,minrad,zminrad
real(kind_phys), parameter :: pfac =2.0, zfmin = 0.01, phifac=8.0
Expand Down Expand Up @@ -7720,23 +7708,23 @@ SUBROUTINE topdown_cloudrad(kts,kte,dz1,zw,xland,kpbl,PBLH, &
bfx0 = max(radsum/rho1(k)/cp,0.)
else ! LAND
radsum=MIN(0.25*radsum,30.0)!practically turn off over land
bfx0 = max(radsum/rho1(k)/cp - max(sflux,0.0),0.)
bfx0 = max(radsum/rho1(k)/cp - max(fltv,0.0),0.)
endif

!entrainment from PBL top thermals
wm3 = grav/thetav(k)*bfx0*MIN(pblh,1500.) ! this is wstar3(i)
wm2 = wm2 + wm3**h2
wm2 = wm2 + wm3**twothirds
bfxpbl = - ent_eff * bfx0
dthvx = max(thetav(k+1)-thetav(k),0.1)
we = max(bfxpbl/dthvx,-sqrt(wm3**h2))
we = max(bfxpbl/dthvx,-sqrt(wm3**twothirds))

DO kk = kts,kpbl+3
!Analytic vertical profile
zfac(kk) = min(max((1.-(zw(kk+1)-zl1)/(zminrad-zl1)),zfmin),1.)
zfacent(kk) = 10.*MAX((zminrad-zw(kk+1))/zminrad,0.0)*(1.-zfac(kk))**3

!Calculate an eddy diffusivity profile (not used at the moment)
wscalek2(kk) = (phifac*karman*wm3*(zfac(kk)))**h1
wscalek2(kk) = (phifac*karman*wm3*(zfac(kk)))**onethird
!Modify shape of Kh to be similar to Lock et al (2000): use pfac = 3.0
KHtopdown(kk) = wscalek2(kk)*karman*(zminrad-zw(kk+1))*(1.-zfac(kk))**3 !pfac
KHtopdown(kk) = MAX(KHtopdown(kk),0.0)
Expand Down