lsoda when setting up parcel_profile function for cape_cin function

[jwynnsmith]

Hello

I am attempting to use the metpy cape_cin function but I cannot get past an error with the parcel_profile function. Here is the error below:

lsoda-- at t (=r1), too much accuracy requested
for precision of machine.. see tolsf (=r2)
in above, r1 = 0.3073318824321D+02 r2 = NaN

I set up the necessary variables and dataframe for the parcel_profile function as follows

import xarray as xr
import numpy as np
from numpy import genfromtxt
from metpy.calc import cape_cin, dewpoint_from_relative_humidity, parcel_profile, surface_based_cape_cin, lcl, add_pressure_to_height, get_layer_heights, pressure_to_height_std
import pandas as pd
from metpy.units import units, pandas_dataframe_to_unit_arrays
from metpy.interpolate import log_interpolate_1d
import math
nc_file = xr.open_dataset('/work/Jonathan.Smith/ENV_VAR/MERRA2_400.inst3_3d_asm_Np.20190629.nc4')
press = nc_file.lev.values
print(press)

temp = nc_file.T.values - 273.15
temp = temp[0,:,300,400]
print(temp)

rel_hum = nc_file.RH.values
rel_hum = rel_hum[0,:,300,400]
print(rel_hum)

Dewpoint = dewpoint_from_relative_humidity(temp * units.degC, rel_hum * units.percent)
print(Dewpoint)

Td0 = Dewpoint[0] * units.degC
print('Surface Td', Td0)
print(np.shape(Td0))For the surface_based_cape_cin function

data = {'P': press,
'T': temp,
'RH': rel_hum,
'Td': Dewpoint}

df = pd.DataFrame(data)
print(df)

profile = parcel_profile(df['P'].values * units.hPa, df['T'][0] * units.degC , df['Td'][0] * units.degC )
print(profile)

Am I applying the units correctly or not setting up the function properly to get the above error?

Thanks
Jonathan

[dopplershift]

So on the current version of MetPy, that output is accompanied by:

ValueError: ODE Integration failed. This is likely due to trying to calculate at too small values of pressure.

Essentially the ODE integrator we use to calculate the moist adiabats doesn't like that you're trying to calculate down to 0.1 hPa.

You're also doing a lot of fighting with the units above. The original netCDF file has unit metadata, so let's just let things flow naturally. This is how I'd write that code above, which slices off the last 15 values of the pressure array so we're only dealing with data up to 50 hPa:

import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr

nc_file = xr.open_dataset('/Users/rmay/Downloads/MERRA2_400.inst3_3d_asm_Np.20190629.nc4')

# Create a slice from 0 to 15 indices from the end of the array
p_slice = slice(0, -15)

# Use this slice on all the arrays to keep them in agreement
press = nc_file.lev[p_slice]
temp = nc_file.T[0, p_slice, 300, 400]
rel_hum = nc_file.RH[0, p_slice, 300, 400]

dewpoint = mpcalc.dewpoint_from_relative_humidity(temp, rel_hum)
Td0 = dewpoint[0]

profile = mpcalc.parcel_profile(press, temp[0], Td0)
sb_cape, sb_cin = mpcalc.surface_based_cape_cin(press, temp, dewpoint)

[jwynnsmith]

Thank you for this solution. I has some help to come to this same conclusion.