#include "macrodef.fh"

module bssn_rhs_opt_mod
  implicit none
  
  ! Block Size Configuration
  ! 8^3 * 8 bytes * ~50 arrays ~= 200KB, fitting comfortably in L2 Cache
  integer, parameter :: BLK = 8
  
  ! Physics Constants
  real*8, parameter :: ZEO = 0.d0, ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0
  real*8, parameter :: EIGHT = 8.D0, HALF = 0.5D0, THR = 3.d0
  real*8, parameter :: PI = 3.14159265358979323846d0
  real*8, parameter :: F1o3 = 1.D0/3.D0, F2o3 = 2.D0/3.D0, F3o2=1.5d0
  real*8, parameter :: F16=1.6d1, F8=8.d0, F1o12=1.0d0/12.0d0, F1o144=1.0d0/144.0d0
  real*8, parameter :: F1o4=0.25d0
  real*8, parameter :: F30=3.d1, F18=1.8d1, F10=1.d1, F6=6.d0, F3=3.d0
  real*8, parameter :: F64=6.4d1, FIT=1.5d1, TWT=2.d1, SIX=6.d0
  
  real*8, parameter :: FF = 0.75d0, eta=2.d0
  
contains

  ! Main Entry Point for Optimized Computation
  subroutine compute_rhs_bssn_opt(ex, T,X, Y, Z, &
               chi, trK, dxx, gxy, gxz, dyy, gyz, dzz, &
               Axx, Axy, Axz, Ayy, Ayz, Azz, &
               Gamx, Gamy, Gamz, Lap, betax, betay, betaz, &
               dtSfx, dtSfy, dtSfz, &
               chi_rhs, trK_rhs, &
               gxx_rhs, gxy_rhs, gxz_rhs, gyy_rhs, gyz_rhs, gzz_rhs, &
               Axx_rhs, Axy_rhs, Axz_rhs, Ayy_rhs, Ayz_rhs, Azz_rhs, &
               Gamx_rhs, Gamy_rhs, Gamz_rhs, &
               Lap_rhs, betax_rhs, betay_rhs, betaz_rhs, &
               dtSfx_rhs, dtSfy_rhs, dtSfz_rhs, &
               rho, Sx, Sy, Sz, Sxx, Sxy, Sxz, Syy, Syz, Szz, &
               Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz, &
               Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz, &
               Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz, &
               Rxx, Rxy, Rxz, Ryy, Ryz, Rzz, &
               ham_Res, movx_Res, movy_Res, movz_Res, &
               Gmx_Res, Gmy_Res, Gmz_Res, &
               Symmetry, Lev, eps, co, gont)

    implicit none
    ! Arguments
    integer,intent(in ):: ex(1:3), Symmetry,Lev,co
    real*8, intent(in ):: T
    real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3))
    real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: chi,dxx,dyy,dzz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: trK,gxy,gxz,gyz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Axx,Axy,Axz,Ayy,Ayz,Azz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Gamx,Gamy,Gamz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Lap, betax, betay, betaz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: dtSfx,  dtSfy,  dtSfz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: chi_rhs,trK_rhs
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: gxx_rhs,gxy_rhs,gxz_rhs
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: gyy_rhs,gyz_rhs,gzz_rhs
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Axx_rhs,Axy_rhs,Axz_rhs
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Ayy_rhs,Ayz_rhs,Azz_rhs
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamx_rhs,Gamy_rhs,Gamz_rhs
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Lap_rhs, betax_rhs, betay_rhs, betaz_rhs
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: dtSfx_rhs,dtSfy_rhs,dtSfz_rhs
    real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: rho,Sx,Sy,Sz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Sxx,Sxy,Sxz,Syy,Syz,Szz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz
    real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Rxx,Rxy,Rxz,Ryy,Ryz,Rzz
    real*8,intent(in) :: eps
    real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: ham_Res, movx_Res, movy_Res, movz_Res
    real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gmx_Res, Gmy_Res, Gmz_Res
    integer, intent(out) :: gont

    ! Local variables
    integer :: kb, jb, ib, ie, je, ke
    real*8 :: dX_val, dY_val, dZ_val
    real*8 :: idx, idy, idz
    real*8 :: dX_check
    
    real*8, parameter :: SYM = 1.D0, ANTI= -1.D0

    dX_val = X(2) - X(1)
    dY_val = Y(2) - Y(1)
    dZ_val = Z(2) - Z(1)
    idx = ONE/dX_val
    idy = ONE/dY_val
    idz = ONE/dZ_val

    ! Optimization 1: 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
       gont = 1
       return
    endif
