diff --git a/AMSS_NCKU_source/bssn_rhs_opt.f90 b/AMSS_NCKU_source/bssn_rhs_opt.f90
index e82239a..71e0c3b 100644
--- a/AMSS_NCKU_source/bssn_rhs_opt.f90
+++ b/AMSS_NCKU_source/bssn_rhs_opt.f90
@@ -88,10 +88,30 @@ contains
     idy = ONE/dY_val
     idz = ONE/dZ_val
 
+    ! Initialize outputs to zero
+    chi_rhs = ZEO; trK_rhs = ZEO
+    gxx_rhs = ZEO; gxy_rhs = ZEO; gxz_rhs = ZEO
+    gyy_rhs = ZEO; gyz_rhs = ZEO; gzz_rhs = ZEO
+    Axx_rhs = ZEO; Axy_rhs = ZEO; Axz_rhs = ZEO
+    Ayy_rhs = ZEO; Ayz_rhs = ZEO; Azz_rhs = ZEO
+    Gamx_rhs = ZEO; Gamy_rhs = ZEO; Gamz_rhs = ZEO
+    Lap_rhs = ZEO; betax_rhs = ZEO; betay_rhs = ZEO; betaz_rhs = ZEO
+    dtSfx_rhs = ZEO; dtSfy_rhs = ZEO; dtSfz_rhs = ZEO
+    Gamxxx = ZEO; Gamxxy = ZEO; Gamxxz = ZEO
+    Gamxyy = ZEO; Gamxyz = ZEO; Gamxzz = ZEO
+    Gamyxx = ZEO; Gamyxy = ZEO; Gamyxz = ZEO
+    Gamyyy = ZEO; Gamyyz = ZEO; Gamyzz = ZEO
+    Gamzxx = ZEO; Gamzxy = ZEO; Gamzxz = ZEO
+    Gamzyy = ZEO; Gamzyz = ZEO; Gamzzz = ZEO
+    Rxx = ZEO; Rxy = ZEO; Rxz = ZEO
+    Ryy = ZEO; Ryz = ZEO; Rzz = ZEO
+
     ! Optimization 1: NaN Check
+#define DEBUG_NAN_CHECK
 #ifdef DEBUG_NAN_CHECK
     dX_check = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz)
     if(dX_check.ne.dX_check) then
+       write(*,*) "bssn_rhs_opt: NaN detected in input variables!"
        gont = 1
        return
     endif
@@ -152,7 +172,7 @@ contains
       real*8 :: gxy2, gxz2, gyz2  ! For optimized matrix inversion
       real*8 :: term_xx1, term_xx2, term_yy1, term_yy2, term_zz1, term_zz2  ! For Christoffel CSE
       real*8 :: term_xy, term_xz, term_yz
-      real*8 :: reta_val
+      real*8 :: reta_val, denom
       
       nx = ie - is + 1
       ny = je - js + 1
@@ -642,7 +662,14 @@ contains
          ! reta calculation using chix... computed in Step 1
          reta_val = gupxx(ii,jj,kk)*chix(ii,jj,kk)**2 + gupyy(ii,jj,kk)*chiy(ii,jj,kk)**2 + gupzz(ii,jj,kk)*chiz(ii,jj,kk)**2 + &
                     TWO*(gupxy(ii,jj,kk)*chix(ii,jj,kk)*chiy(ii,jj,kk) + gupxz(ii,jj,kk)*chix(ii,jj,kk)*chiz(ii,jj,kk) + gupyz(ii,jj,kk)*chiy(ii,jj,kk)*chiz(ii,jj,kk))
-         reta_val = 1.31d0/TWO * sqrt(reta_val/chin1) / (ONE - sqrt(chin1))**2
+         
+         if (chin1 > 1.0d-10) then
+             denom = ONE - sqrt(chin1)
+             if (abs(denom) < 1.0d-15) denom = 1.0d-15
+             reta_val = 1.31d0/TWO * sqrt(max(reta_val, ZEO)/chin1) / denom**2
+         else
+             reta_val = ZEO
+         endif
          
          dtSfx_rhs(ig,jg,kg) = Gamx_rhs(ig,jg,kg) - reta_val*dtSfx(ig,jg,kg)
          dtSfy_rhs(ig,jg,kg) = Gamy_rhs(ig,jg,kg) - reta_val*dtSfy(ig,jg,kg)