anchor.diff:          do i = lo(1), hi(1)
   
              r = position(i,j,k) - center
  -
  -           ! We account for symmetric boundaries in this sum as usual,
  -           ! by adding to the position the locations that would exist
  -           ! on the opposite side of the symmetric boundary.
  -
  -           r = merge(r + (problo - r), r, physbc_lo(:) .eq. Symmetry)
  -           r = merge(r + (r - probhi), r, physbc_hi(:) .eq. Symmetry)
   
              if (rho(i,j,k) > ZERO) then
                 rhoInv = ONE / rho(i,j,k)
  
positive.code:            vel(2) = ymom(i,j,k) * rhoInv
             vel(3) = zmom(i,j,k) * rhoInv
  
             ! Account for symmetric boundaries.
  
             vel = merge(ZERO, vel, physbc_lo(:) .eq. Symmetry)
             vel = merge(ZERO, vel, physbc_hi(:) .eq. Symmetry)
  
             vel = inertial_velocity(pos, vel, time)
  
             g(1) = gravx(i,j,k)
  
negative.code:                                         vol, vo_lo, vo_hi, &
                                          lo, hi, dx, time, Qtt) bind(C,name='quadrupole_tensor_double_dot')
  
    use bl_constants_module, only: ZERO, THIRD, HALF, ONE, TWO, M_PI
    use prob_params_module, only: center, problo, probhi, physbc_lo, physbc_hi, Symmetry, dim
    use wdmerger_util_module, only: inertial_rotation, inertial_velocity
    use castro_util_module, only: position
  
  