anchor.diff:      do i           = 1 , nb_it
          temp        = (t_max-t_min)/(nb_it-1)*(i-1) + t_min
          sigma       = dichotomy(Smin, Smax, eps, temp, omega)
  +       Smin=1.0e-3
  +       Smax=1.0d0
          
          write(18,*)temp," ",sigma
       enddo
  
positive.code:   end if
  
    select case(optical_depth)
    case(1:)
       Fz = 4._x_precision * c**2 * T**4/(27. * sqrt(3._x_precision) * (K_ff + K_e) * Sigma * Sigma_0)
    case default
       Fz = 6.22d20 * 2._x_precision / state_0%Omega_0 * state_0%rho_0 *  H * rho**2 * sqrt(T*state_0%T_0)
  
negative.code:   
  contains
    !-------------------------------------------------------------------------
    ! Subroutine in order to find the result of the equation Q+ - Q- = 0
    ! output : datafile with the temperature T and the surface density Sigma
    ! Range of surface density [Smin, Smax] and temperature [T_min, T_max]
    ! values should be changed in order to find more coherent results
    !-------------------------------------------------------------------------
    subroutine curve()  
      implicit none
      integer::i = 0
      real(kind = x_precision)                               :: temp = 0.0d0
      real(kind = x_precision), parameter                    :: t_min = 1.0d-3
      real(kind = x_precision), parameter                    :: t_max = 1.0d0
      integer,                  parameter                    :: nb_it = 10
      real(kind = x_precision)                               :: sigma = 0.0d0
      real(kind = x_precision)                               :: Smin = 1.0d-3
      real(kind = x_precision)                               :: Smax = 1.0d0
      real(kind = x_precision)                               :: eps = 1.0d-5
      real(kind = x_precision)                               :: omega = 0.0d0
      real(kind = x_precision)                               :: r = 0.0d0
      type(parameters)                                       :: param
  