anchor.diff:  
       implicit none
   
  -    real(kind = x_precision),dimension(n_cell),intent(in) :: sigma_c_thin
  -    real(kind = x_precision),dimension(n_cell),intent(in) :: temp_c_thin
  -    real(kind = x_precision),dimension(n_cell),intent(in) :: tau_thin
  -    real(kind = x_precision),dimension(n_cell),intent(in) :: sigma_c_thick
  -    real(kind = x_precision),dimension(n_cell),intent(in) :: temp_c_thick
  -    real(kind = x_precision),dimension(n_cell),intent(in) :: tau_thick
  +    real(kind = x_precision), dimension(n_cell), intent(in) :: sigma_c_thin
  +    real(kind = x_precision), dimension(n_cell), intent(in) :: temp_c_thin
  +    real(kind = x_precision), dimension(n_cell), intent(in) :: tau_thin
  +    real(kind = x_precision), dimension(n_cell), intent(in) :: sigma_c_thick
  +    real(kind = x_precision), dimension(n_cell), intent(in) :: temp_c_thick
  +    real(kind = x_precision), dimension(n_cell), intent(in) :: tau_thick
       !-----------------------------------------------------------------------
       integer             :: i
  -    integer             :: fid_4 = 11
  -    character(len = 64) :: fname_4
  -
  -    !-----------------------------------------------------------------------
  -
  -    fname_4 = 'critical_points/file.dat'
  -
  -    open(fid_4,file  = fname_4, status='unknown',action='readwrite')
  -    write(fid_4,'(7(A16))') 'Radius','Temp_thin','Sigma_thin','Tau_thin',&
  +    integer             :: fid_c = 11
  +    character(len = 64) :: fname_c
  +
  +    !-----------------------------------------------------------------------
  +
  +    fname_c = 'critical_points/file.dat'
  +
  +    open(fid_c,file  = fname_c, status='unknown',action='readwrite')
  +    write(fid_c,'(7(A16))') 'Radius','Temp_thin','Sigma_thin','Tau_thin',&
                                        'Temp_thick','Sigma_thick','Tau_thick'
   
  -    do i=1, n_cell
  -      write(fid_4,'(7(e16.6e2))') r_state%r(i), temp_c_thin(i), sigma_c_thin(i), tau_thin(i),&
  +    do i = 1, n_cell
  +      write(fid_c,'(7(e16.6e2))') r_state%r(i), temp_c_thin(i), sigma_c_thin(i), tau_thin(i),&
                                                   temp_c_thick(i), sigma_c_thick(i), tau_thick(i)
       end do
   
  -    close(fid_4)
  +    close(fid_c)
   
     end subroutine write_critical_points
   
  
positive.code:   !-------------------------------------------------------------------------
    ! Subroutine in order to display the critical points
    !-------------------------------------------------------------------------
    subroutine display_critical_points(sigma_c_thin, temp_c_thin, sigma_c_thick, temp_c_thick,k)
      implicit none
  
      real(kind = x_precision),intent(in) :: sigma_c_thin
      real(kind = x_precision),intent(in) :: temp_c_thin
      real(kind = x_precision),intent(in) :: sigma_c_thick
      real(kind = x_precision),intent(in) :: temp_c_thick
      integer                 ,intent(in) :: k
  
      !------------------------------------------------------------------------
  
      write(*,*)'**** Critical Point',k,'********'
      write(*,*)'****************************************'
      write(*,"(' Optically thin (T,sigma) :',1p,E12.4,4x,1p,E12.4)")temp_c_thin, sigma_c_thin
      write(*,"(' Optically thick (T,sigma):',1p,E12.4,4x,1p,E12.4)")temp_c_thick, sigma_c_thick
      write(*,*)'****************************************'
  
    end subroutine display_critical_points
  
negative.code:   !-------------------------------------------------------------------------
    ! Subroutine in order to build the S curve
    !-------------------------------------------------------------------------
    subroutine build_s_curve(sigma_real_thick, sigma_real_thin, index_scp, sigma_real)
      implicit none
  
      real(kind = x_precision), dimension(nb_it), intent(in)  :: sigma_real_thick
      real(kind = x_precision), dimension(nb_it), intent(in)  :: sigma_real_thin
      integer,                                    intent(in)  :: index_scp
  
      real(kind = x_precision), dimension(nb_it), intent(out) :: sigma_real
      !-----------------------------------------------------------------------
  
      integer :: i
  
      do i = 1, index_scp
        sigma_real(i) = sigma_real_thick(i)
      enddo
      do i = index_scp + 1, nb_it
        sigma_real(i) = sigma_real_thin(i)
      enddo
  
    end subroutine build_s_curve
  