Added more OpenMP Directives. Fixed serious bug with cloud fraction and implemented new behavior allowing users to select the vertical coordinate type and select their own cloud thresholds. Fixes #25 .

fortran/calc_uh.f90

@@ -61,21 +61,26 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, &
 
     twodx = 2.0*dx
     twody = 2.0*dy
+    !$OMP PARALLEL DO COLLAPSE(3)
     DO k=2,nz-2
         DO j=2,ny-1
             DO i=2,nx-1
-                wavg = 0.5*(w(i,j,k)+w(i,j,k+1))
-                tem1(i,j,k) = wavg*((vs(i+1,j,k) - vs(i-1,j,k))/(twodx*mapfct(i,j)) - &
-                                    (us(i,j+1,k) - us(i,j-1,k))/(twody*mapfct(i,j)))
+                !wavg = 0.5*(w(i,j,k)+w(i,j,k+1))
+                tem1(i,j,k) = (0.5*(w(i,j,k)+w(i,j,k+1)))*((vs(i+1,j,k) - &
+                              vs(i-1,j,k))/(twodx*mapfct(i,j)) - &
+                              (us(i,j+1,k) - us(i,j-1,k))/(twody*mapfct(i,j)))
                 tem2(i,j,k) = 0.5*(zp(i,j,k) + zp(i,j,k+1))
             END DO
         END DO
     END DO
+    !$OMP END PARALLEL DO
 
     ! Integrate over depth uhminhgt to uhmxhgt AGL
     !
     !  WRITE(6,'(a,f12.1,a,f12.1,a)') &
     !        'Calculating UH from ',uhmnhgt,' to ',uhmxhgt,' m AGL'
+    !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i, j, k, zbot, ztop, kbot, ktop, &
+    !$OMP wgtlw, wbot, wtop, wsum, wmean, sum, helbot, heltop)
     DO j=2,ny-2
         DO i=2,nx-2
             zbot = zp(i,j,2) + uhmnhgt
@@ -142,6 +147,7 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, &
             END IF
         END DO
     END DO
+    !$OMP END PARALLEL DO
 
     uh = uh*1000.   ! Scale according to Kain et al. (2008)
 

fortran/eqthecalc.f90

@@ -32,6 +32,8 @@ SUBROUTINE DEQTHECALC(qvp, tmk, prs, eth, miy, mjx, mkzh)
     REAL(KIND=8) :: tlcl
     INTEGER :: i, j, k
 
+    ! Note:  removing temporaries did not improve performance for this algorithm.
+    !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i,j,k,q,t,p,e,tlcl)
     DO k = 1,mkzh
         DO j = 1,mjx
             DO i = 1,miy
@@ -40,11 +42,12 @@ SUBROUTINE DEQTHECALC(qvp, tmk, prs, eth, miy, mjx, mkzh)
                 p = prs(i,j,k)/100.
                 e = q*p/(EPS + q)
                 tlcl = TLCLC1/(LOG(t**TLCLC2/e) - TLCLC3) + TLCLC4
-                eth(i,j,k) = t*(1000.D0/p)**(GAMMA*(1.D0 + GAMMAMD*q))* &
+                eth(i,j,k) = tmk(i,j,k)*(1000.D0/p)**(GAMMA*(1.D0 + GAMMAMD*q))* &
                         EXP((THTECON1/tlcl - THTECON2)*q*(1.D0 + THTECON3*q))
             END DO
         END DO
     END DO
+    !$OMP END PARALLEL DO
 
     RETURN
 

fortran/rip_cape.f90

@@ -353,7 +353,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
 
     !CALL cpu_time(t1)
     !CALL OMP_SET_NUM_THREADS(16)
-!$OMP PARALLEL DO
+
+    !$OMP PARALLEL DO
     DO j = 1,mjy
        DO i = 1,mix
           DO k = 1,mkzh
@@ -364,7 +365,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
           END DO
        END DO
     END DO
-!$OMP END PARALLEL DO
+    !$OMP END PARALLEL DO
     
     CALL DPFCALC(prs_new, sfp, prsf, mix, mjy, mkzh, ter_follow)
 
@@ -377,12 +378,12 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
         RETURN
     END IF
 
-!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(tlcl, ethpari, &
-!$OMP zlcl, kk, ilcl, klcl, tmklift, tvenv, tvlift, ghtlift, & 
-!$OMP facden, tmkenv, qvpenv, eslift, qvplift, buoy, benamin, &
-!$OMP benaccum, zrel, kmax, dz, elfound, &
-!$OMP kel, klfc, &
-!$OMP i,j,k,kpar)
+    !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(tlcl, ethpari, &
+    !$OMP zlcl, kk, ilcl, klcl, tmklift, tvenv, tvlift, ghtlift, &
+    !$OMP facden, tmkenv, qvpenv, eslift, qvplift, buoy, benamin, &
+    !$OMP benaccum, zrel, kmax, dz, elfound, &
+    !$OMP kel, klfc, &
+    !$OMP i,j,k,kpar)
     DO j = 1,mjy
       DO i = 1,mix
           cape(i,j,1) = 0.d0
