function [Profile,Bedrock] = p1dpl(x,nx,dx,dt,t,totaltime,LBC_TYPE,LBC, ...
   RBC_TYPE,RBC,vx,vz,kappa,Ba,Bb,ProfileInitCond,BedrockInitCond,waithandle)

simtime = dt;

x1 = x(1);
xend = x(end);

Profile = ProfileInitCond;
Bedrock = BedrockInitCond;

X = x;

indicies = 1:length(x);

while (simtime <= t)

  % First, produce regolith for time step:

  Bedrock = Bedrock - (Ba*exp(-Bb.* ((Profile-Bedrock) > 0) .* ...
     (Profile-Bedrock)))*dt;

  % Transport (Linear Diffusion):

  Qt = -kappa.*diff(Profile)./dx;
  
  % Apply Boundary Conditions:
 
  if(LBC_TYPE)
    LBCn = LBC;
  else
    LBCn = Qt(1);
    Profile(1) = LBC;
  end

  if(RBC_TYPE)
    RBCn = RBC;
  else
    RBCn = Qt(nx-1);
    Profile(nx) = RBC;
  end

  Qt = [LBCn,Qt,RBCn];

  % Mass conservation:

  dH = dt.*(diff(Qt)./dx);

  % Do production-limited conditions exist?
  PLLocation = ((Profile - dH) < Bedrock);
  ProdLimited = sum(PLLocation);

  if(ProdLimited > 0)

    % For computational efficiency, we must adjust flux and dH values from 
    %   the highest point in the profile to the lowest.  If we do not 
    %   adjust flux and dH values in this order, we must converge the fluxes
    %   and dH values, since the fluxes and dHs are affected by all upslope
    %   fluxes in a production-limited model.  First, we sort the points
    %   from highest to lowest:

    [ProfSort,j] = sort(Profile);
    j = fliplr(j);

    % Remove the boundary conditions from the profile:

    Qt = Qt(2:end-1);

    % Now, we work down the profile from the highest to the lowest points,
    %   adjusting fluxes and dH values accordingly. 

    for(i=1:(length(Profile)-2))
      if((Profile(j(i))-dH(j(i))) < Bedrock(j(i)))
        pos = j(i);
        rpos = pos + 1;
        lpos = pos - 1;

        % If we are at a production-limited topographic divide, divide up
        %   the remaining regolith between the two directions according to
        %   their slopes.  Then, adjust dH at the surrounding locations
        %   accordingly:

        if((rpos < length(Profile)) & (lpos > 1) & (Profile(lpos) ...
           < Profile(pos)) & (Profile(rpos) < Profile(pos)))
          scale = abs(Profile(pos)-Profile(lpos))./(abs(Profile(pos) - ...
            Profile(lpos)) + abs(Profile(pos) - Profile(rpos)));
          Qt(lpos) = -(dx.*(Profile(pos)-Bedrock(pos))./dt).*scale;
          Qt(pos) = -(dx.*(Profile(pos)-Bedrock(pos))./dt).*(1-scale);
          dH(pos) = Profile(pos) - Bedrock(pos);
          dH(lpos) = dt.*((Qt(lpos) - Qt(lpos-1))./dx);
          dH(rpos) = dt.*((Qt(rpos) - Qt(pos))./dx);

        % If we are at a production-limited point that is not a topographic
        %   divide and the downslope direction is to the right, push regolith
        %   downhill to the right.  Adjust the current location's dH and the
        %   node to the right:

        elseif((rpos < length(Profile)) & (lpos > 1))
          if(Profile(pos) > Profile(rpos))
            Qt(pos) = (dx.*(Profile(pos)-Bedrock(pos))./dt) + Qt(lpos);
            dH(pos) = Profile(pos) - Bedrock(pos);
            dH(rpos) = dt.*((Qt(rpos) - Qt(pos))./dx);
          end

        % Finally, if we are at a production limited point and downslope is
        %    to the left, and we are not at a divide, adjust flux to the left
        %    and adjust dH at current location and the node to the left:

        elseif((lpos > 1) & (rpos < length(Profile)))
          if(Profile(pos) > Profile(lpos))
            Qt(lpos) = -(dx.*(Profile(pos)-Bedrock(pos))./dt) + Qt(pos);
            dH(pos) = Profile(pos) - Bedrock(pos);
            dH(lpos) = dt.*((Qt(lpos) - Qt(lpos-1))./dx);
          end
        end
      end
    end

    % Re-Apply Boundary Conditions:
 
    if(LBC_TYPE)
      LBCn = LBC;
    else
      LBCn = Qt(1);
      Profile(1) = LBC;
    end

    if(RBC_TYPE)
      RBCn = RBC;
    else
      RBCn = Qt(end);
      Profile(nx) = RBC;
    end

    Qt = [LBCn,Qt,RBCn];

    % Mass Conservation (after flux adjustment):

    dH = dt.*(diff(Qt)./dx);

  end

  %  Add Elevation Change to Profile:
  Profile = Profile - dH;


  % Add tectonic displacements to profile, x, and bedrock:

  x = x + vx.*dt;
  Profile = Profile + vz.*dt;
  Bedrock = Bedrock + vz.*dt;

  %  Interpolate new profile:
  X(1) = x1;
  x(end) = xend;

  if(sum(x == X) ~= 0)
    Profile = interp1(x,Profile,X,'linear');
    x = X;
  end

  %  Increment time:

  if(((t-simtime) < dt) & ((t-simtime) > 0))
    dt = t-simtime;
  end

  simtime = simtime + dt;

  %  Ensure that Profile is above bedrock (prevents roundoff errors):
  Profile = (Profile >= Bedrock).*Profile + (Profile < Bedrock).*Bedrock;

  %  Update waitbox handle:
  waitbar(simtime./totaltime,waithandle);

end











