nilyin
/
wrf-python
forked from 3rdparty/wrf-python
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Added vertical interpolation routine. Fixed isssues related to it

main
Bill Ladwig 9 years ago
parent
a44dd47892
commit
ad66cd6c07
9 changed files with 501 additions and 106 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 22
      wrf_open/var/src/python/wrf/var/decorators.py
  2. 65
      wrf_open/var/src/python/wrf/var/extension.py
  3. 7
      wrf_open/var/src/python/wrf/var/geoht.py
  4. 131
      wrf_open/var/src/python/wrf/var/interp.py
  5. 2
      wrf_open/var/src/python/wrf/var/temp.py
  6. 65
      wrf_open/var/src/python/wrf/var/util.py
  7. 175
      wrf_open/var/src/python/wrf/var/wrfext.f90
  8. 51
      wrf_open/var/src/python/wrf/var/wrfext.pyf
  9. 21
      wrf_open/var/test/utests.py

22
wrf_open/var/src/python/wrf/var/decorators.py
Unescape View File

@ -1,12 +1,12 @@ @@ -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): @@ -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, @@ -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(): @@ -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,:])

65
wrf_open/var/src/python/wrf/var/extension.py
Unescape View File

@ -8,7 +8,8 @@ from wrf.var._wrfext import (f_interpz3d, f_interp2dxy,f_interp1d, @@ -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): @@ -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

7
wrf_open/var/src/python/wrf/var/geoht.py
Unescape View File

@ -35,6 +35,13 @@ def _get_geoht(wrfnc, timeidx, height=True, msl=True): @@ -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

131
wrf_open/var/src/python/wrf/var/interp.py
Unescape View File

@ -3,11 +3,17 @@ from math import floor, ceil @@ -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, @@ -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)

2
wrf_open/var/src/python/wrf/var/temp.py
Unescape View File

@ -43,7 +43,7 @@ def get_eth(wrfnc, timeidx=0, units="k"): @@ -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

65
wrf_open/var/src/python/wrf/var/util.py
Unescape View File

@ -1,10 +1,13 @@ @@ -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): @@ -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
@ -114,6 +120,63 @@ def is_standard_wrf_var(wrfnc, var): @@ -114,6 +120,63 @@ def is_standard_wrf_var(wrfnc, var):
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

175
wrf_open/var/src/python/wrf/var/wrfext.f90
Unescape View File

@ -1,6 +1,8 @@ @@ -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) @@ -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) @@ -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) @@ -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 @@ -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,& @@ -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) @@ -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) @@ -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) @@ -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) @@ -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, @@ -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) @@ -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,& @@ -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) @@ -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) @@ -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 @@ -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) @@ -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 @@ -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) @@ -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 @@ -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) @@ -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,& @@ -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,& @@ -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,& @@ -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,& @@ -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,& @@ -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) @@ -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,15 +1882,17 @@ SUBROUTINE f_computectt(prs,tk,qci,qcw,qvp,ght,ter,haveqci,ctt,ew,ns,nz) @@ -1817,15 +1882,17 @@ 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
@ -1840,7 +1907,8 @@ SUBROUTINE f_filter2d(a, b, missing, it, nx, ny) @@ -1840,7 +1907,8 @@ SUBROUTINE f_filter2d(a, b, missing, it, nx, ny)
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
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))
@ -1850,34 +1918,39 @@ SUBROUTINE f_filter2d(a, b, missing, it, nx, ny) @@ -1850,34 +1918,39 @@ SUBROUTINE f_filter2d(a, b, missing, it, nx, ny)
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
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=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,11 +1960,13 @@ SUBROUTINE f_monotonic(out,in,lvprs,cor,idir,delta,ew,ns,nz,icorsw) @@ -1887,11 +1960,13 @@ 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
@ -1900,24 +1975,24 @@ SUBROUTINE f_monotonic(out,in,lvprs,cor,idir,delta,ew,ns,nz,icorsw) @@ -1900,24 +1975,24 @@ SUBROUTINE f_monotonic(out,in,lvprs,cor,idir,delta,ew,ns,nz,icorsw)
DO k = 1,nz
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 @@ -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,& @@ -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,& @@ -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,& @@ -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,& @@ -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,& @@ -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
@ -2221,7 +2305,7 @@ SUBROUTINE f_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,& @@ -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,& @@ -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,& @@ -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

51
wrf_open/var/src/python/wrf/var/wrfext.pyf
Unescape View File

@ -327,6 +327,57 @@ python module _wrfext ! in @@ -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

21
wrf_open/var/test/utests.py
Unescape View File

@ -5,7 +5,7 @@ import numpy.ma as ma @@ -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, @@ -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__": @@ -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

Write Preview
Loading…
Cancel
Save