function RTM1(varargin)
disp('Model reactive transport in one dimension')
nx = 100; % number of instances along x
length = 100.0; % m
deltax = length / nx; % m
deltay = 1.0; % m
deltaz = 1.0;

time_end = 10.0; % years
delta_years = time_end / 5.0;
next_output = 0.0;

% Configure and initialize the inlet and interior instances
cp_inlet = ChemPlugin();
% Cell Array of interior instances
cp = cell(nx,1);
for i = 1:nx
    cp{i} = ChemPlugin();
end
Console(cp{1}, 'RTM_stdout.txt');
Config(cp{1}, 'pluses = banner');

cmd = sprintf('volume = %f m3; time end = %f years;', ...
    deltax*deltay*deltaz, time_end);

for i = 1:nx
    Config(cp{i}, cmd);
end

if (nargin)
    file_name = varargin{1};
else
    file_name = input('Please enter RTM input script: ');
end

fprintf('Reading from file %s', file_name);
script = fopen(file_name, 'r');

scope = 0;
line = fgets(script);
while(ischar(line))
    line1 = strcat(line);
   if strcmp(line1, 'go')
        break;
    elseif strcmp(line1, 'scope inlet')
        scope = 1;
    elseif strcmp(line1, 'scope initial')
        scope = 2;
    elseif scope == 1
        error = Config(cp_inlet, line1);
    elseif scope == 2
        for i = 1:nx 
            Config(cp{i}, line1);
        end
    end
    line = fgets(script);
end
fclose(script);

if Initialize(cp_inlet) 
    error('Inlet failed to initialize');
end
PrintOutput(cp_inlet, 'RTM.txt', 'Inlet fluid');

for i = 1:nx
    if Initialize(cp{i})
        error('Instance %i failed to initialize', i);
    end
end
file_name = 'RTM.txt';
next_output = write_results(cp, nx, delta_years, next_output, file_name);

% Query user for velocity. 
if (nargin)
    if ischar(varargin{2})
        veloc_in = str2double(varargin{2});
    else
        veloc_in = varargin{2};
    end
else
    veloc_in = input('\n Please Enter Fluid velocity in m/yr: ');
end

velocity = veloc_in / 31557600;
porosity = Report1(cp{1}, 'porosity');
flow = deltay * deltaz * porosity * velocity;

diffcoef = 1e-10; % m2/s
dispersivity = 1.0; % m
dispcoef = velocity * dispersivity + diffcoef; % m2/s
trans = deltay * deltaz * porosity * dispcoef / deltax;

tcond = 2.0; % W/m/K
ttrans = deltay * deltaz * tcond / deltax;

% Link the instances
link1 = Link(cp{1}, cp_inlet);
Transmissivity(link1, trans, 'm3/s');
error = HeatTrans(link1, ttrans, 'W/K');
FlowRate(link1, flow, 'm3/s');

for i = 2:nx
    link1 = Link(cp{i}, cp{i-1});
    Transmissivity(link1, trans, 'm3/s');
    HeatTrans(link1, ttrans, 'W/K');
    FlowRate(link1, flow, 'm3/s');
end

% free outlet link
link1 = Link(cp{nx});
FlowRate(link1,-flow, 'm3/s'); 

NotDone = 1;
while NotDone
    deltat = 1e99;
    for i = 1:nx
        deltat = min(deltat, ReportTimeStep(cp{i}));
    end
    for i = 1:nx
       if (AdvanceTimeStep(cp{i}, deltat))
           NotDone = 0;
       end
    end
    for i = 1:nx
        if (AdvanceTransport(cp{i}))
            error('Advance Transport failed for instance %i',i);
        end
    end
     for i = 1:nx
        if (AdvanceHeatTransport(cp{i}))
            error('Advance Heat Transport failed for instance %i',i);
        end
    end
    for i = 1:nx
        if (AdvanceChemical(cp{i}))
            error('Advance Transport failed for instance %i',i);
        end
    end
    next_output = write_results(cp, nx, delta_years, next_output, file_name);

end
end
