anchor.diff:      Delta = b1**2 - (4._x_precision * a1 * c1)
       !--------------------------------------------------------
       do i=1,n_cell
  -      state_out%H(i)     = - 0.5_x_precision * (b1(i) + sign(sqrt(Delta(i)),b1(i))) / a1(i) 
  -      state_out%P_rad(i) = T_in(i)**4
  -      state_out%cs(i)    = Omega(i) * state_out%H(i)
  -      state_out%rho(i)   = S_in(i) / (state_out%H(i) * x(i))
  -      state_out%nu(i)    = para%alpha * state_out%cs(i) * state_out%H(i)
  -      state_out%P_gaz(i) = state_out%rho(i) * (T_in(i)**4)
  -      state_out%beta(i)  = state_out%P_gaz(i) / (state_out%P_gaz(i) + state_out%P_rad(i))
  -      state_out%Cv(i)    = (kmp/mu) * ((12._x_precision * gammag -1._x_precision ) &
  -                           * (1._x_precision - state_out%beta(i)) + state_out%beta(i)) &
  -                           / ( state_out%beta(i) * (gammag - 1._x_precision))
  -      !------------limit condition to compute v-------------
  -      if (i .eq. 1) then
  -        state_out%v(i)  = 0._x_precision
  -      else if (i .eq. n_cell) then
  -        state_out%v(i)  = - 1._x_precision / S_in(i) / x(i) 
  -      else
  -        state_out%v(i)  = - 1._x_precision / S_in(i) / x(i) * ( ((state_out%nu(i+1) &
  -                          * S_in(i+1)) - (state_out%nu(i-1) * S_in(i-1))) / (2._x_precision * (x(i+1) - x(i-1))))
  -      endif
  -      state_out%M_dot(i) = - state_out%v(i) * S_in(i) * x(i)
  -      !--------------Compute varaibles for Fz---------------
  -      kappa_ff(i) = 6.13e22_x_precision * state_0%rho_0 * state_out%rho(i) * &
  -                    (state_0%T_0 * T_in(i))**(-7._x_precision/2._x_precision) 
  -      tau(i)      = 0.5_x_precision * sqrt(0.2_x_precision * (1._x_precision * para%X) * 6.13e22_x_precision * state_0%rho_0 &
  -                    * state_out%rho(i) * (state_0%T_0 * T_in(i))**(-7._x_precision/2._x_precision)) * state_0%S_0 &
  -                    * S_in(i) / x(i)
  -      epsil(i)    = 6.22e20_x_precision * (state_0%rho_0 * state_out%rho(i))**2  &
  -                    * sqrt((state_0%T_0 * T_in(i))) 
  -
  -      if (tau(i) .ge. 1.0) then
  -        state_out%Fz(i) = (4._x_precision * c * c * T_in(i)**4) / (27._x_precision * sqrt(3.0) * &
  -                          (kappa_ff(i) + kappa_e) * (S_in(i)/x(i) * state_0%S_0)) 
  -      else
  -        state_out%Fz(i) = epsil(i) * state_out%H(i) * state_0%temps_0 / state_0%rho_0
  -      endif
  +       state_out%x(i)     = x(i)
  +       state_out%Omega(i) = Omega(i)
  +       state_out%H(i)     = - 0.5_x_precision * (b1(i) + sign(sqrt(Delta(i)),b1(i))) &
  +            / a1(i) 
  +       state_out%P_rad(i) = T_in(i)**4
  +       state_out%cs(i)    = Omega(i) * state_out%H(i)
  +       state_out%rho(i)   = S_in(i) / (state_out%H(i) * x(i))
  +       state_out%nu(i)    = para%alpha * state_out%cs(i) * state_out%H(i)
  +       state_out%P_gaz(i) = state_out%rho(i) * (T_in(i)**4)
  +       state_out%beta(i)  = state_out%P_gaz(i) / (state_out%P_gaz(i) + state_out%P_rad(i))
  +       state_out%Cv(i)    = (kmp/mu * state_0%T_0) * ((12._x_precision * gammag -1._x_precision ) &
  +            * (1._x_precision - state_out%beta(i)) + state_out%beta(i)) &
  +            / ( state_out%beta(i) * (gammag - 1._x_precision))
  +       !------------limit condition to compute v-------------
  +       if (i .eq. 1) then
  +          state_out%v(i)  = 0._x_precision
  +       else if (i .eq. n_cell) then
  +          state_out%v(i)  = - 1._x_precision / S_in(i) / x(i) 
  +       else
  +          state_out%v(i)  = - 1._x_precision / S_in(i) / x(i) * ( ((state_out%nu(i+1) &
  +               * S_in(i+1)) - (state_out%nu(i-1) * S_in(i-1))) /  &
  +               (2._x_precision * (x(i+1) - x(i-1))))
  +       endif
  +       state_out%M_dot(i) = - state_out%v(i) * S_in(i) * x(i)
  +       !--------------Compute varaibles for Fz---------------
  +       kappa_ff(i) = 6.13e22_x_precision * state_0%rho_0 * state_out%rho(i) * &
  +            (state_0%T_0 * T_in(i))**(-7._x_precision/2._x_precision)
  +       
  +       tau(i)      = 0.5_x_precision * sqrt(0.2_x_precision * (1._x_precision &
  +            * para%X) * 6.13e22_x_precision * state_0%rho_0 * state_out%rho(i) &
  +            * (state_0%T_0 * T_in(i))**(-7._x_precision/2._x_precision)) &
  +            * state_0%S_0 * S_in(i) / x(i)
  +       
  +       epsil(i)    = 6.22e20_x_precision * (state_0%rho_0 * state_out%rho(i))**2  &
  +            * sqrt((state_0%T_0 * T_in(i))) 
  +       
  +       if (tau(i) .ge. 1.0) then
  +          state_out%Fz(i) = (4._x_precision * c * c * T_in(i)**4) / (27._x_precision &
  +               * sqrt(3.0) * (kappa_ff(i) + kappa_e) * (S_in(i)/x(i) * state_0%S_0)) 
  +       else
  +          state_out%Fz(i) = epsil(i) * state_out%H(i) * state_0%temps_0 / state_0%rho_0
  +       endif
       enddo
     end subroutine compute_variables
     !--------------------------------------------------------------
  
positive.code:     state_0%H_0     = rs
      state_0%M_dot_0 = para%Mdot
      state_0%rho_0   = state_0%S_0 / (2._x_precision*rs)
      state_0%T_0     = (1._x_precision/sqrt(27.0) * 1._x_precision/48._x_precision * para%Mdot * c2 &
                        / ( pi * rs * rs * stefan ))**0.25_x_precision
      state_0%Fz_0    = state_0%S_0 / 2._x_precision / state_0%temps_0
      state_0%Cv_0    = 1._x_precision / state_0%T_0
      state_0%P_gaz_0 = state_0%rho_0 * kmp * state_0%T_0 / mu
  
negative.code:     c2 = c**2
      
      call get_parameters(para)
      rs = 2._x_precision*G*para%M/c2
      
      state_0%omega_0 = sqrt( G*para%M/(rs)**3 )
      state_0%temps_0 = 2._x_precision / state_0%omega_0
  