This is a very nice paper describing the UBER 1 code, its potential applications and in particular, the advantages of the stochastic approach to solving convection-diffusion equations including mixed diffusion terms. I have mostly minor comments that should be easy for authors to address. I would also recommend adding an example that would show how the code behaves in case of strong gradients and not smooth initial distributions.

• I would clarify around line 20 that Js are the canonical momentums; otherwise alpha would go from 1 to 6.
• I thought that we neglect the phases of J simply because we assume isotropy in all phis . I could not precisely follow the argument of electric potential energy. Please elaborate.
• While the paper discusses in detail the advantages of this approach, it omits the disadvantages. I would suggest clearly saying how long the presented runs take to run on a regular PC and on a supercomputer. I would also specify how much wall clock computing would take to calculate one day charging along a given satellite orbit. Also, specify how much slower this approach than more traditional approaches for the multi-dimensional diffusion equations.
• The presented tests all show examples that are initiated with smooth initial conditions. Please provide examples similar to Aseev et al., 2016 with strong gradients in initial conditions.

Aseev N. A., Y. Y. Shprits, A. Y. Drozdov, A. C. Kellerman, (2016), Numerical applications of the advective-diffusive codes for the inner magnetosphere, Space Weather, 14, 993-1010, doi:10.1002/2016SW001484

This paper presents a very nice universal solver for transport equations with examples in various situations, utilizing the SDE method.  It is always happy to see that significant improvements have been made to the SDE method in terms of speed and handling of boundary conditions, since its first adoption to the radiation belt problem in 2008. The SDE method clearly has its advantage in terms of handling cross-diffusion terms and efficiency the case when only a few solutions in (t,x) are needed. In principle, I think this paper is publishable. I only have a few minor comments and suggestions for the authors.

1, The main problem I have with the paper is that it failed to discuss the main concern of the  the SDE method: its speed, especially compared with finite difference methods etc.  Clearly, the SDE code with newer techniques is much much more efficient than the traditional one. But how does it compare with simple finite difference methods in the three cases discussed?  I think the authors could discuss the comparison in different situations; e.g., when one cares about only a few snapshots, and when needs to know the whole history of evolution.  Having a fair discussion about the disadvantage of the method does not make the paper less significant, but instead show future directions where improvements can be made.

2,  Lines 45-50:  These listed numerical codes use either finite difference type methods or SDE/layer methods.  None of the methods are limited by the choice of coordinates. Those authors chose a particular choice of coordinates probably because they did not intend to build a general library or because they tended to demonstrate a new method. I think it would nice for the authors to take this into consideration when discussing previous models.

3, Lines 53-55: The complicated geometry is a problem for some of the models, mainly because radiation belt people seem to have a preference for finite difference methods. There exist general powerful finite volume/elements methods that can handle complicated boundary geometry.  So here “… would be challending for numerical methods” should really be “… would be challenging for finite difference methods.”  

4,  In Table 1, the UBER library input,  can the current code handle time-dependent D? I know the method can, but not sure if it has been implemented by the code.

5,  Lines 158: No, the grids in the layer method need NOT to be uniform.  It was simply chosen for simplicity for estimating error and demonstration purposes.  

6, Lines 155-159: Yes, in SDE methods, one can design this kind of method to obtain the global solutions of f and its history. However, the key to implementing this in SDE is actually about finding an appropriate interpolation method.  “The only operation … would be interpolation” sounds like finding such an interpolation method is easy, while in fact it is NOT.  For the described choice of irregular domain or nonuniform domain, the interpolation method still needs to be of high order to reduce systematic error from interpolation, and to preserve positivity of the solution.  For example, if a simple 4^th order polynomial interpolation method is used, one might introduce oscillation of f, and hence negative value s of f,  from interpolation, and hence violates one of the key advantage of the SDE method.That is why Tao et al., 2016 (doi:10.1002/2015JA022064) introduced those higher-order positivity-preserving methods to be used with layer method. 

7:  Lines 310-315: There actually exists flux-limited Lax-Wendroff methods that can avoid introducing unphysical negative solutions.  And in the case of forming steep gradients, it is well-known that one should add flux-limiters to Lax-Wendroff type methods.

Overall, I think this is a nice paper, and it would be an important contribution to the community of radiation belt modeling. I simply hope that the authors can give a more balanced description of the UBER library, especially when comparing with other methods.

The paper presents a new numerical solver, UBER, which can solve the general form of radiation belt Fokker-Planck equation and Boltzmann equation in arbitrarily provided coordinate systems, and with user-specified boundary geometry, boundary conditions, and equation terms. The mathematical theory and numerical techniques of UBER are clearly described, and three example problems are presented which well demonstrates the capability and reliability of the solver. Overall, this is a nicely written paper with convincing results showing that UBER could make important contribution to radiation belt modeling. The reviewer only has a few minor comments before recommending the paper for publication in GMD.

1. Could the authors clarify or discuss further if the UBER solver can handle boundary conditions at varying locations or not? If we look at the radial diffusion model as a simple example, in the traditional solvers based on finite different method, it is often driven by data-driven outer boundary conditions at a fixed L*. However, the L* location of the satellite data providing the outer boundary is actually varying in time, which could lead to data gaps in the outer boundary condition and uncertainties in the model results. Discussions on if and how this type of boundary conditions can be implemented in UBER will further increase the significance of the work.
2. The reviewer is also curious how the coordinate conversion between adiabatic invariants and energy and pitch angle can be implemented in UBER. As discussed in Subbotin and Shprits [JGR, 2012], several 3D diffusion radiation belt models “utilized two grids to solve the Fokker-Planck equation; one grid, which keeps the first and second adiabatic invariants constant, was used for the computation of the radial diffusion, and the other grid, orthogonal in energy and pitch angle at each fixed radial distance, was used for the computation of energy diffusion, pitch angle diffusion, and mixed energy and pitch angle diffusion.” At each time step, the results were converted and interpolated between the two grids, which could lead to uncertainties in the model results depending on how the conversion and interpolation are performed. Is this coordinate conversion still needed in UBER? If not, why?
3. It will greatly enhance the significance of the paper if discussions are included in comparing the efficiency, stability, and accuracy between the UBER code and traditional finite difference codes. It is demonstrated the UBER is more efficient than previous SDE codes, but how does it compare with the traditional solvers based on the finite difference method? Is it still much less efficient if global distributions of radiation belt electrons in L, energy, and pitch angle are targeted? How about if we only need to solve for the distribution locally at certain L, energy, and pitch angle? The authors could perhaps use example problem 3 to compare the efficiency, stability, and the accuracy between the two different types of solvers and then expand the discussion to higher-dimension models such as 3D diffusion models.