anchor.diff:        ! write(*,*)sigma_real_0(i), sigma_real_1(i)
        !  write(*,*)(dabs(sigma_real_1(i) - sigma_real_0(i+1)))
       end do
  -    
  +
           do i = index_fcp, nb_it - 1
   
              if(dabs(sigma_real_1(i) - sigma_real_0(i+1)) .lt. eps2)then
  
positive.code:     if (coeff_a == 0.0) then
        write(*,*)'Coefficient a in the quadratic equation is nought.'
        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
  
negative.code: 
         open(fid_4,file  = fname_4, status='unknown',action='write')
         !write(fid_4)'nb_cell,temp_c_thin,sigma_c_thin,temp_c_thick,sigma_thin'
  
            do i=1, n
               write(fid_4,'(E12.6,1p,E12.6,4x,1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sqrt(radius(i)/rs),&
                          temp_c_thin(i),sigma_c_thin(i),temp_c_thick(i),sigma_c_thick(i)
            end do
  
         close(fid_4)
  
      end subroutine write_critical_points
  
  
  
  
    !-------------------------------------------------------------------------
    !Subroutine for resolving a quadratic equation as long as the solutions are
    !a set of real numbers
  