Segmentation Fault when using EmanuelConvection

[HenryDane]

While working with climlab for a class, I noticed some odd behavior when the EmanuelConvection process is used in a system with DailyInsolation and that the EmanuelConvection code seems to trigger a segmentation fault.

I'm using climlab version 0.8.1 and I'm on Ubuntu 20.04 (64-bit). This problem also occurred on the Jupyter hub I was using which is a CentOS Linux 8 system.

I've attached the full code needed to reproduce this crash below in repro.txt (GitHub won't let me upload a .py file).

The segmentation fault is usually preceded by several instances of a RuntimeWarning related to a divide-by-zero or an overflow. Looking through the source code, it seemed very likely (at least to me) that the issue was probably coming from the convect routine but I am less sure of it now. I am not sure how else python could be triggering a segmentation fault aside from the Fortran code but, again, I am not confident.

I'm opening this issue less because of the specifics of the problem (I may have set something up wrong with the model, although I doubt it) and more to draw attention to the fact that a numeric instability is causing a segmentation fault (which seems very odd).

I spoke to @nfeldl about this in class and she recommended that I open an issue here. I'm happy to help fix this and I've already begun translating the Fortran code (convect, driver, etc) into Python/numpy code and I'm happy to share it with you if that would be useful.

[brian-rose]

Thanks @HenryDane for this report!

Copying and pasting your code here for reference:

import climlab
from climlab import constants as const
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

# Slightly modified example from https://climlab.readthedocs.io/en/latest/api/climlab.convection.EmanuelConvection.html

# State information
full_state = climlab.column_state(num_lev=40, num_lat=40, water_depth=2.5)
temperature_state = {'Tatm':full_state.Tatm,'Ts':full_state.Ts}
q = np.ones_like(full_state.Tatm) * 5.E-6
full_state['q'] = q

#  The top-level model
model = climlab.TimeDependentProcess(state=full_state,
                              timestep=const.seconds_per_hour)
#  Daily insolation as a function of latitude and time of year
sun = climlab.radiation.DailyInsolation(name='Insolation', domains=temperature_state['Ts'].domain) # changing domain to Tatm causes ValueError (failed in converting ... coszen ... to C/Fortran array)
#  Radiation coupled to water vapor
rad = climlab.radiation.RRTMG(state=temperature_state,
                              specific_humidity=full_state.q,
                              albedo=0.3,
                              timestep=const.seconds_per_day,
                              insolation=sun.insolation,
                              coszen=sun.coszen
                              )
#  Convection scheme -- water vapor is a state variable
conv = climlab.convection.EmanuelConvection(state=full_state,
                              timestep=const.seconds_per_hour)
#  Surface heat flux processes
shf = climlab.surface.SensibleHeatFlux(state=temperature_state, Cd=0.5E-3,
                              timestep=const.seconds_per_hour)
lhf = climlab.surface.LatentHeatFlux(state=full_state, Cd=0.5E-3,
                              timestep=const.seconds_per_hour)
#  Couple all the submodels together
model.add_subprocess('Radiation', rad)
model.add_subprocess('Convection', conv)
model.add_subprocess('SHF', shf)
model.add_subprocess('LHF', lhf)
model.add_subprocess('Sun', sun)

# Debug model info
print(model)

model.integrate_years(1)

plt.plot(model.lat, model.precipitation)
plt.xlabel('Lat')
plt.ylabel('Precip')
plt.show()