Skip to content

Commit

Permalink
Merge pull request #40 from csyhuang/fix_issue_39
Browse files Browse the repository at this point in the history 
Fix issue 39
  • Loading branch information
@csyhuang
csyhuang authored Oct 19, 2021
2 parents 432755d + 71e1a0a commit 2de8e28
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 98 deletions.
34 changes: 24 additions & 10 deletions hn2016_falwa/oopinterface.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import numpy as np
import math
import warnings
import numpy as np
from scipy.interpolate import interp1d

from hn2016_falwa import utilities
Expand Down Expand Up @@ -44,8 +45,6 @@ class QGField(object):
Number of iteration by the Successive over-relaxation (SOR) solver to compute the reference states.
dz : float, optional
Size of uniform pseudoheight grids (in meters).
prefactor : float, optional
Vertical air density summed over height.
npart : int, optional
Number of partitions used to compute equivalent latitude.
If not initialized, it will be set to nlat.
Expand All @@ -72,12 +71,9 @@ class QGField(object):
"""

def __init__(self, xlon, ylat, plev,
u_field, v_field, t_field,
kmax=49, maxit=100000, dz=1000., prefactor=6500.,
npart=None, tol=1.e-5, rjac=0.95,
scale_height=SCALE_HEIGHT, cp=CP, dry_gas_constant=DRY_GAS_CONSTANT,
omega=EARTH_OMEGA, planet_radius=EARTH_RADIUS):
def __init__(self, xlon, ylat, plev, u_field, v_field, t_field, kmax=49, maxit=100000, dz=1000., npart=None,
tol=1.e-5, rjac=0.95, scale_height=SCALE_HEIGHT, cp=CP, dry_gas_constant=DRY_GAS_CONSTANT,
omega=EARTH_OMEGA, planet_radius=EARTH_RADIUS, prefactor=None):

"""
Create a QGField object.
Expand Down Expand Up @@ -155,7 +151,6 @@ def __init__(self, xlon, ylat, plev,
self.kmax = kmax
self.maxit = maxit
self.dz = dz
self.prefactor = prefactor
self.tol = tol
self.rjac = rjac

Expand All @@ -165,6 +160,15 @@ def __init__(self, xlon, ylat, plev,
self.dry_gas_constant = dry_gas_constant
self.omega = omega
self.planet_radius = planet_radius
self._compute_prefactor() # Compute normalization prefactor

# Modification on Oct 19, 2021 - deprecation of prefactor (it should be computed from kmax and dz)
if prefactor is not None:
warnings.warn(
"The optional input prefactor will be deprecated since it can be determined directly from " +
"kmax and dz. Given your input kmax = {kmax} and dz = {dz}, ".format(kmax=self.kmax, dz=self.dz) +
"the computed normalization prefactor is {prefactor}. ".format(prefactor=self.prefactor) +
"Your input value {input} would be ignored.".format(input=prefactor))

# === Variables that will be computed in methods ===
self._qgpv_temp = None
Expand Down Expand Up @@ -200,6 +204,16 @@ def __init__(self, xlon, ylat, plev,
self._lwa = None
self._divergence_eddy_momentum_flux = None

def _compute_prefactor(self):
"""
Private function. Compute prefactor for normalization by evaluating
int^{z=kmax*dz}_0 e^{-z/H} dz
using rectangular rule consistent with the integral evaluation in compute_lwa_and_barotropic_fluxes.f90.
TODO: evaluate numerical integration scheme used in the fortran module.
"""
self.prefactor = sum([math.exp(-k * self.dz / self.scale_height) * self.dz for k in range(1, self.kmax-1)])
return self.prefactor

def _check_valid_plev(self, plev, scale_height, kmax, dz):
"""
Private function. Check the validity of plev to see
Expand Down
104 changes: 16 additions & 88 deletions tests/test_oopinterface.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,26 +28,11 @@ def test_qgfield():
# Create a QGField object out of some masked array for testing. The test below is to ensure
# warning is raised for masked array.
with pytest.warns(UserWarning):
qgfield = QGField(
xlon=xlon,
ylat=np.ma.masked_equal(ylat, 0.0),
plev=plev,
u_field=np.ma.masked_equal(u_field, 0.0),
v_field=np.ma.masked_equal(v_field, 0.0),
t_field=np.ma.masked_equal(t_field, 0.0),
kmax=kmax,
maxit=100000,
dz=1000.,
prefactor=6500.,
npart=None,
tol=1.e-5,
rjac=0.95,
scale_height=SCALE_HEIGHT,
cp=CP,
dry_gas_constant=DRY_GAS_CONSTANT,
omega=EARTH_OMEGA,
planet_radius=EARTH_RADIUS
)
qgfield = QGField(xlon=xlon, ylat=np.ma.masked_equal(ylat, 0.0), plev=plev,
u_field=np.ma.masked_equal(u_field, 0.0), v_field=np.ma.masked_equal(v_field, 0.0),
t_field=np.ma.masked_equal(t_field, 0.0), kmax=kmax, maxit=100000, dz=1000., npart=None,
tol=1.e-5, rjac=0.95, scale_height=SCALE_HEIGHT, cp=CP, dry_gas_constant=DRY_GAS_CONSTANT,
omega=EARTH_OMEGA, planet_radius=EARTH_RADIUS)

# Check that the input fields are interpolated onto a grid of correct dimension
# and the interpolated values are bounded.
Expand Down Expand Up @@ -105,26 +90,9 @@ def test_qgfield():
def test_qgfield_full_globe():

# Create a QGField object for testing
qgfield = QGField(
xlon=xlon,
ylat=ylat,
plev=plev,
u_field=u_field,
v_field=v_field,
t_field=t_field,
kmax=kmax,
maxit=100000,
dz=1000.,
prefactor=6500.,
npart=None,
tol=1.e-5,
rjac=0.95,
scale_height=SCALE_HEIGHT,
cp=CP,
dry_gas_constant=DRY_GAS_CONSTANT,
omega=EARTH_OMEGA,
planet_radius=EARTH_RADIUS
)
qgfield = QGField(xlon=xlon, ylat=ylat, plev=plev, u_field=u_field, v_field=v_field, t_field=t_field, kmax=kmax,
maxit=100000, dz=1000., npart=None, tol=1.e-5, rjac=0.95, scale_height=SCALE_HEIGHT, cp=CP,
dry_gas_constant=DRY_GAS_CONSTANT, omega=EARTH_OMEGA, planet_radius=EARTH_RADIUS)

# Check that the input fields are interpolated onto a grid of correct dimension
# and the interpolated values are bounded.
Expand Down Expand Up @@ -161,26 +129,9 @@ def test_raise_error_for_unrealistic_fields():
"""
Error shall be raised if the SOR algorithm for computing reference state does not converge.
"""
qgfield = QGField(
xlon=xlon,
ylat=ylat,
plev=plev,
u_field=u_field,
v_field=u_field,
t_field=u_field,
kmax=kmax,
maxit=100000,
dz=1000.,
prefactor=6500.,
npart=None,
tol=1.e-5,
rjac=0.95,
scale_height=SCALE_HEIGHT,
cp=CP,
dry_gas_constant=DRY_GAS_CONSTANT,
omega=EARTH_OMEGA,
planet_radius=EARTH_RADIUS
)
qgfield = QGField(xlon=xlon, ylat=ylat, plev=plev, u_field=u_field, v_field=u_field, t_field=u_field, kmax=kmax,
maxit=100000, dz=1000., npart=None, tol=1.e-5, rjac=0.95, scale_height=SCALE_HEIGHT, cp=CP,
dry_gas_constant=DRY_GAS_CONSTANT, omega=EARTH_OMEGA, planet_radius=EARTH_RADIUS)
qgfield.interpolate_fields()
with pytest.raises(ValueError):
qgfield.compute_reference_states()
Expand All @@ -194,14 +145,7 @@ def test_raise_error_for_unrealistic_kmax():
"""
too_large_kmax = 1000
with pytest.raises(ValueError):
QGField(
xlon=xlon,
ylat=ylat,
plev=plev,
u_field=u_field,
v_field=v_field,
t_field=t_field,
kmax=too_large_kmax)
QGField(xlon=xlon, ylat=ylat, plev=plev, u_field=u_field, v_field=v_field, t_field=t_field, kmax=too_large_kmax)


def test_interpolate_fields_even_grids():
Expand All @@ -220,26 +164,10 @@ def test_interpolate_fields_even_grids():
fill_value="extrapolate")(ylat_even)

# Create a QGField object for testing
qgfield_even = QGField(
xlon=xlon,
ylat=ylat_even,
plev=plev,
u_field=u_field_even,
v_field=v_field_even,
t_field=t_field_even,
kmax=kmax,
maxit=100000,
dz=1000.,
prefactor=6500.,
npart=None,
tol=1.e-5,
rjac=0.95,
scale_height=SCALE_HEIGHT,
cp=CP,
dry_gas_constant=DRY_GAS_CONSTANT,
omega=EARTH_OMEGA,
planet_radius=EARTH_RADIUS
)
qgfield_even = QGField(xlon=xlon, ylat=ylat_even, plev=plev, u_field=u_field_even, v_field=v_field_even,
t_field=t_field_even, kmax=kmax, maxit=100000, dz=1000., npart=None, tol=1.e-5, rjac=0.95,
scale_height=SCALE_HEIGHT, cp=CP, dry_gas_constant=DRY_GAS_CONSTANT, omega=EARTH_OMEGA,
planet_radius=EARTH_RADIUS)

qgpv, interpolated_u, interpolated_v, interpolated_theta, static_stability = \
qgfield_even.interpolate_fields()
Expand Down

0 comments on commit 2de8e28

Please sign in to comment.