anchor.diff:  
       delta          = coeff_b**2 - 4_x_precision * coeff_a * coeff_c
   
  +    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
  
positive.code:       !------------------------------------------------------------------------
        write(*,*)'**** Critical Point',k,'********'
        write(*,*)'****************************************'
        write(*,"(' Opticaly thin (T,sigma) :',1p,E12.4,4x,1p,E12.4)")temp_c_thin,sigma_c_thin
        write(*,"(' Opticaly thick (T,sigma):',1p,E12.4,4x,1p,E12.4)")temp_c_thick,sigma_c_thick
  
        write(*,*)'****************************************'
  
      end subroutine display_critical_points
  
  
  
      
  
  
  
negative.code:     ! eps -> Precision
      ! T-> Fixed variable
      !-------------------------------------------------------------------------
      S_center             = (Smin+Smax)/2.
      j = 0
      call variables(T, Smin, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e,&
           tau_eff, P_rad, P_gaz,E_ff,Fz,f_min,Sigma_0, Omega_0,rs,T_0,rho_0,optical_depth)
  