anchor.diff:  		real(dp) :: eta, g14
           integer ::  size_tab ! = 4*size(Tb)
           real(dp), dimension(:), allocatable :: tabTb9, tabTeff, tabTeff6_4
  +        integer :: i
           
           ! make a very dense table of Tb(Teff); then interpolate to get Teff(Tb)        
           size_tab = 4*size(Tb)
  
positive.code:         real(dp) :: root_ph,rho,dfdrho
          integer, pointer :: ipar(:) => null() ! (lipar)
          real(dp), pointer :: rpar(:) => null()  ! (lrpar)
          integer :: i
          real(dp), parameter :: sigma_Th = 8.0_dp*onethird*pi*(electroncharge**2/Melectron/clight2)**2
          real(dp) :: fallback_Pphoto
          real(dp) :: rho_guess, kappa_Th, A, Z, Pgas
          real(dp) :: rho1, rho3, drho, y1, y3
  
negative.code:         ! compute dense table
          do i = 1, size_tab
              ! get Pph(Teff)
  !            call find_photospheric_pressure(Teff,grav,tau,Pphoto,eos_handle,ierr) 
  !		write(*,*) tabTeff_4(i), Pphoto           
          end do
          
          ! interpolate from dense table to get finished product
  