function [ U ] = isostaticmassflux1D(x,t,varargin)
%   [ U ] = isostaticmassflux1D(x,t)
%   [ U ] = isostaticmassflux1D(x,t,water)
%
%   Numerically determine the mass flux in the asthenosphere for an across
%axis profile deriving from isostasy driven flow with differential thermal
%subsidence
%
%   x = location vector (assume km, but units are actually arbirtrary)
%   t = thermal age of lithosphere (governs subsidence rate)
%   water = optional binary vector defining which of locations are
%   submarine (1) and which are subaerial (0). Default is all submarine.
%
%   U = mass flux (km^3/Myr/km of ridge length)
%
%   Generalization of the mathematics in Conder (2012), Lithosphere
%   Developed by James Conder, March 2021
%
%   Included as supplemental materials to
%   Conder, J.A., Oceanic isostasy as a trigger to the rift to drift
%   transition, submitted to Geology, 2021

nx = length(x);
cstar = 0.11;   % km/sqrt(Myr)
t = max(0.02,t);

massgain = cstar./(2*sqrt(t));  %mass gain per area (km/Myr)
if nargin > 2
    submarine = varargin{1};
    if length(submarine) == nx
        massgain = massgain.*submarine;  %no ocean = no mass gain
    end
end

dx = (x(3:end) - x(1:end-2))/2;
dx = [(x(2)-x(1)) dx (x(end)-x(end-1))];
meanmassgain = sum(dx.*massgain)/(x(end)-x(1));

dUdx = massgain - meanmassgain;
U = 0*x;
U(1) = dUdx(1)*dx(1);

for i = 2:length(x)
    U(i) = U(i-1) + dUdx(i)*dx(i);
end

end