anchor.diff:    use fundamental_constants_module, only: Gconst, c_light, parsec
     use probdata_module, only: gw_dist, axis_1, axis_2, axis_3
     use meth_params_module, only: rot_axis
  +  use prob_params_module, only: dim
   
     implicit none
   
  
positive.code: 
       h(:,:) = h(:,:) * TWO * Gconst / (c_light**4 * r)
  
       ! If rot_axis == 3, then h_+ = h_{11} = -h_{22} and h_x = h_{12} = h_{21}.
       ! Analogous statements hold along the other axes.
  
       ! We are adding here so that this calculation makes sense on multiple levels.
  
       if (dir .eq. axis_1) then
  
          h_plus_1  = h_plus_1  + h(axis_2,axis_2)
          h_cross_1 = h_cross_1 + h(axis_2,axis_3)
  
       else if (dir .eq. axis_2) then
  
          h_plus_2  = h_plus_2  + h(axis_3,axis_3)
          h_cross_2 = h_cross_2 + h(axis_3,axis_1)
  
       else if (dir .eq. axis_3) then
  
          h_plus_3  = h_plus_3  + h(axis_1,axis_1)
          h_cross_3 = h_cross_3 + h(axis_1,axis_2)
  
       endif
  
  
negative.code:         do i = lo(1), hi(1)
  
             r = position(i,j,k) - center
  
             if (rho(i,j,k) > ZERO) then
                rhoInv = ONE / rho(i,j,k)
  