fixed divide by zero in Shin-Hong PBL

phys/module_bl_shinhong.F

@@ -840,9 +840,13 @@ subroutine shinhong2d(j,ux,vx,tx,qx,p2d,p2di,pi2d,                          &
        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.)
+       if (wstar(i) .ne. 0) then
+          uwst  = abs(ust(i)/wstar(i)-0.5)
+          uwstx = -80.*uwst+14.
+          csfac = 0.5*(tanh(uwstx)+3.)
+       else
+          csfac = 1
+       endif
        cslen(i) = csfac*hpbl(i)
      endif
    enddo
@@ -962,7 +966,11 @@ subroutine shinhong2d(j,ux,vx,tx,qx,p2d,p2di,pi2d,                          &
        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.)
+       if ((dux .ne. 0) .or. (dvx .ne. 0)) then
+         rigs(i) = govrth(i)*dthvx(i)*deltaoh(i)/(dux**2.+dvx**2.)
+       else
+         rigs(i) = rigsmax
+       endif
        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
@@ -1086,9 +1094,15 @@ subroutine shinhong2d(j,ux,vx,tx,qx,p2d,p2di,pi2d,                          &
      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)
+     if (pblflg(i)) then
+        amf2        = -snlflux0/(mlfrac-sfcfracn)
+        bmf2        = -mlfrac*amf2
+     endif
+     if ((deltaoh(i) .eq. 0) .and. (enlfrac2(i) .eq. 0)) then
+        amf3       = 0.
+     else
+        amf3       = snlflux0*enlfrac2(i)/deltaoh(i)
+     endif
      bmf3        = -amf3*mlfrac
      hfxpbl(i)   = amf3+bmf3
      pth1=pthnl(delxy,cslen(i))
@@ -2283,10 +2297,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)
 !
@@ -2305,10 +2323,14 @@ function pq(d,h)
    real,parameter :: b1 = 2.0
    real :: d,h,doh,num,den
 !
-   doh=d/h
-   num=a1*(doh)**b1+a2
-   den=a3*(doh)**b1+a4
-   pq=a5*num/den+(1.-a5)
+   if (h .ne. 0) then
+      doh=d/h
+      num=a1*(doh)**b1+a2
+      den=a3*(doh)**b1+a4
+      pq=a5*num/den+(1.-a5)
+   else
+      pq=1.
+   endif
    pq=max(pq,pmin)
    pq=min(pq,pmax)
 !
@@ -2328,10 +2350,14 @@ 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)
 !
@@ -2351,10 +2377,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)
 !
@@ -2374,10 +2404,14 @@ function ptke(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
-   ptke=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
+      ptke=num/den
+   else
+      ptke=1.
+   endif
    ptke=max(ptke,pmin)
    ptke=min(ptke,pmax)
 !