From ad66cd6c0751f871ef5dfc475a278c1179042a9f Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Thu, 14 Jan 2016 09:16:37 -0700 Subject: [PATCH] Added vertical interpolation routine. Fixed isssues related to it --- wrf_open/var/src/python/wrf/var/decorators.py | 22 +- wrf_open/var/src/python/wrf/var/extension.py | 65 ++++- wrf_open/var/src/python/wrf/var/geoht.py | 7 + wrf_open/var/src/python/wrf/var/interp.py | 131 +++++++++- wrf_open/var/src/python/wrf/var/temp.py | 2 +- wrf_open/var/src/python/wrf/var/util.py | 65 ++++- wrf_open/var/src/python/wrf/var/wrfext.f90 | 243 ++++++++++++------ wrf_open/var/src/python/wrf/var/wrfext.pyf | 51 ++++ wrf_open/var/test/utests.py | 21 +- 9 files changed, 501 insertions(+), 106 deletions(-) diff --git a/wrf_open/var/src/python/wrf/var/decorators.py b/wrf_open/var/src/python/wrf/var/decorators.py index 824f610..bb52605 100644 --- a/wrf_open/var/src/python/wrf/var/decorators.py +++ b/wrf_open/var/src/python/wrf/var/decorators.py @@ -1,12 +1,12 @@ from functools import wraps from inspect import getargspec -from itertools import product import numpy as n import numpy.ma as ma from wrf.var.units import do_conversion, check_units from wrf.var.destag import destagger +from wrf.var.util import iter_left_indexes __all__ = ["convert_units", "handle_left_iter", "uvmet_left_iter", "handle_casting"] @@ -43,22 +43,6 @@ def convert_units(unit_type, alg_unit): return convert_decorator -def _left_indexes(dims): - """A generator which yields the iteration tuples for a sequence of - dimensions sizes. - - For example, if an array shape is (3,3), then this will yield: - - (0,0), (0,1), (1,0), (1,1) - - Arguments: - - - dims - a sequence of dimensions sizes (e.g. ndarry.shape) - - """ - arg = [xrange(dim) for dim in dims] - for idxs in product(*arg): - yield idxs def _calc_out_dims(outvar, left_dims): left_dims = [x for x in left_dims] @@ -109,7 +93,7 @@ def handle_left_iter(ref_var_expected_dims, ref_var_idx=-1, extra_dims = [ref_var_shape[x] for x in xrange(extra_dim_num)] out_inited = False - for left_idxs in _left_indexes(extra_dims): + for left_idxs in iter_left_indexes(extra_dims): # Make the left indexes plus a single slice object # The single slice will handle all the dimensions to # the right (e.g. [1,1,:]) @@ -232,7 +216,7 @@ def uvmet_left_iter(): output = n.zeros(outdims, u.dtype) - for left_idxs in _left_indexes(extra_dims): + for left_idxs in iter_left_indexes(extra_dims): # Make the left indexes plus a single slice object # The single slice will handle all the dimensions to # the right (e.g. [1,1,:]) diff --git a/wrf_open/var/src/python/wrf/var/extension.py b/wrf_open/var/src/python/wrf/var/extension.py index acc7e5a..cfa1c13 100755 --- a/wrf_open/var/src/python/wrf/var/extension.py +++ b/wrf_open/var/src/python/wrf/var/extension.py @@ -8,7 +8,8 @@ from wrf.var._wrfext import (f_interpz3d, f_interp2dxy,f_interp1d, f_computeuvmet, f_computeomega, f_computetv, f_computewetbulb, f_computesrh, f_computeuh, f_computepw, f_computedbz, - f_lltoij, f_ijtoll, f_converteta, f_computectt) + f_lltoij, f_ijtoll, f_converteta, f_computectt, + f_monotonic, f_filter2d, f_vintrp) from wrf.var._wrfcape import f_computecape from wrf.var.decorators import (handle_left_iter, uvmet_left_iter, handle_casting) @@ -361,6 +362,68 @@ def computectt(p_hpa,tk,qice,qcld,qv,ght,ter,haveqci): return res.T +@handle_left_iter(2,0,ignore_args=(1,)) +@handle_casting(arg_idxs=(0,)) +def smooth2d(field, passes): + # Unlike NCL, this routine will not modify the values in place, but + # copies the original data before modifying it. This allows the decorator + # to work properly and also behaves like the other methods. + + if isinstance(field, n.ma.MaskedArray): + missing = field.fill_value + else: + missing = Constants.DEFAULT_FILL + + field_copy = field.copy() + field_tmp = n.zeros(field_copy.shape, field_copy.dtype) + + f_filter2d(field_copy.T, + field_tmp.T, + missing, + passes) + + # Don't transpose here since the fortran routine is not returning an + # array. It's only modifying the existing array. + return field_copy + +@handle_left_iter(3,0,ignore_args=(3,4,5)) +@handle_casting(arg_idxs=(0,1,2)) +def monotonic(var,lvprs,coriolis,idir,delta,icorsw): + res = f_monotonic(var.T, + lvprs.T, + coriolis.T, + idir, + delta, + icorsw) + + return res.T + +# TODO: Check those decorators, they might be wrong +@handle_left_iter(3,0,ignore_args=(9,10,11,12,13,14)) +@handle_casting(arg_idxs=(0,1,2,3,4,5,6,7,8,9)) +def vintrp(field,pres,tk,qvp,ght,terrain,sfp,smsfp, + vcarray,interp_levels,icase,extrap,vcor,logp, + missing): + + res = f_vintrp(field.T, + pres.T, + tk.T, + qvp.T, + ght.T, + terrain.T, + sfp.T, + smsfp.T, + vcarray.T, + interp_levels, + icase, + extrap, + vcor, + logp, + missing) + + return res.T + + diff --git a/wrf_open/var/src/python/wrf/var/geoht.py b/wrf_open/var/src/python/wrf/var/geoht.py index 21e2529..e8396c8 100755 --- a/wrf_open/var/src/python/wrf/var/geoht.py +++ b/wrf_open/var/src/python/wrf/var/geoht.py @@ -35,6 +35,13 @@ def _get_geoht(wrfnc, timeidx, height=True, msl=True): if msl: return geopt_unstag / Constants.G else: + # Due to broadcasting with multifile/multitime, the 2D terrain + # array needs to be reshaped to a 3D array so the right dims + # line up + new_dims = list(hgt.shape) + new_dims.insert(-2,1) + hgt = hgt.reshape(new_dims) + return (geopt_unstag / Constants.G) - hgt else: return geopt_unstag diff --git a/wrf_open/var/src/python/wrf/var/interp.py b/wrf_open/var/src/python/wrf/var/interp.py index 973192c..1a3cec2 100755 --- a/wrf_open/var/src/python/wrf/var/interp.py +++ b/wrf_open/var/src/python/wrf/var/interp.py @@ -3,11 +3,17 @@ from math import floor, ceil import numpy as n import numpy.ma as ma -from wrf.var.extension import interpz3d,interp2dxy,interp1d +from wrf.var.extension import (interpz3d,interp2dxy,interp1d, + smooth2d,monotonic,vintrp) from wrf.var.decorators import handle_left_iter, handle_casting -from wrf.var.constants import Constants +from wrf.var.util import extract_vars, is_staggered +from wrf.var.constants import Constants, ConversionFactors +from wrf.var.terrain import get_terrain +from wrf.var.geoht import get_height +from wrf.var.temp import get_theta, get_temp, get_eth +from wrf.var.pressure import get_pressure -__all__ = ["interplevel", "vertcross", "interpline"] +__all__ = ["interplevel", "vertcross", "interpline", "vinterp"] # Note: Extension decorator is good enough to handle left dims def interplevel(data3d,zdata,desiredloc,missingval=Constants.DEFAULT_FILL): @@ -258,6 +264,125 @@ def interpline(data2d, pivot_point=None, return var1dtmp[0,:] +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} + + # These constants match what's in the fortran code. + rgas = 287.04 #J/K/kg + ussalr = .0065 # deg C per m, avg lapse rate + sclht = rgas*256./9.81 + + # interp_levels might be a list or tuple, make a numpy array + if not isinstance(interp_levels, n.ndarray): + interp_levels = n.array(interp_levels, "float64") + + # TODO: Check if field is staggered + if is_staggered(field, wrfnc): + raise RuntimeError("Please unstagger field in the vertical") + + # Check for valid coord + if vert_coord not in valid_coords: + raise RuntimeError("'%s' is not a valid vertical " + "coordinate type" % vert_coord) + + # Check for valid field type + if field_type not in valid_field_types: + raise RuntimeError("'%s' is not a valid field type" % field_type) + + log_p_int = 1 if log_p else 0 + + icase = 0 + extrap = 0 + + if extrapolate: + extrap = 1 + icase = icase_lookup[field_type.lower()] + + # Extract vriables + timeidx = -1 # Should this be an argument? + ncvars = extract_vars(wrfnc, timeidx, vars=("PSFC", "QVAPOR", "F")) + + sfp = ncvars["PSFC"] * ConversionFactors.PA_TO_HPA + qv = ncvars["QVAPOR"] + coriolis = ncvars["F"] + + terht = get_terrain(wrfnc, timeidx) + t = get_theta(wrfnc, timeidx) + tk = get_temp(wrfnc, timeidx, units="k") + p = get_pressure(wrfnc, timeidx) + ght = get_height(wrfnc, timeidx, msl=True) + ht_agl = get_height(wrfnc, timeidx, msl=False) + + smsfp = smooth2d(sfp, 3) + + # Vertical coordinate type + vcor = 0 + + if vert_coord in ("pressure", "pres"): + vcor = 1 + vcord_array = p * ConversionFactors.PA_TO_HPA + + elif vert_coord == "ght_msl": + vcor = 2 + vcord_array = n.exp(-ght/sclht) + + elif vert_coord == "ght_agl": + vcor = 3 + vcord_array = n.exp(-ht_agl/sclht) + + elif vert_coord == "theta": + vcor = 4 + idir = 1 + icorsw = 0 + delta = 0.01 + + p_hpa = p * ConversionFactors.PA_TO_HPA + + vcord_array = monotonic(t,p_hpa,coriolis,idir,delta,icorsw) + + elif vert_coord == "theta-e": + vcor = 5 + icorsw = 0 + idir = 1 + delta = 0.01 + + eth = get_eth(wrfnc, timeidx) + + p_hpa = p * ConversionFactors.PA_TO_HPA + + vcord_array = monotonic(eth,p_hpa,coriolis,idir,delta,icorsw) + # We only extrapolate temperature fields below ground if we are + # interpolating to pressure or height vertical surfaces + icase = 0 + + # Set the missing value + if isinstance(field, n.ma.MaskedArray): + missing = field.fill_value + else: + 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) + + return ma.masked_values(res,missing) + diff --git a/wrf_open/var/src/python/wrf/var/temp.py b/wrf_open/var/src/python/wrf/var/temp.py index 10f27f4..70999a3 100755 --- a/wrf_open/var/src/python/wrf/var/temp.py +++ b/wrf_open/var/src/python/wrf/var/temp.py @@ -43,7 +43,7 @@ def get_eth(wrfnc, timeidx=0, units="k"): full_p = p + pb tk = computetk(full_p, full_t) - eth = computeeth ( qv, tk, full_p ) + eth = computeeth(qv, tk, full_p) return eth diff --git a/wrf_open/var/src/python/wrf/var/util.py b/wrf_open/var/src/python/wrf/var/util.py index 9f00990..97f55c1 100644 --- a/wrf_open/var/src/python/wrf/var/util.py +++ b/wrf_open/var/src/python/wrf/var/util.py @@ -1,10 +1,13 @@ from collections import Iterable +from itertools import product import datetime as dt import numpy as n __all__ = ["extract_vars", "extract_global_attrs", "extract_dim", - "combine_files", "is_standard_wrf_var", "extract_times"] + "combine_files", "is_standard_wrf_var", "extract_times", + "iter_left_indexes", "get_left_indexes", "get_right_slices", + "is_staggered"] def _is_multi_time(timeidx): if timeidx == -1: @@ -39,6 +42,9 @@ def extract_global_attrs(wrfnc, attrs): return {attr:_get_attr(wrfnc, attr) for attr in attrlist} def extract_dim(wrfnc, dim): + if _is_multi_file(wrfnc): + wrfnc = wrfnc[0] + d = wrfnc.dimensions[dim] if not isinstance(d, int): return len(d) #netCDF4 @@ -113,6 +119,63 @@ def is_standard_wrf_var(wrfnc, var): if multifile: wrfnc = wrfnc[0] return var in wrfnc.variables + +def is_staggered(var, wrfnc): + we = extract_dim(wrfnc, "west_east") + sn = extract_dim(wrfnc, "south_north") + bt = extract_dim(wrfnc, "bottom_top") + + if (var.shape[-1] != we or var.shape[-2] != sn or var.shape[-3] != bt): + return True + + return False + +def get_left_indexes(ref_var, expected_dims): + """Returns the extra left side dimensions for a variable with an expected + shape. + + For example, if a 2D variable contains an additional left side dimension + for time, this will return the time dimension size. + + """ + extra_dim_num = ref_var.ndim - expected_dims + + if (extra_dim_num == 0): + return [] + + return [ref_var.shape[x] for x in xrange(extra_dim_num)] + +def iter_left_indexes(dims): + """A generator which yields the iteration tuples for a sequence of + dimensions sizes. + + For example, if an array shape is (3,3), then this will yield: + + (0,0), (0,1), (1,0), (1,1) + + Arguments: + + - dims - a sequence of dimensions sizes (e.g. ndarry.shape) + + """ + arg = [xrange(dim) for dim in dims] + for idxs in product(*arg): + yield idxs + +def get_right_slices(var, right_ndims, fixed_val=0): + extra_dim_num = var.ndim - right_ndims + if extra_dim_num == 0: + return [slice(None,None,None)] * right_ndims + + return [fixed_val]*extra_dim_num + [slice(None,None,None)]*right_ndims + + + + + + + + diff --git a/wrf_open/var/src/python/wrf/var/wrfext.f90 b/wrf_open/var/src/python/wrf/var/wrfext.f90 index d68fae3..4837418 100755 --- a/wrf_open/var/src/python/wrf/var/wrfext.f90 +++ b/wrf_open/var/src/python/wrf/var/wrfext.f90 @@ -1,6 +1,8 @@ SUBROUTINE f_interpz3d(data3d,zdata,desiredloc,missingval,out2d,nx,ny,nz) + IMPLICIT NONE + INTEGER, INTENT(IN) :: nx,ny,nz REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: data3d REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: out2d @@ -45,10 +47,13 @@ SUBROUTINE f_interpz3d(data3d,zdata,desiredloc,missingval,out2d,nx,ny,nz) END DO RETURN + END SUBROUTINE f_interpz3d SUBROUTINE f_interp2dxy(v3d,xy,v2d,nx,ny,nz,nxy) + IMPLICIT NONE + INTEGER,INTENT(IN) :: nx,ny,nz,nxy REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: v3d REAL(KIND=8),DIMENSION(nxy,nz),INTENT(OUT) :: v2d @@ -73,10 +78,13 @@ SUBROUTINE f_interp2dxy(v3d,xy,v2d,nx,ny,nz,nxy) END DO RETURN + END SUBROUTINE f_interp2dxy SUBROUTINE f_interp1d(v_in,z_in,z_out,vmsg,v_out,nz_in,nz_out) + IMPLICIT NONE + INTEGER,INTENT(IN) :: NZ_IN,NZ_OUT REAL(KIND=8),DIMENSION(nz_in),INTENT(IN) :: v_in,z_in REAL(KIND=8),DIMENSION(nz_out),INTENT(IN) :: z_out @@ -116,6 +124,7 @@ SUBROUTINE f_interp1d(v_in,z_in,z_out,vmsg,v_out,nz_in,nz_out) END DO RETURN + END SUBROUTINE f_interp1d ! This routine assumes @@ -132,7 +141,9 @@ END SUBROUTINE f_interp1d SUBROUTINE f_computeslp(z,t,p,q,t_sea_level,t_surf,level,throw_exception,& sea_level_pressure,nx,ny,nz) + IMPLICIT NONE + EXTERNAL throw_exception ! Estimate sea level pressure. INTEGER, INTENT(IN) :: nx,ny,nz @@ -269,11 +280,14 @@ SUBROUTINE f_computeslp(z,t,p,q,t_sea_level,t_surf,level,throw_exception,& ! * sea_level_pressure(20,2),sea_level_pressure(20,3) RETURN + END SUBROUTINE f_computeslp ! Temperature from potential temperature in kelvin. SUBROUTINE f_computetk(pressure,theta,tk,nx) + IMPLICIT NONE + INTEGER, INTENT(IN) :: nx REAL(KIND=8) :: pi REAL(KIND=8), DIMENSION(nx), INTENT(IN) :: pressure @@ -289,10 +303,12 @@ SUBROUTINE f_computetk(pressure,theta,tk,nx) END DO RETURN + END SUBROUTINE f_computetk ! Dewpoint. Note: 1D array arguments. SUBROUTINE f_computetd(pressure,qv_in,td,nx) + IMPLICIT NONE INTEGER, INTENT(IN) :: nx @@ -315,10 +331,12 @@ SUBROUTINE f_computetd(pressure,qv_in,td,nx) END DO RETURN + END SUBROUTINE f_computetd ! Relative Humidity. Note: 1D array arguments SUBROUTINE f_computerh(qv,p,t,rh,nx) + IMPLICIT NONE INTEGER, INTENT(IN) :: nx @@ -345,10 +363,13 @@ SUBROUTINE f_computerh(qv,p,t,rh,nx) END DO RETURN + END SUBROUTINE f_computerh SUBROUTINE f_computeabsvort(u,v,msfu,msfv,msft,cor,dx,dy,av,nx,ny,nz,nxp1,nyp1) + IMPLICIT NONE + INTEGER, INTENT(IN) :: nx,ny,nz,nxp1,nyp1 REAL(KIND=8),DIMENSION(nxp1,ny,nz),INTENT(IN) :: u REAL(KIND=8),DIMENSION(nx,nyp1,nz),INTENT(IN) :: v @@ -392,11 +413,14 @@ SUBROUTINE f_computeabsvort(u,v,msfu,msfv,msft,cor,dx,dy,av,nx,ny,nz,nxp1,nyp1) END DO RETURN + END SUBROUTINE f_computeabsvort SUBROUTINE f_computepvo(u,v,theta,prs,msfu,msfv,msft,cor,dx,dy,pv,nx,ny,nz,nxp1,nyp1) + IMPLICIT NONE + INTEGER,INTENT(IN) :: nx,ny,nz,nxp1,nyp1 REAL(KIND=8),DIMENSION(nxp1,ny,nz),INTENT(IN) :: u REAL(KIND=8),DIMENSION(nx,nyp1,nz),INTENT(IN) :: v @@ -454,12 +478,16 @@ SUBROUTINE f_computepvo(u,v,theta,prs,msfu,msfv,msft,cor,dx,dy,pv,nx,ny,nz,nxp1, END DO END DO END DO + RETURN + END SUBROUTINE f_computepvo ! Theta-e SUBROUTINE f_computeeth(qvp,tmk,prs,eth,miy,mjx,mkzh) + IMPLICIT NONE + ! Input variables ! Qvapor [g/kg] REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: qvp @@ -513,11 +541,13 @@ SUBROUTINE f_computeeth(qvp,tmk,prs,eth,miy,mjx,mkzh) END DO RETURN + END SUBROUTINE f_computeeth SUBROUTINE f_computeuvmet(u,v,longca,longcb,flong,flat,& cen_long,cone,rpd,istag,is_msg_val,umsg,vmsg,uvmetmsg,& uvmet,nx,ny,nxp1,nyp1,nz) + IMPLICIT NONE ! ISTAG should be 0 if the U,V grids are not staggered. @@ -594,11 +624,14 @@ SUBROUTINE f_computeuvmet(u,v,longca,longcb,flong,flat,& END DO RETURN + END SUBROUTINE f_computeuvmet SUBROUTINE f_computeomega(qvp,tmk,www,prs,omg,mx,my,mz) + IMPLICIT NONE + INTEGER,INTENT(IN) :: mx, my, mz REAL(KIND=8),INTENT(IN),DIMENSION(mx,my,mz) :: qvp REAL(KIND=8),INTENT(IN),DIMENSION(mx,my,mz) :: tmk @@ -622,10 +655,13 @@ SUBROUTINE f_computeomega(qvp,tmk,www,prs,omg,mx,my,mz) END DO ! RETURN + END SUBROUTINE f_computeomega SUBROUTINE f_computetv(temp,qv,tv,nx,ny,nz) + IMPLICIT NONE + INTEGER, INTENT(IN) :: nx,ny,nz REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: temp REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: qv @@ -641,12 +677,16 @@ SUBROUTINE f_computetv(temp,qv,tv,nx,ny,nz) END DO END DO END DO + RETURN + END SUBROUTINE f_computetv ! Need to deal with the fortran stop statements SUBROUTINE f_computewetbulb(prs,tmk,qvp,PSADITHTE,PSADIPRS,PSADITMK,throw_exception,twb,nx,ny,nz) + IMPLICIT NONE + EXTERNAL throw_exception INTEGER,INTENT(IN) :: nx, ny, nz REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: prs @@ -751,10 +791,13 @@ SUBROUTINE f_computewetbulb(prs,tmk,qvp,PSADITHTE,PSADIPRS,PSADITMK,throw_except END DO RETURN + END SUBROUTINE f_computewetbulb SUBROUTINE f_computesrh(u, v, ght, ter, top, sreh, miy, mjx, mkzh) + IMPLICIT NONE + INTEGER,INTENT(IN) :: miy, mjx, mkzh REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: u, v, ght REAL(KIND=8),INTENT(IN) :: top @@ -825,10 +868,13 @@ SUBROUTINE f_computesrh(u, v, ght, ter, top, sreh, miy, mjx, mkzh) END DO RETURN + END SUBROUTINE f_computesrh SUBROUTINE f_computeuh(zp,mapfct,dx,dy,uhmnhgt,uhmxhgt,us,vs,w,tem1,tem2,uh,nx,ny,nz,nzp1) + IMPLICIT NONE + INTEGER, INTENT(IN) :: nx,ny,nz,nzp1 REAL(KIND=8),DIMENSION(nx,ny,nzp1),INTENT(IN) :: zp REAL(KIND=8),DIMENSION(nx,ny),INTENT(IN) :: mapfct @@ -954,10 +1000,13 @@ SUBROUTINE f_computeuh(zp,mapfct,dx,dy,uhmnhgt,uhmxhgt,us,vs,w,tem1,tem2,uh,nx,n uh = uh * 1000. ! Scale according to Kain et al. (2008) RETURN + END SUBROUTINE f_computeuh SUBROUTINE f_computepw(p,tv,qv,ht,zdiff,pw,nx,ny,nz,nzh) + IMPLICIT NONE + REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: p,tv,qv REAL(KIND=8),DIMENSION(nx,ny,nzh),INTENT(IN) :: ht REAL(KIND=8),DIMENSION(nx,ny),INTENT(OUT) :: pw @@ -979,11 +1028,14 @@ SUBROUTINE f_computepw(p,tv,qv,ht,zdiff,pw,nx,ny,nz,nzh) RETURN + END SUBROUTINE f_computepw SUBROUTINE f_computedbz(prs,tmk,qvp,qra,qsn,qgr,sn0,ivarint,iliqskin,dbz,nx,ny,nz) + IMPLICIT NONE + ! Arguments INTEGER, INTENT(IN) :: nx,ny,nz INTEGER, INTENT(IN) :: sn0,ivarint,iliqskin @@ -1135,10 +1187,13 @@ SUBROUTINE f_computedbz(prs,tmk,qvp,qra,qsn,qgr,sn0,ivarint,iliqskin,dbz,nx,ny,n END DO RETURN + END SUBROUTINE f_computedbz SUBROUTINE rotatecoords(ilat,ilon,olat,olon,lat_np,lon_np,lon_0,direction) + IMPLICIT NONE + REAL(KIND=8),INTENT(IN) :: ilat,ilon REAL(KIND=8),INTENT(OUT) :: olat,olon REAL(KIND=8),INTENT(IN) :: lat_np,lon_np,lon_0 @@ -1183,12 +1238,15 @@ SUBROUTINE rotatecoords(ilat,ilon,olat,olon,lat_np,lon_np,lon_0,direction) olon = DEG_PER_RAD*(ATAN2(sinlam,coslam)-dlam-lam_0+lam_np) RETURN + END SUBROUTINE rotatecoords SUBROUTINE f_lltoij(map_proj,truelat1,truelat2,stdlon,lat1,lon1,& pole_lat,pole_lon,knowni,knownj,dx,latinc,& loninc,lat,lon,throw_exception,loc) + IMPLICIT NONE + EXTERNAL throw_exception ! Converts input lat/lon values to the cartesian (i,j) value ! for the given projection. @@ -1386,6 +1444,7 @@ SUBROUTINE f_lltoij(map_proj,truelat1,truelat2,stdlon,lat1,lon1,& loc(2) = i RETURN + END SUBROUTINE f_lltoij @@ -1396,6 +1455,7 @@ SUBROUTINE f_ijtoll(map_proj,truelat1,truelat2,stdlon,lat1,lon1,& ! converts input lat/lon values to the cartesian (i,j) value ! for the given projection. IMPLICIT NONE + EXTERNAL throw_exception INTEGER,INTENT(IN) :: map_proj REAL(KIND=8),INTENT(IN) :: stdlon @@ -1626,6 +1686,7 @@ SUBROUTINE f_ijtoll(map_proj,truelat1,truelat2,stdlon,lat1,lon1,& loc(2) = lon RETURN + END SUBROUTINE f_ijtoll ! Eta = (p - ptop) / (psfc - ptop) @@ -1638,6 +1699,7 @@ SUBROUTINE f_converteta(full_t, znu, psfc, ptop, pcalc, mean_t, temp_t,& z, nx,ny,nz) IMPLICIT NONE + REAL(KIND=8),DIMENSION(nx,ny,nz),INTENT(IN) :: full_t REAL(KIND=8),DIMENSION(nz),INTENT(IN) :: znu REAL(KIND=8),DIMENSION(nx,ny),INTENT(IN) :: psfc @@ -1692,6 +1754,7 @@ SUBROUTINE f_converteta(full_t, znu, psfc, ptop, pcalc, mean_t, temp_t,& END DO RETURN + END SUBROUTINE f_converteta @@ -1730,6 +1793,8 @@ SUBROUTINE f_computectt(prs,tk,qci,qcw,qvp,ght,ter,haveqci,ctt,ew,ns,nz) !miy = ns !mkzh = nz + prsctt = 0 ! removes the warning + ! Calculate the surface pressure DO j=1,ns DO i=1,ew @@ -1817,67 +1882,75 @@ SUBROUTINE f_computectt(prs,tk,qci,qcw,qvp,ght,ter,haveqci,ctt,ew,ns,nz) ! 40 CONTINUE ! 190 CONTINUE RETURN + END SUBROUTINE f_computectt SUBROUTINE f_filter2d(a, b, missing, it, nx, ny) + IMPLICIT NONE INTEGER, INTENT(IN) :: nx, ny, it - REAL(KIND=8), DIMENSION(nx, ny), INTENT(IN) :: a + REAL(KIND=8), DIMENSION(nx, ny), INTENT(INOUT) :: a REAL(KIND=8), INTENT(IN) :: missing - REAL(KIND=8), DIMENSION(nx, ny), INTENT(OUT) :: b + REAL(KIND=8), DIMENSION(nx, ny), INTENT(INOUT) :: b REAL(KIND=8), PARAMETER :: COEF=0.25D0 INTEGER :: i, j, iter DO iter=1,it - DO j=1,ny - DO i = 1,nx - b(i,j) = a(i,j) - END DO - END DO + DO j=1,ny + DO i = 1,nx + b(i,j) = a(i,j) + END DO + END DO - DO j=2,ny-1 - DO i=1,nx - IF (b(i,j-1).EQ.missing .OR. b(i,j).EQ.missing .OR. b(i,j+1).EQ.missing) THEN - a(i,j) = a(i,j) - ELSE - a(i,j) = a(i,j) + COEF*(b(i,j-1)-2*b(i,j) + b(i,j+1)) - END IF - END DO - END DO + DO j=2,ny-1 + DO i=1,nx + IF (b(i,j-1) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & + b(i,j+1) .EQ. missing) THEN + a(i,j) = a(i,j) + ELSE + a(i,j) = a(i,j) + COEF*(b(i,j-1)-2*b(i,j) + b(i,j+1)) + END IF + END DO + END DO - DO j=1,ny - DO i=2,nx-1 - IF (b(i-1,j).EQ.missing .OR. b(i,j).EQ.missing .OR. b(i+1,j).EQ.missing) THEN - a(i,j) = a(i,j) - ELSE - a(i,j) = a(i,j) + COEF*(b(i-1,j)-2*b(i,j)+b(i+1,j)) - END IF - END DO - END DO - c do j=1,ny - c do i=1,nx - c b(i,j) = a(i,j) - c enddo - c enddo - c do j=2,ny-1 - c do i=1,nx - c a(i,j) = a(i,j) - .99*coef*(b(i,j-1)-2*b(i,j)+b(i,j+1)) - c enddo - c enddo - c do j=1,ny - c do i=2,nx-1 - c a(i,j) = a(i,j) - .99*coef*(b(i-1,j)-2*b(i,j)+b(i+1,j)) - c enddo - c enddo + DO j=1,ny + DO i=2,nx-1 + IF (b(i-1,j) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & + b(i+1,j) .EQ. missing) THEN + a(i,j) = a(i,j) + ELSE + a(i,j) = a(i,j) + COEF*(b(i-1,j)-2*b(i,j)+b(i+1,j)) + END IF + END DO + END DO + ! do j=1,ny + ! do i=1,nx + ! b(i,j) = a(i,j) + ! enddo + ! enddo + ! do j=2,ny-1 + ! do i=1,nx + ! a(i,j) = a(i,j) - .99*coef*(b(i,j-1)-2*b(i,j)+b(i,j+1)) + ! enddo + ! enddo + ! do j=1,ny + ! do i=2,nx-1 + ! a(i,j) = a(i,j) - .99*coef*(b(i-1,j)-2*b(i,j)+b(i+1,j)) + ! enddo + ! enddo END DO + RETURN -END + +END SUBROUTINE f_filter2d SUBROUTINE f_monotonic(out,in,lvprs,cor,idir,delta,ew,ns,nz,icorsw) + IMPLICIT NONE + INTEGER, INTENT(IN) :: idir,ew,ns,nz,icorsw REAL(KIND=8), INTENT(IN) :: delta REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: in @@ -1887,37 +1960,39 @@ SUBROUTINE f_monotonic(out,in,lvprs,cor,idir,delta,ew,ns,nz,icorsw) INTEGER :: i,j,k,ripk,k300 + k300 = 0 ! removes the warning + DO j=1,ns DO i=1,ew - IF (icorsw.EQ.1.AND.cor(i,j).LT.0.) THEN + IF (icorsw .EQ. 1 .AND. cor(i,j) .LT. 0.) THEN DO k=1,nz - in(i,j,k)=-in(i,j,k) + in(i,j,k) = -in(i,j,k) END DO END IF ! First find k index that is at or below (height-wise) ! the 300 hPa level. DO k = 1,nz - ripk = nz-k+1 + ripk = nz-k+1 IF (lvprs(i,j,k) .LE. 300.d0) THEN - k300=k + k300 = k EXIT END IF END DO DO k = k300, 1,-1 - IF (idir.EQ.1) THEN - out(i,j,k)=MIN(in(i,j,k),in(i,j,k+1)+delta) - ELSE IF (idir.EQ.-1) THEN - out(i,j,k)=MAX(in(i,j,k),in(i,j,k+1)-delta) + IF (idir .EQ. 1) THEN + out(i,j,k) = MIN(in(i,j,k), in(i,j,k+1)+delta) + ELSE IF (idir .EQ. -1) THEN + out(i,j,k) = MAX(in(i,j,k), in(i,j,k+1)-delta) END IF END DO DO k = k300+1, nz - IF (idir.EQ.1) THEN - out(i,j,k)=MAX(in(i,j,k),in(i,j,k-1)-delta) - ELSE IF (idir.EQ.-1) THEN - out(i,j,k)=MIN(in(i,j,k),in(i,j,k-1)+delta) + IF (idir .EQ. 1) THEN + out(i,j,k) = MAX(in(i,j,k), in(i,j,k-1)-delta) + ELSE IF (idir .EQ. -1) THEN + out(i,j,k) = MIN(in(i,j,k), in(i,j,k-1)+delta) END IF END DO END DO @@ -1929,6 +2004,7 @@ END SUBROUTINE f_monotonic FUNCTION f_intrpvalue (wvalp0,wvalp1,vlev,vcp0,vcp1,icase) + IMPLICIT NONE INTEGER, INTENT(IN) :: icase @@ -1976,25 +2052,26 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& icase,ew,ns,nz,extrap,vcor,logp,rmsg) IMPLICIT NONE + INTEGER, INTENT(IN) :: ew,ns,nz,icase,extrap INTEGER, INTENT(IN) :: vcor,numlevels,logp REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: datain,pres,tk,qvp REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: ght REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: terrain,sfp,smsfp REAL(KIND=8), DIMENSION(ew,ns,numlevels), INTENT(OUT) :: dataout - REAL(KIND-8), DIMENSION(ew,ns,nz), INTENT(IN) :: vcarray(ew,ns,nz) - REAL(KIND=8), DIMENSION(numlevels), INTENT(IN) :: interp_levels(numlevels) + REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: vcarray + REAL(KIND=8), DIMENSION(numlevels), INTENT(IN) :: interp_levels REAL(KIND=8), INTENT(IN) :: rmsg - INTEGER :: njx,niy,nreqlvs,ripk - INTEGER :: i,j,k,itriv,kupper - INTEGER :: ifound,miy,mjx,isign + INTEGER :: nreqlvs,ripk !njx,niy + INTEGER :: i,j,k,kupper !itriv + INTEGER :: ifound,isign !miy,mjx REAL(KIND=8), DIMENSION(ew,ns) :: tempout REAL(KIND=8) :: rlevel,vlev,diff REAL(KIND=8) :: tmpvlev REAL(KIND=8) :: vcp1,vcp0,valp0,valp1 - REAL(KIND=8) :: cvc - REAL(KIND=8) :: qvlhsl,ttlhsl,vclhsl,vctophsl +! REAL(KIND=8) :: cvc + REAL(KIND=8) :: vclhsl,vctophsl !qvlhsl,ttlhsl REAL(KIND=8) :: f_intrpvalue REAL(KIND=8) :: plhsl,zlhsl,ezlhsl,tlhsl,psurf,pratio,tlev REAL(KIND=8) :: ezsurf,psurfsm,zsurf,qvapor,vt @@ -2023,6 +2100,12 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& REAL(KIND=8), PARAMETER :: GAMMA = RGAS/CP REAL(KIND=8), PARAMETER :: GAMMAMD = RGASMD-CPMD + ! Removes the warnings for uninitialized variables + cvcord = '' + plev = 0 + zlev = 0 + vlev = 0 + IF (vcor .EQ. 1) THEN cvcord = 'p' ELSE IF ((vcor .EQ. 2) .OR. (vcor .EQ. 3)) THEN @@ -2044,7 +2127,7 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& DO nreqlvs = 1,numlevels IF (cvcord .EQ. 'z') THEN - !Convert rlevel to meters from km + ! Convert rlevel to meters from km rlevel = interp_levels(nreqlvs) * 1000.d0 vlev = EXP(-rlevel/SCLHT) ELSE IF (cvcord .EQ. 'p') THEN @@ -2056,7 +2139,7 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& DO j=1,ns DO i=1,ew - !Get the interpolated value that is within the model domain + ! Get the interpolated value that is within the model domain ifound = 0 DO k = 1,nz-1 ripk = nz-k+1 @@ -2094,7 +2177,8 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& !115 CONTINUE IF (ifound .EQ. 1) THEN !Grid point is in the model domain - GOTO 333 ! CYCLE + !GOTO 333 ! CYCLE + CYCLE END IF !If the user has requested no extrapolatin then just assign @@ -2122,9 +2206,9 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& CYCLE END IF - ! Only remaining possibility is that the specified level is below - ! lowest model level. If lowest model level value is missing, - ! set interpolated value to missing. + ! Only remaining possibility is that the specified level is below + ! lowest model level. If lowest model level value is missing, + ! set interpolated value to missing. IF (datain(i,j,1) .EQ. rmsg) THEN tempout(i,j) = rmsg @@ -2133,18 +2217,18 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& 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. + ! 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. tempout(i,j) = datain(i,j,1) - ! For the special cases of pressure on height levels or height on - ! pressure levels, or temperature-related variables on pressure or - ! height levels, perform a special extrapolation based on - ! US Standard Atmosphere. Here we calcualate the surface pressure - ! with the altimeter equation. This is how RIP calculates the - ! surface pressure. + ! For the special cases of pressure on height levels or height on + ! pressure levels, or temperature-related variables on pressure or + ! height levels, perform a special extrapolation based on + ! US Standard Atmosphere. Here we calcualate the surface pressure + ! with the altimeter equation. This is how RIP calculates the + ! surface pressure. IF (icase .GT. 0) THEN plhsl = pres(i,j,1) * 0.01d0 !pressure at lowest model level zlhsl = ght(i,j,1) !grid point height a lowest model level @@ -2221,7 +2305,7 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& END IF ELSE IF (cvcord .EQ. 'z') THEN zlev = -sclht*LOG(vlev) - plev = pbot*(1. + USSALR/vt*(zbot-zlev))**/EXPONI + plev = pbot*(1. + USSALR/vt*(zbot-zlev))**EXPONI IF (icase .EQ. 1) THEN tempout(i,j) = plev !GOTO 333 ! CYCLE @@ -2245,9 +2329,9 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& ! extraolation for equivalent potential temperature ELSE IF (icase .EQ. 6) THEN e = qvapor*plev/(EPS + qvapor) - tlcl = TLCLC1/(LOG(tlev**TLCLC2/e)-TLCLC3) + TLCLC4 + tlcl = TLCLC1/(LOG(tlev**TLCLC2/e) - TLCLC3) + TLCLC4 tempout(i,j)=tlev*(1000.d0/plev)**(gammam)*& - EXP((THTECON1/tlcl-THTECON2)*& + EXP((THTECON1/tlcl - THTECON2)*& qvapor*(1. + THTECON3*qvapor)) END IF END IF @@ -2267,7 +2351,8 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& END DO !end for the nreqlvs RETURN -END SUBROUTINE f_vinterp + +END SUBROUTINE f_vintrp diff --git a/wrf_open/var/src/python/wrf/var/wrfext.pyf b/wrf_open/var/src/python/wrf/var/wrfext.pyf index 8c54b8f..174f79b 100644 --- a/wrf_open/var/src/python/wrf/var/wrfext.pyf +++ b/wrf_open/var/src/python/wrf/var/wrfext.pyf @@ -327,6 +327,57 @@ python module _wrfext ! in integer, optional,intent(in),check(shape(prs,1)==ns),depend(prs) :: ns=shape(prs,1) integer, optional,intent(in),check(shape(prs,2)==nz),depend(prs) :: nz=shape(prs,2) end subroutine f_computectt + subroutine f_filter2d(a,b,missing,it,nx,ny) ! in :_wrfext:wrfext.f90 + real(kind=8) dimension(nx,ny),intent(inout) :: a + real(kind=8) dimension(nx,ny),intent(inout),depend(nx,ny) :: b + real(kind=8) intent(in) :: missing + integer intent(in) :: it + integer, optional,intent(in),check(shape(a,0)==nx),depend(a) :: nx=shape(a,0) + integer, optional,intent(in),check(shape(a,1)==ny),depend(a) :: ny=shape(a,1) + end subroutine f_filter2d + subroutine f_monotonic(out,in,lvprs,cor,idir,delta,ew,ns,nz,icorsw) ! in :_wrfext:wrfext.f90 + real(kind=8) dimension(ew,ns,nz),intent(out),depend(ew,ns,nz) :: out + real(kind=8) dimension(ew,ns,nz),intent(inout) :: in + real(kind=8) dimension(ew,ns,nz),depend(ew,ns,nz) :: lvprs + real(kind=8) dimension(ew,ns),depend(ew,ns) :: cor + integer intent(in) :: idir + real(kind=8) intent(in) :: delta + integer, optional,intent(in),check(shape(in,0)==ew),depend(in) :: ew=shape(in,0) + integer, optional,intent(in),check(shape(in,1)==ns),depend(in) :: ns=shape(in,1) + integer, optional,intent(in),check(shape(in,2)==nz),depend(in) :: nz=shape(in,2) + integer intent(in) :: icorsw + end subroutine f_monotonic + function f_intrpvalue(wvalp0,wvalp1,vlev,vcp0,vcp1,icase) ! in :_wrfext:wrfext.f90 + real(kind=8) intent(in) :: wvalp0 + real(kind=8) intent(in) :: wvalp1 + real(kind=8) intent(in) :: vlev + real(kind=8) intent(in) :: vcp0 + real(kind=8) intent(in) :: vcp1 + integer intent(in) :: icase + real(kind=8) :: f_intrpvalue + end function f_intrpvalue + subroutine f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,sfp,smsfp,vcarray,interp_levels,numlevels,icase,ew,ns,nz,extrap,vcor,logp,rmsg) ! in :_wrfext:wrfext.f90 + real(kind=8) dimension(ew,ns,nz),intent(in) :: datain + real(kind=8) dimension(ew,ns,numlevels),intent(out),depend(ew,ns,numlevels) :: dataout + real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: pres + real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: tk + real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: qvp + real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: ght + real(kind=8) dimension(ew,ns),intent(in),depend(ew,ns) :: terrain + real(kind=8) dimension(ew,ns),intent(in),depend(ew,ns) :: sfp + real(kind=8) dimension(ew,ns),intent(in),depend(ew,ns) :: smsfp + real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: vcarray + real(kind=8) dimension(numlevels),intent(in) :: interp_levels + integer, optional,intent(in),check(len(interp_levels)>=numlevels),depend(interp_levels) :: numlevels=len(interp_levels) + integer intent(in) :: icase + integer, optional,intent(in),check(shape(datain,0)==ew),depend(datain) :: ew=shape(datain,0) + integer, optional,intent(in),check(shape(datain,1)==ns),depend(datain) :: ns=shape(datain,1) + integer, optional,intent(in),check(shape(datain,2)==nz),depend(datain) :: nz=shape(datain,2) + integer intent(in) :: extrap + integer intent(in) :: vcor + integer intent(in) :: logp + real(kind=8) intent(in) :: rmsg + end subroutine f_vintrp end interface end python module _wrfext diff --git a/wrf_open/var/test/utests.py b/wrf_open/var/test/utests.py index 77205cb..091f47b 100644 --- a/wrf_open/var/test/utests.py +++ b/wrf_open/var/test/utests.py @@ -5,7 +5,7 @@ import numpy.ma as ma import os, sys import subprocess -from wrf.var import getvar, interplevel, interpline, vertcross +from wrf.var import getvar, interplevel, interpline, vertcross, vinterp NCL_EXE = "/Users/ladwig/nclbuild/6.3.0/bin/ncl" TEST_FILE = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00" @@ -289,6 +289,23 @@ def make_interp_test(varname, wrf_in, referent, multi=False, end_point=end_point) nt.assert_allclose(t2_line1, t2_line2) + elif (varname == "vinterp"): + tk = getvar(in_wrfnc, "temp", units="k") + + interp_levels = [200,1000,50] + + field = vinterp(in_wrfnc, + field=tk, + vert_coord="theta", + interp_levels=interp_levels, + extrapolate=True, + field_type="tk", + log_p=True) + + print field.shape + print field + + return test @@ -306,7 +323,7 @@ if __name__ == "__main__": "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", "wa", "uvmet10", "uvmet", "z", "ctt", "cape_2d", "cape_3d"] - interp_methods = ["interplevel", "vertcross", "interpline"] + interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"] try: import netCDF4