diff --git a/AMSS_NCKU_source/bssn_rhs_c.C b/AMSS_NCKU_source/bssn_rhs_c.C
index ba0f615..9bba3ca 100644
--- a/AMSS_NCKU_source/bssn_rhs_c.C
+++ b/AMSS_NCKU_source/bssn_rhs_c.C
@@ -50,7 +50,7 @@ int f_compute_rhs_bssn(int *ex, double &T,
     double Gamxx[all],Gamxy[all],Gamxz[all];
     double Gamyx[all],Gamyy[all],Gamyz[all];
     double Gamzx[all],Gamzy[all],Gamzz[all];
-    double Kx[all], Ky[all], Kz[all], div_beta[all], S[all];
+    double Kx[all], Ky[all], Kz[all], S[all];
     double f[all], fxx[all], fxy[all], fxz[all], fyy[all], fyz[all], fzz[all];
     double Gamxa[all], Gamya[all], Gamza[all], alpn1[all], chin1[all];
     double gupxx[all], gupxy[all], gupxz[all];
@@ -123,21 +123,21 @@ int f_compute_rhs_bssn(int *ex, double &T,
 
         // 3ms //
         for(int i=0;i<all;i+=1){
-            div_beta[i] = betaxx[i] + betayy[i] + betazz[i];
-            chi_rhs[i] = F2o3 * chin1[i] * (alpn1[i] * trK[i] - div_beta[i]); 
-            gxx_rhs[i] = -TWO * alpn1[i] * Axx[i] - F2o3 * (dxx[i] + ONE) * div_beta[i] + 
+            const double divb = betaxx[i] + betayy[i] + betazz[i];
+            chi_rhs[i] = F2o3 * chin1[i] * (alpn1[i] * trK[i] - divb); 
+            gxx_rhs[i] = -TWO * alpn1[i] * Axx[i] - F2o3 * (dxx[i] + ONE) * divb + 
                     TWO * ((dxx[i] + ONE) * betaxx[i] + gxy[i] * betayx[i] + gxz[i] * betazx[i]);
-            gyy_rhs[i] = -TWO * alpn1[i] * Ayy[i] - F2o3 * (dyy[i] + ONE) * div_beta[i] + 
+            gyy_rhs[i] = -TWO * alpn1[i] * Ayy[i] - F2o3 * (dyy[i] + ONE) * divb + 
                     TWO * (gxy[i] * betaxy[i] + (dyy[i] + ONE) * betayy[i] + gyz[i] * betazy[i]);
-            gzz_rhs[i] = -TWO * alpn1[i] * Azz[i] - F2o3 * (dzz[i] + ONE) * div_beta[i] + 
+            gzz_rhs[i] = -TWO * alpn1[i] * Azz[i] - F2o3 * (dzz[i] + ONE) * divb + 
                         TWO * (gxz[i] * betaxz[i] + gyz[i] * betayz[i] + (dzz[i] + ONE) * betazz[i]);
-            gxy_rhs[i] = -TWO * alpn1[i] * Axy[i] + F1o3 * gxy[i] * div_beta[i] + 
+            gxy_rhs[i] = -TWO * alpn1[i] * Axy[i] + F1o3 * gxy[i] * divb + 
                         (dxx[i] + ONE) * betaxy[i] + gxz[i] * betazy[i] + (dyy[i] + ONE) * betayx[i] 
                         + gyz[i] * betazx[i] - gxy[i] * betazz[i];        
-            gyz_rhs[i] = -TWO * alpn1[i] * Ayz[i] + F1o3 * gyz[i] * div_beta[i] + 
+            gyz_rhs[i] = -TWO * alpn1[i] * Ayz[i] + F1o3 * gyz[i] * divb + 
                         gxy[i] * betaxz[i] + (dyy[i] + ONE) * betayz[i] + gxz[i] * betaxy[i] 
                         + (dzz[i] + ONE) * betazy[i] - gyz[i] * betaxx[i];
-            gxz_rhs[i] = -TWO * alpn1[i] * Axz[i] + F1o3 * gxz[i] * div_beta[i] + 
+            gxz_rhs[i] = -TWO * alpn1[i] * Axz[i] + F1o3 * gxz[i] * divb + 
                         (dxx[i] + ONE) * betaxz[i] + gxy[i] * betayz[i] + gyz[i] * betayx[i] 
                         + (dzz[i] + ONE) * betazx[i] - gxz[i] * betayy[i];                  
         }
@@ -328,6 +328,7 @@ int f_compute_rhs_bssn(int *ex, double &T,
 
         // Fused: fxx/Gamxa + Gamma_rhs part 2 (2 loops -> 1)
         for(int i=0;i<all;i+=1){
+            const double divb = betaxx[i] + betayy[i] + betazz[i];
             double lfxx = gxxx[i] + gxyy[i] + gxzz[i];
             double lfxy = gxyx[i] + gyyy[i] + gyzz[i];
             double lfxz = gxzx[i] + gyzy[i] + gzzz[i];
@@ -342,21 +343,21 @@ int f_compute_rhs_bssn(int *ex, double &T,
             Gamxa[i] = gxa; Gamya[i] = gya; Gamza[i] = gza;
 
             Gamx_rhs[i] = Gamx_rhs[i]
-                        + F2o3 * gxa * div_beta[i]
+                        + F2o3 * gxa * divb
                         - gxa * betaxx[i] - gya * betaxy[i] - gza * betaxz[i]
                         + F1o3 * ( gupxx[i] * lfxx + gupxy[i] * lfxy + gupxz[i] * lfxz )
                         + gupxx[i] * gxxx[i] + gupyy[i] * gyyx[i] + gupzz[i] * gzzx[i]
                         + TWO * ( gupxy[i] * gxyx[i] + gupxz[i] * gxzx[i] + gupyz[i] * gyzx[i] );
 
             Gamy_rhs[i] = Gamy_rhs[i]
-                        + F2o3 * gya * div_beta[i]
+                        + F2o3 * gya * divb
                         - gxa * betayx[i] - gya * betayy[i] - gza * betayz[i]
                         + F1o3 * ( gupxy[i] * lfxx + gupyy[i] * lfxy + gupyz[i] * lfxz )
                         + gupxx[i] * gxxy[i] + gupyy[i] * gyyy[i] + gupzz[i] * gzzy[i]
                         + TWO * ( gupxy[i] * gxyy[i] + gupxz[i] * gxzy[i] + gupyz[i] * gyzy[i] );
 
             Gamz_rhs[i] = Gamz_rhs[i]
-                        + F2o3 * gza * div_beta[i]
+                        + F2o3 * gza * divb
                         - gxa * betazx[i] - gya * betazy[i] - gza * betazz[i]
                         + F1o3 * ( gupxz[i] * lfxx + gupyz[i] * lfxy + gupzz[i] * lfxz )
                         + gupxx[i] * gxxz[i] + gupyy[i] * gyyz[i] + gupzz[i] * gzzz[i]
@@ -758,6 +759,7 @@ int f_compute_rhs_bssn(int *ex, double &T,
         }
         // 2.5ms //
         for (int i=0;i<all;i+=1) {
+            const double divb = betaxx[i] + betayy[i] + betazz[i];
 
             S[i] = chin1[i] * (
                     gupxx[i] * Sxx[i] + gupyy[i] * Syy[i] + gupzz[i] * Szz[i]
@@ -883,19 +885,19 @@ int f_compute_rhs_bssn(int *ex, double &T,
                 f[i] * Axx_rhs[i]
                 + alpn1[i] * ( trK[i] * Axx[i] - TWO * fxx[i] )
                 + TWO * ( Axx[i] * betaxx[i] + Axy[i] * betayx[i] + Axz[i] * betazx[i] )
-                - F2o3 * Axx[i] * div_beta[i];
+                - F2o3 * Axx[i] * divb;
 
             Ayy_rhs[i] =
                 f[i] * Ayy_rhs[i]
                 + alpn1[i] * ( trK[i] * Ayy[i] - TWO * fyy[i] )
                 + TWO * ( Axy[i] * betaxy[i] + Ayy[i] * betayy[i] + Ayz[i] * betazy[i] )
-                - F2o3 * Ayy[i] * div_beta[i];
+                - F2o3 * Ayy[i] * divb;
 
             Azz_rhs[i] =
                 f[i] * Azz_rhs[i]
                 + alpn1[i] * ( trK[i] * Azz[i] - TWO * fzz[i] )
                 + TWO * ( Axz[i] * betaxz[i] + Ayz[i] * betayz[i] + Azz[i] * betazz[i] )
-                - F2o3 * Azz[i] * div_beta[i];
+                - F2o3 * Azz[i] * divb;
 
             Axy_rhs[i] =
                 f[i] * Axy_rhs[i]
@@ -904,7 +906,7 @@ int f_compute_rhs_bssn(int *ex, double &T,
                 + Axz[i] * betazy[i]
                 + Ayy[i] * betayx[i]
                 + Ayz[i] * betazx[i]
-                + F1o3 * Axy[i] * div_beta[i]
+                + F1o3 * Axy[i] * divb
                 - Axy[i] * betazz[i];
 
             Ayz_rhs[i] =
@@ -914,7 +916,7 @@ int f_compute_rhs_bssn(int *ex, double &T,
                 + Ayy[i] * betayz[i]
                 + Axz[i] * betaxy[i]
                 + Azz[i] * betazy[i]
-                + F1o3 * Ayz[i] * div_beta[i]
+                + F1o3 * Ayz[i] * divb
                 - Ayz[i] * betaxx[i];
 
             Axz_rhs[i] =
@@ -924,7 +926,7 @@ int f_compute_rhs_bssn(int *ex, double &T,
                 + Axy[i] * betayz[i]
                 + Ayz[i] * betayx[i]
                 + Azz[i] * betazx[i]
-                + F1o3 * Axz[i] * div_beta[i]
+                + F1o3 * Axz[i] * divb
                 - Axz[i] * betayy[i];
 
             /* Compute trace of S_ij */