anchor.diff:      
       subroutine find_photospheric_pressure(Teff,grav,tau,Pphoto,eos_handle,ierr)
           use constants_def
  -        use nucchem_def, only : nuclide_not_found
  +        use nucchem_def, only : nuclide_not_found, nuclib
       	use nucchem_lib, only : get_nuclide_index
       	use num_lib
           integer, parameter :: default_maximum_iterations_photosphere = 20
  
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: 		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)
  