anchor.diff:  100 format(1x, 5(g20.10))
   
     write (*, 100) pf%tm, kinetic_energy, internal_energy, -potential_energy, &
  -       kinetic_energy + internal_energy - potential_energy
  +       kinetic_energy + internal_energy - potential_energy, max_mach
   
   
     call destroy(pf)
  
positive.code:                          p(ii,jj,kk,rhoe_comp) * &
                           (dx(1)/rr)*(dx(2)/rr)*(dx(3)/rr)
  
                      imask(ii*r1:(ii+1)*r1-1, &
                            jj*r1:(jj+1)*r1-1, &
                            kk*r1:(kk+1)*r1-1) = .false.
                   end if
  
                end do
  
negative.code:   kinetic_energy = ZERO
    internal_energy = ZERO
    potential_energy = ZERO
    mach_number = ZERO
    max_mach = ZERO
    !Radius of the star. Model radius for secondary is 7.60156250e8, Model radius for primary is 9.69531250e8
    rad_star = 7.60156250e8
  
    imask(:,:,:) = .true.
  
  