anchor.diff:            do i = rloc(1,3), rloc(1,3) + 1
                if (i .ge. lo(1) .and. j .ge. lo(2) .and. k .ge. lo(3) .and. &
                    i .le. hi(1) .and. j .le. hi(2) .and. k .le. hi(3)) then
  +
                   omegasq = omegasq + c(i-rloc(1,3),j-rloc(2,3),k-rloc(3,3),3) &
  -                        * TWO * phi(i,j,k) / (sum(d_vector(:,3)**2)-sum(d_vector(:,2)**2))
  +                        * TWO * phi(i,j,k) / (sin(theta3)**2 * sum(d_vector(:,3)**2) - sin(theta2)**2 * sum(d_vector(:,2)**2))
  +
                endif
             enddo
          enddo
  
positive.code:     double precision :: r(3)
      double precision :: old_rho, drho
      double precision :: dV
      double precision :: omega
      double precision :: max_dist
  
      type (eos_t) :: eos_state, ambient_state
  
negative.code:              if (i .ge. lo(1) .and. j .ge. lo(2) .and. k .ge. lo(3) .and. &
                   i .le. hi(1) .and. j .le. hi(2) .and. k .le. hi(3)) then
                  bernoulli_2 = bernoulli_2 + c(i-rloc(1,2),j-rloc(2,2),k-rloc(3,2),2) &
                              * (phi(i,j,k) - HALF * sum(cross_product(omega, d_vector(:,2))**2))
               endif
            enddo
         enddo
  