function [z,b] =hillslope_profile3(z,b)
%
% hillslope_profile is a generic 2-D hillslope profile evolution solver
% hillslope_profile solves the eqns: qs = kx^ms^n ; dz/dt = -1/(1-porosity)*(dqs/dx)
% soil production as well as transport may be modeled
%
% note  that s = -dz/dx (flow from left to right).
% no flow and no sediment flux is assumed at the left boundary, i=1
%
% bcflag specifies variable boundary conditions at the right margin (toe of hill, i=Maxi)
% bcflag 1: z(Maxi) = constant (all sediment carried away)
% bcflag 2: z(Maxi) lowered at a constant rate (stream incision)
% bcflag 3: zero slope at i=Maxi; zero flux (all sediment accumulates)
% wflag 1: do not compute weathering (transport limited)
% wflag 2: compute weathering
%
% time step dt is determined from std advection-diffusion eqn stability (FTCS explicit)
% possible stability problems when n > 1,
% use stabil < 1 if problems arise.  stabil = 2.25 is max for n=1; m=1;
% stabil = 3 is max for n=1, m=2.
% for best stability with diffusion and non-linear diffusion (m=0) use CS in calc qs
% and dqs/dx
% for best stability with advection (m>0) use FS in calc qs but CS for
% dqs/dx
%
% units: z [m], dx [m], timeyr [yr], k [m^(2-n)/yr], E [m/yr]; all time and rates in yr
%
	hold off
%
% specify initial conditions if not specified in command line (e.g., [z2,b2]=hillslope_profile3(z1,b1))
% where z1 and b1 are vectors in your workspace.
% alternatively you can read in a saved file with "load"
%
   z=[15:-.5:1];
	b=[13:-.5:-1];
