anchor.diff:  
       real(kind = x_precision),dimension(n_cell),intent(out) :: temperature    ! Temperature for the critical point
       real(kind = x_precision),dimension(n_cell),intent(out) :: s              ! Surface density for the critical point
  - !  __________________________________________________________________________
  -    integer                                                :: i     = 0
  -    integer                                                :: j     = 0
  -    integer                                                :: k     = 0
  -    integer                                                :: file_length = 1       ! Number of lines of the data saving file
  -
  - ! For the temperature and s
  +    !------------------------------------------------------------------------
  +    integer                                                :: i,j,k
  +    integer                                                :: file_length    ! Number of lines of the data saving file
  +
  +    ! For the temperature and s
       real(kind = x_precision)                               :: temp
       real(kind = x_precision)                               :: Smin
       real(kind = x_precision)                               :: Smax
  -    real(kind = x_precision)                               :: Sigma                 ! S
  -
  -    real(kind = x_precision)                               :: omega                 ! Angular velocity
  -    real(kind = x_precision)                               :: r                     ! Radius
  -    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)                               :: Sigma          ! S
  +
  +    real(kind = x_precision)                               :: omega          ! Angular velocity
  +    real(kind = x_precision)                               :: r              ! Radius
  +    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)              :: temp_t_thick
       real(kind = x_precision),dimension(nb_it)              :: sigma_t_thick
       real(kind = x_precision),dimension(nb_it)              :: temp_t_thin
       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
  +    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
   
  - ! For the critical points
  +    ! 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
  +    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)
  - !____________________________________________________________________________________
  +    real(kind = x_precision),dimension(nb_it)              :: tau_thin       ! Optical thickness for the second critical point
  +
  +    !------------------------------------------------------------------------
  +    ! Loop for n_cell values of x (radius)
  +    !------------------------------------------------------------------------
   
       do k = 1, n_cell
         file_length = 0
   
  -       r     = r_state%r(k)
  -       omega = x_state%Omega(k)
  -
  -       ! For each point of the S curve
  -       do i = 1, nb_it
  -          temp          = dt * (i-1) + t_min   ! Chosing the temperature value
  -
  -          ! Optical thick case (tau >= 1)
  -          optical_depth = 1
  -
  -          ! S range
  -          Smin          = S_min
  -          Smax          = S_max
  -
  -          ! S found with the dichotomy approach
  -          sigma         = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
  -
  -          temp_t_thick(i)  = temp
  -          sigma_t_thick(i) = sigma
  -
  -
  -          ! Optical thin case (tau < 1)
  -          optical_depth = 0
  -
  -          ! S range
  -          Smin          = S_min
  -          Smax          = S_max
  -
  -          ! S found with the dichotomy approach
  -          sigma         = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
  -
  -          temp_t_thin(i)   =  temp
  -          sigma_t_thin(i)  =  sigma
  -
  -        enddo
  -
  -       ! For the first critical point (the one one the right of the S shape)
  -       call first_critical_point(sigma_t_thick, temp_t_thick, index_fcp,sigma_c_thick, temp_c_thick, nb_it)
  -
  -       ! For the second critical point (the one one the left of the S shape)
  -       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)
  -
  -       ! Combining the two separate curves into one that forms the S shape
  -       call build_s_curve(sigma_t_thick, sigma_t_thin, temp_t_thick, nb_it, index_scp, sigma_real, temp_real)
  -
  -       call display_critical_points(sigma_c_thin, temp_c_thin,sigma_c_thick, temp_c_thick, k)
  +      r     = r_state%r(k)
  +      omega = x_state%Omega(k)
  +
  +      ! For each point of the S curve
  +      do i = 1, nb_it
  +        temp          = dt * (i-1) + t_min   ! Chosing the temperature value
  +
  +        ! Optical thick case (tau >= 1)
  +        optical_depth = 1
  +
  +        ! S range
  +        Smin          = S_min
  +        Smax          = S_max
  +
  +        ! S found with the dichotomy approach
  +        sigma         = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
  +
  +        temp_t_thick(i)  = temp
  +        sigma_t_thick(i) = sigma
  +
  +
  +        ! Optical thin case (tau < 1)
  +        optical_depth = 0
  +
  +        ! S range
  +        Smin          = S_min
  +        Smax          = S_max
  +
  +        ! S found with the dichotomy approach
  +        sigma         = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
  +
  +        temp_t_thin(i)   =  temp
  +        sigma_t_thin(i)  =  sigma
  +
  +      enddo
  +
  +      ! For the first critical point (the one one the right of the S shape)
  +      call first_critical_point(sigma_t_thick, temp_t_thick, index_fcp,sigma_c_thick, temp_c_thick, nb_it)
  +
  +      ! For the second critical point (the one one the left of the S shape)
  +      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)
  +
  +      ! Combining the two separate curves into one that forms the S shape
  +      call build_s_curve(sigma_t_thick, sigma_t_thin, temp_t_thick, nb_it, index_scp, sigma_real, temp_real)
  +
  +      call display_critical_points(sigma_c_thin, temp_c_thin,sigma_c_thick, temp_c_thick, k)
   
   
         ! Saving the DIMENSIONED values of S and T for the S curve in a file
  
positive.code: 
      real(kind = x_precision),dimension(nb_it),intent(out) :: sigma_real
      real(kind = x_precision),dimension(nb_it),intent(out) :: temp_real
   !____________________________________________________________________________________
      integer::i = 0
  
      do i = 1,index_scp
        sigma_real(i) = sigma_real_thick(i)
  
negative.code: 
          endif
  
          S_center = (Smin + Smax) / 2._x_precision
          j = j + 1
  
        end do iteration
  