anchor.diff:  
   
       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)
  +         tau_eff, P_rad, P_gaz,E_ff,Fz,f_min,Sigma_0, Omega_0,rs,T_0,rho_0, optical_depth)
   
       call variables(T, Smax, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e, &
  -         tau_eff, P_rad, P_gaz,E_ff,Fz,f_max,Sigma_0, Omega_0,rs,T_0, rho_0)
  +         tau_eff, P_rad, P_gaz,E_ff,Fz,f_max,Sigma_0, Omega_0,rs,T_0, rho_0, optical_depth)
   
       call variables(T, S_center, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
  -         K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f_center,Sigma_0, Omega_0,rs,T_0, rho_0)
  +         K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f_center,Sigma_0, Omega_0,rs,T_0, rho_0, optical_depth)
   
   
             if(f_min * f_center .gt. 0.) then
  
positive.code: 
      !Saving all the data
      fname = 'Temperature_Sigma_'//trim(number_of_cell)//'_r.dat'
         open(fid,file  = fname, status='unknown',action='readwrite')
         !  write(fid,*)"T     &         Sigma"
         do i          = 1, nb_it
            do j          = 1, nb_it
               do while( abs( temp_real_1(i) - temp_real_0(j) ) .gt. eps .AND. &
                          abs( sigma_real_1(i) - sigma_real_0(j) ) .gt. eps) then
  
                  write(fid,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real_1(i),temp_real_1(i)
  
                enddo
             enddo
          enddo
  
          do i          = j, nb_it
  
               write(fid,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real_0(i),temp_real_0(i)
  
          enddo
         close(fid)
  
      enddo
  
  
negative.code:   ! Dichotomic function in order to determine the change of sign in a given
    ! interval [Smin,Smax] with an epsilon precision
    !-------------------------------------------------------------------------
    real(kind=x_precision) function dichotomy(Smin, Smax, eps, T, omega, sigma_0, Omega_0,rs, T_0, rho_0, optical_depth)
      use mod_read_parameters
      use mod_constants
      use mod_variables
  