@@ -395,10 +396,10 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
               ! (note, qvppari and tmkpari already calculated above for 2d case.)
              
               tlcl = TLCLC1/(LOG(tmk_new(kpar,i,j)**TLCLC2/(MAX(1.D-20,qvp_new(kpar,i,j)*prs_new(kpar,i,j)/ &
-              (EPS + qvp_new(kpar,i,j))))) - TLCLC3) + TLCLC4
+                     (EPS + qvp_new(kpar,i,j))))) - TLCLC3) + TLCLC4
              
               ethpari = tmk_new(kpar,i,j)*(1000.D0/prs_new(kpar,i,j))**(GAMMA*(1.D0 + GAMMAMD*qvp_new(kpar,i,j)))* &
-                  EXP((THTECON1/tlcl - THTECON2)*qvp_new(kpar,i,j)*(1.D0 + THTECON3*qvp_new(kpar,i,j)))
+                        EXP((THTECON1/tlcl - THTECON2)*qvp_new(kpar,i,j)*(1.D0 + THTECON3*qvp_new(kpar,i,j)))
                 
               zlcl = ght_new(kpar,i,j) + (tmk_new(kpar,i,j) - tlcl)/(G/CP * (1.D0 + CPMD*qvp_new(kpar,i,j)))
 

fortran/wrf_cloud_fracf.f90

@@ -18,7 +18,11 @@ SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew)
     kchi = 0
     kcmi = 0
     kclo = 0
+    lowc = 0
+    midc = 0
+    highc = 0
 
+    !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i, j, k, kchi, kcmi, kclo)
     DO j = 1,ns
         DO i = 1,ew
             DO k = 1,nz-1
@@ -27,30 +31,124 @@ SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew)
                 IF ( pres(i,j,k) .GT. 45000. ) kchi=k
             END DO
 
-        DO k = 1,nz-1
-            IF (k .GE. kclo .AND. k .LT. kcmi) THEN
-                lowc(i,j) = MAX(rh(i,j,k), lowc(i,j))
-            ELSE IF (k .GE. kcmi .AND. k .LT. kchi) THEN ! mid cloud
-                midc(i,j) = MAX(rh(i,j,k), midc(i,j))
-            ELSE if (k .GE. kchi) THEN                  ! high cloud
-                highc(i,j) = MAX(rh(i,j,k), highc(i,j))
-            END IF
-        END DO
+            DO k = 1,nz-1
+                IF (k .GE. kclo .AND. k .LT. kcmi) THEN
+                    lowc(i,j) = MAX(rh(i,j,k), lowc(i,j))
+                ELSE IF (k .GE. kcmi .AND. k .LT. kchi) THEN ! mid cloud
+                    midc(i,j) = MAX(rh(i,j,k), midc(i,j))
+                ELSE if (k .GE. kchi) THEN                  ! high cloud
+                    highc(i,j) = MAX(rh(i,j,k), highc(i,j))
+                END IF
+            END DO
 
-        lowc(i,j)  = 4.0*lowc(i,j)/100. - 3.0
-        midc(i,j)  = 4.0*midc(i,j)/100. - 3.0
-        highc(i,j) = 2.5*highc(i,j)/100. - 1.5
+            lowc(i,j)  = 4.0*lowc(i,j)/100. - 3.0
+            midc(i,j)  = 4.0*midc(i,j)/100. - 3.0
+            highc(i,j) = 2.5*highc(i,j)/100. - 1.5
 
-        lowc(i,j)  = MIN(lowc(i,j), 1.0)
-        lowc(i,j)  = MAX(lowc(i,j), 0.0)
-        midc(i,j)  = MIN(midc(i,j), 1.0)
-        midc(i,j)  = MAX(midc(i,j), 0.0)
-        highc(i,j) = MIN(highc(i,j), 1.0)
-        highc(i,j) = MAX(highc(i,j), 0.0)
+            lowc(i,j)  = MIN(lowc(i,j), 1.0)
+            lowc(i,j)  = MAX(lowc(i,j), 0.0)
+            midc(i,j)  = MIN(midc(i,j), 1.0)
+            midc(i,j)  = MAX(midc(i,j), 0.0)
+            highc(i,j) = MIN(highc(i,j), 1.0)
+            highc(i,j) = MAX(highc(i,j), 0.0)
 
        END DO
     END DO
+    !$OMP END PARALLEL DO
 
     RETURN
 
 END SUBROUTINE DCLOUDFRAC
