anchor.diff:      close(fid_c)
   
     end subroutine write_critical_points
  +
  +  !-------------------------------------------------------------------------
  +  ! Subroutine to save in a file the S curve
  +  !-------------------------------------------------------------------------
  +  subroutine save_data(k, sigma)
  +    implicit none
  +
  +    integer,                                    intent(in) :: k
  +    real(kind = x_precision), dimension(nb_it), intent(in) :: sigma
  +
  +    !------------------------------------------------------------------------
  +    ! For the output file
  +    character(len = 8)                         :: number_of_cell
  +    character(len = 64)                        :: fname_tot
  +    integer                                    :: fid_tot
  +    integer                                    :: ios
  +
  +    real(kind = x_precision), dimension(nb_it) :: sigma_real ! Dimensioned surface density
  +    real(kind = x_precision)                   :: f
  +    real(kind = x_precision)                   :: tau_eff
  +    integer                                    :: optical_depth
  +    integer                                    :: j
  +    !------------------------------------------------------------------------
  +
  +    write(number_of_cell,'(I5.5)') k
  +    fid_tot = 22 + k
  +    fname_tot = 's_curves/Temperature_Sigma_'//trim(number_of_cell)//'_tot.dat'
  +
  +    open(fid_tot,file  = fname_tot, action='write', status = 'replace', iostat = ios)
  +    if (ios /= 0) then
  +      write(*,*)"Error while opening the ", fname_tot," file."
  +      stop
  +    endif
  +
  +    write(fid_tot,'(3(A16))') 'Surface_density', 'Temperature', 'Optical_depth'
  +
  +    do j = 1, nb_it
  +      call variables(k, temperature(j), sigma(j), f, optical_depth, tau_eff)
  +
  +      sigma_real(j) = sigma(j) * state_0%S_0
  +
  +      write(fid_tot,fmt = '(3(e16.6e2))') sigma_real(j), temp_real(j), tau_eff
  +    enddo
  +
  +    close(fid_tot)
  +
  +  end subroutine save_data
   
   
     !-------------------------------------------------------------------------
  
positive.code:   implicit none
  
    real(kind = x_precision), dimension(nb_it) :: temperature ! Temperatures for which we need to solve the dichotomy
    real(kind = x_precision)                   :: eps         ! Precision required for the dichotomy
    real(kind = x_precision)                   :: S_min       ! Minimum surface density limit
    real(kind = x_precision)                   :: S_max       ! Maximum surface density limit
  
negative.code:       temperature(i) = dt * (i-1) + Tmin
      enddo
  
      temp_real  = temperature * state_0%T_0
  
    end subroutine make_temperature
  
  
  