From 0c270c7e028b5475d8ec0bea901fb3e9867371c1 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Mon, 30 Oct 2017 16:52:05 -0600 Subject: [PATCH] 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 | 12 +++- fortran/eqthecalc.f90 | 5 +- fortran/rip_cape.f90 | 21 +++--- fortran/wrf_cloud_fracf.f90 | 134 +++++++++++++++++++++++++++++++----- fortran/wrffortran.pyf | 23 +++++-- src/wrf/cloudfrac.py | 84 ++++++++++++++++++++-- src/wrf/computation.py | 58 ++++++++++++++-- src/wrf/extension.py | 36 ++++++---- src/wrf/metadecorators.py | 70 +++++++++++++++---- src/wrf/routines.py | 3 +- src/wrf/specialdec.py | 46 ++++++------- test/comp_utest.py | 4 +- test/listBug.ncl | 12 ++++ test/utests.py | 6 +- 14 files changed, 412 insertions(+), 102 deletions(-) diff --git a/fortran/calc_uh.f90 b/fortran/calc_uh.f90 index a893189..b99a48e 100644 --- a/fortran/calc_uh.f90 +++ b/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) diff --git a/fortran/eqthecalc.f90 b/fortran/eqthecalc.f90 index a4cb6ba..d805be9 100644 --- a/fortran/eqthecalc.f90 +++ b/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 diff --git a/fortran/rip_cape.f90 b/fortran/rip_cape.f90 index 51c3916..c4429b4 100644 --- a/fortran/rip_cape.f90 +++ b/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))) diff --git a/fortran/wrf_cloud_fracf.f90 b/fortran/wrf_cloud_fracf.f90 index 61cbe4c..c925704 100644 --- a/fortran/wrf_cloud_fracf.f90 +++ b/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 diff --git a/fortran/wrffortran.pyf b/fortran/wrffortran.pyf index 5a01ba5..41eecbb 100644 --- a/fortran/wrffortran.pyf +++ b/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 diff --git a/src/wrf/cloudfrac.py b/src/wrf/cloudfrac.py index 75b500f..38805c6 100644 --- a/src/wrf/cloudfrac.py +++ b/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. @@ -59,7 +91,27 @@ 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, @@ -84,5 +136,29 @@ def get_cloudfrac(wrfin, timeidx=0, method="cat", squeeze=True, tk = _tk(full_p, full_t) rh = _rh(qv, full_p, tk) + + 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 _cloudfrac(full_p, rh) + return ma.masked_values(cfrac, missing) + diff --git a/src/wrf/computation.py b/src/wrf/computation.py index 81e37dd..1a0e1c3 100644 --- a/src/wrf/computation.py +++ b/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 @@ -972,6 +973,33 @@ def cloudfrac(pres, relh, meta=True): - return_val[0,...] will contain LOW level cloud fraction - 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 @@ -979,8 +1007,8 @@ def cloudfrac(pres, relh, meta=True): 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, diff --git a/src/wrf/extension.py b/src/wrf/extension.py index 90bb8cc..777e29d 100755 --- a/src/wrf/extension.py +++ b/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 diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index e34b71d..0fe6b35 100644 --- a/src/wrf/metadecorators.py +++ b/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. @@ -613,7 +642,8 @@ def set_latlon_metadata(xy=False): return da 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. @@ -798,7 +829,8 @@ 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. @@ -1012,7 +1044,6 @@ def _set_cross_meta(wrapped, instance, args, kwargs): return DataArray(result, name=outname, dims=outdimnames, coords=outcoords, attrs=outattrs) - def _set_line_meta(wrapped, instance, args, kwargs): """A decorator implementation that sets the metadata for a wrapped line @@ -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 = {} diff --git a/src/wrf/routines.py b/src/wrf/routines.py index b72423b..4aed3d8 100644 --- a/src/wrf/routines.py +++ b/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" : [] } diff --git a/src/wrf/specialdec.py b/src/wrf/specialdec.py index c093ca5..8d5efc7 100644 --- a/src/wrf/specialdec.py +++ b/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)) diff --git a/test/comp_utest.py b/test/comp_utest.py index 846466f..6af64d2 100644 --- a/test/comp_utest.py +++ b/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): @@ -538,7 +538,7 @@ def make_func(varname, wrfnc, timeidx, method, squeeze, meta): if meta: self.assertEqual(result.dims, ref.dims) - + return func diff --git a/test/listBug.ncl b/test/listBug.ncl index 49b81f0..8e9745a 100644 --- a/test/listBug.ncl +++ b/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 diff --git a/test/utests.py b/test/utests.py index 234239a..b7b26f1 100644 --- a/test/utests.py +++ b/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