+
+
+! NCLFORTSTART
+SUBROUTINE DCLOUDFRAC2(vert, rh, vert_inc_w_height, low_thresh, mid_thresh, &
+                       high_thresh, msg, lowc, midc, highc, nz, ns, ew)
+    IMPLICIT NONE
+
+    !f2py threadsafe
+    !f2py intent(in,out) :: lowc, midc, highc
+
+    INTEGER  nz, ns, ew
+    REAL(KIND=8), DIMENSION(ew, ns, nz), INTENT(IN) :: rh, vert
+    REAL(KIND=8), INTENT(IN) :: low_thresh, mid_thresh, high_thresh, msg
+    INTEGER, INTENT(IN) :: vert_inc_w_height
+    REAL(KIND=8), DIMENSION(ew, ns), INTENT(OUT) :: lowc, midc, highc
+
+! NCLEND
+
+    INTEGER i, j, k, kstart, kend
+    INTEGER kchi, kcmi, kclo
+
+    ! Initialize the output
+    lowc = 0
+    midc = 0
+    highc = 0
+
+    !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i, j, k, kchi, kcmi, kclo)
+    DO j = 1,ns
+        DO i = 1,ew
+            ! A value of -1 means 'not found'.  This is needed to handle
+            ! the mountains, where the level thresholds are below the lowest
+            ! model level.
+            kchi = -1
+            kcmi = -1
+            kclo = -1
+
+            IF (vert_inc_w_height .NE. 0) THEN ! Vert coord increase with height
+                DO k = 1,nz
+                    IF (vert(i,j,k) .LT. low_thresh) kclo=k
+                    IF (vert(i,j,k) .LT. mid_thresh) kcmi=k
+                    IF (vert(i,j,k) .LT. high_thresh) kchi=k
+                END DO
+            ELSE ! Vert coord decrease with height
+                DO k = 1,nz
+                    IF (vert(i,j,k) .GT. low_thresh) kclo=k
+                    IF (vert(i,j,k) .GT. mid_thresh) kcmi=k
+                    IF (vert(i,j,k) .GT. high_thresh) kchi=k
+                END DO
+            ENDIF
+
+            DO k = 1,nz
+                IF (k .GE. kclo .AND. k .LT. kcmi) THEN
+                    lowc(i,j) = MAX(rh(i,j,k), lowc(i,j))
+                ELSE IF (k .GE. kcmi .AND. k .LT. kchi) THEN ! mid cloud
+                    midc(i,j) = MAX(rh(i,j,k), midc(i,j))
+                ELSE if (k .GE. kchi) THEN                  ! high cloud
+                    highc(i,j) = MAX(rh(i,j,k), highc(i,j))
+                END IF
+            END DO
+
+            ! Only do this when a cloud threshold is in the model vertical
+            ! domain, otherwise it will be set to missing
+            IF (kclo .GE. 1) THEN
+                lowc(i,j)  = 4.0*lowc(i,j)/100. - 3.0
+                lowc(i,j)  = MIN(lowc(i,j), 1.0)
+                lowc(i,j)  = MAX(lowc(i,j), 0.0)
+            ELSE
+                lowc(i,j) = msg
+            END IF
+
+            IF (kcmi .GE. 1) THEN
+                midc(i,j)  = 4.0*midc(i,j)/100. - 3.0
+                midc(i,j)  = MIN(midc(i,j), 1.0)
+                midc(i,j)  = MAX(midc(i,j), 0.0)
+            ELSE
+                midc(i,j) = msg
+            END IF
+
+            IF (kchi .GE. 1) THEN
+                highc(i,j) = 2.5*highc(i,j)/100. - 1.5
+                highc(i,j) = MIN(highc(i,j), 1.0)
+                highc(i,j) = MAX(highc(i,j), 0.0)
+            ELSE
+                highc(i,j) = msg
+            END IF
+
+       END DO
+    END DO
+    !$OMP END PARALLEL DO
+
+    RETURN
+
+END SUBROUTINE DCLOUDFRAC2

fortran/wrffortran.pyf

@@ -123,6 +123,22 @@ python module _wrffortran ! in
             integer, optional,check(shape(pres,1)==ns),depend(pres) :: ns=shape(pres,1)
             integer, optional,check(shape(pres,0)==ew),depend(pres) :: ew=shape(pres,0)
         end subroutine dcloudfrac
+        subroutine dcloudfrac2(vert,rh,vert_inc_w_height,low_thresh,mid_thresh,high_thresh,msg,lowc,midc,highc,nz,ns,ew) ! in :_wrffortran:wrf_cloud_fracf.f90
+            threadsafe 
+            real(kind=8) dimension(ew,ns,nz),intent(in) :: vert
+            real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: rh
+            integer intent(in) :: vert_inc_w_height
+            real(kind=8) intent(in) :: low_thresh
+            real(kind=8) intent(in) :: mid_thresh
+            real(kind=8) intent(in) :: high_thresh
+            real(kind=8) intent(in) :: msg
+            real(kind=8) dimension(ew,ns),intent(out,in),depend(ew,ns) :: lowc
+            real(kind=8) dimension(ew,ns),intent(out,in),depend(ew,ns) :: midc
+            real(kind=8) dimension(ew,ns),intent(out,in),depend(ew,ns) :: highc
+            integer, optional,check(shape(vert,2)==nz),depend(vert) :: nz=shape(vert,2)
+            integer, optional,check(shape(vert,1)==ns),depend(vert) :: ns=shape(vert,1)
+            integer, optional,check(shape(vert,0)==ew),depend(vert) :: ew=shape(vert,0)
+        end subroutine dcloudfrac2
         module wrf_constants ! in :_wrffortran:wrf_constants.f90
             real(kind=8), parameter,optional :: wrf_earth_radius=6370000.d0
             real(kind=8), parameter,optional :: rhowat=1000.d0
@@ -527,8 +543,8 @@ python module _wrffortran ! in
             threadsafe 
             real(kind=8) dimension(ew,ns,nz),intent(out,in) :: out
             real(kind=8) dimension(ew,ns,nz),intent(inout),depend(ew,ns,nz) :: in
-            real(kind=8) dimension(ew,ns,nz),depend(ew,ns,nz) :: lvprs
-            real(kind=8) dimension(ew,ns),depend(ew,ns) :: cor
+            real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: lvprs
+            real(kind=8) dimension(ew,ns),intent(in),depend(ew,ns) :: cor
             integer intent(in) :: idir
             real(kind=8) intent(in) :: delta
             integer, optional,intent(in),check(shape(out,0)==ew),depend(out) :: ew=shape(out,0)
@@ -536,7 +552,7 @@ python module _wrffortran ! in
             integer, optional,intent(in),check(shape(out,2)==nz),depend(out) :: nz=shape(out,2)
             integer intent(in) :: icorsw
         end subroutine wrf_monotonic
