The atomization process of turbulent liquid jets is not well understood. Detailed numerical simulations can help study the fundamental mechanisms in regions where experimental access and analysis is virtually impossible. However, simulating atomization accurately is a huge numerical challenge since time and length scales vary over several orders of magnitude, the phase interface is a material discontinuity, and surface tension forces are singular. The Refined Level Set Grid (RLSG) method presented here is one numerical approach to simulate the primary breakup process of liquid jets and sheets in detail. With it, the liquid/gas phase interface is tracked by a level set method using an auxiliary high resolution equidistant Cartesian grid. This not only allows for application of higher-order WENO schemes retaining their full order of accuracy both for advecting and reinitializing the level set scalar, but it also provides the necessary high resolution of the phase interface geometry during topology change events in an efficient manner. When atomization occurs in non-isothermal environments, such as in combustion devices, thermal fluctuations can be significant on length scales associated with the liquid atomization process. Since the surface tension force is a function of local temperature, these thermal fluctuations may result in large local variations of the surface tension force, thereby potentially impacting the atomization process. Here, we present a numerical technique to incorporate these thermal Marangoni forces into the balanced force algorithm, together with verification and validation test cases geared towards testing the applicability of the proposed methods to the case of atomization.