%	z=[12:-.8:0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%	b=[12:-.8:0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%    z=[14:-.05:12.05,12:-.8:0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%	b=[14:-.05:12.05,12:-.8:0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%    z=[12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12:-.8:0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%	b=[12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12:-.8:0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%    z=[12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12:-.8:2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2];
%	b=[10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10:-.8:0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%    z=[20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20:-.8:0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%	b=[20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20:-.8:0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
    
%
% define parameters
%
	n=1;
	m=0;
	timeyr=2000;
	dx=1;
	bcflag=1;
	wflag=1;
    rflag=1;        % rflag=2 is roering non-linear diffusion model (m=0,n=1) (only for wflag=1 case for now)
    Sc=.85;
	k=.01;
	porosity=0.0;
	E=.0001;
	Wd=2.3;
	Ws=.001;
	h = z - b;
	Maxi=length(z);
	ax=(0:(Maxi-1)).*dx;
	x=Maxi.*dx;
    nplot=5;
    scarp_h = max(z) - min(z)
%
% compute dt for stability (first guess)
%
	stabil = 2.25;
    if rflag == 2
        stabil2 = .02;          % reduce from 0.95 to 0.02 for roering model if approaching angle of repose
    else 
        stabil2 = .95;
    end
	if m==0
		B=dx.*dx.*(1-porosity)./(2.*k);
		dt=stabil2.*B             
	else
		A=2.*dx.*(1-porosity)./(k.*m.*x.^(m-1));
		B=dx.*dx.*(1-porosity)./(2.*k.*x.^m);
		dt=stabil.*min(A,B)
	end
%
% compute constants
%
	k1 = -1.*k./(dx.^n);
	k3 = dt./((1-porosity).*dx);
    Scdx = 2.*dx.*Sc;
%
% determine # timesteps for the simulation and specify number of times 
% to output data
%
	nsteps = timeyr./dt
	tplot=round(nsteps./nplot);
%
% initialize variables and compute x^n outside timeloops for maximized efficiency
%
	qs = zeros(size(z));
	xm = zeros(size(z));
	k1xm = zeros(size(z));
%
	for i = 1:Maxi
		xm(i) = ((i-1).*dx).^m;
	end
	xm(1)=xm(1)+(.2.*dx);
%
	k1xm1 = k1.*xm(1);
	k1xmm = k1.*xm(Maxi);
	k1xm = k1*xm;
%
% plot z initial for reference and then loop over time t = 1:nsteps
%
	hold off
	figure(1)
	clf
	plot(ax,z,'m',ax,b,'c')
	hold on
	set(gca,'xlabel',text(0,0,'distance (m)'))
	set(gca,'ylabel',text(0,0,'elevation (m)'))
%
% BEGIN BIG TIME LOOP
%
	for t = 1:nsteps
%
% first calc qs(i); BC -- no flow over ridge (i=1)
%
		if wflag==1
			i=1;
			if m==0
				qs(i) = k1.*(z(i+1) - z(i)).^n;
                Smax = abs((z(i+1)-z(i))/dx);
                if rflag==2
                    qs(i) = k1.*(z(i+1) - z(i)).^n*(1/(1-(abs(z(i+1) - z(i))/(Scdx./2))^2));
                    Smax = abs((z(i+1)-z(i))/dx);
                end
			else
				qs(i) = k1xm1.*(z(i+1) - z(i)).^n;
                Smax = abs((z(i+1)-z(i))/dx);
			end
            for i = 2:Maxi-1
                if m==0
                    if rflag==2
                        qs(i) = k1xm(i)./2.*(z(i+1) - z(i-1)).^n*(1/(1-(abs(z(i+1) - z(i-1))/Scdx)^2));
                        Smax = max(Smax,abs((z(i+1)-z(i-1))/(2*dx)));
                    else
				        qs(i) = k1xm(i)./2.*(z(i+1) - z(i-1)).^n;
                        Smax = max(Smax,abs((z(i+1)-z(i-1))/(2*dx)));
                    end
                else
                    qs(i) = k1xm(i).*(z(i+1) - z(i)).^n;
                    Smax = max(Smax,abs((z(i+1)-z(i-1))/(2*dx)));
                end
			end
			i = Maxi;
            if rflag==2
                qs(i) = k1xmm.*(z(i) - z(i-1)).^n*(1/(1-(abs(z(i) - z(i-1))/(Scdx./2))^2));
                Smax = max(Smax,abs((z(i)-z(i-1))/dx));
            else
      		    qs(i) = k1xmm.*(z(i) - z(i-1)).^n;
                Smax = max(Smax,abs((z(i)-z(i-1))/dx));
            end
%		
		elseif wflag==2
			i=1;
				W(i)=Ws.*exp(-1.*Wd.*h(i));
				if m==0
					qs(i) = k1.*(z(i+1) - z(i)).^n;
                    qs(i) = min(qs(i),(W(i) + h(i)./dt).*dx);
				else
					qs(i) = k1xm1.*(z(i+1) - z(i)).^n;
                    qs(i) = min(qs(i),(W(i) + h(i)./dt).*dx);
				end
%
			for i = 2:Maxi-1
					W(i)=Ws.*exp(-1.*Wd.*h(i));
					qs(i) = k1xm(i)./2.*(z(i+1) - z(i-1)).^n;
                    qs(i) = min(qs(i),(W(i) + h(i)./dt).*dx+qs(i-1));
			end
%
			i = Maxi;	
					W(i)=Ws.*exp(-1.*Wd.*h(i));	
					qs(i) = k1xmm.*(z(i) - z(i-1)).^n;
                    qs(i) = min(qs(i),(W(i) + h(i)./dt).*dx+qs(i-1));	
		end
% 		
			qs=abs(qs);
%
% now calc dz/dt based on qs(i), again with a no-flow BC at ridge (i=1)
%
		i=1;
			tempz = z(i);
			z(i) = z(i) - k3.*(qs(i));
			if wflag==2
				b(i) = -1.*dt.*W(i) + b(i);
				if z(i)<b(i)
					i
                    z(i)
                    b(i)
				end
				h(i)=z(i)-b(i);
			end
			dzdt(t)= (tempz-z(i))/dt;
            dzdti(i) = (tempz-z(i))/dt;
			n_dzdt(t) = dzdt(t)/E;
			time(t)=t*dt;
%			
		for i = 2:Maxi-2
            tempz=z(i);
			z(i) = z(i) - (k3./2).*(qs(i+1) - qs(i-1));
%            z(i) = z(i) - (k3).*(qs(i+1) - qs(i));
            dzdti(i) = (tempz-z(i))/dt;
			if wflag==2
				b(i) = -1.*dt.*W(i) + b(i);
				if z(i)<b(i)
					i
                    z(i)
                    b(i)
				end
				h(i)=z(i)-b(i);
			end
		end
%
		i=Maxi-1;
            tempz=z(i);
			z(i) = z(i) - k3.*(qs(i) - qs(i-1));
            dzdti(i) = (tempz-z(i))/dt;
			if wflag==2
				b(i) = -1.*dt.*W(i) + b(i);
				if z(i)<b(i)
					i
                    z(i)
                    b(i)
				end
				h(i)=z(i)-b(i);
			end
%
% calc last point applying appropriate BC as specified by bcflag
% note if bcflag = 1, z(Maxi) is constant, so do nothing
%
		i=Maxi;
		if bcflag==2
            tempz=z(i);
			z(i) = z(i) - E.*dt;
            dzdti(i) = (tempz-z(i))/dt;
			if wflag==2
				b(i) = -1.*dt.*W(i) + b(i);
				if z(i)<b(i)
					i
                    z(i)
                    b(i)
				end
				h(i)=z(i)-b(i);
			end
		elseif bcflag==1
            dzdti(i) = 0;
			if wflag==2
				b(i) = -1.*dt.*W(i) + b(i);
				if z(i)<b(i)
					i
                    z(i)
                    b(i)
				end
				h(i)=z(i)-b(i);
			end
		elseif bcflag==3
            tempz=z(i);
			z(i) = z(i) + (k3./2).*(qs(i-1));
            dzdti(i) = (tempz-z(i))/dt;
			if wflag==2
				b(i) = -1.*dt.*W(i) + b(i);
				if z(i)<b(i)
					i
                    z(i)
                    b(i)
				end
				h(i)=z(i)-b(i);
			end
      	end
%
% plot z at appropriate times.  use "hold on" to avoid erasing previous plots
%
		if rem(t,tplot)==0
%            t
			if wflag==1
				figure(1)
                plot (ax,z,'g')
			elseif wflag==2
				figure(1)
                plot (ax,z,'r',ax,b,'b')
				title (['Ws=',num2str(Ws),' Wd=',num2str(Wd),' timeyr=',num2str(timeyr),' m=',num2str(m),' n=',num2str(n),' k=',num2str(k),' bc=',num2str(bcflag),' E=',num2str(E)])
                figure(3)
                plot(h,dzdti,'m')
                hold on
                xlabel('soil thickness (m)')
                ylabel('dzdt')
                title('Soil thickness and erosion rate (curvature)')
            end
		end
%
% end time loop -- repeat for all nsteps
%
	end
%
% for Pierce and Colman test, calc effective diffusivity constant from Smax
% -- which holds the maximum value from the last time step
% using analytical result in Pierce and Colman paper (assuming initial
% scarp gradient = 0.8 = tan_alpha)
%
years = max(time)
max_slope = Smax
Keff = (scarp_h/(4*sqrt(max(time))*0.8*erf(Smax/0.8)))^2
%
% plot erosion rate at divide as a function of time in a second figure
%
	figure(2)
	plot(time,n_dzdt,'r')
	xlabel('time in years (t * dt)')
	ylabel('dzdt/E ')
	title([ 'time to steady state: E =' num2str(E)])
    figure(3)
    plot(h,dzdti,'r')
    xlabel('soil thickness (m)')
    ylabel('dzdt')
    title('Soil thickness and erosion rate (curvature)')
	hold off
	z=z;
 % write final value dzdt/E to screen
    norm_erosion_rate = n_dzdt(t)