-        function wrf_intrp_value(wvalp0,wvalp1,vlev,vcp0,vcp1,icase,errstat,errmsg) ! in :_wrffortran:wrf_vinterp.f90
+        function wrf_intrp_value(wvalp0,wvalp1,vlev,vcp0,vcp1,icase,errstat) ! in :_wrffortran:wrf_vinterp.f90
             threadsafe 
             use wrf_constants, only: sclht,algerr
             real(kind=8) intent(in) :: wvalp0
@@ -546,7 +562,6 @@ python module _wrffortran ! in
             real(kind=8) intent(in) :: vcp1
             integer intent(in) :: icase
             integer intent(inout) :: errstat
-            character*(*) intent(inout) :: errmsg
             real(kind=8) :: wrf_intrp_value
         end function wrf_intrp_value
         subroutine wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,sfp,smsfp,vcarray,interp_levels,numlevels,icase,ew,ns,nz,extrap,vcor,logp,rmsg,errstat,errmsg) ! in :_wrffortran:wrf_vinterp.f90

src/wrf/cloudfrac.py

@@ -5,11 +5,16 @@ from .constants import Constants
 from .extension import _tk, _rh, _cloudfrac
 from .metadecorators import set_cloudfrac_metadata
 from .util import extract_vars
+from .geoht import _get_geoht
+
+import numpy.ma as ma
 
 @set_cloudfrac_metadata()
 def get_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True, 
-                 cache=None, meta=True, _key=None):
-    """Return the cloud fraction.
+                 cache=None, meta=True, _key=None,
+                 vert_type="pres", low_thresh=None, mid_thresh=None, 
+                 high_thresh=None, missing=Constants.DEFAULT_FILL):
+    """Return the cloud fraction for low, mid, and high level clouds.
     
     The leftmost dimension of the returned array represents three different 
     quantities:
@@ -18,6 +23,33 @@ def get_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
         - return_val[1,...] will contain MID level cloud fraction
         - return_val[2,...] will contain HIGH level cloud fraction
     
+    For backwards compatibility, the default vertical coordinate type is 
+    pressure, with default cloud levels defined as: 
+    
+        97000 Pa <= low_cloud < 80000 Pa
+        80000 Pa <= mid_cloud < 45000 Pa
+        45000 Pa <= high_cloud
+        
+    If the vertical coordinate type is 'height_agl' or 'height_msl', the 
+    default cloud levels are defined as:
+    
+        300 m <= low_cloud < 2000 m
+        2000 m <= mid_cloud < 6000 m
+        6000 m <= high_cloud
+        
+    Note that the default low cloud levels are chosen to
+    exclude clouds near the surface (fog). If you want fog included, set 
+    *low_thresh* to ~99500 Pa if *vert_type* is set to 'pres', or 15 m if using 
+    'height_msl' or 'height_agl'. Keep in mind that the lowest mass grid points 
+    are slightly above the ground, and in order to find clouds, the 
+    *low_thresh* needs to be set to values that are slightly greater than 
+    (less than) the lowest height (pressure) values.
+    
+    When using 'pres' or 'height_agl' for *vert_type*, there is a possibility 
+    that the lowest WRF level will be higher than the low_cloud or mid_cloud 
+    threshold, particularly for mountainous regions.  When this happens, a 
+    fill value will be used in the output.
+    
     This functions extracts the necessary variables from the NetCDF file 
     object in order to perform the calculation.
     
@@ -60,6 +92,26 @@ def get_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
         _key (:obj:`int`, optional): A caching key. This is used for internal 
             purposes only.  Default is None.
             
