Guard THINC density ratio computation against divide-by-zero

When advxb approaches 0 or 1 at an interface, rho_b and rho_e
could divide by zero. Guard with sgm_eps (1e-16) using the same
pattern already used for C = (aC - qmin + sgm_eps)/(qmax + sgm_eps).

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: sbryngelson

src/simulation/m_muscl.fpp

@@ -280,8 +280,10 @@ contains
                                 A = (B/cosh(ic_beta) - 1._wp)/tanh(ic_beta)
 
                                 ! Save original density ratios before THINC overwrites them
-                                rho_b = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/vL_rs_vf_${XYZ}$ (j, k, l, advxb)
-                                rho_e = vL_rs_vf_${XYZ}$ (j, k, l, contxe)/(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb))
+                                rho_b = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/ &
+                                        max(vL_rs_vf_${XYZ}$ (j, k, l, advxb), sgm_eps)
+                                rho_e = vL_rs_vf_${XYZ}$ (j, k, l, contxe)/ &
+                                        max(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb), sgm_eps)
 
                                 ! Left reconstruction
                                 aTHINC = qmin + 5e-1_wp*qmax*(1._wp + sign*A)