anchor.diff:  ! and given a density cutoff, computes the total volume of all zones
   ! whose density is greater or equal to that density cutoff.
   ! We also impose a distance requirement so that we only look 
  -! at zones that are within the effective Roche radius of the white dwarf.
  +! at zones within the Roche lobe of the white dwarf.
   
   subroutine ca_volumeindensityboundary(rho,r_l1,r_l2,r_l3,r_h1,r_h2,r_h3,lo,hi,dx,&
                                         volp,vols,rho_cutoff)
   
     use bl_constants_module
  -  use probdata_module, only: mass_P, mass_S, com_P, com_S, roche_rad_P, roche_rad_S, single_star
  +  use probdata_module, only: mass_P, mass_S, com_P, com_S, single_star
     use prob_params_module, only: problo
   
     implicit none
  
positive.code: 
    use bl_constants_module, only: HALF, ZERO, ONE
    use prob_params_module, only: problo
    use probdata_module, only: mass_P, mass_S, com_P, com_S, single_star, roche_rad_P, roche_rad_S
  
    implicit none
  
  
negative.code:   double precision :: rho(r_l1:r_h1,r_l2:r_h2,r_l3:r_h3)
  
    integer          :: i, j, k
    double precision :: dV
    double precision :: r(3), r_P, r_S, grav_force_P, grav_force_S
  
    dV  = dx(1)*dx(2)*dx(3)
    volp = ZERO
    vols = ZERO
  
    do k = lo(3), hi(3)
       r(3) = problo(3) + (dble(k) + HALF) * dx(3)
       do j = lo(2), hi(2)
          r(2) = problo(2) + (dble(j) + HALF) * dx(2)
          do i = lo(1), hi(1)
             r(1) = problo(1) + (dble(i) + HALF) * dx(1)
  
             if (rho(i,j,k) > rho_cutoff) then
  
                r_P = sqrt(sum((r - com_p)**2))
                r_S = sqrt(sum((r - com_s)**2))
  
                grav_force_P = mass_P / r_P
                if (single_star .or. mass_S < 1.d-12 * mass_P) then
                   grav_force_S = ZERO
                else
                   grav_force_S = mass_S / r_S
                endif
                
                if (grav_force_P > grav_force_S) then
                   volp = volp + dV
                else
                   vols = vols + dV
                endif
               
  