anchor.diff:    subroutine curve(eps, temperature, S_min, S_max, temperature_c, sigma_c)
       implicit none
   
  -    real(kind = x_precision)                  ,intent(in)  :: eps            ! Precision required for the dichotomy
  -    real(kind = x_precision)                  ,intent(in)  :: S_min          ! Minimum surface density limit
  -    real(kind = x_precision)                  ,intent(in)  :: S_max          ! Maximum surface density limit
  -
  -    real(kind = x_precision),dimension(nb_it) ,intent(out) :: temperature    ! Temperatures for which we need a solution
  -    real(kind = x_precision),dimension(n_cell),intent(out) :: temperature_c  ! Temperature for the critical point
  -    real(kind = x_precision),dimension(n_cell),intent(out) :: sigma_c        ! Surface density for the critical point
  +    real(kind = x_precision),                    intent(in)  :: eps            ! Precision required for the dichotomy
  +    real(kind = x_precision),                    intent(in)  :: S_min          ! Minimum surface density limit
  +    real(kind = x_precision),                    intent(in)  :: S_max          ! Maximum surface density limit
  +
  +    real(kind = x_precision), dimension(nb_it) , intent(out) :: temperature    ! Temperatures for which we need a solution
  +    real(kind = x_precision), dimension(n_cell), intent(out) :: temperature_c  ! Temperature for the critical point
  +    real(kind = x_precision), dimension(n_cell), intent(out) :: sigma_c        ! Surface density for the critical point
       !------------------------------------------------------------------------
  -    integer                                                :: i,j,k
  +    integer                                     :: i,j,k
   
       ! For variables
  -    real(kind = x_precision)                               :: f              ! Q+ - Q-
  -    real(kind = x_precision)                               :: tau_eff        ! Effective optical depth
  -    integer                                                :: optical_depth  ! Indicator for the optical thickness
  +    real(kind = x_precision)                    :: f              ! Q+ - Q-
  +    real(kind = x_precision)                    :: tau_eff        ! Effective optical depth
  +    integer                                     :: optical_depth  ! Indicator for the optical thickness
   
       ! For the separate curves given the thickness
  -    real(kind = x_precision),dimension(nb_it)              :: sigma_t_thick
  -    real(kind = x_precision),dimension(nb_it)              :: sigma_t_thin
  -    real(kind = x_precision),dimension(nb_it)              :: temp_real      ! Temperature for the S curve
  -    real(kind = x_precision),dimension(nb_it)              :: sigma_real     ! Surface density for the S curve
  +    real(kind = x_precision), dimension(nb_it)  :: sigma_t_thick
  +    real(kind = x_precision), dimension(nb_it)  :: sigma_t_thin
  +    real(kind = x_precision), dimension(nb_it)  :: temp_real      ! Temperature for the S curve
  +    real(kind = x_precision), dimension(nb_it)  :: sigma_real     ! Surface density for the S curve
   
       ! For the output file
  -    character(len = 8)                                     :: number_of_cell
  -    character(len = 64)                                    :: fname_tot
  -    integer                                                :: fid_tot
  -    integer                                                :: ios
  +    character(len = 8)                          :: number_of_cell
  +    character(len = 64)                         :: fname_tot
  +    integer                                     :: fid_tot
  +    integer                                     :: ios
   
       ! For the critical points
  -    integer                                                :: index_fcp
  -    real(kind = x_precision)                               :: sigma_c_thick
  -    real(kind = x_precision)                               :: temp_c_thick
  -    real(kind = x_precision),dimension(n_cell)             :: sigma_thick
  -    real(kind = x_precision),dimension(n_cell)             :: temp_thick
  -    real(kind = x_precision),dimension(nb_it)              :: tau_thick      ! Optical thickness for the first critical point
  -    integer                                                :: index_scp
  -    real(kind = x_precision)                               :: sigma_c_thin
  -    real(kind = x_precision)                               :: temp_c_thin
  -    real(kind = x_precision),dimension(n_cell)             :: sigma_thin
  -    real(kind = x_precision),dimension(n_cell)             :: temp_thin
  -    real(kind = x_precision),dimension(nb_it)              :: tau_thin       ! Optical thickness for the second critical point
  +    integer                                     :: index_fcp
  +    real(kind = x_precision)                    :: sigma_c_thick
  +    real(kind = x_precision)                    :: temp_c_thick
  +    real(kind = x_precision), dimension(n_cell) :: sigma_thick
  +    real(kind = x_precision), dimension(n_cell) :: temp_thick
  +    real(kind = x_precision), dimension(nb_it)  :: tau_thick      ! Optical thickness for the first critical point
  +    integer                                     :: index_scp
  +    real(kind = x_precision)                    :: sigma_c_thin
  +    real(kind = x_precision)                    :: temp_c_thin
  +    real(kind = x_precision), dimension(n_cell) :: sigma_thin
  +    real(kind = x_precision), dimension(n_cell) :: temp_thin
  +    real(kind = x_precision), dimension(nb_it)  :: tau_thin       ! Optical thickness for the second critical point
   
       !------------------------------------------------------------------------
  -    ! Loop for n_cell values of x (radius)
  +    ! Loop over values of r/x (radius)
       !------------------------------------------------------------------------
  -
       do k = 1, n_cell
         ! For each point of the S curve
         do i = 1, nb_it
  
positive.code:        index_fcp, index_scp, sigma_c_thin, temp_c_thin)
  
        ! Combining the two separate curves into one that forms the S shape
        call build_s_curve(sigma_t_thick, sigma_t_thin, temperature, index_scp, sigma_real, temp_real)
  
        call display_critical_points(sigma_c_thin, temp_c_thin, sigma_c_thick, temp_c_thick, k)
  
  
negative.code:   !-------------------------------------------------------------------------
    ! Subroutine in order to build the S curve
    !-------------------------------------------------------------------------
    subroutine build_s_curve(sigma_real_thick, sigma_real_thin, index_scp, sigma_real)
      implicit none
  
      real(kind = x_precision), dimension(nb_it), intent(in)  :: sigma_real_thick
      real(kind = x_precision), dimension(nb_it), intent(in)  :: sigma_real_thin
      integer,                                    intent(in)  :: index_scp
  
      real(kind = x_precision), dimension(nb_it), intent(out) :: sigma_real
      !-----------------------------------------------------------------------
  
      integer :: i
  
      do i = 1, index_scp
        sigma_real(i) = sigma_real_thick(i)
      enddo
      do i = index_scp + 1, nb_it
        sigma_real(i) = sigma_real_thin(i)
      enddo
  
    end subroutine build_s_curve
  