+        vert_type (:obj:`str`, optional):  The type of vertical coordinate used 
+            to determine cloud type thresholds.  Must be 'pres', 'height_msl',
+            or 'height_agl'.  For backwards compatibility, the default 
+            is 'pres'.
+            
+        low_thresh (:obj:`float`, optional): The lower bound for what is 
+            considered a low cloud.  If *vert_type* is 'pres', the default is 
+            97000 Pa.  If *vert_type* is 'height_agl' or 'height_msl', then the 
+            default is 300 m. 
+            
+        mid_thresh (:obj:`float`, optional): The lower bound for what is 
+            considered a mid level cloud.  If *vert_type* is 'pres', the 
+            default is 80000 Pa.  If *vert_type* is 'height_agl' or 
+            'height_msl', then the default is 2000 m. 
+            
+        high_thresh (:obj:`float`, optional): The lower bound for what is 
+            considered a high level cloud.  If *vert_type* is 'pres', the 
+            default is 45000 Pa.  If *vert_type* is 'height_agl' or 
+            'height_msl', then the default is 6000 m.  
+            
     Returns:
         :class:`xarray.DataArray` or :class:`numpy.ndarray`: The 
         cloud fraction array whose leftmost dimension is 3 (LOW=0, MID=1, 
@@ -85,4 +137,28 @@ def get_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True,
     tk = _tk(full_p, full_t)
     rh = _rh(qv, full_p, tk)
                 
-    return _cloudfrac(full_p, rh)
+    if vert_type.lower() == "pres" or vert_type.lower() == "pressure":
+        v_coord = full_p
+        _low_thresh = 97000. if low_thresh is None else low_thresh
+        _mid_thresh = 80000. if mid_thresh is None else mid_thresh
+        _high_thresh = 45000. if high_thresh is None else high_thresh
+        vert_inc_w_height = 0
+    elif (vert_type.lower() == "height_msl" 
+          or vert_type.lower() == "height_agl"):
+        is_msl = vert_type.lower() == "height_msl"
+        v_coord = _get_geoht(wrfin, timeidx, method, squeeze, 
+                             cache, meta=False, _key=_key, height=True, 
+                             msl=is_msl)
+        _low_thresh = 300. if low_thresh is None else low_thresh
+        _mid_thresh = 2000. if mid_thresh is None else mid_thresh
+        _high_thresh = 6000. if high_thresh is None else high_thresh
+        vert_inc_w_height = 1
+    else:
+        raise ValueError("'vert_type' must be 'pres', 'height_msl', "
+                         "or 'height_agl'")
+    
+    cfrac = _cloudfrac(v_coord, rh, vert_inc_w_height, 
+                      _low_thresh, _mid_thresh, _high_thresh, missing)
+    
+    return ma.masked_values(cfrac, missing)
+

src/wrf/computation.py

@@ -962,8 +962,9 @@ def cape_3d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow,
     return ma.masked_values(cape_cin, missing)
 
 
-@set_cloudfrac_alg_metadata(copyarg="pres")
-def cloudfrac(pres, relh, meta=True):
+@set_cloudfrac_alg_metadata(copyarg="vert")
+def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh,
+              high_thresh, missing=Constants.DEFAULT_FILL, meta=True):
     """Return the cloud fraction.
     
     The leftmost dimension of the returned array represents three different 
@@ -973,14 +974,41 @@ def cloudfrac(pres, relh, meta=True):
         - return_val[1,...] will contain MID level cloud fraction
         - return_val[2,...] will contain HIGH level cloud fraction
         
+    For backwards compatibility, the default vertical coordinate type is 
+    pressure, with default cloud levels defined as: 
+    
+        97000 Pa <= low_cloud < 80000 Pa
+        80000 Pa <= mid_cloud < 45000 Pa
+        45000 Pa <= high_cloud
+        
+    If the vertical coordinate type is 'height_agl' or 'height_msl', the 
+    default cloud levels are defined as:
+    
+        300 m <= low_cloud < 2000 m
+        2000 m <= mid_cloud < 6000 m
+        6000 m <= high_cloud
+        
+    Note that the default low cloud levels are chosen to
+    exclude clouds near the surface (fog). If you want fog included, set 
+    *low_thresh* to ~99500 Pa if *vert_type* is set to 'pres', or 15 m if using 
+    'height_msl' or 'height_agl'. Keep in mind that the lowest mass grid points 
+    are slightly above the ground, and in order to find clouds, the 
+    *low_thresh* needs to be set to values that are slightly greater than 
+    (less than) the lowest height (pressure) values.
+    
+    When using 'pres' or 'height_agl' for *vert_type*, there is a possibility 
+    that the lowest WRF level will be higher than the low_cloud or mid_cloud 
+    threshold, particularly for mountainous regions.  When this happens, a 
+    fill value will be used in the output.
+    
     This is the raw computational algorithm and does not extract any variables 
     from WRF output files.  Use :meth:`wrf.getvar` to both extract and compute
     diagnostic variables.
     
     Args:  
             
-        pres (:class:`xarray.DataArray` or :class:`numpy.ndarray`): Full 
-            pressure (perturbation + base state pressure) in [Pa], with the 
+        vert (:class:`xarray.DataArray` or :class:`numpy.ndarray`): The  
+            vertical coordinate variable (usually pressure or height) with the 
             rightmost dimensions as bottom_top x south_north x west_east
             
             Note:
@@ -990,8 +1018,21 @@ def cloudfrac(pres, relh, meta=True):
                 dimension names to the output.  Otherwise, default names will
                 be used.
             
-        relh (:class:`xarray.DataArray` or :class:`numpy.ndarray`) Relative 
-            humidity with the same dimensionality as *pres*
+        relh (:class:`xarray.DataArray` or :class:`numpy.ndarray`): Relative 
+            humidity with the same dimensionality as *vert*
+            
+        vert_inc_w_height (:obj:`int`): Set to 1 if the vertical coordinate 
+            values increase with height (height values). Set to 0 if the 
+            vertical coordinate values decrease with height (pressure values).
+            
+        low_thresh (:obj:`float`): The bottom vertical threshold for what is 
+            considered a low cloud.
+            
+        mid_thresh (:obj:`float`): The bottom vertical threshold for what is 
+            considered a mid level cloud.
+            
+        high_thresh (:obj:`float`): The bottom vertical threshold for what is 
+            considered a high cloud.
 
         meta (:obj:`bool`): Set to False to disable metadata and return 
             :class:`numpy.ndarray` instead of 
@@ -1016,7 +1057,10 @@ def cloudfrac(pres, relh, meta=True):
         :meth:`wrf.getvar`, :meth:`wrf.rh`
     
     """
