anchor.diff:  
       real(kind = x_precision),dimension(nb_it)              :: temp_t_thick  =0.0d0
       real(kind = x_precision),dimension(nb_it)              :: sigma_t_thick =0.0d0
  -    real(kind = x_precision),dimension(nb_it)              :: temp_t_thin  =0.0d0
  -    real(kind = x_precision),dimension(nb_it)              :: sigma_t_thin =0.0d0
  -    real(kind = x_precision),dimension(nb_it)              :: temp_real =0.0d0
  -    real(kind = x_precision),dimension(nb_it)              :: sigma_real=0.0d0
  +    real(kind = x_precision),dimension(nb_it)              :: temp_t_thin   =0.0d0
  +    real(kind = x_precision),dimension(nb_it)              :: sigma_t_thin  =0.0d0
  +    real(kind = x_precision),dimension(nb_it)              :: temp_real     =0.0d0
  +    real(kind = x_precision),dimension(nb_it)              :: sigma_real    =0.0d0
   
       character(len = 8)                                     :: number_of_cell
  -    character(len = 64)                                    :: fname_1
  -    character(len = 64)                                    :: fname_2
  -    character(len = 64)                                    :: fname_3
  -
  -    integer                                                :: fid_1
  -    integer                                                :: fid_2
  -    integer                                                :: fid_3
  +    !character(len = 64)                                    :: fname_thick
  +    !character(len = 64)                                    :: fname_thin
  +    character(len = 64)                                    :: fname_tot
  +
  +    !integer                                                :: fid_thick
  +    !integer                                                :: fid_thin
  +    integer                                                :: fid_tot
   
       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(n_cell)             :: sigma_thick
  +    !real(kind = x_precision),dimension(n_cell)             :: temp_thick
       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(n_cell)             :: sigma_thin
  +    !real(kind = x_precision),dimension(n_cell)             :: temp_thin
   
       !-------------------------------------------------------------------------
       ! Test for 1 value of r
       !-------------------------------------------------------------------------
       do k              = 1 , n_cell
   
  -       !r = (rmax-rmin)/(n_cell-1)*(k-1) + rmin
  -       !r              = 10._x_precision*G*params%M/(c**2)
  -       r   = r_state%r(k)
  -
  -       omega          = x_state%Omega(k)
  -
  -       write(number_of_cell,'(I5.5)') k
  -       fid_1 = 20 + k
  -       fid_2 = 21 + k
  -       fid_3 = 22 + k
  -
  -
  -       fname_1 = 's_curves/Temperature_Sigma_'//trim(number_of_cell)//'_1.dat'
  -       fname_2 = 's_curves/Temperature_Sigma_'//trim(number_of_cell)//'_0.dat'
  -       fname_3 = 's_curves/Temperature_Sigma_'//trim(number_of_cell)//'_tot.dat'
  -
  -
  -       open(fid_1,file  = fname_1, status='unknown',action='readwrite')
  -       open(fid_2,file  = fname_2, status='unknown',action='readwrite')
  -       open(fid_3,file  = fname_3, status='unknown',action='readwrite')
  -
  -       do i          = 1, nb_it
  -
  -          Smin        = 1d1
  -          Smax        = 1d4
  -
  -          temp        = dt * (i-1) + t_min
  -          optical_depth = 1
  -
  -          sigma       = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
  -
  -          call variables(temp, sigma, omega, f, optical_depth)
  -
  -          temp_t_thick(i)   =  temp
  -          sigma_t_thick(i)  = sigma
  -
  -          !do j          = 1, nb_it
  -          !   write(fid_2,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real_1(j),temp_real_1(j)
  -          !enddo
  -
  -          Smin        = 1d1
  -          Smax        = 1d4
  -
  -          optical_depth = 0
  -
  -          sigma       = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
  -
  -          call variables(temp, sigma, omega, f, optical_depth)
  -
  -          temp_t_thin(i)   =  temp
  -          sigma_t_thin(i)  =  sigma
  -
  -
  -         ! do l          = 1, nb_it
  -         !    write(fid_1,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real_0(l),temp_real_0(l)
  -         ! enddo
  -
  -       enddo
  -
  -
  -       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)
  -
  -       call display_critical_points(sigma_c_thin, temp_c_thin,sigma_c_thick, temp_c_thick, k)
  -
  -
  -
  -          do l          = 1, nb_it
  -
  -             temp_real(l)  = log10( temp_real(l) * state_0%T_0 )
  -             sigma_real(l) = log10( sigma_real(l) * state_0%S_0 )
  -
  -             write(fid_3,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real(l),temp_real(l)
  -          enddo
  -
  -       close(fid_1)
  -       close(fid_2)
  -       close(fid_3)
  -
  -          temperature(k) = temp_c_thick
  -          s(k)           = sigma_c_thick
  -
  -
  -       !  temp_thick(k)  = log10( temp_c_thick * state_0%T_0 )
  -       !  sigma_thick(k) = log10( sigma_c_thick * state_0%S_0 )
  -       !  temp_thin(k)  = log10( temp_c_thin * state_0%T_0 )
  -       !  sigma_thin(k) = log10( sigma_c_thin * state_0%S_0 ) !PROBLEMATIC
  +      !r = (rmax-rmin)/(n_cell-1)*(k-1) + rmin
  +      !r              = 10._x_precision*G*params%M/(c**2)
  +      r   = r_state%r(k)
  +
  +      omega          = x_state%Omega(k)
  +
  +      write(number_of_cell,'(I5.5)') k
  +      !fid_thick = 20 + k
  +      !fid_thin = 21 + k
  +      fid_tot = 22 + k
  +
  +
  +      !fname_thick = 's_curves/Temperature_Sigma_'//trim(number_of_cell)//'_thick.dat'
  +      !fname_thin = 's_curves/Temperature_Sigma_'//trim(number_of_cell)//'_thin.dat'
  +      fname_tot = 's_curves/Temperature_Sigma_'//trim(number_of_cell)//'_tot.dat'
  +
  +
  +      !open(fid_thick,file  = fname_thick, status='unknown',action='readwrite')
  +      !open(fid_thin,file  = fname_thin, status='unknown',action='readwrite')
  +      open(fid_tot,file  = fname_tot, status='unknown',action='readwrite')
  +
  +      do i = 1, nb_it
  +
  +        temp          = dt * (i-1) + t_min
  +
  +        ! Optical thick case
  +
  +        Smin          = 1d1
  +        Smax          = 1d4
  +
  +        optical_depth = 1
  +
  +        sigma         = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
  +
  +        call variables(temp, sigma, omega, f, optical_depth)
  +
  +        temp_t_thick(i)  = temp
  +        sigma_t_thick(i) = sigma
  +
  +        ! Optical thin case
  +
  +        Smin          = 1d1
  +        Smax          = 1d4
  +
  +        optical_depth = 0
  +
  +        sigma         = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
  +
  +        call variables(temp, sigma, omega, f, optical_depth)
  +
  +        temp_t_thin(i)   =  temp
  +        sigma_t_thin(i)  =  sigma
  +
  +
  +        !do j = 1, nb_it
  +        !  write(fid_thick,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)') sigma_real_thick(j), temp_real_thick(j)
  +        !  write(fid_thin,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)') sigma_real_thin(j),temp_real_thin(j)
  +        !enddo
  +
  +      enddo
  +
  +
  +      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)
  +
  +      call display_critical_points(sigma_c_thin, temp_c_thin,sigma_c_thick, temp_c_thick, k)
  +
  +
  +
  +      do j = 1, nb_it
  +
  +        temp_real(l)  = log10( temp_real(l) * state_0%T_0 )
  +        sigma_real(l) = log10( sigma_real(l) * state_0%S_0 )
  +
  +        write(fid_tot,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)') sigma_real(l), temp_real(l)
  +      enddo
  +
  +      !close(fid_thick)
  +      !close(fid_thin)
  +      close(fid_tot)
  +
  +      temperature(k) = temp_c_thick
  +      s(k)           = sigma_c_thick
  +
  +
  +      !temp_thick(k)  = log10( temp_c_thick * state_0%T_0 )
  +      !sigma_thick(k) = log10( sigma_c_thick * state_0%S_0 )
  +      !temp_thin(k)   = log10( temp_c_thin * state_0%T_0 )
  +      !sigma_thin(k)  = log10( sigma_c_thin * state_0%S_0 ) !PROBLEMATIC
   
   
       enddo
   
  -     !  call write_critical_points(n_cell, r_state%r, rs, sigma_thin, temp_thin, sigma_thick, temp_thick)
  -
  +    !call write_critical_points(n_cell, r_state%r, rs, sigma_thin, temp_thin, sigma_thick, temp_thick)
   
     end subroutine curve
  -
   
   
     !-------------------------------------------------------------------------
  
positive.code:       continue
      else
  
  
      if (delta .lt. 0.) then
        write(*,*)'No solutions in the R field.'
  
     else
        sol_1         = -0.5_x_precision *(coeff_b + sign(sqrt(delta),coeff_b))/coeff_a
        if (sol_1 .lt. 0.d-12) then
           write(*,*)'Problem'
           stop
        endif
        sol_2 = (coeff_c / (coeff_a * sol_1))
        sol = max(sol_1,sol_2)
  
      end if
  
      end if
    end subroutine quadratic
  
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
  
  