diff --git a/wrf_open/var/src/python/wrf/var/__init__.py b/wrf_open/var/src/python/wrf/var/__init__.py index bf8fa90..b6b58bc 100755 --- a/wrf_open/var/src/python/wrf/var/__init__.py +++ b/wrf_open/var/src/python/wrf/var/__init__.py @@ -103,6 +103,8 @@ _FUNC_MAP = {"cape2d" : get_2dcape, "slp" : get_slp, "theta" : get_theta, "temp" : get_temp, + "tk" : get_tk, + "tc" : get_tc, "theta_e" : get_eth, "tv" : get_tv, "twb" : get_tw, @@ -117,7 +119,8 @@ _FUNC_MAP = {"cape2d" : get_2dcape, "wa" : get_w_destag, "lat" : get_lat, "lon" : get_lon, - "pressure" : get_pressure, + "pressure" : get_pressure_hpa, + "pres" : get_pressure, "wspddir" : get_destag_wspd_wdir, "wspddir_uvmet" : get_uvmet_wspd_wdir, "wspddir_uvmet10" : get_uvmet10_wspd_wdir, @@ -140,6 +143,8 @@ _VALID_KARGS = {"cape2d" : ["missing"], "rh2m" : [], "slp" : ["units"], "temp" : ["units"], + "tk" : [], + "tc" : [], "theta" : ["units"], "theta_e" : ["units"], "tv" : ["units"], @@ -155,6 +160,7 @@ _VALID_KARGS = {"cape2d" : ["missing"], "wa" : ["units"], "lat" : [], "lon" : [], + "pres" : ["units"], "pressure" : ["units"], "wspddir" : ["units"], "wspddir_uvmet" : ["units"], @@ -172,8 +178,7 @@ _ALIASES = {"cape_2d" : "cape2d", "latitude" : "lat", "longitude" : "lon", "omg" : "omega", - "pres" : "pressure", - "p" : "pressure", + "p" : "pres", "rh2" : "rh2m", "z": "height", "ter" : "terrain", diff --git a/wrf_open/var/src/python/wrf/var/decorators.py b/wrf_open/var/src/python/wrf/var/decorators.py index bb52605..ea3975d 100644 --- a/wrf_open/var/src/python/wrf/var/decorators.py +++ b/wrf_open/var/src/python/wrf/var/decorators.py @@ -26,10 +26,18 @@ def convert_units(unit_type, alg_unit): # If units are provided to the method call, use them. # Otherwise, need to parse the argspec to find what the default # value is since wraps does not preserve this. + argspec = getargspec(func) + try: + unit_idx = argspec.args.index("units") + except ValueError: + unit_idx = None + if ("units" in kargs): desired_units = kargs["units"] + elif (len(args) > unit_idx and unit_idx is not None): + desired_units = args[unit_idx] else: - argspec = getargspec(func) + arg_idx_from_right = (len(argspec.args) - argspec.args.index("units")) desired_units = argspec.defaults[-arg_idx_from_right] diff --git a/wrf_open/var/src/python/wrf/var/interp.py b/wrf_open/var/src/python/wrf/var/interp.py index 1a3cec2..3f5c63c 100755 --- a/wrf_open/var/src/python/wrf/var/interp.py +++ b/wrf_open/var/src/python/wrf/var/interp.py @@ -266,22 +266,30 @@ def interpline(data2d, pivot_point=None, def vinterp(wrfnc, field, vert_coord, interp_levels, extrapolate=False, field_type=None, log_p=False): - valid_coords = ("pressure", "pres", "ght_msl", - "ght_agl", "theta", "theta-e") - - valid_field_types = (None,"none", "pressure","pres","p","z", - "tc","tk", "theta","theta-e", "ght") - - icase_lookup = { None : 0, - "p" : 1, - "pres" : 1, - "pressure" : 1, - "z" : 2, - "ght" : 2, - "tc" : 3, - "tk" : 4, - "theta" : 5, - "theta-e" : 6} + # Remove case sensitivity + field_type = field_type.lower() if field_type is not None else "none" + vert_coord = vert_coord.lower() if vert_coord is not None else "none" + + valid_coords = ("pressure", "pres", "p", "ght_msl", + "ght_agl", "theta", "th", "theta-e", "thetae", "eth") + + valid_field_types = ("none", "pressure", "pres", "p", "z", + "tc", "tk", "theta", "th", "theta-e", "thetae", + "eth", "ght") + + icase_lookup = {"none" : 0, + "p" : 1, + "pres" : 1, + "pressure" : 1, + "z" : 2, + "ght" : 2, + "tc" : 3, + "tk" : 4, + "theta" : 5, + "th" : 5, + "theta-e" : 6, + "thetae" : 6, + "eth" : 6} # These constants match what's in the fortran code. rgas = 287.04 #J/K/kg @@ -312,7 +320,7 @@ def vinterp(wrfnc, field, vert_coord, interp_levels, extrapolate=False, if extrapolate: extrap = 1 - icase = icase_lookup[field_type.lower()] + icase = icase_lookup[field_type] # Extract vriables timeidx = -1 # Should this be an argument? @@ -325,7 +333,7 @@ def vinterp(wrfnc, field, vert_coord, interp_levels, extrapolate=False, terht = get_terrain(wrfnc, timeidx) t = get_theta(wrfnc, timeidx) tk = get_temp(wrfnc, timeidx, units="k") - p = get_pressure(wrfnc, timeidx) + p = get_pressure(wrfnc, timeidx, units="pa") ght = get_height(wrfnc, timeidx, msl=True) ht_agl = get_height(wrfnc, timeidx, msl=False) @@ -334,7 +342,7 @@ def vinterp(wrfnc, field, vert_coord, interp_levels, extrapolate=False, # Vertical coordinate type vcor = 0 - if vert_coord in ("pressure", "pres"): + if vert_coord in ("pressure", "pres", "p"): vcor = 1 vcord_array = p * ConversionFactors.PA_TO_HPA @@ -346,7 +354,7 @@ def vinterp(wrfnc, field, vert_coord, interp_levels, extrapolate=False, vcor = 3 vcord_array = n.exp(-ht_agl/sclht) - elif vert_coord == "theta": + elif vert_coord in ("theta", "th"): vcor = 4 idir = 1 icorsw = 0 @@ -356,7 +364,12 @@ def vinterp(wrfnc, field, vert_coord, interp_levels, extrapolate=False, vcord_array = monotonic(t,p_hpa,coriolis,idir,delta,icorsw) - elif vert_coord == "theta-e": + # We only extrapolate temperature fields below ground + # if we are interpolating to pressure or height vertical surfaces. + + icase = 0 + + elif vert_coord in ("theta-e", "thetae", "eth"): vcor = 5 icorsw = 0 idir = 1 @@ -378,8 +391,8 @@ def vinterp(wrfnc, field, vert_coord, interp_levels, extrapolate=False, missing = Constants.DEFAULT_FILL res = vintrp(field,p,tk,qv,ght,terht,sfp,smsfp, - vcord_array,interp_levels, - icase,extrap,vcor,log_p_int,missing) + vcord_array,interp_levels, + icase,extrap,vcor,log_p_int,missing) return ma.masked_values(res,missing) diff --git a/wrf_open/var/src/python/wrf/var/pressure.py b/wrf_open/var/src/python/wrf/var/pressure.py index 399f439..ff9e64a 100755 --- a/wrf_open/var/src/python/wrf/var/pressure.py +++ b/wrf_open/var/src/python/wrf/var/pressure.py @@ -1,12 +1,11 @@ -from wrf.var.constants import Constants from wrf.var.decorators import convert_units from wrf.var.util import extract_vars -__all__ = ["get_pressure"] +__all__ = ["get_pressure", "get_pressure_hpa"] @convert_units("pressure", "pa") -def get_pressure(wrfnc, timeidx=0, units="hpa"): +def get_pressure(wrfnc, timeidx=0, units="pa"): try: p_vars = extract_vars(wrfnc, timeidx, vars=("P", "PB")) @@ -24,5 +23,9 @@ def get_pressure(wrfnc, timeidx=0, units="hpa"): return pres +def get_pressure_hpa(wrfnc, timeidx=0, units="hpa"): + return get_pressure(wrfnc, timeidx, units=units) + + \ No newline at end of file diff --git a/wrf_open/var/src/python/wrf/var/psadlookup.py b/wrf_open/var/src/python/wrf/var/psadlookup.py index 1edceb4..730144f 100755 --- a/wrf_open/var/src/python/wrf/var/psadlookup.py +++ b/wrf_open/var/src/python/wrf/var/psadlookup.py @@ -4568,13 +4568,15 @@ _LOOKUP_TABLE = ( """) _LOOKUP_LIST = _LOOKUP_TABLE.split() -_PSADITHTE = n.array([float(x) for x in _LOOKUP_LIST[0:_THETA_DIM]],n.double) +_PSADITHTE = n.array([float(x) for x in _LOOKUP_LIST[0:_THETA_DIM]], + "float64") _PSADIPRS = n.array([float(x) - for x in _LOOKUP_LIST[_THETA_DIM:_THETA_DIM+_PRS_DIM]], n.double) + for x in _LOOKUP_LIST[_THETA_DIM:_THETA_DIM+_PRS_DIM]], + "float64") # Will keep C-indexing at the python level, be sure to transpose this _PSADITMK = n.array([float(x) for x in _LOOKUP_LIST[_THETA_DIM+_PRS_DIM:]], - n.double).reshape((_PRS_DIM,_THETA_DIM)) + "float64").reshape((_PRS_DIM,_THETA_DIM)) def get_lookup_tables(): return _PSADITHTE,_PSADIPRS,_PSADITMK diff --git a/wrf_open/var/src/python/wrf/var/temp.py b/wrf_open/var/src/python/wrf/var/temp.py index 70999a3..d94b3d3 100755 --- a/wrf_open/var/src/python/wrf/var/temp.py +++ b/wrf_open/var/src/python/wrf/var/temp.py @@ -4,7 +4,8 @@ from wrf.var.extension import computetk, computeeth, computetv, computewetbulb from wrf.var.decorators import convert_units from wrf.var.util import extract_vars -__all__ = ["get_theta", "get_temp", "get_eth", "get_tv", "get_tw"] +__all__ = ["get_theta", "get_temp", "get_eth", "get_tv", "get_tw", + "get_tk", "get_tc"] @convert_units("temp", "k") def get_theta(wrfnc, timeidx=0, units="k"): @@ -84,6 +85,12 @@ def get_tw(wrfnc, timeidx=0, units="k"): tw = computewetbulb(full_p,tk,qv) return tw + +def get_tk(wrfnc, timeidx=0): + return get_temp(wrfnc, timeidx, units="k") + +def get_tc(wrfnc, timeidx=0): + return get_temp(wrfnc, timeidx, units="c") diff --git a/wrf_open/var/src/python/wrf/var/units.py b/wrf_open/var/src/python/wrf/var/units.py index e22fbba..4399723 100755 --- a/wrf_open/var/src/python/wrf/var/units.py +++ b/wrf_open/var/src/python/wrf/var/units.py @@ -29,7 +29,7 @@ def _c_to_k(var): return var + Constants.TCK0 def _c_to_f(var): - return ((9.0/5.0)*var) + 32.0 + return 1.8 * var + 32.0 # Temperature is a more complicated operation so requres functions def _apply_temp_conv(var, var_unit, dest_unit): @@ -115,9 +115,10 @@ def check_units(unit, unit_type): def do_conversion(var, vartype, var_unit, dest_unit): if vartype != "temp": - return _apply_conv_fact(var, vartype, var_unit, dest_unit) + return _apply_conv_fact(var, vartype, + var_unit.lower(), dest_unit.lower()) else: - return _apply_temp_conv(var, var_unit, dest_unit) + return _apply_temp_conv(var, var_unit.lower(), dest_unit.lower()) diff --git a/wrf_open/var/src/python/wrf/var/wrfcape.f90 b/wrf_open/var/src/python/wrf/var/wrfcape.f90 index 73e748c..60567ec 100755 --- a/wrf_open/var/src/python/wrf/var/wrfcape.f90 +++ b/wrf_open/var/src/python/wrf/var/wrfcape.f90 @@ -248,6 +248,7 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,& klev = 0 klcl = 0 + ! the comments were taken from a mark stoelinga email, 23 apr 2007, ! in response to a user getting the "outside of lookup table bounds" ! error message. @@ -410,8 +411,10 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,& tvlift = tvirtual(tmklift,qvplift) ghtlift = zlcl ilcl = 1 + ELSE tmklift = tonpsadiabat(ethpari,prs(i,j,k),PSADITHTE,PSADIPRS,PSADITMK,GAMMA,throw_exception) + eslift = EZERO*exp(ESLCON1* (tmklift-CELKEL)/(tmklift-ESLCON2)) qvplift = EPS*eslift/ (prs(i,j,k)-eslift) tvenv = tvirtual(tmk(i,j,k),qvp(i,j,k)) @@ -421,6 +424,7 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,& ! buoyancy buoy(kk) = GRAV* (tvlift-tvenv)/tvenv zrel(kk) = ghtlift - ghtpari + IF ((kk.gt.1).and.(buoy(kk)*buoy(kk-1).lt.0.0d0)) THEN ! parcel ascent curve crosses sounding curve, so create a new level @@ -432,7 +436,9 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,& buoy(kk-1) = 0.d0 zrel(kk-1) = zrel(kk-2) +buoy(kk-2)/& (buoy(kk-2)-buoy(kk))*(zrel(kk)-zrel(kk-2)) + END IF + IF (ilcl.eq.1) THEN klcl = kk ilcl = 2 @@ -547,8 +553,8 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,& cin(i,j,mkzh-1) = zrel(klcl) + ghtpari - ter(i,j) ! meters agl cin(i,j,mkzh-2) = zrel(klfc) + ghtpari - ter(i,j) - ENDIF + ENDIF END DO END DO diff --git a/wrf_open/var/src/python/wrf/var/wrfext.f90 b/wrf_open/var/src/python/wrf/var/wrfext.f90 index 4837418..1964602 100755 --- a/wrf_open/var/src/python/wrf/var/wrfext.f90 +++ b/wrf_open/var/src/python/wrf/var/wrfext.f90 @@ -295,7 +295,7 @@ SUBROUTINE f_computetk(pressure,theta,tk,nx) REAL(KIND=8), DIMENSION(nx), INTENT(OUT) :: tk INTEGER :: i - REAL(KIND=8), PARAMETER :: P1000MB=100000.D0, R_D=287.06D0, CP=7.D0*R_D/2.D0 + REAL(KIND=8), PARAMETER :: P1000MB=100000.D0, R_D=287.D0, CP=7.D0*R_D/2.D0 DO i = 1,nx pi = (pressure(i)/P1000MB)**(R_D/CP) @@ -639,7 +639,6 @@ SUBROUTINE f_computeomega(qvp,tmk,www,prs,omg,mx,my,mz) REAL(KIND=8),INTENT(IN),DIMENSION(mx,my,mz) :: prs REAL(KIND=8),INTENT(OUT),DIMENSION(mx,my,mz) :: omg - ! Local variables INTEGER :: i, j, k REAL(KIND=8),PARAMETER :: GRAV=9.81, RGAS=287.04, EPS=0.622 @@ -648,8 +647,8 @@ SUBROUTINE f_computeomega(qvp,tmk,www,prs,omg,mx,my,mz) DO j=1,my DO i=1,mx omg(i,j,k)=-GRAV*prs(i,j,k) / & - (RGAS*((tmk(i,j,k)*(EPS+qvp(i,j,k))) / & - (EPS*(1.+qvp(i,j,k)))))*www(i,j,k) + (RGAS*((tmk(i,j,k)*(EPS+qvp(i,j,k))) / & + (EPS*(1.+qvp(i,j,k)))))*www(i,j,k) END DO END DO END DO @@ -2034,7 +2033,7 @@ FUNCTION f_intrpvalue (wvalp0,wvalp1,vlev,vcp0,vcp1,icase) rvalue = (vlev-vcp0)*(valp1-valp0)/(vcp1-vcp0)+valp0 IF (icase .EQ. 2) THEN !GHT - f_intrpvalue = -SCLHT*log(rvalue) + f_intrpvalue = -SCLHT*LOG(rvalue) ELSE f_intrpvalue = rvalue END IF @@ -2099,6 +2098,7 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& REAL(KIND=8), PARAMETER :: CPMD = 0.887d0 REAL(KIND=8), PARAMETER :: GAMMA = RGAS/CP REAL(KIND=8), PARAMETER :: GAMMAMD = RGASMD-CPMD + REAL(KIND=8), PARAMETER :: CELKEL = 273.16d0 ! Removes the warnings for uninitialized variables cvcord = '' @@ -2136,7 +2136,6 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& vlev = interp_levels(nreqlvs) END IF - DO j=1,ns DO i=1,ew ! Get the interpolated value that is within the model domain @@ -2147,6 +2146,7 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& vcp0 = vcarray(i,j,ripk) valp0 = datain(i,j,ripk) valp1 = datain(i,j,ripk-1) + IF ((vlev .GE. vcp0 .AND. vlev .LE. vcp1) .OR. & (vlev .LE. vcp0 .AND. vlev .GE. vcp1)) THEN ! print *,i,j,valp0,valp1 @@ -2167,6 +2167,7 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& tmpvlev = vlev END IF tempout(i,j) = f_intrpvalue(valp0,valp1,tmpvlev,vcp0,vcp1,icase) + ! print *,"one ",i,j,tempout(i,j) ifound = 1 END IF @@ -2189,7 +2190,6 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& CYCLE END IF - ! The grid point is either above or below the model domain ! First we will check to see if the grid point is above the ! model domain. @@ -2216,7 +2216,6 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& CYCLE END IF - ! If the field comming in is not a pressure,temperature or height ! field we can set the output array to the value at the lowest ! model level. @@ -2244,7 +2243,7 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& ezsurf = EXP(-zsurf/SCLHT) ! The if for checking above ground - if ((cvcord .EQ. 'z' .AND. vlev .LT. ezsurf) .OR. & + IF ((cvcord .EQ. 'z' .AND. vlev .LT. ezsurf) .OR. & (cvcord .EQ. 'p' .AND. vlev .LT. psurf)) THEN ! We are below the lowest data level but above the ground. @@ -2320,7 +2319,7 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& qvapor = MAX(qvp(i,j,1),1.e-15) gammam = GAMMA*(1. + GAMMAMD*qvapor) IF (icase .EQ. 3) THEN - tempout(i,j) = tlev - 273.16d0 + tempout(i,j) = tlev - CELKEL ELSE IF (icase .EQ. 4) THEN tempout(i,j) = tlev ! Potential temperature - theta diff --git a/wrf_open/var/test/ncl_get_var.ncl b/wrf_open/var/test/ncl_get_var.ncl index e5037aa..17a38e4 100644 --- a/wrf_open/var/test/ncl_get_var.ncl +++ b/wrf_open/var/test/ncl_get_var.ncl @@ -145,6 +145,7 @@ z_500 = wrf_user_intrp3d(z,p,"h",plev,0.,False) z_500_dims = dimsizes(z_500) + ; 2D interpolation along line t2 = wrf_user_getvar(input_file, "T2", 0) dimst2 = dimsizes(t2) @@ -177,6 +178,155 @@ fout->z_500 = (/z_500/) fout->t2_line = (/t2_line/) + + ; 3D interpolation to new vertical coordinates + ; interp t to theta + fld1 = wrf_user_getvar(input_file, "tk", -1) + vert_coord = "theta" + interp_levels = (/200,300,500,1000/) + + opts = True + opts@extrapolate = True + opts@field_type = "T" + opts@logP = True + + fld1_intrp = wrf_user_vert_interp(input_file,fld1,vert_coord,interp_levels,opts) + + fld1_dims = dimsizes(fld1_intrp) + + filedimdef(fout, (/"fld1_time", "fld1_levs", "fld1_sn", "fld1_we"/), \ + (/fld1_dims(0), fld1_dims(1), fld1_dims(2), fld1_dims(3)/), \ + (/False,False,False,False/)) + + filevardef(fout, "fld_tk_theta", typeof(fld1_intrp), (/"fld1_time", "fld1_levs", "fld1_sn", "fld1_we"/)) + filevarattdef(fout,"fld_tk_theta", fld1_intrp) + fout->fld_tk_theta = (/fld1_intrp/) + + ; interp t to theta-e + fld2 = fld1 + vert_coord := "theta-e" + fld2_intrp = wrf_user_vert_interp(input_file,fld2,vert_coord,interp_levels,opts) + + + fld2_dims = dimsizes(fld2_intrp) + + filedimdef(fout, (/"fld2_time", "fld2_levs", "fld2_sn", "fld2_we"/), \ + (/fld2_dims(0), fld2_dims(1), fld2_dims(2), fld2_dims(3)/), \ + (/False,False,False,False/)) + + filevardef(fout, "fld_tk_theta_e", typeof(fld2_intrp), (/"fld2_time", "fld2_levs", "fld2_sn", "fld2_we"/)) + filevarattdef(fout,"fld_tk_theta_e", fld2_intrp) + fout->fld_tk_theta_e = (/fld2_intrp/) + + + ; interp t to pressure + fld3 = fld1 + vert_coord := "pressure" + interp_levels := (/850,500/) + fld3_intrp = wrf_user_vert_interp(input_file,fld3,vert_coord,interp_levels,opts) + + + fld3_dims = dimsizes(fld3_intrp) + + filedimdef(fout, (/"fld3_time", "fld3_levs", "fld3_sn", "fld3_we"/), \ + (/fld3_dims(0), fld3_dims(1), fld3_dims(2), fld3_dims(3)/), \ + (/False,False,False,False/)) + + filevardef(fout, "fld_tk_pres", typeof(fld3_intrp), (/"fld3_time", "fld3_levs", "fld3_sn", "fld3_we"/)) + filevarattdef(fout,"fld_tk_pres", fld3_intrp) + fout->fld_tk_pres = (/fld3_intrp/) + + ; interp t to ght_msl + fld4 = fld1 + vert_coord := "ght_msl" + interp_levels := (/1,2/) + fld4_intrp = wrf_user_vert_interp(input_file,fld4,vert_coord,interp_levels,opts) + + + fld4_dims = dimsizes(fld4_intrp) + + filedimdef(fout, (/"fld4_time", "fld4_levs", "fld4_sn", "fld4_we"/), \ + (/fld4_dims(0), fld4_dims(1), fld4_dims(2), fld4_dims(3)/), \ + (/False,False,False,False/)) + + filevardef(fout, "fld_tk_ght_msl", typeof(fld4_intrp), (/"fld4_time", "fld4_levs", "fld4_sn", "fld4_we"/)) + filevarattdef(fout,"fld_tk_ght_msl", fld4_intrp) + fout->fld_tk_ght_msl = (/fld4_intrp/) + + + ; interp t to ght_agl + fld5 = fld1 + vert_coord := "ght_agl" + interp_levels := (/1,2/) + fld5_intrp = wrf_user_vert_interp(input_file,fld1,vert_coord,interp_levels,opts) + + + fld5_dims = dimsizes(fld5_intrp) + + filedimdef(fout, (/"fld5_time", "fld5_levs", "fld5_sn", "fld5_we"/), \ + (/fld5_dims(0), fld5_dims(1), fld5_dims(2), fld5_dims(3)/), \ + (/False,False,False,False/)) + + filevardef(fout, "fld_tk_ght_agl", typeof(fld5_intrp), (/"fld5_time", "fld5_levs", "fld5_sn", "fld5_we"/)) + filevarattdef(fout,"fld_tk_ght_agl", fld5_intrp) + fout->fld_tk_ght_agl = (/fld5_intrp/) + + ; interp ht to pres + fld6 = wrf_user_getvar(input_file, "height", -1) + vert_coord := "pressure" + opts@field_type = "ght" + interp_levels := (/500,50/) + fld6_intrp = wrf_user_vert_interp(input_file,fld6,vert_coord,interp_levels,opts) + + + fld6_dims = dimsizes(fld6_intrp) + + filedimdef(fout, (/"fld6_time", "fld6_levs", "fld6_sn", "fld6_we"/), \ + (/fld6_dims(0), fld6_dims(1), fld6_dims(2), fld6_dims(3)/), \ + (/False,False,False,False/)) + + filevardef(fout, "fld_ht_pres", typeof(fld6_intrp), (/"fld6_time", "fld6_levs", "fld6_sn", "fld6_we"/)) + filevarattdef(fout,"fld_ht_pres", fld6_intrp) + fout->fld_ht_pres = (/fld6_intrp/) + + + ; interp pres to theta + fld7 = wrf_user_getvar(input_file, "pressure", -1) + vert_coord := "theta" + opts@field_type = "pressure" + interp_levels := (/200,300,500,1000/) + fld7_intrp = wrf_user_vert_interp(input_file,fld7,vert_coord,interp_levels,opts) + + + fld7_dims = dimsizes(fld7_intrp) + + filedimdef(fout, (/"fld7_time", "fld7_levs", "fld7_sn", "fld7_we"/), \ + (/fld7_dims(0), fld7_dims(1), fld7_dims(2), fld7_dims(3)/), \ + (/False,False,False,False/)) + + filevardef(fout, "fld_pres_theta", typeof(fld7_intrp), (/"fld7_time", "fld7_levs", "fld7_sn", "fld7_we"/)) + filevarattdef(fout,"fld_pres_theta", fld7_intrp) + fout->fld_pres_theta = (/fld7_intrp/) + + ; interp theta-e to pressure + fld8 = wrf_user_getvar(input_file, "eth", -1) + vert_coord := "pressure" + opts@field_type = "T" + interp_levels := (/850,500,5/) + fld8_intrp = wrf_user_vert_interp(input_file,fld8,vert_coord,interp_levels,opts) + + + fld8_dims = dimsizes(fld8_intrp) + + filedimdef(fout, (/"fld8_time", "fld8_levs", "fld8_sn", "fld8_we"/), \ + (/fld8_dims(0), fld8_dims(1), fld8_dims(2), fld8_dims(3)/), \ + (/False,False,False,False/)) + + filevardef(fout, "fld_thetae_pres", typeof(fld8_intrp), (/"fld8_time", "fld8_levs", "fld8_sn", "fld8_we"/)) + filevarattdef(fout,"fld_thetae_pres", fld8_intrp) + fout->fld_thetae_pres = (/fld8_intrp/) + + delete(fout) diff --git a/wrf_open/var/test/utests.py b/wrf_open/var/test/utests.py index 091f47b..1905dc1 100644 --- a/wrf_open/var/test/utests.py +++ b/wrf_open/var/test/utests.py @@ -100,71 +100,28 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False): if (varname == "tc"): my_vals = getvar(in_wrfnc, "temp", units="c") tol = 0 - atol = .1 - nt.assert_allclose(my_vals, ref_vals, tol, atol) - elif (varname == "tk"): - my_vals = getvar(in_wrfnc, "temp", units="k") - tol = 0 - atol = .1 - nt.assert_allclose(my_vals, ref_vals, tol, atol) - elif (varname == "td"): - my_vals = getvar(in_wrfnc, "td", units="c") - tol = 0 - atol = .1 - nt.assert_allclose(my_vals, ref_vals, tol, atol) - elif (varname == "pressure"): - my_vals = getvar(in_wrfnc, varname, units="hpa") - tol = 1/100. - atol = 0 - nt.assert_allclose(my_vals, ref_vals, tol, atol) - elif (varname == "p"): - my_vals = getvar(in_wrfnc, varname, units="pa") - tol = 1/100. - atol = 0 + atol = .1 # Note: NCL uses 273.16 as conversion for some reason nt.assert_allclose(my_vals, ref_vals, tol, atol) - elif (varname == "slp"): - my_vals = getvar(in_wrfnc, varname, units="hpa") - tol = 1/100. - atol = 0 - nt.assert_allclose(my_vals, ref_vals, tol, atol) - - elif (varname == "uvmet"): - my_vals = getvar(in_wrfnc, varname) - tol = 0/100. - atol = 0.0001 - nt.assert_allclose(my_vals, ref_vals, tol, atol) - elif (varname == "uvmet10"): - my_vals = getvar(in_wrfnc, varname) - tol = 1/100. - atol = .5 - nt.assert_allclose(my_vals, ref_vals, tol, atol) - - elif (varname == "omg"): - my_vals = getvar(in_wrfnc, varname) - tol = 2/100. - atol = 0 - nt.assert_allclose(my_vals, ref_vals, tol, atol) - elif (varname == "ctt"): - my_vals = getvar(in_wrfnc, varname) - tol = 1/100. - atol = 0 + elif (varname == "pw"): + my_vals = getvar(in_wrfnc, "pw") + tol = .5/100.0 + atol = 0 # NCL uses different constants and doesn't use same + # handrolled virtual temp in method nt.assert_allclose(my_vals, ref_vals, tol, atol) elif (varname == "cape_2d"): cape_2d = getvar(in_wrfnc, varname) - tol = 1/100. - atol = 0 + tol = 0/100. # Not sure why different, F77 vs F90? + atol = .1 nt.assert_allclose(cape_2d, ref_vals, tol, atol) elif (varname == "cape_3d"): cape_3d = getvar(in_wrfnc, varname) - tol = 1/100. - atol = 0 - + tol = 0/100. # Not sure why different, F77 vs F90? + atol = .01 nt.assert_allclose(cape_3d, ref_vals, tol, atol) - else: my_vals = getvar(in_wrfnc, varname) - tol = 1/100. - atol = 0 + tol = 0/100. + atol = 0.0001 nt.assert_allclose(my_vals, ref_vals, tol, atol) @@ -239,7 +196,7 @@ def make_interp_test(varname, wrf_in, referent, multi=False, if (varname == "interplevel"): ref_ht_500 = _get_refvals(referent, "z_500", repeat, multi) hts = getvar(in_wrfnc, "z") - p = getvar(in_wrfnc, "pressure", units="hpa") + p = getvar(in_wrfnc, "pressure") hts_500 = interplevel(hts, p, 500) nt.assert_allclose(hts_500, ref_ht_500) @@ -249,7 +206,7 @@ def make_interp_test(varname, wrf_in, referent, multi=False, ref_p_cross = _get_refvals(referent, "p_cross", repeat, multi) hts = getvar(in_wrfnc, "z") - p = getvar(in_wrfnc, "pressure", units="hpa") + p = getvar(in_wrfnc, "pressure") pivot_point = (hts.shape[-2] / 2, hts.shape[-1] / 2) ht_cross = vertcross(hts, p, pivot_point=pivot_point,angle=90.) @@ -290,9 +247,13 @@ def make_interp_test(varname, wrf_in, referent, multi=False, nt.assert_allclose(t2_line1, t2_line2) elif (varname == "vinterp"): + # Tk to theta + fld_tk_theta = _get_refvals(referent, "fld_tk_theta", repeat, multi) + fld_tk_theta = n.squeeze(fld_tk_theta) + tk = getvar(in_wrfnc, "temp", units="k") - interp_levels = [200,1000,50] + interp_levels = [200,300,500,1000] field = vinterp(in_wrfnc, field=tk, @@ -302,10 +263,131 @@ def make_interp_test(varname, wrf_in, referent, multi=False, field_type="tk", log_p=True) - print field.shape - print field + tol = 0/100. + atol = 0.0001 + + field = n.squeeze(field) + nt.assert_allclose(field, fld_tk_theta, tol, atol) + + # Tk to theta-e + fld_tk_theta_e = _get_refvals(referent, "fld_tk_theta_e", repeat, multi) + fld_tk_theta_e = n.squeeze(fld_tk_theta_e) + + interp_levels = [200,300,500,1000] + field = vinterp(in_wrfnc, + field=tk, + vert_coord="theta-e", + interp_levels=interp_levels, + extrapolate=True, + field_type="tk", + log_p=True) + + tol = 0/100. + atol = 0.0001 + + field = n.squeeze(field) + nt.assert_allclose(field, fld_tk_theta_e, tol, atol) + + # Tk to pressure + fld_tk_pres = _get_refvals(referent, "fld_tk_pres", repeat, multi) + fld_tk_pres = n.squeeze(fld_tk_pres) + + interp_levels = [850,500] + + field = vinterp(in_wrfnc, + field=tk, + vert_coord="pressure", + interp_levels=interp_levels, + extrapolate=True, + field_type="tk", + log_p=True) + + field = n.squeeze(field) + nt.assert_allclose(field, fld_tk_pres, tol, atol) + + # Tk to geoht_msl + fld_tk_ght_msl = _get_refvals(referent, "fld_tk_ght_msl", repeat, multi) + fld_tk_ght_msl = n.squeeze(fld_tk_ght_msl) + interp_levels = [1,2] + + field = vinterp(in_wrfnc, + field=tk, + vert_coord="ght_msl", + interp_levels=interp_levels, + extrapolate=True, + field_type="tk", + log_p=True) + + field = n.squeeze(field) + nt.assert_allclose(field, fld_tk_ght_msl, tol, atol) + + # Tk to geoht_agl + fld_tk_ght_agl = _get_refvals(referent, "fld_tk_ght_agl", repeat, multi) + fld_tk_ght_agl = n.squeeze(fld_tk_ght_agl) + interp_levels = [1,2] + + field = vinterp(in_wrfnc, + field=tk, + vert_coord="ght_agl", + interp_levels=interp_levels, + extrapolate=True, + field_type="tk", + log_p=True) + + field = n.squeeze(field) + nt.assert_allclose(field, fld_tk_ght_agl, tol, atol) + + # Hgt to pressure + fld_ht_pres = _get_refvals(referent, "fld_ht_pres", repeat, multi) + fld_ht_pres = n.squeeze(fld_ht_pres) + + z = getvar(in_wrfnc, "height", units="m") + interp_levels = [500,50] + field = vinterp(in_wrfnc, + field=z, + vert_coord="pressure", + interp_levels=interp_levels, + extrapolate=True, + field_type="ght", + log_p=True) + + field = n.squeeze(field) + nt.assert_allclose(field, fld_ht_pres, tol, atol) + + # Pressure to theta + fld_pres_theta = _get_refvals(referent, "fld_pres_theta", repeat, multi) + fld_pres_theta = n.squeeze(fld_pres_theta) + + p = getvar(in_wrfnc, "pressure") + interp_levels = [200,300,500,1000] + field = vinterp(in_wrfnc, + field=p, + vert_coord="theta", + interp_levels=interp_levels, + extrapolate=True, + field_type="pressure", + log_p=True) + + field = n.squeeze(field) + nt.assert_allclose(field, fld_pres_theta, tol, atol) + + # Theta-e to pres + fld_thetae_pres = _get_refvals(referent, "fld_thetae_pres", repeat, multi) + fld_thetae_pres = n.squeeze(fld_thetae_pres) + + eth = getvar(in_wrfnc, "eth") + interp_levels = [850,500,5] + field = vinterp(in_wrfnc, + field=eth, + vert_coord="pressure", + interp_levels=interp_levels, + extrapolate=True, + field_type="theta-e", + log_p=True) + field = n.squeeze(field) + nt.assert_allclose(field, fld_thetae_pres, tol, atol) return test