-    return _cloudfrac(pres, relh)
+    cfrac = _cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh, 
+                      high_thresh, missing)
+    
+    return ma.masked_values(cfrac, missing)
 
 
 @set_alg_metadata(2, "pres_hpa", refvarndims=3, 

src/wrf/extension.py

@@ -7,7 +7,7 @@ from .constants import Constants
 
 from ._wrffortran import (dcomputetk, dinterp3dz, dinterp2dxy, dinterp1d,
                           dcomputeseaprs, dfilter2d, dcomputerh, dcomputeuvmet,
-                          dcomputetd, dcapecalc2d, dcapecalc3d, dcloudfrac, 
+                          dcomputetd, dcapecalc2d, dcapecalc3d, dcloudfrac2, 
                           wrfcttcalc, calcdbz, dcalrelhl, dcalcuh, dcomputepv, 
                           dcomputeabsvort, dlltoij, dijtoll, deqthecalc,
                           omgcalc, virtual_temp, wetbulbcalc, dcomputepw,
@@ -624,28 +624,34 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow,
 
 @check_args(0, 3, (3,3))
 @cloudfrac_left_iter()
-@cast_type(arg_idxs=(0, 1), outviews=("lowview", "medview", "highview"))
-@extract_and_transpose(outviews=("lowview", "medview", "hightview"))
-def _cloudfrac(p, rh, lowview=None, medview=None, highview=None):
-    """Wrapper for dcloudfrace.
+@cast_type(arg_idxs=(0, 1), outviews=("lowview", "midview", "highview"))
+@extract_and_transpose(outviews=("lowview", "midview", "highview"))
+def _cloudfrac(vert, rh, vert_inc_w_height, low_thresh, mid_thresh,
+               high_thresh, missing, lowview=None, midview=None, highview=None):
+    """Wrapper for dcloudfrac2.
     
     Located in wrf_cloud_fracf.f90.
     
     """
     if lowview is None:
-        lowview = np.zeros(p.shape[0:2], p.dtype, order="F")
+        lowview = np.zeros(vert.shape[0:2], vert.dtype, order="F")
         
-    if medview is None:
-        medview = np.zeros(p.shape[0:2], p.dtype, order="F")
+    if midview is None:
+        midview = np.zeros(vert.shape[0:2], vert.dtype, order="F")
         
     if highview is None:
-        highview = np.zeros(p.shape[0:2], p.dtype, order="F")
-    
-    result = dcloudfrac(p, 
-                        rh, 
-                        lowview, 
-                        medview, 
-                        highview)
+        highview = np.zeros(vert.shape[0:2], vert.dtype, order="F")
+    
+    result = dcloudfrac2(vert, 
+                         rh,
+                         vert_inc_w_height,
+                         low_thresh,
+                         mid_thresh,
+                         high_thresh,
+                         missing,
+                         lowview, 
+                         midview, 
+                         highview)
     
     return result
 

src/wrf/metadecorators.py

@@ -322,6 +322,7 @@ def set_wind_metadata(copy_varname, name, description,
         
     return func_wrapper
 
+
 def set_cape_metadata(is2d):
     """A decorator that sets the metadata for a wrapped CAPE function's output.
     
@@ -468,14 +469,22 @@ def set_cloudfrac_metadata():
             return wrapped(*args, **kwargs)
         
         argvars = from_args(wrapped, ("wrfin", "timeidx", "method", "squeeze", 
-                                      "cache", "_key"), 
-                          *args, **kwargs)
+                                      "cache", "_key", "vert_type", 
+                                      "low_thresh", "mid_thresh", 
+                                      "high_thresh", "missing"), 
+                            *args, **kwargs)
         wrfin = argvars["wrfin"]
         timeidx = argvars["timeidx"]
         method = argvars["method"]
         squeeze = argvars["squeeze"]
         cache = argvars["cache"]
         _key = argvars["_key"]
+        vert_type = argvars["vert_type"]
+        low_thresh = argvars["low_thresh"]
+        mid_thresh = argvars["mid_thresh"]
+        high_thresh = argvars["high_thresh"]
+        missing = argvars["missing"]
+        
         if cache is None:
             cache = {}
         
@@ -499,6 +508,20 @@ def set_cloudfrac_metadata():
         outattrs.update(copy_var.attrs)
         outdimnames = [None] * result.ndim
         
+        # For printing units
+        unitstr = ("Pa" if vert_type.lower() == "pres" 
+                   or vert_type.lower() == "pressure" else "m")
+        
+        # For setting the threholds in metdata
+        if vert_type.lower() == "pres" or vert_type.lower() == "pressure":
+            _low_thresh = 97000. if low_thresh is None else low_thresh
+            _mid_thresh = 80000. if mid_thresh is None else mid_thresh
+            _high_thresh = 45000. if high_thresh is None else high_thresh
+        elif vert_type.lower() == "height_msl" or "height_agl":
+            _low_thresh = 300. if low_thresh is None else low_thresh
+            _mid_thresh = 2000. if mid_thresh is None else mid_thresh
+            _high_thresh = 6000. if high_thresh is None else high_thresh
+        
         # Right dims
         outdimnames[-2:] = copy_var.dims[-2:]
         # Left dims
@@ -507,6 +530,11 @@ def set_cloudfrac_metadata():
         outattrs["description"] = "low, mid, high clouds"
         outattrs["MemoryOrder"] = "XY"
         outattrs["units"] = "%"
+        outattrs["low_thresh"] = "{} {}".format(_low_thresh, unitstr)
+        outattrs["mid_thresh"] = "{} {}".format(_mid_thresh, unitstr)
+        outattrs["high_thresh"] = "{} {}".format(_high_thresh, unitstr)
+        outattrs["_FillValue"] = missing
+        outattrs["missing_value"] = missing
         outname = "cloudfrac"
 
         # xarray doesn't line up coordinate dimensions based on 
@@ -529,6 +557,7 @@ def set_cloudfrac_metadata():
         
     return func_wrapper
 
