I've spent the last week doing radiative transfer calculations with Molpop-CEP in JPL with Moshe Elitzur, Paul Goldsmith and Karen Willacy in photodissociation regions (PDR). I have managed to compute CO, H2O, OH, C+, C and O2. The only molecule that really gave us some problems was water. We found that, in some models, the temperature goes down to 10 K, while part of the PDR is still with temperature above 350 K. This poses a problem to any radiative transfer method that solves the full non-local problem. The reason is that in some parts of the atmosphere we have many levels properly pumped but the upper energy levels are completely unpopulated in the low temperature regions. Then, the convergence is dominated by these levels and the iteration never converges.

We have an idea to update the code to deal with an adaptive control of the number of levels in each zone. The idea is to delete energy levels whose population is below a certain threshold and only solve for the remaining low-energy levels. This would produce that the mean intensity in these regions is not affected. If these region happen close to the boundary, they will simply see the boundary condition. If we have accumulated mean intensity in the radiative transfer calculation, these regions where the line has disappeared will be plainly transparent. So, the mean intensity will be still computed for all transitions, but it will be updated only in the regions where both levels of the transition are active. I will think about the practical implementation of this idea and see if it works.