anchor.diff:            M_enclosed(i) = M_enclosed(i-1) + &
                  FOUR3RD*M_PI*(xznr(i) - xznl(i))* &
                  (xznr(i)**2 +xznl(i)*xznr(i) + xznl(i)**2)*model_hse(i,idens_model)
  -
  -          cs_hse(i) = eos_state%cs
   
          enddo  ! End loop over zones
   
  
positive.code: 
      integer, parameter :: MAX_ITER = 250
  
      integer :: iter, iter_mass, num_iters
  
      integer :: icutoff
  
  
negative.code: 
      double precision :: temp_base, delta
  
      double precision :: xzn_hse(nx), xznl(nx), xznr(nx), M_enclosed(nx)
  
      integer :: i, n
  
  