+
 def set_latlon_metadata(xy=False):
     """A decorator that sets the metadata for a wrapped latlon function's 
     output.
@@ -614,6 +643,7 @@ def set_latlon_metadata(xy=False):
     
     return func_wrapper
 
+
 def set_height_metadata(geopt=False):
     """A decorator that sets the metadata for a wrapped height function's 
     output.
@@ -704,6 +734,7 @@ def set_height_metadata(geopt=False):
                          dims=outdimnames, coords=outcoords, attrs=outattrs)
     return func_wrapper
 
+
 def _set_horiz_meta(wrapped, instance, args, kwargs):  
     """A decorator implementation that sets the metadata for a wrapped 
     horizontal interpolation function.
@@ -799,6 +830,7 @@ def _set_horiz_meta(wrapped, instance, args, kwargs):
     return DataArray(result, name=outname, dims=outdimnames, 
                      coords=outcoords, attrs=outattrs)
 
+
 def _set_cross_meta(wrapped, instance, args, kwargs):
     """A decorator implementation that sets the metadata for a wrapped cross \
     section interpolation function.
@@ -1013,7 +1045,6 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
                      coords=outcoords, attrs=outattrs)  
     
 
-
 def _set_line_meta(wrapped, instance, args, kwargs):
     """A decorator implementation that sets the metadata for a wrapped line 
     interpolation function.
@@ -1770,6 +1801,7 @@ def set_alg_metadata(alg_ndims, refvarname,
     
     return func_wrapper
 
+
 def set_smooth_metdata():
     @wrapt.decorator
     def func_wrapper(wrapped, instance, args, kwargs): 
@@ -1994,14 +2026,14 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"):
     return func_wrapper
 
 
-def set_cloudfrac_alg_metadata(copyarg="pres"):
+def set_cloudfrac_alg_metadata(copyarg="vert"):
     """A decorator that sets the metadata for the wrapped raw cloud fraction 
     diagnostic function.
     
     Args:
             
         copyarg (:obj:`str`): The wrapped function argument to use for 
-            copying dimension names.  Default is 'pres'.
+            copying dimension names.  Default is 'vert'.
     
     Returns:
     
@@ -2025,6 +2057,16 @@ def set_cloudfrac_alg_metadata(copyarg="pres"):
             
         result = wrapped(*args, **kwargs)
         
+        argvals = from_args(wrapped, (copyarg, "low_thresh", 
+                                 "mid_thresh", "high_thresh", 
+                                 "missing"), 
+                       *args, **kwargs)
+        cp = argvals[copyarg]
+        low_thresh = argvals["low_thresh"]
+        mid_thresh = argvals["mid_thresh"]
+        high_thresh = argvals["high_thresh"]
+        missing = argvals["missing"]
+        
         # Default dimension names
         outdims = ["dim_{}".format(i) for i in py3range(result.ndim)]
         
@@ -2034,15 +2076,17 @@ def set_cloudfrac_alg_metadata(copyarg="pres"):
         outattrs["description"] = "low, mid, high clouds"
         outattrs["units"] = "%"
         outattrs["MemoryOrder"] = "XY"
+        outattrs["low_thresh"] = low_thresh
+        outattrs["mid_thresh"] = mid_thresh
+        outattrs["high_thresh"] = high_thresh
+        outattrs["_FillValue"] = missing
+        outattrs["missing_value"] = missing
             
-        
-        p = from_args(wrapped, copyarg, *args, **kwargs)[copyarg]
-            
-        if isinstance(p, DataArray):
+        if isinstance(cp, DataArray):
             # Right dims
-            outdims[-2:] = p.dims[-2:]
+            outdims[-2:] = cp.dims[-2:]
             # Left dims
-            outdims[1:-2] = p.dims[0:-3]
+            outdims[1:-2] = cp.dims[0:-3]
                 
         
         outcoords = {}     

src/wrf/routines.py

@@ -111,7 +111,8 @@ _VALID_KARGS = {"cape2d" : ["missing"],
              "uvmet_wspd_wdir" : ["units"],
              "uvmet10_wspd_wdir" : ["units"],  
              "ctt" : [],
-             "cloudfrac" : [],
+             "cloudfrac" : ["vert_type", "low_thresh", 
+                            "mid_thresh", "high_thresh"],
              "default" : []
             }
   

src/wrf/specialdec.py

@@ -383,6 +383,7 @@ def cape_left_iter(alg_dtype=np.float64):
     
     return func_wrapper
 
+
 def cloudfrac_left_iter(alg_dtype=np.float64):
     """A decorator to handle iterating over the leftmost dimensions for the 
     cloud fraction diagnostic.
@@ -409,46 +410,45 @@ def cloudfrac_left_iter(alg_dtype=np.float64):
     """
     @wrapt.decorator
     def func_wrapper(wrapped, instance, args, kwargs):
-        # The cape calculations use an ascending vertical pressure coordinate
         new_args = list(args)
         new_kwargs = dict(kwargs)
         
-        p = args[0]
+        vert = args[0]
         rh = args[1]
             
-        num_left_dims = p.ndim - 3
-        orig_dtype = p.dtype
+        num_left_dims = vert.ndim - 3
+        orig_dtype = vert.dtype
         
-        # No special left side iteration, build the output from the cape,cin
-        # result
+        # No special left side iteration, build the output from the 
+        # low, mid, high results.
         if (num_left_dims == 0):
