anchor.diff:          use constants_def
           use nucchem_def, only : nuclide_not_found
       	use nucchem_lib, only : get_nuclide_index
  +    	use num_lib
           
           real(dp), intent(in) :: Teff,grav,tau
           real(dp), intent(out) :: Pphoto
  
positive.code:         real(dp) :: root_ph,rho,dfdrho
          integer, pointer :: ipar(:) => null() ! (lipar)
          real(dp), pointer :: rpar(:) => null()  ! (lrpar)
          
          ierr = 0
          allocate(ipar(number_photosphere_ipar), rpar(number_photosphere_rpar))
  
negative.code:         rpar(itau) = tau
          ipar(ihandle) = eos_handle
  
       	! set up rootfind
          imax = 20
          epsx = 1.0d-8
          epsy = 1.0d-8
          
          rhoph_guess = 1.d03
  
          !root_ph = photosphere(rho,dfdrho,lrpar,rpar,lipar,ipar,ierr)        
  		root_ph = safe_root_with_initial_guess(photosphere,rhoph_guess,x1,x3,y1,y3 &
              &   imax,epsx,epsy,lrpar,rpar,lipar,ipar,ierr)
              if (ierr /= 0) then
                  write(*,*) 'unable to converge', , x1, x3, y1, y3
                  cycle
              end if
          Pphoto = rpar(iPph)
          
          deallocate(ipar, rpar)
  