diff --git a/dyn_em/module_diffusion_em.F b/dyn_em/module_diffusion_em.F index 57fd621e7a..dc5c12c818 100644 --- a/dyn_em/module_diffusion_em.F +++ b/dyn_em/module_diffusion_em.F @@ -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) @@ -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 @@ -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 @@ -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) @@ -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) @@ -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) @@ -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