anchor.diff:         rpar(iKph) = kappa
   
          dlnT4dlnP(1) = 0.75_dp*kappa*P/grav/exp(lnT4(1))
  +       
  +    contains
  +        subroutine get_rho_from_PT(ierr)
  +            use num_lib, only: safe_root_without_brackets
  +            integer, intent(out) :: ierr
  +            lnrho = safe_root_without_brackets(eval_pressure,rho_guess,dlnrho,max_newt, &
  +            &   imax, epsx,epsy,lrpar,rpar,lipar,ipar,ierr)
  +            
  +        end subroutine get_rho_from_PT
  +
       end subroutine deriv
   
       real(dp) function eval_pressure(lnrho, chi_rho, lrpar, rpar, lipar, ipar, ierr)
  
positive.code:         maximum_iterations = default_maximum_iterations_photosphere
          eps_lnrho = default_tolerance_photosphere_lnrho
          eps_ph = default_tolerance_photosphere_condition
          fallback_Pphoto = 2.0_dp*onethird*arad*Teff**4   
  
          ! scale tolerance to a thomson scaterring atmosphere
          sigma_Th = 8.0_dp*onethird*pi*(electroncharge**2/Melectron/clight2)**2
  
negative.code: module bc09
      use constants_def, only: dp
      
      ! data structure for the parameter arrays
      ! moments of the composition
  