anchor.diff:      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 * state_0%S_0
   
  -
       !-------------------------------------------------------------------------
       !Select the case for the optical depth to compute Fz
       !-------------------------------------------------------------------------
   
  -   ! if (tau_eff .ge. 1.)  then
  -   !    optical_depth     = 1
  -   ! else
  -   !    optical_depth     = 0
  -   ! end if
  -
       select case(optical_depth)
   
       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)
  +       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)
  
positive.code: 
    end function dichotomy
  
  
    !-------------------------------------------------------------------------
    !Subroutine in order to display variables
    !-------------------------------------------------------------------------
    subroutine display_variables(temp,Omega,r,sigma, H, rho, cs, nu, Q_plus, &
         Q_minus, K_ff, K_e, tau_eff, E_ff,Fz,f)
      implicit none
  
      real(kind = x_precision), intent(in)                     :: temp
      real(kind = x_precision), intent(in)                     :: Omega
      real(kind = x_precision), intent(in)                     :: r
      real(kind = x_precision), intent(in)                     :: sigma
      real(kind = x_precision), intent(in)                     :: H
      real(kind = x_precision), intent(in)                     :: rho
      real(kind = x_precision), intent(in)                     :: cs
      real(kind = x_precision), intent(in)                     :: nu
      real(kind = x_precision), intent(in)                     :: Q_plus
      real(kind = x_precision), intent(in)                     :: Q_minus
      real(kind = x_precision), intent(in)                     :: K_ff
      real(kind = x_precision), intent(in)                     :: K_e
      real(kind = x_precision), intent(in)                     :: tau_eff
      real(kind = x_precision), intent(in)                     :: f
      real(kind = x_precision), intent(in)                     :: E_ff
      real(kind = x_precision), intent(in)                     :: Fz
      !------------------------------------------------------------------------
  
      write(*,*)'               Variables                '
      write(*,*)'****************************************'
      write(*,"(' r           =',1p,E12.4)") r
      write(*,"(' Temp        =',1p,E12.4)") temp
      write(*,"(' Sigma       =',1p,E12.4)") sigma
      write(*,"(' Omega       =',1p,E12.4)") Omega
      write(*,"(' H           =',1p,E12.4)") H
      write(*,"(' rho         =',1p,E12.4)") rho
      write(*,"(' cs          =',1p,E12.4)") cs
      write(*,"(' nu          =',1p,E12.4)") nu
      write(*,"(' Q plus      =',1p,E12.4)") Q_plus
      write(*,"(' Q minus     =',1p,E12.4)") Q_minus
      write(*,"(' K_ff        =',1p,E12.4)") K_ff
      write(*,"(' K_e         =',1p,E12.4)") K_e
      write(*,"(' tau_eff     =',1p,E12.4)") tau_eff
      write(*,"(' E_ff        =',1p,E12.4)") E_ff
      write(*,"(' delta Q     =',1p,E12.4)") f
      write(*,"(' Fz          =',1p,E12.4)") Fz
      write(*,*)'****************************************'
    end subroutine display_variables
  
  
  end module mod_S_curve
  
negative.code:     coeff_c              = - params%RTM * T * Sigma * state_0%S_0 / (2._x_precision * state_0%H_0**2)
  
      call quadratic(coeff_a , coeff_b , coeff_c , H)
  
      rho                  = Sigma / H
      cs                   = Omega * H
  