anchor.diff:      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
  +      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,&
  +  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
  -    integer,intent(in)                                   :: nb_it
   
       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
   
  
positive.code:     real(kind = x_precision),dimension(n_cell),intent(out)  :: s
  
      integer                                                :: i     = 0
      integer                                                :: k     = 0
      integer                                                :: l     = 0
  
      real(kind = x_precision)                               :: temp  = 0.0d0
      real(kind = x_precision), parameter                    :: t_min = 5.0d-1
      real(kind = x_precision), parameter                    :: t_max = 4.49d0
  
      integer, parameter                                     :: nb_it = 100
      real(kind = x_precision), parameter                    :: dt    = (t_max-t_min) / (nb_it-1)
  
      real(kind = x_precision)                               :: eps   = 1.0d-4
  
negative.code: 
      case(1)
  
        Fz = 2._x_precision * c**2 * T**4 /(27._x_precision * sqrt(3._x_precision) * (K_ff + K_e) * Sigma * state_0%S_0)
  
      case (0)
  
        Fz = 4._x_precision * state_0%H_0 * E_ff * H / (state_0%Omega_0 * state_0%S_0)
  
      end select
  
      Q_plus  = 3._x_precision  * state_0%H_0**2 * nu * (Omega * state_0%Omega_0)**2
      Q_minus = Fz / Sigma
  
      f = Q_plus - Q_minus
  
    end subroutine variables
  
  