I'm trying to create a partial differential equation that approximates 1-D climate in a rocky planet's atmosphere, which accounts for energy transport via radiation and convection. I am only interested in steady-states. I've come up with a PDE which should have a steady state that satisfies my interests, but I'm having a hard time solving it numerically. Consider the conservation equation, where the conserved quantity is energy density, $E$ (with CGS units ergs cm$^3$)
$$\frac{\partial E}{\partial t} = \sigma_r - \frac{\partial F_c}{\partial z}$$
Here, $t$ is time and $z$ is altitude. $\sigma_r$ represents the energy source/sink from how radiation travels through the atmosphere having a complex dependence on temperature ($T$). $F_c$ is the convective flux of energy. Note that energy density can be related to temperature via,
$$\frac{\partial (\rho c_p T)}{\partial t} = \frac{\partial E}{\partial t}$$
The convective flux has the following form
$$F_c = h \sigma_r \frac{\partial T}{\partial z}$$
Here, $h$ is a sigmoid function that depends strongly on $\frac{\partial T}{\partial z}$. $h$ is very big when the gradient ($\frac{\partial T}{\partial z}$) is unstable to convection, and $h$ is very small when the gradient is convectively stable. Overall, the idea is that the term $\frac{\partial F_c}{\partial z}$ will diffusively transport energy in the atmosphere to account for any convective instability that may arise.
In a finite volume $i$, the convection can be approximated as
$$\frac{\partial F_c}{\partial z} |^i \approx \frac{F_{c}^{i+1/2} - F_{c}^{i-1/2}}{\Delta z^i}$$
I have tried a center difference scheme:
$$F_{c}^{i+1/2} = h^{i+1/2} \sigma_r^{i+1/2} \frac{T^{i+1} - T^{i}}{\Delta z^{i+1/2}}$$
Here, $h^{i+1/2}$ depends on the approximated derivative $\frac{T^{i+1} - T^{i}}{\Delta z^{i+1/2}}$. I then solve for the steady state temperature structure using the hybrj method in MINPACK. The approach seems to work reasonably in some cases when there isn't much "convection" in the atmosphere, but it fails for atmospheres with lots of "convection" (i.e. when $\frac{\partial F_c}{\partial z}$ becomes big over a lot of the atmosphere).
I am wondering if my center differencing scheme could be a problem? What scheme might be more reliable?
Thanks!
Edit:
I've been doing some more experimenting. I've created a new formulation for the convective flux
$$F_c = \epsilon T$$
where $\epsilon$ is given by
$$\log_{10}\left(\frac{\epsilon}{10^{15}}\right) = -\frac{1}{2} + \frac{1}{\pi}\arctan\left(\frac{-\frac{dT}{dz} \frac{1}{\gamma} - 1}{0.05}\right)$$
$\gamma$ is the critical lapse rate. If, $-\frac{dT}{dz} > \frac{1}{\gamma}$ then, $\epsilon$ is big, and if not, then $\epsilon$ is small. In summary, $\epsilon$ can is an approximation of a discontinuity. The values $10^{15}$ and $0.05$ are tunable parameters, that change how the discontinuity is approximated.
When I re-solve the problem with the $F_c$ given above everything works as I want it to. I use 1st order upwind differencing, because there is a wave traveling in the positive $z$ direction.
