anchor.diff:      real(kind = x_precision),dimension(nb_it),intent(in) :: temp_real_thick
       real(kind = x_precision),dimension(nb_it),intent(in) :: sigma_real_thin
   
  -    real(kind = x_precision)                 ,intent(in) :: eps2
  -
       integer                                  ,intent(in) :: index_fcp
       integer                                  ,intent(out):: index_scp
   
       real(kind = x_precision)                 ,intent(out):: sigma_c_thin
       real(kind = x_precision)                 ,intent(out):: temp_c_thin
  -    integer                                              :: i = 0
  -
  -    do i = index_fcp, nb_it - 1
  -
  -       if(dabs((sigma_real_thick(i) - sigma_real_thin(i+1) / sigma_real_thick(i))) .lt. eps2)then
  -          index_scp = i
  -          ! sigma_c_thin = (sigma_real_thick(i) + sigma_real_thin(i+1))/2._x_precision
  -          sigma_c_thin =  sigma_real_thin(i+1)
  -
  -          temp_c_thin = temp_real_thick(i)
  -       end if
  -    enddo
  +    integer                                              :: i
  +    i = max(1, index_fcp)
  +    
  +    ! The change occurs when the difference of sigma changes sign
  +    do while (sigma_real_thick(i) < sigma_real_thin(i) .and. i < nb_it)
  +       i = i + 1
  +    end do
  +
  +    index_scp = i
  +    sigma_c_thin =  sigma_real_thin(i+1)
  +    temp_c_thin = temp_real_thick(i)
   
       endsubroutine second_critical_point
   
   
   
  -    subroutine build_s_curve(sigma_real_thick, sigma_real_thin, temp_real_thick,nb_it, index_scp, sigma_real, temp_real)
  +    subroutine build_s_curve(sigma_real_thick, sigma_real_thin, temp_real_thick, nb_it, index_scp, sigma_real, temp_real)
         implicit none
         integer,intent(in)                                   :: nb_it
   
  
positive.code:     real(kind = x_precision), parameter                    :: dt    = (t_max-t_min) / (nb_it-1)
  
      real(kind = x_precision)                               :: eps   = 1.0d-4
      real(kind = x_precision)                               :: eps2   = 5.0d-1
  
      real(kind = x_precision)                               :: Sigma = 0.0d0
      real(kind = x_precision)                               :: Smin  = 0.0d0
  
negative.code:        call first_critical_point(sigma_t_thick, temp_t_thick, index_fcp,sigma_c_thick, temp_c_thick, nb_it)
  
         call second_critical_point(sigma_t_thick, sigma_t_thin, temp_t_thick,&
           index_fcp, index_scp, sigma_c_thin, temp_c_thin, nb_it)
  
         call build_s_curve(sigma_t_thick, sigma_t_thin, temp_t_thick, nb_it, index_scp, sigma_real, temp_real)
  
  