Skip to content

km_opt=5 causes several divide by 0 #1304

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 1 commit into from
Dec 29, 2020
Merged

km_opt=5 causes several divide by 0 #1304

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
124 changes: 72 additions & 52 deletions dyn_em/module_diffusion_em.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5133,7 +5133,7 @@ SUBROUTINE nonlocal_flux (config_flags,nlflux,gamu,gamv,hpbl,kpbl, &
DO i = i_start, i_end
delxy = sqrt(dx/msftx(i,j)*dy/msfty(i,j))
pu1=pu(delxy,hpbl(i,j))
IF(sfcflg(i,j).and.sflux(i,j).GT.0.0)THEN
IF((sfcflg(i,j).and.sflux(i,j).GT.0.0) .and. (hpbl(i,j) .GT. 0)) THEN
!nonlocal momentum flux based on Brown and Grant (1997)
brint = -sm*ust(i,j)*ust(i,j)*wstar3(i,j)/(hpbl(i,j)*wscale(i,j)**4)
gamu(i,j) = pu1 * brint*u_phy(i,1,j)/wspd(i,j)
Expand Down Expand Up @@ -5165,8 +5165,12 @@ SUBROUTINE nonlocal_flux (config_flags,nlflux,gamu,gamv,hpbl,kpbl, &
deltaoh(i,j) = max(ezfac*deltaoh(i,j),hpbl(i,j)-za(i,kpbl(i,j)-1,j)-1.)
deltaoh(i,j) = min(deltaoh(i,j), hpbl(i,j))

rigs(i,j) = govrth(i,j)*dthv(i,j)*deltaoh(i,j)/(du**2.+dv**2.)
rigs(i,j) = max(min(rigs(i,j), rigsmax),rimin)
if ((du .ne. 0) .or. (dv .ne. 0)) then
rigs(i,j) = govrth(i,j)*dthv(i,j)*deltaoh(i,j)/(du**2.+dv**2.)
rigs(i,j) = max(min(rigs(i,j), rigsmax),rimin)
else
rigs(i,j) = rigsmax
endif
enlfrac2(i,j) = max(min(wm3/wstar3(i,j)/(1.+cpent/rigs(i,j)),entfmax),entfmin)
enlfrac2(i,j) = enlfrac2(i,j)*enlfrac
ENDIF
Expand All @@ -5182,44 +5186,41 @@ SUBROUTINE nonlocal_flux (config_flags,nlflux,gamu,gamv,hpbl,kpbl, &
ENDDO
ENDDO
ENDDO

DO j = j_start, j_end
DO i = i_start, i_end
deltaoh(i,j) = deltaoh(i,j)/hpbl(i,j)
ENDDO
ENDDO

DO j = j_start, j_end
DO i = i_start, i_end
delxy = sqrt(dx/msftx(i,j)*dy/msfty(i,j))
mlfrac = mltop-deltaoh(i,j)
ezfrac = mltop+deltaoh(i,j)
zfacmf(i,1,j) = min(max((zq(i,2,j)/hpbl(i,j)),zfmin),1.)
sfcfracn = max(sfcfracn1,zfacmf(i,1,j))
!
sflux0 = (a11+a12*sfcfracn)*sflux(i,j)
snlflux0 = nlfrac*sflux0
amf1 = snlflux0/sfcfracn
amf2 = -snlflux0/(mlfrac-sfcfracn)
bmf2 = -mlfrac*amf2
amf3 = snlflux0*enlfrac2(i,j)/deltaoh(i,j)
bmf3 = -amf3*mlfrac
hfxpbl(i,j) = amf3+bmf3
pth1 = pthnl(delxy,hpbl(i,j))
hfxpbl(i,j) = hfxpbl(i,j)*pth1
IF (pblflg(i,j)) THEN
deltaoh(i,j) = deltaoh(i,j)/hpbl(i,j)
delxy = sqrt(dx/msftx(i,j)*dy/msfty(i,j))
mlfrac = mltop-deltaoh(i,j)
ezfrac = mltop+deltaoh(i,j)
zfacmf(i,1,j) = min(max((zq(i,2,j)/hpbl(i,j)),zfmin),1.)
sfcfracn = max(sfcfracn1,zfacmf(i,1,j))
!
sflux0 = (a11+a12*sfcfracn)*sflux(i,j)
snlflux0 = nlfrac*sflux0
amf1 = snlflux0/sfcfracn
amf2 = -snlflux0/(mlfrac-sfcfracn)
bmf2 = -mlfrac*amf2
amf3 = snlflux0*enlfrac2(i,j)/deltaoh(i,j)
bmf3 = -amf3*mlfrac
hfxpbl(i,j) = amf3+bmf3
pth1 = pthnl(delxy,hpbl(i,j))
hfxpbl(i,j) = hfxpbl(i,j)*pth1

DO k = kts, ktf
zfacmf(i,k,j) = max((zq(i,k+1,j)/hpbl(i,j)),zfmin)
IF(pblflg(i,j).and.k.LT.kpbl(i,j)) THEN
IF(zfacmf(i,k,j).LE.sfcfracn) THEN
nlflux(i,k,j) = amf1*zfacmf(i,k,j)
ELSE IF (zfacmf(i,k,j).LE.mlfrac) THEN
nlflux(i,k,j) = amf2*zfacmf(i,k,j)+bmf2
ENDIF
nlflux(i,k,j) = nlflux(i,k,j) + hfxpbl(i,j)*exp(-entfacmf(i,k,j))
nlflux(i,k,j) = nlflux(i,k,j)*pth1
ENDIF
ENDDO
DO k = kts, ktf
zfacmf(i,k,j) = max((zq(i,k+1,j)/hpbl(i,j)),zfmin)
IF(k.LT.kpbl(i,j)) THEN
IF(zfacmf(i,k,j).LE.sfcfracn) THEN
nlflux(i,k,j) = amf1*zfacmf(i,k,j)
ELSE IF (zfacmf(i,k,j).LE.mlfrac) THEN
nlflux(i,k,j) = amf2*zfacmf(i,k,j)+bmf2
ENDIF
nlflux(i,k,j) = nlflux(i,k,j) + hfxpbl(i,j)*exp(-entfacmf(i,k,j))
nlflux(i,k,j) = nlflux(i,k,j)*pth1
ENDIF
ENDDO
ENDIF
ENDDO
ENDDO
END SUBROUTINE nonlocal_flux
Expand All @@ -5237,10 +5238,15 @@ FUNCTION pthnl(d,h)
REAL,PARAMETER :: b1 = 2.0, b2 = 0.875
real :: 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)
if (h .ne. 0) then
doh = d/h
num = a1*(doh)**b1 + a2*(doh)**b2+a3
den = a4*(doh)**b1 + a5*(doh)**b2+a6
pthnl = a7*num/den + (1. - a7)
else
pthnl = 1.
endif

pthnl = max(pthnl,pmin)
pthnl = min(pthnl,pmax)

Expand All @@ -5261,10 +5267,14 @@ FUNCTION pthl(d,h)
REAL,PARAMETER :: b1 = 2.0, b2 = 0.5
REAL :: 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)
if (h .ne. 0) then
doh = d/h
num = a1*(doh)**b1 + a2*(doh)**b2+a3
den = a4*(doh)**b1 + a5*(doh)**b2+a6
pthl = a7*num/den+(1. - a7)
else
pthl = 1.
endif
pthl = max(pthl,pmin)
pthl = min(pthl,pmax)

Expand All @@ -5284,10 +5294,14 @@ FUNCTION pu(d,h)
REAL,PARAMETER :: b1 = 2.0, b2 = 0.6666667
REAL :: 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
if (h .ne. 0) then
doh = d/h
num = a1*(doh)**b1 + a2*(doh)**b2
den = a3*(doh)**b1 + a4*(doh)**b2+a5
pu = num/den
else
pu = 1.
endif
pu = max(pu,pmin)
pu = min(pu,pmax)

Expand Down Expand Up @@ -8146,15 +8160,21 @@ SUBROUTINE vertical_diffusion_implicit(ru_tendf, rv_tendf, rw_tendf, rt_tendf,&
DO i = i_start, i_end
DO k = kts, ktf-1
rdz_w = - g/dnw(k)/(c1h(k)*mu(i,j)+c2h(k))
beta = 1.5*sqrt(tke(i,k,j))/l_diss(i,k,j)
IF (l_diss(i,k,j) .ne. 0) THEN
beta = 1.5*sqrt(tke(i,k,j))/l_diss(i,k,j)
ELSE
beta = 0.
ENDIF
a(k) = - 2.0*xkxavg(i,k,j)*dt*rdz_w*rdz(i,k,j)
b(k) = 1.0 + 2.0*dt*rdz_w*(rdz(i,k,j)*xkxavg(i,k,j) &
+ rdz(i,k+1,j)*xkxavg(i,k+1,j)) &
+ dt*beta
c(k) = - 2.0*xkxavg(i,k+1,j)*dt*rdz(i,k+1,j)*rdz_w

d(k) = tke(i,k,j) &
+ 0.5*dt*tke(i,k,j)**1.5/l_diss(i,k,j)
d(k) = tke(i,k,j)
IF (l_diss(i,k,j) .ne. 0) THEN
d(k) = d(k) + 0.5*dt*tke(i,k,j)**1.5/l_diss(i,k,j)
ENDIF
ENDDO

a(ktf) = 0. !-1
Expand Down