anchor.diff:    subroutine curve()  
       implicit none
       integer                                                :: i     = 0
  +    integer                                                :: j     = 0
       real(kind = x_precision)                               :: temp  = 0.0d0
       real(kind = x_precision), parameter                    :: t_min = 1.2d-2
       real(kind = x_precision), parameter                    :: t_max = 1.0d0
  
positive.code:     implicit none
      real(kind=x_precision), intent(in)                     :: coeff_a,coeff_b,coeff_c
      real(kind=x_precision), intent(inout)                  :: sol
      real(kind=x_precision)                                 :: delta=0.
      delta          = coeff_b**2 - 4 * coeff_a * coeff_c
      if (delta .lt. 0.) then
        write(*,*)'No solutions in the R field.'
      else
        sol          = -0.5d0 *(coeff_b+sign(sqrt(delta),coeff_b))/coeff_a
      end if
  end subroutine quadratic
  
  
  
  subroutine initial_variables(rs, rmin, Mdot_0, Sigma_0)
    implicit none
  
    real(kind=x_precision),intent(inout)                      :: rs
    real(kind=x_precision),intent(inout)                      :: rmin
    real(kind=x_precision),intent(inout)                      :: Mdot_0
    real(kind=x_precision),intent(inout)                      :: Sigma_0
    type(parameters)                                          :: param  
    call get_parameters(param)
    rs                 = 2._x_precision * G * param%M/(c**2)
    rmin               = 3._x_precision * rs
    Mdot_0             = param%Mdot 
    Sigma_0            = Mdot_0 /(state_0%Omega_0 * rs**2 * 2 * pi)
    
  end subroutine initial_variables
  
  
  
  subroutine variables(T, Sigma, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e, tau_eff, &
       P_rad, P_gaz,E_ff,Fz,f,Sigma_0)
    implicit none
   
    type(parameters)                                         :: param
    real(kind = x_precision),intent(in)                      :: T,Sigma,Omega
    real(kind = x_precision)                                 :: coeff_a=0.,coeff_b=0.,coeff_c=0.
    real(kind = x_precision),intent(out)                     :: H     
    real(kind = x_precision),intent(out)                     :: rho 
    real(kind = x_precision),intent(out)                     :: cs    
    real(kind = x_precision),intent(out)                     :: nu  
    real(kind = x_precision),intent(out)                     :: Q_plus
    real(kind = x_precision),intent(out)                     :: K_ff 
    real(kind = x_precision),intent(out)                     :: K_e  
    real(kind = x_precision),intent(out)                     :: E_ff  
    real(kind = x_precision),intent(out)                     :: tau_eff
    real(kind = x_precision),intent(out)                     :: Fz  
    real(kind = x_precision),intent(out)                     :: Q_minus
    real(kind = x_precision),intent(out)                     :: P_rad 
    real(kind = x_precision),intent(out)                     :: P_gaz 
    real(kind = x_precision),intent(in)                      :: Sigma_0
    real(kind = x_precision),intent(out)                     :: f 
    integer                                                  :: optical_depth = 0
    
    call get_parameters(param)
    
  
    coeff_a              = (Omega**2 * state_0%Omega_0**2 * Sigma * Sigma_0)/2.
    coeff_b              = (-1._x_precision/3._x_precision) * cst_rad*T**4 * state_0%T_0**4
    coeff_c              = (-1._x_precision * kmp * T * state_0%T_0 * Sigma * Sigma_0)/(2._x_precision * mu)
    
    call quadratic(coeff_a,coeff_b,coeff_c,H)
    
    rho                  = Sigma/H
    cs                   = Omega * H
    nu                   = param%alpha * cs * H
    P_rad                = T**4
    P_gaz                = rho * T
    K_ff                 = 6.13d22 * state_0%rho_0 * rho * (state_0%T_0 * T)**(-3.5)
    K_e                  = 0.2_x_precision * (1_x_precision + param%X)
    E_ff                 = 6.22d20 * (state_0%rho_0 * rho)**2 * sqrt(state_0%T_0 * T)
    tau_eff              = 0.5_x_precision * sqrt(K_e * K_ff) * Sigma * Sigma_0    
    !-------------------------------------------------------------------------
    !Select the case for the opticaly depth
    !-------------------------------------------------------------------------
    if (tau_eff .ge. 1.)  then
       optical_depth     = 1
    else
       optical_depth     = 0
    end if
    !write(*,*)'tau',tau_eff
    select case(optical_depth)
    case(1)
       Fz = 4._x_precision * c**2 * T**4/(27. * sqrt(3._x_precision) * (K_ff + K_e) * Sigma * Sigma_0)
    case default
       Fz = 6.22d20 * 2._x_precision / state_0%Omega_0 * state_0%rho_0 *  H * rho**2 * sqrt(T*state_0%T_0)
    end select
    Q_minus             = 2._x_precision * Fz /Sigma
    Q_plus              = 9._x_precision /4._x_precision * nu * Omega**2
    f                   = Q_plus - Q_minus
    end subroutine variables
  
  
    !-------------------------------------------------------------------------
    ! Function in order to determine the change of sign in a given interval
    ! [Smin,Smax] with an epsilon precision
  
negative.code:     !-------------------------------------------------------------------------
      ! Test for 1 value of r
      !-------------------------------------------------------------------------
      do j              = 1 , 1
         !do j              = 1 , n_cell
         !  r = (param%rmax-rmin)/(n_cell-1)*(j-1) + rmin
         r              = 6*G*param%M/(c**2)     
         omega          = G*param%M/r**3
         write(number_of_cell,'(I5.5)') j
         fid = 20 + j
         fname = 'Temperature_Sigma_'//trim(number_of_cell)//'.dat'
         
         open(fid,file   = fname, status='unknown',action='readwrite') 
         write(fid,*)"T              Sigma"
         do i           = 1 , nb_it
            Smin        = 1.0d-5
            Smax        = 1.0d-0
            temp        = (t_max-t_min)/(nb_it-1)*(i-1) + t_min
            sigma       = dichotomy(Smin, Smax, eps, temp, omega)
            call variables(temp, sigma, omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
                 K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f,Sigma_0)
            call display_variables(temp,Omega,r, sigma, H, rho, cs, nu, Q_plus, Q_minus,&
                 K_ff, K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f)
            
            write(fid,'(1p,E12.6,4x,1p,E12.6)')temp,sigma
         enddo
         close(fid)
      enddo
    end subroutine curve
  
  
  