anchor.diff:  
       call variables(T, Smax, Omega, f_max, optical_depth)
   
  -  !   write(*,*)'fmin = ',f_min
  -  !   write(*,*)'fmax = ',f_max
  -
  -    if ( f_max * f_min .gt. 0.) then
  -   !    write(*,*)'This function image does not switch its sign in this particular interval.'
  -     !  dichotomy = 0
  -
  -    else if( f_max * f_min .lt. 0.) then
  -       iteration:do while ( dabs( Smax - Smin ) .ge. eps .and. j .lt. 10000)
  -
  -
  -    call variables(T, Smin, Omega, f_min, optical_depth)
  -
  -    call variables(T, Smax, Omega, f_max, optical_depth)
  -
  -    call variables(T, S_center, Omega, f_center, optical_depth)
  -
  -
  -          if(f_min * f_center .gt. 0.) then
  -
  -             Smin         = S_center
  -
  -          else if (f_max * f_center .gt. 0.) then
  -             Smax         = S_center
  -
  -          endif
  -
  -          S_center       = (Smin + Smax) * 1._x_precision / 2._x_precision
  -          j               = j + 1
  -
  -       end do iteration
  -         dichotomy = S_center
  +    !write(*,*)'fmin = ',f_min
  +    !write(*,*)'fmax = ',f_max
  +
  +    if ( f_max * f_min > 0.) then
  +      !write(*,*)'This function image does not switch its sign in this particular interval.'
  +      !dichotomy = 0
  +
  +    else if( f_max * f_min < 0.) then
  +      iteration:do while ( dabs( Smax - Smin ) >= eps .and. j < 10000)
  +
  +        call variables(T, Smin, Omega, f_min, optical_depth)
  +
  +        call variables(T, Smax, Omega, f_max, optical_depth)
  +
  +        call variables(T, S_center, Omega, f_center, optical_depth)
  +
  +        if (f_min * f_center > 0.) then
  +
  +          Smin = S_center
  +
  +        else if (f_max * f_center > 0.) then
  +
  +          Smax = S_center
  +
  +        endif
  +
  +        S_center = (Smin + Smax) * 1._x_precision / 2._x_precision
  +        j = j + 1
  +
  +      end do iteration
  +
  +      dichotomy = S_center
   
       endif
   
  -
     end function dichotomy
   
  -end module mod_S_curve
  +end module mod_s_curve
  
positive.code:     real(kind = x_precision)                     :: Q_minus
      !------------------------------------------------------------------------
  
      coeff_a              = (Omega * state_0%Omega_0)**2 * Sigma * state_0%S_0 / 2._x_precision
      coeff_b              = (-1._x_precision/3._x_precision) * cst_rad * (T * state_0%T_0)**4 / state_0%H_0
      coeff_c              = - params%RTM * T * Sigma * state_0%S_0 / (2._x_precision * state_0%H_0**2)
  
      call quadratic(coeff_a , coeff_b , coeff_c , H)
  
      rho                  = Sigma / H
      cs                   = Omega * H
      nu                   = params%alpha * cs * H
      K_ff                 = 6.13d22 * state_0%rho_0 * rho * (state_0%T_0 * T)**(-3.5_x_precision)
      K_e                  = params%kappa_e
      E_ff                 = 6.22d20 * (state_0%rho_0 * rho)**2 * sqrt(state_0%T_0 * T)
      tau_eff              = 0.5_x_precision * sqrt(K_e * K_ff) * Sigma * state_0%S_0
  
      !-------------------------------------------------------------------------
      !Select the case for the optical depth to compute Fz
  
negative.code:     real(kind = x_precision)                 ,intent(out):: sigma_c_thick
      real(kind = x_precision)                 ,intent(out):: temp_c_thick
      integer                                              :: i
  
      i = 1
      do while (sigma_real_thick(i) < sigma_real_thick(i+1) .and. i < nb_it - 1)
        index_fcp     = i
        sigma_c_thick = sigma_real_thick(i)
        i             = i + 1
      enddo
      temp_c_thick  = temp_real_thick(i)
  
    endsubroutine first_critical_point
  
  
    !-------------------------------------------------------------------------
    ! Subroutine in order to find the second critical point
    !-------------------------------------------------------------------------
    subroutine second_critical_point(sigma_real_thick, sigma_real_thin, temp_real_thick,&
           index_fcp, index_scp, sigma_c_thin, temp_c_thin, nb_it)
      implicit none
  
      real(kind = x_precision),dimension(nb_it),intent(in) :: sigma_real_thick
      real(kind = x_precision),dimension(nb_it),intent(in) :: temp_real_thick
      real(kind = x_precision),dimension(nb_it),intent(in) :: sigma_real_thin
  
      integer                                  ,intent(in) :: nb_it
      integer                                  ,intent(in) :: index_fcp
      integer                                  ,intent(out):: index_scp
  
  