-            low, med, high = wrapped(*new_args, **new_kwargs)
+            low, mid, high = wrapped(*new_args, **new_kwargs)
             
             output_dims = (3,)
-            output_dims += p.shape[-2:]
+            output_dims += vert.shape[-2:]
             output = np.empty(output_dims, orig_dtype)
             
             output[0,:] = low[:]
-            output[1,:] = med[:]
+            output[1,:] = mid[:]
             output[2,:] = high[:]
             
             return output
                 
 
-        # Initial output is ...,cape_cin,nz,ny,nx to create contiguous views
-        outdims = p.shape[0:num_left_dims]
+        # Initial output is ...,low_mid_high,nz,ny,nx to create contiguous views
+        outdims = vert.shape[0:num_left_dims]
         extra_dims = tuple(outdims) # Copy the left-most dims for iteration
         
         outdims += (3,) # low_mid_high
         
-        outdims += p.shape[-2:]
+        outdims += vert.shape[-2:]
         
         outview_array = np.empty(outdims, alg_dtype)
         
         # Create the output array where the leftmost dim is the cloud type
         output_dims = (3,)
         output_dims += extra_dims
-        output_dims += p.shape[-2:]
+        output_dims += vert.shape[-2:]
         output = np.empty(output_dims, orig_dtype)
         
         has_missing = False
@@ -456,25 +456,25 @@ def cloudfrac_left_iter(alg_dtype=np.float64):
         for left_idxs in iter_left_indexes(extra_dims):
             left_and_slice_idxs = left_idxs + (slice(None),)
             low_idxs = left_idxs + (0, slice(None))
-            med_idxs = left_idxs + (1, slice(None))
+            mid_idxs = left_idxs + (1, slice(None))
             high_idxs = left_idxs + (2, slice(None))
             
             low_output_idxs = (0,) + left_idxs + (slice(None),)
-            med_output_idxs = (1,) + left_idxs + (slice(None),)
+            mid_output_idxs = (1,) + left_idxs + (slice(None),)
             high_output_idxs = (2,) + left_idxs + (slice(None),)
             
-            new_args[0] = p[left_and_slice_idxs]
+            new_args[0] = vert[left_and_slice_idxs]
             new_args[1] = rh[left_and_slice_idxs]
             
             # Skip the possible empty/missing arrays for the join method
-            # Note: Masking handled by cape.py or computation.py, so only 
+            # Note: Masking handled by cloudfrac.py or computation.py, so only 
             # supply the fill values here.
             skip_missing = False
             for arg in (new_args[0:2]):
                 if isinstance(arg, np.ma.MaskedArray):
                     if arg.mask.all():
                         output[low_output_idxs] = missing
-                        output[med_output_idxs] = missing
+                        output[mid_output_idxs] = missing
                         output[high_output_idxs] = missing
                         
                         skip_missing = True
@@ -484,14 +484,14 @@ def cloudfrac_left_iter(alg_dtype=np.float64):
                 continue
             
             lowview = outview_array[low_idxs]
-            medview = outview_array[med_idxs]
+            midview = outview_array[mid_idxs]
             highview = outview_array[high_idxs]
             
             new_kwargs["lowview"] = lowview
-            new_kwargs["medview"] = medview
+            new_kwargs["midview"] = midview
             new_kwargs["highview"] = highview
             
-            low, med, high = wrapped(*new_args, **new_kwargs)
+            low, mid, high = wrapped(*new_args, **new_kwargs)
             
             # Make sure the result is the same data as what got passed in 
             # Can delete this once everything works
@@ -501,8 +501,8 @@ def cloudfrac_left_iter(alg_dtype=np.float64):
             
             output[low_output_idxs] = (
                             outview_array[low_idxs].astype(orig_dtype))
-            output[med_output_idxs] = (
-                            outview_array[med_idxs].astype(orig_dtype))
+            output[mid_output_idxs] = (
+                            outview_array[mid_idxs].astype(orig_dtype))
             output[high_output_idxs] = (
                             outview_array[high_idxs].astype(orig_dtype))
         

test/comp_utest.py

@@ -512,7 +512,7 @@ def get_args(varname, wrfnc, timeidx, method, squeeze):
         tkel = tk(full_p, full_t)
         relh = rh(qv, full_p, tkel)
         
-        return (full_p, relh)
+        return (full_p, relh, 0, 97000., 80000., 45000.)
         
         
 class WRFVarsTest(ut.TestCase):

test/listBug.ncl

@@ -1,6 +1,18 @@
+; Bug1: This segfaults
 l = NewList("fifo")
 name = "foo"
 ListAppend(l, (/name/))
 print(l)
 print(l[0])
 name = "bar"
+
+; Bug2 Variables disappear
+a = addfile("/Users/ladwig/Documents/wrf_files/wrfout_d02_2010-06-13_21:00:00.nc", "r")
+b := wrf_user_getvar(a, "slp", -1)
+c = NewList("fifo") 
+ListAppend(c, (/b/)) 
+b := wrf_user_getvar(a, "rh", -1) 
+ListAppend(c, (/b/))
+
+print(c[0])
+print(c[1]) ; Variables start disappearing

test/utests.py

@@ -170,7 +170,11 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
             tol = 2/100.
             atol = 0.1
             #print (np.amax(np.abs(to_np(my_vals) - ref_vals)))
-            nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol)
+            try:
+                nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol)
+            except:
+                print (np.amax(np.abs(to_np(my_vals) - ref_vals)))
+                raise
     
     
     return test