function [g, chease_data, nout_file] = find_consistent_CHEASE_equilibrium(RAPTOR_out, RAPTOR_Ip, RAPTOR_config, RAPTOR_model, ...
                                       B0, R0, RZ_contours, nout_file, free_fns, index)
% Update the MHD equilibrium based on imposed radial profiles

persistent CHEASE_updt_file_dir timetag

if index == 1
  CHEASE_updt_file_dir = fullfile(RAPTOR_model.envopts.RAPTOR_path, 'data', 'CHEASE_updt_EQ');
  if ~exist(CHEASE_updt_file_dir, 'dir');
    mkdir(CHEASE_updt_file_dir);
  end
  
  timetag = datestr(now, 30);
  
  if ~exist(fullfile(CHEASE_updt_file_dir, timetag), 'dir');
    mkdir(fullfile(CHEASE_updt_file_dir, timetag));
  end
end

mu0 = 4e-7*pi;

% Create EXPEQ file
expeq_filename = strcat('EXPEQ_iteration', num2str(index));
expeq_file = fullfile(CHEASE_updt_file_dir, timetag, expeq_filename);

expeq_struct = write_expeq(); % default
switch free_fns
  case 'ttprime_pprime_rhopol'
    nsttp = 1;
    TTprime = RAPTOR_out.ttprime * (1/B0);
    expeq_struct.TTprime = TTprime;
    
    sgn_propt = 1;
    nppfun = 4;
    pprime = RAPTOR_out.pprime * (mu0*R0^2/B0);
    expeq_struct.Pprime = pprime;
    
    rhotor = false;
    scal = 2; % 1 could also be implemented
  case 'jpar_p_rhotor'
    nsttp = 4;
    jpar = RAPTOR_out.jpar * (mu0*R0/B0);
    expeq_struct.jdotb_over_b0 = jpar;
    
    sgn_propt = -1;
    nppfun = 8;
    p = RAPTOR_out.p * (mu0/B0^2);
    expeq_struct.pressure = p;
    
    rhotor = true;
    scal = 0;
  case 'jpar_p_rhopol'    
    nsttp = 4;
    jpar = RAPTOR_out.jpar * (mu0*R0/B0);
    expeq_struct.jdotb_over_b0 = jpar;
    
    sgn_propt = -1;
    nppfun = 8;
    p = RAPTOR_out.p * (mu0/B0^2);
    expeq_struct.pressure = p;
    
    rhotor = false;
    scal = 0;
  case 'q_p_rhopol' 
    nsttp = 5;
    expeq_struct.q = RAPTOR_out.q;
    
    sgn_propt = -1;
    nppfun = 8;
    p = RAPTOR_out.p * (mu0/B0^2);
    expeq_struct.pressure = p;
    
    rhotor = false;
    scal = 0;
end
expeq_struct.pedge = RAPTOR_out.p(end) * (mu0/B0^2);

if rhotor
  expeq_struct.rho = RAPTOR_out.rho;
  expeq_struct.nrhotype = 1; % 0 -rho_psi, 1 -rho_tor. Make sure 'datain' has NFUNRHO=1
else
  expeq_struct.rho = RAPTOR_out.rhopol;
  expeq_struct.nrhotype = 0;
end
% Grid of rho, can be rho_psi or rho_tor (set nrhotype accordingly)
expeq_struct.RZ_psi = RZ_contours; % nx2 matrix with R/R0 and Z/R0 coordinates of boundary
expeq_struct.nsttp = nsttp; % 1 -Give TT', 2 -Give I*, 3 -Give I_||, 4 -Give J_||, 5 -Give Q
expeq_struct.nppfun = nppfun; % 4 => pprime, 8 => pressure given (if only pprime given but nppfun=8, then p in pprime)

expeq_struct_wrt = write_expeq(expeq_struct, expeq_file);
if index == 1 
  plot_expeq(expeq_struct_wrt.fnamefull, 'b'); % plot input EXPEQ file
end

% Create namelist file
namelist_filename = strcat('namelist_iteration', num2str(index));
[~,~,namelist_struct] = run_chease(1);
namelist_file = fullfile(CHEASE_updt_file_dir, timetag, namelist_filename);

namelist_struct.b0exp = B0;
namelist_struct.r0exp = R0;
namelist_struct.pedge = expeq_struct.pedge;

if scal == 2
  namelist_struct.ncscal = 2; 
  Ip = RAPTOR_Ip * mu0/(R0*B0);
  namelist_struct.currt = Ip;
elseif scal == 1
  error('ncscal=1 not yet implemented')
elseif scal == 0
  namelist_struct.ncscal = 4; 
end

namelist_struct.nfunrho = expeq_struct.nrhotype; 
namelist_struct.nrhomesh = expeq_struct.nrhotype;
namelist_struct.nsttp = nsttp;
namelist_struct.nppfun = nppfun;
namelist_struct.npropt = sgn_propt * nsttp;
namelist_struct.tensbnd = -1; %-5;
namelist_struct.tensprof = -1; %-5;
namelist_struct.tensbnd = 1e-7;
namelist_struct.tensprof = 1e-7;
namelist_struct.nrbox = 257;
namelist_struct.nzbox = 257;
namelist_struct.relax = 0.5;

if isempty(nout_file)
  namelist_struct.nopt = 0; % don't use NOUT file
else
  namelist_struct.nopt = -2; % use NOUT file as initial condition input
end
write_namelist_chease(namelist_file, namelist_struct);

% Run CHEASE
fname_out = run_chease(namelist_file, expeq_file, [], [], [], [], [], nout_file);

% Find eqdsk file
for ii=1:length(fname_out);
  if ~isempty(strfind(fname_out{ii},'EQDSK'))
    eqdsk_file_out = fname_out{ii}; break
  end
end

% Find cols file
for ii=1:length(fname_out);
  if ~isempty(strfind(fname_out{ii},'.cols'))
    ocols_file_out = fname_out{ii}; break
  end
end

% Find nout file
last_slash = strfind(fname_out{1},'/');
nout_file = [fname_out{1}(1:last_slash(end)),'NOUT'];

profiledata_out = read_ocols(ocols_file_out);
cocos = 1;
eqdskdata_out = read_eqdsk(eqdsk_file_out, cocos);

% Calculate equilibrium profiles
chease_data.eqdsk = eqdskdata_out;
chease_data.profiledata = profiledata_out;

% Calculate g vector
g = update_g_from_profiledata(profiledata_out, B0, R0, RAPTOR_config, RAPTOR_model);
end