#endif
    gont = 0

    ! Optimization 2: Blocked Execution
    do kb = 1, ex(3), BLK
      ke = min(kb+BLK-1, ex(3))
      do jb = 1, ex(2), BLK
        je = min(jb+BLK-1, ex(2))
        do ib = 1, ex(1), BLK
           ie = min(ib+BLK-1, ex(1))
           call compute_block_kernel(ib, ie, jb, je, kb, ke)
        end do
      end do
    end do

  contains

    subroutine compute_block_kernel(is, ie, js, je, ks, ke)
      integer, intent(in) :: is, ie, js, je, ks, ke
      integer :: nx, ny, nz
      integer :: ii, jj, kk, ig, jg, kg
      
      ! Local Arrays (Stack Allocated)
      ! Derivatives
      real*8 :: gxxx(BLK,BLK,BLK), gxxy(BLK,BLK,BLK), gxxz(BLK,BLK,BLK)
      real*8 :: gxyx(BLK,BLK,BLK), gxyy(BLK,BLK,BLK), gxyz(BLK,BLK,BLK)
      real*8 :: gxzx(BLK,BLK,BLK), gxzy(BLK,BLK,BLK), gxzz(BLK,BLK,BLK)
      real*8 :: gyyx(BLK,BLK,BLK), gyyy(BLK,BLK,BLK), gyyz(BLK,BLK,BLK)
      real*8 :: gyzx(BLK,BLK,BLK), gyzy(BLK,BLK,BLK), gyzz(BLK,BLK,BLK)
      real*8 :: gzzx(BLK,BLK,BLK), gzzy(BLK,BLK,BLK), gzzz(BLK,BLK,BLK)
      real*8 :: chix(BLK,BLK,BLK), chiy(BLK,BLK,BLK), chiz(BLK,BLK,BLK)
      real*8 :: betaxx(BLK,BLK,BLK), betaxy(BLK,BLK,BLK), betaxz(BLK,BLK,BLK)
      real*8 :: betayx(BLK,BLK,BLK), betayy(BLK,BLK,BLK), betayz(BLK,BLK,BLK)
      real*8 :: betazx(BLK,BLK,BLK), betazy(BLK,BLK,BLK), betazz(BLK,BLK,BLK)
      real*8 :: Lapx(BLK,BLK,BLK), Lapy(BLK,BLK,BLK), Lapz(BLK,BLK,BLK)
      real*8 :: Kx(BLK,BLK,BLK), Ky(BLK,BLK,BLK), Kz(BLK,BLK,BLK)
      real*8 :: Gamxx(BLK,BLK,BLK), Gamxy(BLK,BLK,BLK), Gamxz(BLK,BLK,BLK)
      real*8 :: Gamyx(BLK,BLK,BLK), Gamyy(BLK,BLK,BLK), Gamyz(BLK,BLK,BLK)
      real*8 :: Gamzx(BLK,BLK,BLK), Gamzy(BLK,BLK,BLK), Gamzz(BLK,BLK,BLK)
      
      ! 2nd Derivatives
      real*8 :: fxx(BLK,BLK,BLK), fxy(BLK,BLK,BLK), fxz(BLK,BLK,BLK)
      real*8 :: fyy(BLK,BLK,BLK), fyz(BLK,BLK,BLK), fzz(BLK,BLK,BLK)
      real*8 :: f(BLK,BLK,BLK)
      
      ! Physics Intermediates
      real*8 :: gupxx(BLK,BLK,BLK), gupxy(BLK,BLK,BLK), gupxz(BLK,BLK,BLK)
      real*8 :: gupyy(BLK,BLK,BLK), gupyz(BLK,BLK,BLK), gupzz(BLK,BLK,BLK)
      real*8 :: Rxx_l(BLK,BLK,BLK), Rxy_l(BLK,BLK,BLK), Rxz_l(BLK,BLK,BLK)
      real*8 :: Ryy_l(BLK,BLK,BLK), Ryz_l(BLK,BLK,BLK), Rzz_l(BLK,BLK,BLK)
      
      real*8 :: alpn1, chin1, gdet, igdet, div_beta
      real*8 :: Gamxa, Gamya, Gamza, S_val
      real*8 :: val_gxx, val_gyy, val_gzz
      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
      
      nx = ie - is + 1
      ny = je - js + 1
      nz = ke - ks + 1
      
      ! === 1. Compute 1st Derivatives ===
      call calc_derivs(betax, betaxx, betaxy, betaxz, ANTI, SYM, SYM, is, ie, js, je, ks, ke)
      call calc_derivs(betay, betayx, betayy, betayz, SYM, ANTI, SYM, is, ie, js, je, ks, ke)
      call calc_derivs(betaz, betazx, betazy, betazz, SYM, SYM, ANTI, is, ie, js, je, ks, ke)
      
      call calc_derivs(chi, chix, chiy, chiz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      
      call calc_derivs(dxx, gxxx, gxxy, gxxz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      call calc_derivs(gxy, gxyx, gxyy, gxyz, ANTI, ANTI, SYM, is, ie, js, je, ks, ke)
      call calc_derivs(gxz, gxzx, gxzy, gxzz, ANTI, SYM, ANTI, is, ie, js, je, ks, ke)
      call calc_derivs(dyy, gyyx, gyyy, gyyz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      call calc_derivs(gyz, gyzx, gyzy, gyzz, SYM, ANTI, ANTI, is, ie, js, je, ks, ke)
      call calc_derivs(dzz, gzzx, gzzy, gzzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      
      call calc_derivs(Lap, Lapx, Lapy, Lapz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      call calc_derivs(trK, Kx, Ky, Kz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      
      call calc_derivs(Gamx, Gamxx, Gamxy, Gamxz, ANTI, SYM, SYM, is, ie, js, je, ks, ke)
      call calc_derivs(Gamy, Gamyx, Gamyy, Gamyz, SYM, ANTI, SYM, is, ie, js, je, ks, ke)
      call calc_derivs(Gamz, Gamzx, Gamzy, Gamzz, SYM, SYM, ANTI, is, ie, js, je, ks, ke)

      ! === 2. Compute Evolution RHS (Part 1) ===
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         
         div_beta = betaxx(ii,jj,kk) + betayy(ii,jj,kk) + betazz(ii,jj,kk)
         alpn1 = Lap(ig,jg,kg) + ONE
         chin1 = chi(ig,jg,kg) + ONE
         
         ! chi_rhs
         chi_rhs(ig,jg,kg) = F2o3 * chin1 * (alpn1 * trK(ig,jg,kg) - div_beta)
         
         ! Metric RHS
         val_gxx = dxx(ig,jg,kg) + ONE
         val_gyy = dyy(ig,jg,kg) + ONE
         val_gzz = dzz(ig,jg,kg) + ONE
         
         gxx_rhs(ig,jg,kg) = -TWO*alpn1*Axx(ig,jg,kg) - F2o3*val_gxx*div_beta + &
             TWO*(val_gxx*betaxx(ii,jj,kk) + gxy(ig,jg,kg)*betayx(ii,jj,kk) + gxz(ig,jg,kg)*betazx(ii,jj,kk))
         gyy_rhs(ig,jg,kg) = -TWO*alpn1*Ayy(ig,jg,kg) - F2o3*val_gyy*div_beta + &
             TWO*(gxy(ig,jg,kg)*betaxy(ii,jj,kk) + val_gyy*betayy(ii,jj,kk) + gyz(ig,jg,kg)*betazy(ii,jj,kk))
         gzz_rhs(ig,jg,kg) = -TWO*alpn1*Azz(ig,jg,kg) - F2o3*val_gzz*div_beta + &
             TWO*(gxz(ig,jg,kg)*betaxz(ii,jj,kk) + gyz(ig,jg,kg)*betayz(ii,jj,kk) + val_gzz*betazz(ii,jj,kk))
             
         gxy_rhs(ig,jg,kg) = -TWO*alpn1*Axy(ig,jg,kg) + F1o3*gxy(ig,jg,kg)*div_beta + &
             val_gxx*betaxy(ii,jj,kk) + gxz(ig,jg,kg)*betazy(ii,jj,kk) + &
             val_gyy*betayx(ii,jj,kk) + gyz(ig,jg,kg)*betazx(ii,jj,kk) - gxy(ig,jg,kg)*betazz(ii,jj,kk)
             
         gyz_rhs(ig,jg,kg) = -TWO*alpn1*Ayz(ig,jg,kg) + F1o3*gyz(ig,jg,kg)*div_beta + &
             gxy(ig,jg,kg)*betaxz(ii,jj,kk) + val_gyy*betayz(ii,jj,kk) + &
             gxz(ig,jg,kg)*betaxy(ii,jj,kk) + val_gzz*betazy(ii,jj,kk) - gyz(ig,jg,kg)*betaxx(ii,jj,kk)
             
         gxz_rhs(ig,jg,kg) = -TWO*alpn1*Axz(ig,jg,kg) + F1o3*gxz(ig,jg,kg)*div_beta + &
             val_gxx*betaxz(ii,jj,kk) + gxy(ig,jg,kg)*betayz(ii,jj,kk) + &
             gyz(ig,jg,kg)*betayx(ii,jj,kk) + val_gzz*betazx(ii,jj,kk) - gxz(ig,jg,kg)*betayy(ii,jj,kk)
             
         ! Inverse Metric (Optimization 4: Optimized 3x3 matrix inversion)
         ! Precompute squared terms to avoid redundant multiplications
         gxy2 = gxy(ig,jg,kg)**2
         gxz2 = gxz(ig,jg,kg)**2
         gyz2 = gyz(ig,jg,kg)**2

         ! Compute determinant using precomputed squares
         gdet = val_gxx*val_gyy*val_gzz + TWO*gxy(ig,jg,kg)*gxz(ig,jg,kg)*gyz(ig,jg,kg) - &
                val_gxx*gyz2 - val_gyy*gxz2 - val_gzz*gxy2
         igdet = ONE / gdet

         ! Compute cofactors and multiply by inverse determinant
         gupxx(ii,jj,kk) = (val_gyy*val_gzz - gyz2) * igdet
         gupyy(ii,jj,kk) = (val_gxx*val_gzz - gxz2) * igdet
         gupzz(ii,jj,kk) = (val_gxx*val_gyy - gxy2) * igdet
         gupxy(ii,jj,kk) = (gxz(ig,jg,kg)*gyz(ig,jg,kg) - gxy(ig,jg,kg)*val_gzz) * igdet
         gupxz(ii,jj,kk) = (gxy(ig,jg,kg)*gyz(ig,jg,kg) - gxz(ig,jg,kg)*val_gyy) * igdet
         gupyz(ii,jj,kk) = (gxy(ig,jg,kg)*gxz(ig,jg,kg) - gyz(ig,jg,kg)*val_gxx) * igdet
         
         ! Compute Gam^i Residual (Constraint)
         if(co == 0) then
            Gmx_Res(ig,jg,kg) = Gamx(ig,jg,kg) - ( &
               gupxx(ii,jj,kk)*(gupxx(ii,jj,kk)*gxxx(ii,jj,kk)+gupxy(ii,jj,kk)*gxyx(ii,jj,kk)+gupxz(ii,jj,kk)*gxzx(ii,jj,kk)) + &
               gupxy(ii,jj,kk)*(gupxx(ii,jj,kk)*gxyx(ii,jj,kk)+gupxy(ii,jj,kk)*gyyx(ii,jj,kk)+gupxz(ii,jj,kk)*gyzx(ii,jj,kk)) + &
               gupxz(ii,jj,kk)*(gupxx(ii,jj,kk)*gxzx(ii,jj,kk)+gupxy(ii,jj,kk)*gyzx(ii,jj,kk)+gupxz(ii,jj,kk)*gzzx(ii,jj,kk)) + &
               gupxx(ii,jj,kk)*(gupxy(ii,jj,kk)*gxxy(ii,jj,kk)+gupyy(ii,jj,kk)*gxyy(ii,jj,kk)+gupyz(ii,jj,kk)*gxzy(ii,jj,kk)) + &
               gupxy(ii,jj,kk)*(gupxy(ii,jj,kk)*gxyy(ii,jj,kk)+gupyy(ii,jj,kk)*gyyy(ii,jj,kk)+gupyz(ii,jj,kk)*gyzy(ii,jj,kk)) + &
               gupxz(ii,jj,kk)*(gupxy(ii,jj,kk)*gxzy(ii,jj,kk)+gupyy(ii,jj,kk)*gyzy(ii,jj,kk)+gupyz(ii,jj,kk)*gzzy(ii,jj,kk)) + &
               gupxx(ii,jj,kk)*(gupxz(ii,jj,kk)*gxxz(ii,jj,kk)+gupyz(ii,jj,kk)*gxyz(ii,jj,kk)+gupzz(ii,jj,kk)*gxzz(ii,jj,kk)) + &
               gupxy(ii,jj,kk)*(gupxz(ii,jj,kk)*gxyz(ii,jj,kk)+gupyz(ii,jj,kk)*gyyz(ii,jj,kk)+gupzz(ii,jj,kk)*gyzz(ii,jj,kk)) + &
               gupxz(ii,jj,kk)*(gupxz(ii,jj,kk)*gxzz(ii,jj,kk)+gupyz(ii,jj,kk)*gyzz(ii,jj,kk)+gupzz(ii,jj,kk)*gzzz(ii,jj,kk)) )
         end if
      end do; end do; end do

      ! === 3. Christoffel Symbols (Optimization 5: CSE) ===
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1

         term_xx1 = TWO*gxyx(ii,jj,kk) - gxxy(ii,jj,kk)
         term_xx2 = TWO*gxzx(ii,jj,kk) - gxxz(ii,jj,kk)
         term_yy1 = TWO*gxyy(ii,jj,kk) - gyyx(ii,jj,kk)
         term_yy2 = TWO*gyzy(ii,jj,kk) - gyyz(ii,jj,kk)
         term_zz1 = TWO*gxzz(ii,jj,kk) - gzzx(ii,jj,kk)
         term_zz2 = TWO*gyzz(ii,jj,kk) - gzzy(ii,jj,kk)
         term_xy = gxzy(ii,jj,kk) + gyzx(ii,jj,kk) - gxyz(ii,jj,kk)
         term_xz = gxyz(ii,jj,kk) + gyzx(ii,jj,kk) - gxzy(ii,jj,kk)
         term_yz = gxyz(ii,jj,kk) + gxzy(ii,jj,kk) - gyzx(ii,jj,kk)

         Gamxxx(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*gxxx(ii,jj,kk) + gupxy(ii,jj,kk)*term_xx1 + gupxz(ii,jj,kk)*term_xx2)
         Gamyxx(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*gxxx(ii,jj,kk) + gupyy(ii,jj,kk)*term_xx1 + gupyz(ii,jj,kk)*term_xx2)
         Gamzxx(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*gxxx(ii,jj,kk) + gupyz(ii,jj,kk)*term_xx1 + gupzz(ii,jj,kk)*term_xx2)

         Gamxyy(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*term_yy1 + gupxy(ii,jj,kk)*gyyy(ii,jj,kk) + gupxz(ii,jj,kk)*term_yy2)
         Gamyyy(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*term_yy1 + gupyy(ii,jj,kk)*gyyy(ii,jj,kk) + gupyz(ii,jj,kk)*term_yy2)
         Gamzyy(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*term_yy1 + gupyz(ii,jj,kk)*gyyy(ii,jj,kk) + gupzz(ii,jj,kk)*term_yy2)

         Gamxzz(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*term_zz1 + gupxy(ii,jj,kk)*term_zz2 + gupxz(ii,jj,kk)*gzzz(ii,jj,kk))
         Gamyzz(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*term_zz1 + gupyy(ii,jj,kk)*term_zz2 + gupyz(ii,jj,kk)*gzzz(ii,jj,kk))
         Gamzzz(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*term_zz1 + gupyz(ii,jj,kk)*term_zz2 + gupzz(ii,jj,kk)*gzzz(ii,jj,kk))

         Gamxxy(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*gxxy(ii,jj,kk) + gupxy(ii,jj,kk)*gyyx(ii,jj,kk) + gupxz(ii,jj,kk)*term_xy)
         Gamyxy(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*gxxy(ii,jj,kk) + gupyy(ii,jj,kk)*gyyx(ii,jj,kk) + gupyz(ii,jj,kk)*term_xy)
         Gamzxy(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*gxxy(ii,jj,kk) + gupyz(ii,jj,kk)*gyyx(ii,jj,kk) + gupzz(ii,jj,kk)*term_xy)

         Gamxxz(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*gxxz(ii,jj,kk) + gupxy(ii,jj,kk)*term_xz + gupxz(ii,jj,kk)*gzzx(ii,jj,kk))
         Gamyxz(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*gxxz(ii,jj,kk) + gupyy(ii,jj,kk)*term_xz + gupyz(ii,jj,kk)*gzzx(ii,jj,kk))
         Gamzxz(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*gxxz(ii,jj,kk) + gupyz(ii,jj,kk)*term_xz + gupzz(ii,jj,kk)*gzzx(ii,jj,kk))

         Gamxyz(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*term_yz + gupxy(ii,jj,kk)*gyyz(ii,jj,kk) + gupxz(ii,jj,kk)*gzzy(ii,jj,kk))
         Gamyyz(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*term_yz + gupyy(ii,jj,kk)*gyyz(ii,jj,kk) + gupyz(ii,jj,kk)*gzzy(ii,jj,kk))
         Gamzyz(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*term_yz + gupyz(ii,jj,kk)*gyyz(ii,jj,kk) + gupzz(ii,jj,kk)*gzzy(ii,jj,kk))
      end do; end do; end do

      ! === 4. Ricci Tensor (Part 1 - Terms from A_ij) ===
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         
         Rxx_l(ii,jj,kk) = gupxx(ii,jj,kk)**2*Axx(ig,jg,kg) + gupxy(ii,jj,kk)**2*Ayy(ig,jg,kg) + gupxz(ii,jj,kk)**2*Azz(ig,jg,kg) + &
                           TWO*(gupxx(ii,jj,kk)*gupxy(ii,jj,kk)*Axy(ig,jg,kg) + gupxx(ii,jj,kk)*gupxz(ii,jj,kk)*Axz(ig,jg,kg) + gupxy(ii,jj,kk)*gupxz(ii,jj,kk)*Ayz(ig,jg,kg))
         Ryy_l(ii,jj,kk) = gupxy(ii,jj,kk)**2*Axx(ig,jg,kg) + gupyy(ii,jj,kk)**2*Ayy(ig,jg,kg) + gupyz(ii,jj,kk)**2*Azz(ig,jg,kg) + &
                           TWO*(gupxy(ii,jj,kk)*gupyy(ii,jj,kk)*Axy(ig,jg,kg) + gupxy(ii,jj,kk)*gupyz(ii,jj,kk)*Axz(ig,jg,kg) + gupyy(ii,jj,kk)*gupyz(ii,jj,kk)*Ayz(ig,jg,kg))
         Rzz_l(ii,jj,kk) = gupxz(ii,jj,kk)**2*Axx(ig,jg,kg) + gupyz(ii,jj,kk)**2*Ayy(ig,jg,kg) + gupzz(ii,jj,kk)**2*Azz(ig,jg,kg) + &
                           TWO*(gupxz(ii,jj,kk)*gupyz(ii,jj,kk)*Axy(ig,jg,kg) + gupxz(ii,jj,kk)*gupzz(ii,jj,kk)*Axz(ig,jg,kg) + gupyz(ii,jj,kk)*gupzz(ii,jj,kk)*Ayz(ig,jg,kg))
                           
         Rxy_l(ii,jj,kk) = gupxx(ii,jj,kk)*gupxy(ii,jj,kk)*Axx(ig,jg,kg) + gupxy(ii,jj,kk)*gupyy(ii,jj,kk)*Ayy(ig,jg,kg) + gupxz(ii,jj,kk)*gupyz(ii,jj,kk)*Azz(ig,jg,kg) + &
                           (gupxx(ii,jj,kk)*gupyy(ii,jj,kk) + gupxy(ii,jj,kk)**2)*Axy(ig,jg,kg) + &
                           (gupxx(ii,jj,kk)*gupyz(ii,jj,kk) + gupxz(ii,jj,kk)*gupxy(ii,jj,kk))*Axz(ig,jg,kg) + &
                           (gupxy(ii,jj,kk)*gupyz(ii,jj,kk) + gupxz(ii,jj,kk)*gupyy(ii,jj,kk))*Ayz(ig,jg,kg)
         Rxz_l(ii,jj,kk) = gupxx(ii,jj,kk)*gupxz(ii,jj,kk)*Axx(ig,jg,kg) + gupxy(ii,jj,kk)*gupyz(ii,jj,kk)*Ayy(ig,jg,kg) + gupxz(ii,jj,kk)*gupzz(ii,jj,kk)*Azz(ig,jg,kg) + &
                           (gupxx(ii,jj,kk)*gupyz(ii,jj,kk) + gupxy(ii,jj,kk)*gupxz(ii,jj,kk))*Axy(ig,jg,kg) + &
                           (gupxx(ii,jj,kk)*gupzz(ii,jj,kk) + gupxz(ii,jj,kk)**2)*Axz(ig,jg,kg) + &
                           (gupxy(ii,jj,kk)*gupzz(ii,jj,kk) + gupxz(ii,jj,kk)*gupyz(ii,jj,kk))*Ayz(ig,jg,kg)
         Ryz_l(ii,jj,kk) = gupxy(ii,jj,kk)*gupxz(ii,jj,kk)*Axx(ig,jg,kg) + gupyy(ii,jj,kk)*gupyz(ii,jj,kk)*Ayy(ig,jg,kg) + gupyz(ii,jj,kk)*gupzz(ii,jj,kk)*Azz(ig,jg,kg) + &
                           (gupxy(ii,jj,kk)*gupyz(ii,jj,kk) + gupyy(ii,jj,kk)*gupxz(ii,jj,kk))*Axy(ig,jg,kg) + &
                           (gupxy(ii,jj,kk)*gupzz(ii,jj,kk) + gupyz(ii,jj,kk)*gupxz(ii,jj,kk))*Axz(ig,jg,kg) + &
                           (gupyy(ii,jj,kk)*gupzz(ii,jj,kk) + gupyz(ii,jj,kk)**2)*Ayz(ig,jg,kg)
      end do; end do; end do

      ! === 5. Gamma^i RHS ===
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         alpn1 = Lap(ig,jg,kg) + ONE
         chin1 = chi(ig,jg,kg) + ONE
         
         Gamx_rhs(ig,jg,kg) = -TWO*(Lapx(ii,jj,kk)*Rxx_l(ii,jj,kk) + Lapy(ii,jj,kk)*Rxy_l(ii,jj,kk) + Lapz(ii,jj,kk)*Rxz_l(ii,jj,kk)) + &
                              TWO*alpn1*( -F3o2/chin1*(chix(ii,jj,kk)*Rxx_l(ii,jj,kk) + chiy(ii,jj,kk)*Rxy_l(ii,jj,kk) + chiz(ii,jj,kk)*Rxz_l(ii,jj,kk)) - &
                              gupxx(ii,jj,kk)*(F2o3*Kx(ii,jj,kk) + EIGHT*PI*Sx(ig,jg,kg)) - &
                              gupxy(ii,jj,kk)*(F2o3*Ky(ii,jj,kk) + EIGHT*PI*Sy(ig,jg,kg)) - &
                              gupxz(ii,jj,kk)*(F2o3*Kz(ii,jj,kk) + EIGHT*PI*Sz(ig,jg,kg)) + &
                              Gamxxx(ig,jg,kg)*Rxx_l(ii,jj,kk) + Gamxyy(ig,jg,kg)*Ryy_l(ii,jj,kk) + Gamxzz(ig,jg,kg)*Rzz_l(ii,jj,kk) + &
                              TWO*(Gamxxy(ig,jg,kg)*Rxy_l(ii,jj,kk) + Gamxxz(ig,jg,kg)*Rxz_l(ii,jj,kk) + Gamxyz(ig,jg,kg)*Ryz_l(ii,jj,kk)) )

         Gamy_rhs(ig,jg,kg) = -TWO*(Lapx(ii,jj,kk)*Rxy_l(ii,jj,kk) + Lapy(ii,jj,kk)*Ryy_l(ii,jj,kk) + Lapz(ii,jj,kk)*Ryz_l(ii,jj,kk)) + &
                              TWO*alpn1*( -F3o2/chin1*(chix(ii,jj,kk)*Rxy_l(ii,jj,kk) + chiy(ii,jj,kk)*Ryy_l(ii,jj,kk) + chiz(ii,jj,kk)*Ryz_l(ii,jj,kk)) - &
                              gupxy(ii,jj,kk)*(F2o3*Kx(ii,jj,kk) + EIGHT*PI*Sx(ig,jg,kg)) - &
                              gupyy(ii,jj,kk)*(F2o3*Ky(ii,jj,kk) + EIGHT*PI*Sy(ig,jg,kg)) - &
                              gupyz(ii,jj,kk)*(F2o3*Kz(ii,jj,kk) + EIGHT*PI*Sz(ig,jg,kg)) + &
                              Gamyxx(ig,jg,kg)*Rxx_l(ii,jj,kk) + Gamyyy(ig,jg,kg)*Ryy_l(ii,jj,kk) + Gamyzz(ig,jg,kg)*Rzz_l(ii,jj,kk) + &
                              TWO*(Gamyxy(ig,jg,kg)*Rxy_l(ii,jj,kk) + Gamyxz(ig,jg,kg)*Rxz_l(ii,jj,kk) + Gamyyz(ig,jg,kg)*Ryz_l(ii,jj,kk)) )

         Gamz_rhs(ig,jg,kg) = -TWO*(Lapx(ii,jj,kk)*Rxz_l(ii,jj,kk) + Lapy(ii,jj,kk)*Ryz_l(ii,jj,kk) + Lapz(ii,jj,kk)*Rzz_l(ii,jj,kk)) + &
                              TWO*alpn1*( -F3o2/chin1*(chix(ii,jj,kk)*Rxz_l(ii,jj,kk) + chiy(ii,jj,kk)*Ryz_l(ii,jj,kk) + chiz(ii,jj,kk)*Rzz_l(ii,jj,kk)) - &
                              gupxz(ii,jj,kk)*(F2o3*Kx(ii,jj,kk) + EIGHT*PI*Sx(ig,jg,kg)) - &
                              gupyz(ii,jj,kk)*(F2o3*Ky(ii,jj,kk) + EIGHT*PI*Sy(ig,jg,kg)) - &
                              gupzz(ii,jj,kk)*(F2o3*Kz(ii,jj,kk) + EIGHT*PI*Sz(ig,jg,kg)) + &
                              Gamzxx(ig,jg,kg)*Rxx_l(ii,jj,kk) + Gamzyy(ig,jg,kg)*Ryy_l(ii,jj,kk) + Gamzzz(ig,jg,kg)*Rzz_l(ii,jj,kk) + &
                              TWO*(Gamzxy(ig,jg,kg)*Rxy_l(ii,jj,kk) + Gamzxz(ig,jg,kg)*Rxz_l(ii,jj,kk) + Gamzyz(ig,jg,kg)*Ryz_l(ii,jj,kk)) )
      end do; end do; end do

      ! === 6. Second Derivatives for Ricci ===
      ! Need 2nd derivatives of gxx etc.
      call calc_dderivs(dxx, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         val_gxx = dxx(ig,jg,kg)+ONE
         val_gyy = dyy(ig,jg,kg)+ONE
         val_gzz = dzz(ig,jg,kg)+ONE
         
         Rxx_l(ii,jj,kk) = -HALF*Rxx_l(ii,jj,kk) + val_gxx*Gamxx(ii,jj,kk) + gxy(ig,jg,kg)*Gamyx(ii,jj,kk) + gxz(ig,jg,kg)*Gamzx(ii,jj,kk) + &
                           (gupxx(ii,jj,kk)*Gamxxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamxyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamxzz(ig,jg,kg) + &
                            TWO*(gupxy(ii,jj,kk)*Gamxxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamxxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamxyz(ig,jg,kg)))*gxxx(ii,jj,kk) + &
                           (gupxx(ii,jj,kk)*Gamyxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamyyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamyzz(ig,jg,kg) + &
                            TWO*(gupxy(ii,jj,kk)*Gamyxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamyxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamyyz(ig,jg,kg)))*gxyx(ii,jj,kk) + &
                           (gupxx(ii,jj,kk)*Gamzxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamzyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamzzz(ig,jg,kg) + &
                            TWO*(gupxy(ii,jj,kk)*Gamzxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamzxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamzyz(ig,jg,kg)))*gxzx(ii,jj,kk) + &
                           gupxx(ii,jj,kk)*(fxx(ii,jj,kk) + TWO*(Gamxxx(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzx(ii,jj,kk)) + Gamxxx(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxxy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxxz(ii,jj,kk)) + &
                           gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + TWO*(Gamxxx(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzx(ii,jj,kk)) + Gamxxy(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxxy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxxz(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxyz(ii,jj,kk)) + &
                           gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + TWO*(Gamxxx(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxx(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzx(ii,jj,kk)) + Gamxxz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzz(ii,jj,kk)) + &
                           gupyy(ii,jj,kk)*(fyy(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzx(ii,jj,kk)) + Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxyz(ii,jj,kk)) + &
                           gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzx(ii,jj,kk)) + Gamxxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk)) + &
                           gupzz(ii,jj,kk)*(fzz(ii,jj,kk) + TWO*(Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzx(ii,jj,kk)) + Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk))
      end do; end do; end do

      call calc_dderivs(dyy, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         val_gyy = dyy(ig,jg,kg)+ONE
         Ryy_l(ii,jj,kk) = -HALF*Ryy_l(ii,jj,kk) + gxy(ig,jg,kg)*Gamxy(ii,jj,kk) + val_gyy*Gamyy(ii,jj,kk) + gyz(ig,jg,kg)*Gamzy(ii,jj,kk) + &
                           (gupxx(ii,jj,kk)*Gamxxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamxyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamxzz(ig,jg,kg) + &
                            TWO*(gupxy(ii,jj,kk)*Gamxxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamxxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamxyz(ig,jg,kg)))*gxyy(ii,jj,kk) + &
                           (gupxx(ii,jj,kk)*Gamyxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamyyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamyzz(ig,jg,kg) + &
                            TWO*(gupxy(ii,jj,kk)*Gamyxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamyxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamyyz(ig,jg,kg)))*gyyy(ii,jj,kk) + &
                           (gupxx(ii,jj,kk)*Gamzxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamzyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamzzz(ig,jg,kg) + &
                            TWO*(gupxy(ii,jj,kk)*Gamzxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamzxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamzyz(ig,jg,kg)))*gyzy(ii,jj,kk) + &
                           gupxx(ii,jj,kk)*(fxx(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzy(ii,jj,kk)) + Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxyz(ii,jj,kk)) + &
                           gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxzy(ii,jj,kk)) + Gamxyy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyyz(ii,jj,kk)) + &
                           gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzy(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupyy(ii,jj,kk)*(fyy(ii,jj,kk) + TWO*(Gamxyy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzy(ii,jj,kk)) + Gamxyy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyyz(ii,jj,kk)) + &
                           gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + TWO*(Gamxyy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzy(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamxyy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupzz(ii,jj,kk)*(fzz(ii,jj,kk) + TWO*(Gamxyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzy(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk))
      end do; end do; end do

      call calc_dderivs(dzz, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         val_gzz = dzz(ig,jg,kg)+ONE
         Rzz_l(ii,jj,kk) = -HALF*Rzz_l(ii,jj,kk) + gxz(ig,jg,kg)*Gamxz(ii,jj,kk) + gyz(ig,jg,kg)*Gamyz(ii,jj,kk) + val_gzz*Gamzz(ii,jj,kk) + &
                           (gupxx(ii,jj,kk)*Gamxxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamxyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamxzz(ig,jg,kg) + &
                            TWO*(gupxy(ii,jj,kk)*Gamxxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamxxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamxyz(ig,jg,kg)))*gxzz(ii,jj,kk) + &
                           (gupxx(ii,jj,kk)*Gamyxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamyyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamyzz(ig,jg,kg) + &
                            TWO*(gupxy(ii,jj,kk)*Gamyxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamyxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamyyz(ig,jg,kg)))*gyzz(ii,jj,kk) + &
                           (gupxx(ii,jj,kk)*Gamzxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamzyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamzzz(ig,jg,kg) + &
                            TWO*(gupxy(ii,jj,kk)*Gamzxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamzxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamzyz(ig,jg,kg)))*gzzz(ii,jj,kk) + &
                           gupxx(ii,jj,kk)*(fxx(ii,jj,kk) + TWO*(Gamxxz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk)) + Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk)) + &
                           gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + TWO*(Gamxxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzz(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + TWO*(Gamxxz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyzz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzzz(ig,jg,kg)*gxzz(ii,jj,kk)) + Gamxzz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzz(ii,jj,kk)) + &
                           gupyy(ii,jj,kk)*(fyy(ii,jj,kk) + TWO*(Gamxyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + TWO*(Gamxyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzzz(ig,jg,kg)*gyzz(ii,jj,kk)) + Gamxzz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzz(ii,jj,kk)) + &
                           gupzz(ii,jj,kk)*(fzz(ii,jj,kk) + TWO*(Gamxzz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzzz(ig,jg,kg)*gzzz(ii,jj,kk)) + Gamxzz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gzzz(ii,jj,kk))
      end do; end do; end do

      call calc_dderivs(gxy, fxx, fxy, fxz, fyy, fyz, fzz, ANTI, ANTI, SYM, is, ie, js, je, ks, ke)
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         
         ! Compute Gam^i at local point
         Gamxa = Gamx(ig,jg,kg)
         Gamya = Gamy(ig,jg,kg)
         Gamza = Gamz(ig,jg,kg)
         val_gxx = dxx(ig,jg,kg)+ONE
         val_gyy = dyy(ig,jg,kg)+ONE

         Rxy_l(ii,jj,kk) = HALF*( -Rxy_l(ii,jj,kk) + val_gxx*Gamxy(ii,jj,kk) + gxy(ig,jg,kg)*Gamyy(ii,jj,kk) + gxz(ig,jg,kg)*Gamzy(ii,jj,kk) + &
                           gxy(ig,jg,kg)*Gamxx(ii,jj,kk) + val_gyy*Gamyx(ii,jj,kk) + gyz(ig,jg,kg)*Gamzx(ii,jj,kk) + &
                           Gamxa*gxyx(ii,jj,kk) + Gamya*gyyx(ii,jj,kk) + Gamza*gyzx(ii,jj,kk) + &
                           Gamxa*gxxy(ii,jj,kk) + Gamya*gxyy(ii,jj,kk) + Gamza*gxzy(ii,jj,kk) ) + &
                           gupxx(ii,jj,kk)*(Gamxxx(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxyz(ii,jj,kk)) + &
                           gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyyz(ii,jj,kk)) + &
                           gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupyy(ii,jj,kk)*(Gamxxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyyz(ii,jj,kk)) + &
                           gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzyy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupzz(ii,jj,kk)*(Gamxxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk))
      end do; end do; end do

      call calc_dderivs(gxz, fxx, fxy, fxz, fyy, fyz, fzz, ANTI, SYM, ANTI, is, ie, js, je, ks, ke)
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         
         ! Compute Gam^i at local point
         Gamxa = Gamx(ig,jg,kg)
         Gamya = Gamy(ig,jg,kg)
         Gamza = Gamz(ig,jg,kg)
         val_gxx = dxx(ig,jg,kg)+ONE
         val_gzz = dzz(ig,jg,kg)+ONE

         Rxz_l(ii,jj,kk) = HALF*( -Rxz_l(ii,jj,kk) + val_gxx*Gamxz(ii,jj,kk) + gxy(ig,jg,kg)*Gamyz(ii,jj,kk) + gxz(ig,jg,kg)*Gamzz(ii,jj,kk) + &
                           gxz(ig,jg,kg)*Gamxx(ii,jj,kk) + gyz(ig,jg,kg)*Gamyx(ii,jj,kk) + val_gzz*Gamzx(ii,jj,kk) + &
                           Gamxa*gxzx(ii,jj,kk) + Gamya*gyzx(ii,jj,kk) + Gamza*gzzx(ii,jj,kk) + &
                           Gamxa*gxxz(ii,jj,kk) + Gamya*gxyz(ii,jj,kk) + Gamza*gxzz(ii,jj,kk) ) + &
                           gupxx(ii,jj,kk)*(fxx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzz(ii,jj,kk)) + &
                           gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxx(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzzz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gzzz(ii,jj,kk)) + &
                           gupyy(ii,jj,kk)*(Gamxxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzzz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzz(ii,jj,kk)) + &
                           gupzz(ii,jj,kk)*(Gamxxz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzzz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzz(ii,jj,kk))
      end do; end do; end do

      call calc_dderivs(gyz, fxx, fxy, fxz, fyy, fyz, fzz, SYM, ANTI, ANTI, is, ie, js, je, ks, ke)
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         
         ! Compute Gam^i at local point
         Gamxa = Gamx(ig,jg,kg)
         Gamya = Gamy(ig,jg,kg)
         Gamza = Gamz(ig,jg,kg)
         val_gyy = dyy(ig,jg,kg)+ONE
         val_gzz = dzz(ig,jg,kg)+ONE

         Ryz_l(ii,jj,kk) = HALF*( -Ryz_l(ii,jj,kk) + gxy(ig,jg,kg)*Gamxz(ii,jj,kk) + val_gyy*Gamyz(ii,jj,kk) + gyz(ig,jg,kg)*Gamzz(ii,jj,kk) + &
                           gxz(ig,jg,kg)*Gamxy(ii,jj,kk) + gyz(ig,jg,kg)*Gamyy(ii,jj,kk) + val_gzz*Gamzy(ii,jj,kk) + &
                           Gamxa*gxzy(ii,jj,kk) + Gamya*gyzy(ii,jj,kk) + Gamza*gzzy(ii,jj,kk) + &
                           Gamxa*gxyz(ii,jj,kk) + Gamya*gyyz(ii,jj,kk) + Gamza*gyzz(ii,jj,kk) ) + &
                           gupxx(ii,jj,kk)*(Gamxxy(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk)) + &
                           gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyzz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzz(ii,jj,kk)) + &
                           gupyy(ii,jj,kk)*(fyy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
                           gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzyy(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gzzz(ii,jj,kk)) + &
                           gupzz(ii,jj,kk)*(fzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzz(ii,jj,kk))
      end do; end do; end do

      ! === 7. Chi Terms for Ricci ===
      call calc_dderivs(chi, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         chin1 = chi(ig,jg,kg) + ONE
         
         ! Compute trace part f (D^2 chi terms)
         f(ii,jj,kk) = gupxx(ii,jj,kk)*(fxx(ii,jj,kk) - F3o2/chin1*chix(ii,jj,kk)**2) + &
                       gupyy(ii,jj,kk)*(fyy(ii,jj,kk) - F3o2/chin1*chiy(ii,jj,kk)**2) + &
                       gupzz(ii,jj,kk)*(fzz(ii,jj,kk) - F3o2/chin1*chiz(ii,jj,kk)**2) + &
                   TWO*(gupxy(ii,jj,kk)*(fxy(ii,jj,kk) - F3o2/chin1*chix(ii,jj,kk)*chiy(ii,jj,kk)) + &
                        gupxz(ii,jj,kk)*(fxz(ii,jj,kk) - F3o2/chin1*chix(ii,jj,kk)*chiz(ii,jj,kk)) + &
                        gupyz(ii,jj,kk)*(fyz(ii,jj,kk) - F3o2/chin1*chiy(ii,jj,kk)*chiz(ii,jj,kk)))
         
         ! Add chi terms to Ricci
         val_gxx = dxx(ig,jg,kg)+ONE; val_gyy = dyy(ig,jg,kg)+ONE; val_gzz = dzz(ig,jg,kg)+ONE
         
         Rxx_l(ii,jj,kk) = Rxx_l(ii,jj,kk) + (fxx(ii,jj,kk) - chix(ii,jj,kk)**2/chin1/TWO + val_gxx*f(ii,jj,kk))/chin1/TWO
         Ryy_l(ii,jj,kk) = Ryy_l(ii,jj,kk) + (fyy(ii,jj,kk) - chiy(ii,jj,kk)**2/chin1/TWO + val_gyy*f(ii,jj,kk))/chin1/TWO
         Rzz_l(ii,jj,kk) = Rzz_l(ii,jj,kk) + (fzz(ii,jj,kk) - chiz(ii,jj,kk)**2/chin1/TWO + val_gzz*f(ii,jj,kk))/chin1/TWO
         
         Rxy_l(ii,jj,kk) = Rxy_l(ii,jj,kk) + (fxy(ii,jj,kk) - chix(ii,jj,kk)*chiy(ii,jj,kk)/chin1/TWO + gxy(ig,jg,kg)*f(ii,jj,kk))/chin1/TWO
         Rxz_l(ii,jj,kk) = Rxz_l(ii,jj,kk) + (fxz(ii,jj,kk) - chix(ii,jj,kk)*chiz(ii,jj,kk)/chin1/TWO + gxz(ig,jg,kg)*f(ii,jj,kk))/chin1/TWO
         Ryz_l(ii,jj,kk) = Ryz_l(ii,jj,kk) + (fyz(ii,jj,kk) - chiy(ii,jj,kk)*chiz(ii,jj,kk)/chin1/TWO + gyz(ig,jg,kg)*f(ii,jj,kk))/chin1/TWO
         
         ! Also need f for A_ij RHS later? No, f is recalculated from Lap later.
      end do; end do; end do

      ! === 8. Lap Terms ===
      call calc_dderivs(Lap, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         
         ! Subtract connection terms from D_i D_j Lap (fxx is just d_i d_j Lap initially)
         fxx(ii,jj,kk) = fxx(ii,jj,kk) - Gamxxx(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyxx(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzxx(ig,jg,kg)*Lapz(ii,jj,kk)
         fyy(ii,jj,kk) = fyy(ii,jj,kk) - Gamxyy(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyyy(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzyy(ig,jg,kg)*Lapz(ii,jj,kk)
         fzz(ii,jj,kk) = fzz(ii,jj,kk) - Gamxzz(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyzz(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzzz(ig,jg,kg)*Lapz(ii,jj,kk)
         fxy(ii,jj,kk) = fxy(ii,jj,kk) - Gamxxy(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyxy(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzxy(ig,jg,kg)*Lapz(ii,jj,kk)
         fxz(ii,jj,kk) = fxz(ii,jj,kk) - Gamxxz(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyxz(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzxz(ig,jg,kg)*Lapz(ii,jj,kk)
         fyz(ii,jj,kk) = fyz(ii,jj,kk) - Gamxyz(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyyz(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzyz(ig,jg,kg)*Lapz(ii,jj,kk)
         
         ! trK_rhs partial update: g^ij D_i D_j Lap
         trK_rhs(ig,jg,kg) = chin1 * (gupxx(ii,jj,kk)*fxx(ii,jj,kk) + gupyy(ii,jj,kk)*fyy(ii,jj,kk) + gupzz(ii,jj,kk)*fzz(ii,jj,kk) + &
                             TWO*(gupxy(ii,jj,kk)*fxy(ii,jj,kk) + gupxz(ii,jj,kk)*fxz(ii,jj,kk) + gupyz(ii,jj,kk)*fyz(ii,jj,kk)))
                             
         ! Prepare fxx etc for A_ij RHS (Add R_ij and S_ij parts)
         ! Need Lap value
         alpn1 = Lap(ig,jg,kg) + ONE
         
         ! Store fxx_new = alpn1 * (Rxx - 8PI Sxx) - D_i D_j Lap
         fxx(ii,jj,kk) = alpn1 * (Rxx_l(ii,jj,kk) - EIGHT*PI*Sxx(ig,jg,kg)) - fxx(ii,jj,kk)
         fyy(ii,jj,kk) = alpn1 * (Ryy_l(ii,jj,kk) - EIGHT*PI*Syy(ig,jg,kg)) - fyy(ii,jj,kk)
         fzz(ii,jj,kk) = alpn1 * (Rzz_l(ii,jj,kk) - EIGHT*PI*Szz(ig,jg,kg)) - fzz(ii,jj,kk)
         fxy(ii,jj,kk) = alpn1 * (Rxy_l(ii,jj,kk) - EIGHT*PI*Sxy(ig,jg,kg)) - fxy(ii,jj,kk)
         fxz(ii,jj,kk) = alpn1 * (Rxz_l(ii,jj,kk) - EIGHT*PI*Sxz(ig,jg,kg)) - fxz(ii,jj,kk)
         fyz(ii,jj,kk) = alpn1 * (Ryz_l(ii,jj,kk) - EIGHT*PI*Syz(ig,jg,kg)) - fyz(ii,jj,kk)
      end do; end do; end do

      ! === 9. A_ij RHS and Final trK RHS ===
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         
         val_gxx = dxx(ig,jg,kg)+ONE; val_gyy = dyy(ig,jg,kg)+ONE; val_gzz = dzz(ig,jg,kg)+ONE
         alpn1 = Lap(ig,jg,kg) + ONE
         chin1 = chi(ig,jg,kg) + ONE
         div_beta = betaxx(ii,jj,kk) + betayy(ii,jj,kk) + betazz(ii,jj,kk)
         
         ! Trace-free part f
         f(ii,jj,kk) = F1o3 * (gupxx(ii,jj,kk)*fxx(ii,jj,kk) + gupyy(ii,jj,kk)*fyy(ii,jj,kk) + gupzz(ii,jj,kk)*fzz(ii,jj,kk) + &
                       TWO*(gupxy(ii,jj,kk)*fxy(ii,jj,kk) + gupxz(ii,jj,kk)*fxz(ii,jj,kk) + gupyz(ii,jj,kk)*fyz(ii,jj,kk)))
         
         ! Initial A_ij RHS (TF(R_ij part))
         Axx_rhs(ig,jg,kg) = fxx(ii,jj,kk) - val_gxx * f(ii,jj,kk)
         Ayy_rhs(ig,jg,kg) = fyy(ii,jj,kk) - val_gyy * f(ii,jj,kk)
         Azz_rhs(ig,jg,kg) = fzz(ii,jj,kk) - val_gzz * f(ii,jj,kk)
         Axy_rhs(ig,jg,kg) = fxy(ii,jj,kk) - gxy(ig,jg,kg) * f(ii,jj,kk)
         Axz_rhs(ig,jg,kg) = fxz(ii,jj,kk) - gxz(ig,jg,kg) * f(ii,jj,kk)
         Ayz_rhs(ig,jg,kg) = fyz(ii,jj,kk) - gyz(ig,jg,kg) * f(ii,jj,kk)
         
         ! Add A_il A^l_j terms (Calculate fxx as A_ik A^k_j)
         ! fxx = A_xx^2 ...
         ! Re-calculate fxx... reusing local arrays
         ! fxx = g^lm A_xl A_mj
         fxx(ii,jj,kk) = gupxx(ii,jj,kk)*Axx(ig,jg,kg)**2 + gupyy(ii,jj,kk)*Axy(ig,jg,kg)**2 + gupzz(ii,jj,kk)*Axz(ig,jg,kg)**2 + &
                         TWO*(gupxy(ii,jj,kk)*Axx(ig,jg,kg)*Axy(ig,jg,kg) + gupxz(ii,jj,kk)*Axx(ig,jg,kg)*Axz(ig,jg,kg) + gupyz(ii,jj,kk)*Axy(ig,jg,kg)*Axz(ig,jg,kg))
         fyy(ii,jj,kk) = gupxx(ii,jj,kk)*Axy(ig,jg,kg)**2 + gupyy(ii,jj,kk)*Ayy(ig,jg,kg)**2 + gupzz(ii,jj,kk)*Ayz(ig,jg,kg)**2 + &
                         TWO*(gupxy(ii,jj,kk)*Axy(ig,jg,kg)*Ayy(ig,jg,kg) + gupxz(ii,jj,kk)*Axy(ig,jg,kg)*Ayz(ig,jg,kg) + gupyz(ii,jj,kk)*Ayy(ig,jg,kg)*Ayz(ig,jg,kg))
         fzz(ii,jj,kk) = gupxx(ii,jj,kk)*Axz(ig,jg,kg)**2 + gupyy(ii,jj,kk)*Ayz(ig,jg,kg)**2 + gupzz(ii,jj,kk)*Azz(ig,jg,kg)**2 + &
                         TWO*(gupxy(ii,jj,kk)*Axz(ig,jg,kg)*Ayz(ig,jg,kg) + gupxz(ii,jj,kk)*Axz(ig,jg,kg)*Azz(ig,jg,kg) + gupyz(ii,jj,kk)*Ayz(ig,jg,kg)*Azz(ig,jg,kg))
         fxy(ii,jj,kk) = gupxx(ii,jj,kk)*Axx(ig,jg,kg)*Axy(ig,jg,kg) + gupyy(ii,jj,kk)*Axy(ig,jg,kg)*Ayy(ig,jg,kg) + gupzz(ii,jj,kk)*Axz(ig,jg,kg)*Ayz(ig,jg,kg) + &
                         gupxy(ii,jj,kk)*(Axx(ig,jg,kg)*Ayy(ig,jg,kg) + Axy(ig,jg,kg)**2) + &
                         gupxz(ii,jj,kk)*(Axx(ig,jg,kg)*Ayz(ig,jg,kg) + Axz(ig,jg,kg)*Axy(ig,jg,kg)) + &
                         gupyz(ii,jj,kk)*(Axy(ig,jg,kg)*Ayz(ig,jg,kg) + Axz(ig,jg,kg)*Ayy(ig,jg,kg))
         fxz(ii,jj,kk) = gupxx(ii,jj,kk)*Axx(ig,jg,kg)*Axz(ig,jg,kg) + gupyy(ii,jj,kk)*Axy(ig,jg,kg)*Ayz(ig,jg,kg) + gupzz(ii,jj,kk)*Axz(ig,jg,kg)*Azz(ig,jg,kg) + &
                         gupxy(ii,jj,kk)*(Axx(ig,jg,kg)*Ayz(ig,jg,kg) + Axy(ig,jg,kg)*Axz(ig,jg,kg)) + &
                         gupxz(ii,jj,kk)*(Axx(ig,jg,kg)*Azz(ig,jg,kg) + Axz(ig,jg,kg)**2) + &
                         gupyz(ii,jj,kk)*(Axy(ig,jg,kg)*Azz(ig,jg,kg) + Axz(ig,jg,kg)*Ayz(ig,jg,kg))
         fyz(ii,jj,kk) = gupxx(ii,jj,kk)*Axy(ig,jg,kg)*Axz(ig,jg,kg) + gupyy(ii,jj,kk)*Ayy(ig,jg,kg)*Ayz(ig,jg,kg) + gupzz(ii,jj,kk)*Ayz(ig,jg,kg)*Azz(ig,jg,kg) + &
                         gupxy(ii,jj,kk)*(Axy(ig,jg,kg)*Ayz(ig,jg,kg) + Ayy(ig,jg,kg)*Axz(ig,jg,kg)) + &
                         gupxz(ii,jj,kk)*(Axy(ig,jg,kg)*Azz(ig,jg,kg) + Ayz(ig,jg,kg)*Axz(ig,jg,kg)) + &
                         gupyz(ii,jj,kk)*(Ayy(ig,jg,kg)*Azz(ig,jg,kg) + Ayz(ig,jg,kg)**2)
         
         ! Update A_ij RHS with trK and A^2 terms and Lie derivative terms (shift part)
         Axx_rhs(ig,jg,kg) = chin1*Axx_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Axx(ig,jg,kg) - TWO*fxx(ii,jj,kk)) + &
                             TWO*(Axx(ig,jg,kg)*betaxx(ii,jj,kk) + Axy(ig,jg,kg)*betayx(ii,jj,kk) + Axz(ig,jg,kg)*betazx(ii,jj,kk)) - F2o3*Axx(ig,jg,kg)*div_beta
         Ayy_rhs(ig,jg,kg) = chin1*Ayy_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Ayy(ig,jg,kg) - TWO*fyy(ii,jj,kk)) + &
                             TWO*(Axy(ig,jg,kg)*betaxy(ii,jj,kk) + Ayy(ig,jg,kg)*betayy(ii,jj,kk) + Ayz(ig,jg,kg)*betazy(ii,jj,kk)) - F2o3*Ayy(ig,jg,kg)*div_beta
         Azz_rhs(ig,jg,kg) = chin1*Azz_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Azz(ig,jg,kg) - TWO*fzz(ii,jj,kk)) + &
                             TWO*(Axz(ig,jg,kg)*betaxz(ii,jj,kk) + Ayz(ig,jg,kg)*betayz(ii,jj,kk) + Azz(ig,jg,kg)*betazz(ii,jj,kk)) - F2o3*Azz(ig,jg,kg)*div_beta
         Axy_rhs(ig,jg,kg) = chin1*Axy_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Axy(ig,jg,kg) - TWO*fxy(ii,jj,kk)) + &
                             Axx(ig,jg,kg)*betaxy(ii,jj,kk) + Axz(ig,jg,kg)*betazy(ii,jj,kk) + &
                             Ayy(ig,jg,kg)*betayx(ii,jj,kk) + Ayz(ig,jg,kg)*betazx(ii,jj,kk) + &
                             F1o3*Axy(ig,jg,kg)*div_beta - Axy(ig,jg,kg)*betazz(ii,jj,kk)
         Axz_rhs(ig,jg,kg) = chin1*Axz_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Axz(ig,jg,kg) - TWO*fxz(ii,jj,kk)) + &
                             Axx(ig,jg,kg)*betaxz(ii,jj,kk) + Axy(ig,jg,kg)*betayz(ii,jj,kk) + &
                             Ayz(ig,jg,kg)*betayx(ii,jj,kk) + Azz(ig,jg,kg)*betazx(ii,jj,kk) + &
                             F1o3*Axz(ig,jg,kg)*div_beta - Axz(ig,jg,kg)*betayy(ii,jj,kk)
         Ayz_rhs(ig,jg,kg) = chin1*Ayz_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Ayz(ig,jg,kg) - TWO*fyz(ii,jj,kk)) + &
                             Axy(ig,jg,kg)*betaxz(ii,jj,kk) + Ayy(ig,jg,kg)*betayz(ii,jj,kk) + &
                             Axz(ig,jg,kg)*betaxy(ii,jj,kk) + Azz(ig,jg,kg)*betazy(ii,jj,kk) + &
                             F1o3*Ayz(ig,jg,kg)*div_beta - Ayz(ig,jg,kg)*betaxx(ii,jj,kk)
                             
         ! Compute trace of S_ij (f is chin1 here)
         S_val = chin1 * (gupxx(ii,jj,kk)*Sxx(ig,jg,kg) + gupyy(ii,jj,kk)*Syy(ig,jg,kg) + gupzz(ii,jj,kk)*Szz(ig,jg,kg) + &
                 TWO*(gupxy(ii,jj,kk)*Sxy(ig,jg,kg) + gupxz(ii,jj,kk)*Sxz(ig,jg,kg) + gupyz(ii,jj,kk)*Syz(ig,jg,kg)))
         
         ! Final trK RHS (Reuse fxx as A^2 trace)
         ! fxx (trace of A^2) = g^ij g^kl A_ik A_jl = f (computed earlier)
         ! Wait, I computed fxx...fzz as A^2 components.
         ! Trace is g^ij f_ij.
         f(ii,jj,kk) = gupxx(ii,jj,kk)*fxx(ii,jj,kk) + gupyy(ii,jj,kk)*fyy(ii,jj,kk) + gupzz(ii,jj,kk)*fzz(ii,jj,kk) + &
                       TWO*(gupxy(ii,jj,kk)*fxy(ii,jj,kk) + gupxz(ii,jj,kk)*fxz(ii,jj,kk) + gupyz(ii,jj,kk)*fyz(ii,jj,kk))
         
         trK_rhs(ig,jg,kg) = -trK_rhs(ig,jg,kg) + alpn1*(F1o3*trK(ig,jg,kg)**2 + f(ii,jj,kk) + FOUR*PI*(rho(ig,jg,kg)+S_val))
         
      end do; end do; end do

      ! === 10. Gauge Conditions and Output Copy ===
      !DIR$ SIMD
      do kk=1,nz; do jj=1,ny; do ii=1,nx
         ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
         
         alpn1 = Lap(ig,jg,kg) + ONE
         chin1 = chi(ig,jg,kg) + ONE
         Lap_rhs(ig,jg,kg) = -TWO*alpn1*trK(ig,jg,kg)
         
         ! Betax RHS (GAUGE == 2)
         betax_rhs(ig,jg,kg) = FF*dtSfx(ig,jg,kg)
         betay_rhs(ig,jg,kg) = FF*dtSfy(ig,jg,kg)
         betaz_rhs(ig,jg,kg) = FF*dtSfz(ig,jg,kg)
         
         ! dtSfx RHS (GAUGE == 2)
         ! 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
         
         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)
         dtSfz_rhs(ig,jg,kg) = Gamz_rhs(ig,jg,kg) - reta_val*dtSfz(ig,jg,kg)
         
         ! Copy Rxx_l back to global Rxx (required by interface)
         Rxx(ig,jg,kg) = Rxx_l(ii,jj,kk)
         Ryy(ig,jg,kg) = Ryy_l(ii,jj,kk)
         Rzz(ig,jg,kg) = Rzz_l(ii,jj,kk)
         Rxy(ig,jg,kg) = Rxy_l(ii,jj,kk)
         Rxz(ig,jg,kg) = Rxz_l(ii,jj,kk)
         Ryz(ig,jg,kg) = Ryz_l(ii,jj,kk)
      end do; end do; end do
      
      call apply_lopsided(dtSfx, dtSfx_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
      call apply_lopsided(dtSfy, dtSfy_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
      call apply_lopsided(dtSfz, dtSfz_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
      
      call apply_lopsided(betax, betax_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
      call apply_lopsided(betay, betay_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
      call apply_lopsided(betaz, betaz_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
      
      ! Apply dissipation to gauge vars
      if (eps > 0.0d0) then
         call apply_kodiss(dtSfx, dtSfx_rhs, is, ie, js, je, ks, ke)
         call apply_kodiss(dtSfy, dtSfy_rhs, is, ie, js, je, ks, ke)
         call apply_kodiss(dtSfz, dtSfz_rhs, is, ie, js, je, ks, ke)
         call apply_kodiss(betax, betax_rhs, is, ie, js, je, ks, ke)
         call apply_kodiss(betay, betay_rhs, is, ie, js, je, ks, ke)
         call apply_kodiss(betaz, betaz_rhs, is, ie, js, je, ks, ke)
      end if

    end subroutine compute_block_kernel

    !-------------------------------------------------------------------------
    ! Helper: Safe 1st Derivative (4th Order)
    !-------------------------------------------------------------------------
    subroutine calc_derivs(f, fx, fy, fz, sym1, sym2, sym3, is, ie, js, je, ks, ke)
      real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f
      real*8, dimension(BLK,BLK,BLK), intent(out) :: fx, fy, fz
      real*8, intent(in) :: sym1, sym2, sym3
      integer, intent(in) :: is, ie, js, je, ks, ke
      
      integer :: ii, jj, kk, ig, jg, kg
      real*8 :: d12dx, d12dy, d12dz
      real*8 :: vM2, vM1, vP1, vP2
      
      d12dx = ONE / (12.d0 * (X(2)-X(1)))
      d12dy = ONE / (12.d0 * (Y(2)-Y(1)))
      d12dz = ONE / (12.d0 * (Z(2)-Z(1)))
      
      !DIR$ SIMD
      do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1
         ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1
         
         ! X-Derivative
         if (ig > 2 .and. ig < ex(1)-1) then
            fx(ii,jj,kk) = d12dx * (f(ig-2,jg,kg) - EIGHT*f(ig-1,jg,kg) + EIGHT*f(ig+1,jg,kg) - f(ig+2,jg,kg))
         else
            vM2 = get_val(f, ig-2, jg, kg, sym1, 1); vM1 = get_val(f, ig-1, jg, kg, sym1, 1)
            vP1 = get_val(f, ig+1, jg, kg, sym1, 1); vP2 = get_val(f, ig+2, jg, kg, sym1, 1)
            fx(ii,jj,kk) = d12dx * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2)
         end if
         
         ! Y-Derivative
         if (jg > 2 .and. jg < ex(2)-1) then
            fy(ii,jj,kk) = d12dy * (f(ig,jg-2,kg) - EIGHT*f(ig,jg-1,kg) + EIGHT*f(ig,jg+1,kg) - f(ig,jg+2,kg))
         else
            vM2 = get_val(f, ig, jg-2, kg, sym2, 2); vM1 = get_val(f, ig, jg-1, kg, sym2, 2)
            vP1 = get_val(f, ig, jg+1, kg, sym2, 2); vP2 = get_val(f, ig, jg+2, kg, sym2, 2)
            fy(ii,jj,kk) = d12dy * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2)
         end if
         
         ! Z-Derivative
         if (kg > 2 .and. kg < ex(3)-1) then
            fz(ii,jj,kk) = d12dz * (f(ig,jg,kg-2) - EIGHT*f(ig,jg,kg-1) + EIGHT*f(ig,jg,kg+1) - f(ig,jg,kg+2))
         else
            vM2 = get_val(f, ig, jg, kg-2, sym3, 3); vM1 = get_val(f, ig, jg, kg-1, sym3, 3)
            vP1 = get_val(f, ig, jg, kg+1, sym3, 3); vP2 = get_val(f, ig, jg, kg+2, sym3, 3)
            fz(ii,jj,kk) = d12dz * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2)
         end if
      end do; end do; end do
    end subroutine calc_derivs

    !-------------------------------------------------------------------------
    ! Helper: Safe 2nd Derivative (Mixed & Diag)
    !-------------------------------------------------------------------------
    subroutine calc_dderivs(f, fxx, fxy, fxz, fyy, fyz, fzz, sym1, sym2, sym3, is, ie, js, je, ks, ke)
      real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f
      real*8, dimension(BLK,BLK,BLK), intent(out) :: fxx, fxy, fxz, fyy, fyz, fzz
      real*8, intent(in) :: sym1, sym2, sym3
      integer, intent(in) :: is, ie, js, je, ks, ke
      
      integer :: ii, jj, kk, ig, jg, kg
      real*8 :: Fdxdx, Fdydy, Fdzdz, Fdxdy, Fdxdz, Fdydz
      real*8 :: vM2, vM1, vP1, vP2
      
      Fdxdx = ONE / (12.d0 * (X(2)-X(1))**2)
      Fdydy = ONE / (12.d0 * (Y(2)-Y(1))**2)
      Fdzdz = ONE / (12.d0 * (Z(2)-Z(1))**2)
      Fdxdy = ONE / (144.d0 * (X(2)-X(1))*(Y(2)-Y(1)))
      Fdxdz = ONE / (144.d0 * (X(2)-X(1))*(Z(2)-Z(1)))
      Fdydz = ONE / (144.d0 * (Y(2)-Y(1))*(Z(2)-Z(1)))
      
      !DIR$ SIMD
      do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1
         ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1
         
         ! FXX
         if (ig > 2 .and. ig < ex(1)-1) then
            fxx(ii,jj,kk) = Fdxdx * (-f(ig-2,jg,kg) + F16*f(ig-1,jg,kg) - F30*f(ig,jg,kg) + F16*f(ig+1,jg,kg) - f(ig+2,jg,kg))
         else
            vM2 = get_val(f, ig-2, jg, kg, sym1, 1); vM1 = get_val(f, ig-1, jg, kg, sym1, 1)
            vP1 = get_val(f, ig+1, jg, kg, sym1, 1); vP2 = get_val(f, ig+2, jg, kg, sym1, 1)
            fxx(ii,jj,kk) = Fdxdx * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2)
         endif
         
         ! FYY
         if (jg > 2 .and. jg < ex(2)-1) then
            fyy(ii,jj,kk) = Fdydy * (-f(ig,jg-2,kg) + F16*f(ig,jg-1,kg) - F30*f(ig,jg,kg) + F16*f(ig,jg+1,kg) - f(ig,jg+2,kg))
         else
            vM2 = get_val(f, ig, jg-2, kg, sym2, 2); vM1 = get_val(f, ig, jg-1, kg, sym2, 2)
            vP1 = get_val(f, ig, jg+1, kg, sym2, 2); vP2 = get_val(f, ig, jg+2, kg, sym2, 2)
            fyy(ii,jj,kk) = Fdydy * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2)
         endif
         
         ! FZZ
         if (kg > 2 .and. kg < ex(3)-1) then
            fzz(ii,jj,kk) = Fdzdz * (-f(ig,jg,kg-2) + F16*f(ig,jg,kg-1) - F30*f(ig,jg,kg) + F16*f(ig,jg,kg+1) - f(ig,jg,kg+2))
         else
            vM2 = get_val(f, ig, jg, kg-2, sym3, 3); vM1 = get_val(f, ig, jg, kg-1, sym3, 3)
            vP1 = get_val(f, ig, jg, kg+1, sym3, 3); vP2 = get_val(f, ig, jg, kg+2, sym3, 3)
            fzz(ii,jj,kk) = Fdzdz * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2)
         endif
         
         ! FXY
         if (ig > 2 .and. ig < ex(1)-1 .and. jg > 2 .and. jg < ex(2)-1) then
            fxy(ii,jj,kk) = Fdxdy * ( (f(ig-2,jg-2,kg)-EIGHT*f(ig-1,jg-2,kg)+EIGHT*f(ig+1,jg-2,kg)-f(ig+2,jg-2,kg)) &
                                -EIGHT*(f(ig-2,jg-1,kg)-EIGHT*f(ig-1,jg-1,kg)+EIGHT*f(ig+1,jg-1,kg)-f(ig+2,jg-1,kg)) &
                                +EIGHT*(f(ig-2,jg+1,kg)-EIGHT*f(ig-1,jg+1,kg)+EIGHT*f(ig+1,jg+1,kg)-f(ig+2,jg+1,kg)) &
                                -      (f(ig-2,jg+2,kg)-EIGHT*f(ig-1,jg+2,kg)+EIGHT*f(ig+1,jg+2,kg)-f(ig+2,jg+2,kg)) )
         else
            ! Simplify boundary mixed terms to 2nd order for robustness
            fxy(ii,jj,kk) = ZEO 
         endif
         
         ! FXZ
         if (ig > 2 .and. ig < ex(1)-1 .and. kg > 2 .and. kg < ex(3)-1) then
            fxz(ii,jj,kk) = Fdxdz * ( (f(ig-2,jg,kg-2)-EIGHT*f(ig-1,jg,kg-2)+EIGHT*f(ig+1,jg,kg-2)-f(ig+2,jg,kg-2)) &
                                -EIGHT*(f(ig-2,jg,kg-1)-EIGHT*f(ig-1,jg,kg-1)+EIGHT*f(ig+1,jg,kg-1)-f(ig+2,jg,kg-1)) &
                                +EIGHT*(f(ig-2,jg,kg+1)-EIGHT*f(ig-1,jg,kg+1)+EIGHT*f(ig+1,jg,kg+1)-f(ig+2,jg,kg+1)) &
                                -      (f(ig-2,jg,kg+2)-EIGHT*f(ig-1,jg,kg+2)+EIGHT*f(ig+1,jg,kg+2)-f(ig+2,jg,kg+2)) )
         else
            fxz(ii,jj,kk) = ZEO
         endif
         
         ! FYZ
         if (jg > 2 .and. jg < ex(2)-1 .and. kg > 2 .and. kg < ex(3)-1) then
            fyz(ii,jj,kk) = Fdydz * ( (f(ig,jg-2,kg-2)-EIGHT*f(ig,jg-1,kg-2)+EIGHT*f(ig,jg+1,kg-2)-f(ig,jg+2,kg-2)) &
                                -EIGHT*(f(ig,jg-2,kg-1)-EIGHT*f(ig,jg-1,kg-1)+EIGHT*f(ig,jg+1,kg-1)-f(ig,jg+2,kg-1)) &
                                +EIGHT*(f(ig,jg-2,kg+1)-EIGHT*f(ig,jg-1,kg+1)+EIGHT*f(ig,jg+1,kg+1)-f(ig,jg+2,kg+1)) &
                                -      (f(ig,jg-2,kg+2)-EIGHT*f(ig,jg-1,kg+2)+EIGHT*f(ig,jg+1,kg+2)-f(ig,jg+2,kg+2)) )
         else
            fyz(ii,jj,kk) = ZEO
         endif
      end do; end do; end do
    end subroutine calc_dderivs

    !-------------------------------------------------------------------------
    ! Helper: Lopsided Derivative
    !-------------------------------------------------------------------------
    subroutine apply_lopsided(f, f_rhs, bx, by, bz, is, ie, js, je, ks, ke)
      real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f, bx, by, bz
      real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: f_rhs
      integer, intent(in) :: is, ie, js, je, ks, ke
      
      integer :: ii, jj, kk, ig, jg, kg
      real*8 :: d12dx, d12dy, d12dz
      
      d12dx = ONE / (12.d0 * (X(2)-X(1)))
      d12dy = ONE / (12.d0 * (Y(2)-Y(1)))
      d12dz = ONE / (12.d0 * (Z(2)-Z(1)))
      
      !DIR$ SIMD
      do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1
         ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1
         
         ! X-Direction
         if (bx(ii,jj,kk) > ZEO) then
            if (ig > 3 .and. ig < ex(1)) then
               f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + bx(ii,jj,kk)*d12dx*(-F3*f(ig-1,jg,kg)-F10*f(ig,jg,kg)+F18*f(ig+1,jg,kg)-F6*f(ig+2,jg,kg)+f(ig+3,jg,kg))
            endif
         elseif (bx(ii,jj,kk) < ZEO) then
            if (ig > 0 .and. ig < ex(1)-3) then
               f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) - bx(ii,jj,kk)*d12dx*(-F3*f(ig+1,jg,kg)-F10*f(ig,jg,kg)+F18*f(ig-1,jg,kg)-F6*f(ig-2,jg,kg)+f(ig-3,jg,kg))
            endif
         endif
         
         ! Y-Direction
         if (by(ii,jj,kk) > ZEO) then
            if (jg > 3 .and. jg < ex(2)) then
               f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + by(ii,jj,kk)*d12dy*(-F3*f(ig,jg-1,kg)-F10*f(ig,jg,kg)+F18*f(ig,jg+1,kg)-F6*f(ig,jg+2,kg)+f(ig,jg+3,kg))
            endif
         elseif (by(ii,jj,kk) < ZEO) then
            if (jg > 0 .and. jg < ex(2)-3) then
               f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) - by(ii,jj,kk)*d12dy*(-F3*f(ig,jg+1,kg)-F10*f(ig,jg,kg)+F18*f(ig,jg-1,kg)-F6*f(ig,jg-2,kg)+f(ig,jg-3,kg))
            endif
         endif
         
         ! Z-Direction
         if (bz(ii,jj,kk) > ZEO) then
            if (kg > 3 .and. kg < ex(3)) then
               f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + bz(ii,jj,kk)*d12dz*(-F3*f(ig,jg,kg-1)-F10*f(ig,jg,kg)+F18*f(ig,jg,kg+1)-F6*f(ig,jg,kg+2)+f(ig,jg,kg+3))
            endif
         elseif (bz(ii,jj,kk) < ZEO) then
            if (kg > 0 .and. kg < ex(3)-3) then
               f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) - bz(ii,jj,kk)*d12dz*(-F3*f(ig,jg,kg+1)-F10*f(ig,jg,kg)+F18*f(ig,jg,kg-1)-F6*f(ig,jg,kg-2)+f(ig,jg,kg-3))
            endif
         endif
      end do; end do; end do
    end subroutine apply_lopsided

    !-------------------------------------------------------------------------
    ! Helper: Dissipation
    !-------------------------------------------------------------------------
    subroutine apply_kodiss(f, f_rhs, is, ie, js, je, ks, ke)
      real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f
      real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: f_rhs
      integer, intent(in) :: is, ie, js, je, ks, ke
      
      integer :: ii, jj, kk, ig, jg, kg
      real*8 :: dX, dY, dZ
      
      dX=X(2)-X(1); dY=Y(2)-Y(1); dZ=Z(2)-Z(1)
      
      !DIR$ SIMD
      do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1
         ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1
         
         if (ig > 3 .and. ig < ex(1)-3 .and. jg > 3 .and. jg < ex(2)-3 .and. kg > 3 .and. kg < ex(3)-3) then
            ! X
            f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + eps/dX/F64 * ((f(ig-3,jg,kg)+f(ig+3,jg,kg)) - SIX*(f(ig-2,jg,kg)+f(ig+2,jg,kg)) + FIT*(f(ig-1,jg,kg)+f(ig+1,jg,kg)) - TWT*f(ig,jg,kg))
            ! Y
            f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + eps/dY/F64 * ((f(ig,jg-3,kg)+f(ig,jg+3,kg)) - SIX*(f(ig,jg-2,kg)+f(ig,jg+2,kg)) + FIT*(f(ig,jg-1,kg)+f(ig,jg+1,kg)) - TWT*f(ig,jg,kg))
            ! Z
            f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + eps/dZ/F64 * ((f(ig,jg,kg-3)+f(ig,jg,kg+3)) - SIX*(f(ig,jg,kg-2)+f(ig,jg,kg+2)) + FIT*(f(ig,jg,kg-1)+f(ig,jg,kg+1)) - TWT*f(ig,jg,kg))
         endif
      end do; end do; end do
    end subroutine apply_kodiss

    !-------------------------------------------------------------------------
    ! Helper: Boundary Value Getter
    !-------------------------------------------------------------------------
    function get_val(f, i, j, k, sym, dir) result(val)
      real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f
      integer, intent(in) :: i, j, k, dir
      real*8, intent(in) :: sym
      real*8 :: val
      
      integer :: ii, jj, kk
      
      ii = i; jj = j; kk = k
      val = 0.0d0
      
      ! Handle Lower Boundary
      if (dir == 1 .and. i < 1) then
         ii = 1 + (1 - i); val = sym
      elseif (dir == 2 .and. j < 1) then
         jj = 1 + (1 - j); val = sym
      elseif (dir == 3 .and. k < 1) then
         kk = 1 + (1 - k); val = sym
      else
         val = 1.0d0 ! No symmetry flip needed
      endif
      
      ! Basic bounds check
      if (ii > 0 .and. ii <= ex(1) .and. jj > 0 .and. jj <= ex(2) .and. kk > 0 .and. kk <= ex(3)) then
         val = val * f(ii, jj, kk)
      else
         val = 0.0d0 ! Fallback
      endif
    end function get_val

  end subroutine compute_rhs_bssn_opt

end module bssn_rhs_opt_mod