anchor.diff:         use prob_params_module,  only: center
          use bl_constants_module, only: ZERO, HALF, ONE, TWO
          use probdata_module,     only: do_initial_relaxation, relaxation_timescale
  -
  +       
          implicit none
   
          integer          :: lo(3),hi(3)
  
positive.code: 
         if (do_initial_relaxation) then
  
            damping_factor = ONE / relaxation_timescale
  
         else
  
            damping_factor = ZERO
            
         endif
  
         do k = lo(3), hi(3)
            r(3) = problo(3) + (dble(k) + HALF) * dx(3) - center(3)
  
            do j = lo(2), hi(2)
               r(2) = problo(2) + (dble(j) + HALF) * dx(2) - center(2)
  
               do i = lo(1), hi(1)
                  r(1) = problo(1) + (dble(i) + HALF) * dx(1) - center(1)
  
                     src(i,j,k,UMX:UMZ) = -new_state(i,j,k,UMX:UMZ) * damping_factor
  
                     ! Kinetic energy source term is v . momentum source term
  
                     src(i,j,k,UEDEN) = dot_product(new_state(i,j,k,UMX:UMZ) / new_state(i,j,k,URHO), src(i,j,k,UMX:UMZ))
  
               enddo
            enddo
         enddo
  
negative.code:        ! Local variables
  
         double precision :: damping_factor
         integer          :: i, j, k
         double precision :: new_mom(3), old_mom(3), new_ke, old_ke, rhoInv, dtInv
  
         ! Note that this function exists in a tiling region so we should only 
         ! modify the zones between lo and hi. 
  