% seacliff generation  code
% written by r.s. anderson in 1998, with
% modifications since then
% Note: This vs code simplified by Merritts
%for GSA workshop on Tectonics & Topograhy, 2001

% requires the file: delo18.file4
% based on oxygen isotope data from deep-sea cores
% for generation of a sealevel curve

clear
clc
clf

%DEVELOP SEA LEVEL CURVE FROM 0-18 data--
% use sin function to generate shape of sea level curve;
% seamp is the half amplitude of sealevel change, and is 60 m;
% the half period for sealevel oscillations is 10,000 yrs
seamean=-60;
seaamp=60;
pi=3.14159;
% sea=seamean+(seaamp*sin(2*pi*time/period));
%time/period = fraction of time along sin curve;
%sea level varies from -120 m to 0 m about mean of -60 m

% establish time frame, from -740,000 yrs to present (0 yrs);
% use time interval of 100 yrs
tmax=740000;
dt=100;
time=0:dt:tmax;
ttmax=max(time);
ttmin=min(time);
jmax=(tmax/dt)+1;

%Open/Load the file that contains the sea level curve data
load delo18.file3;
jj=floor(tmax/2000)+1
a=delo18(1:jj,:);
aflip=flipud(a);
ti=(tmax/1000):-(dt/1000):0;
del=interp1(aflip(:,1),aflip(:,2),ti);
delscale=2*seaamp/(4.04)
% sea=seamean-(del*delscale)-2.0792;
% this sets present sea level to zero;
%NOTE:   Matlab command FLOOR(X) rounds the elements of X to the  
%nearest integers towards minus infinity.

sea = -32.841 - (33.747*del)-(8.5605*del.*del);
% This is the nonlinear (polynomial) predictor of sealevel 
% derived from scatter plots of sea level elevation vs del 0-18 
% using Huon Peninsula and Barbados work;
%see Fig. 3, Anderson et al, 1999
%Refs: Chappell & Shackleton, 1986; Fairbanks, 1989; and Bard, 1990. 

%SET UP INPUT PARAMETERS
% uplift rate = uprate (m/yr);
% sea cliff erosion rate = edot0 (m/yr);
% slope of the local shelf = slope2;
% topographic diffusivity in m^2/year = kappa0.
%All of these input parameters can be varied.
uprate = 0.0005; 
edot0 = 0.25;
slope2 = 0.017; 
kappa0 = 0.015;

%ESTABLISH MIN & MAX VALUES FOR SEALEVEL-TIME PLOT
slope = 0.035; 
zmin=(seamean-seaamp)-(tmax*uprate)-40
zmax=(seamean+seaamp)+200
xmax=(zmax-zmin)/slope;
dx=5;
x=1:dx:xmax;
imax=length(x);
z0=zmin+(slope*x);
z=z0;

% CONTROLS ON SEA-CLIFF EROSION:  WAVE-ENERGY DISSIPATION
%Far-field wave energy is reduced by energy dissipation rate 
%along travel path.
zstar=4;
wavedepth=3*zstar;
dissscale=3*dx*wavedepth;

%SET # IMAGE FRAMES TO SEE DURING TERRACE %GENERATION/DEGRADATION
nprint=40;
jprint=floor(jmax/nprint);
numprint = 0;

%CREATE PLOT
%h=plot(x,z);  
xplot = x/1000;
xmaxplot = max(xplot);

% FIGURE(H) makes H the current figure, forces it to become visible,
%and raises it above all other figures on the screen.  If Figure H
%does not exist, and H is an integer, a new figure is created with
%handle H.

%Plot Time (tplot) versus Sea Level (sea)

figure(1) 
subplot(2,1,1)
subplot('position',[.15 .70 .75 .25])
tplot = (time-tmax)/1000;
plot(tplot,sea,'b')
xlabel(['Time (kyr)'],'fontsize',18,'fontname','times')
ylabel(['Sealevel (m)'],'fontsize',18,'fontname','times')
axis([-ttmax/1000 0 -130 20])
set(gca,'fontsize',14,'fontname','times')

% SUBPLOT Create axes in tiled positions.
%SUBPLOT(mnp), breaks the Figure window
% into an m-by-n matrix of small axes, and selects the p-th axes 
%for the current plot.  The axes are counted along the top row 
%of the Figure window, then the second row, etc.  For example,
 
      %  SUBPLOT(2,1,1)
      %  SUBPLOT(2,1,2)

% SUBPLOT('position',[left bottom width height]) creates an
% axis at the specified position in normalized coordinates 
% (in the range from 0.0 to 1.0).
%NOTE:   GCA is Matlab command to get handle to current axis.
% SET(H,'PropertyName',PropertyValue) sets the value of the
    %specified property for the graphics object with handle H.


subplot(2,1,2)
subplot('position',[.15 .10 .75 .50])
plot(xplot,z)
axis([0 xmaxplot min(z0) max(z0+(uprate*tmax))]);
xlabel(['Distance (m)'],'fontsize',18,'fontname','times')
ylabel(['Elevation (m)'],'fontsize',18,'fontname','times')

zsave=[x;z];

M = moviein(nprint,gcf);

%NOTE:   GCF is Matlab command to get handle to current figure.
    %H = GCF returns the handle of the current figure. The current
    %figure is the window into which graphics commands like PLOT,
    %TITLE, SURF, etc. will draw.


% ********** start the time clock and run program *************

for j=1:jmax

% Type 'Help for' (w/o quotes) to see defn of statement.
% Recall that jmax=(tmax/dt)+1;
% tmax = 740,000, and dt = 100, so jmax = 7401.
% The uplift pattern is taken to be uniform in x.

uplift=uprate*dt*ones(1,length(x));
%Ones creates M-by-N matrix of ones.
z=z+uplift;
% Find the intersection of the sea surface and the topography
surface=sea(j)*ones(1,length(x));
k=find(z>surface);
xsea1(j)=x(k(1));

% Here we insert a dependence of the cliff
% erosion rate on the width of the shelf.
% Wave energy is dissipated as waves approach shore,
% as a function of the weighted water depth.

gg=find(((sea(j)-wavedepth)<z)&(sea(j)>z));
dissrate=exp(-(sea(j)-z(1,gg))/zstar)*dx;
dissipate=sum(dissrate);
edot=edot0*(dissscale/dissipate);
ero(j)=edot;

% Recall that edot0 is the cliff erosion rate in m/yr, 
% set initially at 0.25 m/yr.

ksea=k(1)+round(edot*dt/dx);
xsea(j)=x(ksea);
shelf=z(k(1))+(slope2*(x-x(k(1))));

% Recall that slope2 = slope of local shelf, 
% set initially at 0.017.

surface(1:(k(1)-1))=9999*ones(size(1:k(1)-1));
surface(ksea+1:imax)=9999*ones(size(ksea+1:imax));
z=min(z,surface);

% Now take into account erosion in this reach of
% nearshore area where water depths are < "wavedepth".

ggnew=find(((sea(j)-wavedepth)<z)&(sea(j)>z));
nearshore=wavedepth/(x(max(ggnew))-x(ggnew(1)));
znew=z(1,ggnew(1))+(nearshore*(x(1,ggnew)-x(1,ggnew(1))));
z(1,ggnew)=min(z(1,ggnew),znew);

% Now diffuse wave action over the terrestrial 
% portion of the shelf space.

kappa=zeros(1,length(x));
kkk=find(x>xsea(j));
kappa1=kappa0*ones(1,length(kkk));
kappa(kkk(1):imax)=kappa1;

dzdt = 4 * kappa .* ((del2(z)./(dx*dx)));
z=z+(dzdt*dt);

% Now determine whether or not we want to take a snapshot.

p=rem(j,jprint);
if p == 0
numprint = numprint+1;

figure (1)

%subplot(2,1,1) is Time versus Sea Level

subplot('position',[.15 .70 .75 .25])
tplot = (time-tmax)/1000;
plot(tplot,sea,'b')
xlabel(['Time (kyr)'],'fontsize',18,'fontname','times')
ylabel(['Sealevel (m)'],'fontsize',18,'fontname','times')
axis([-ttmax/1000 0 -130 20])
hold on
tnow = (time(j)-tmax)/1000;
plot(tnow,sea(j),'rd','linewidth',3.0)
set(gca,'fontsize',14,'fontname','times')
hold off

	xp1 = xplot;
	xpatch = [0 xp1 xp1(imax:-1:1)];
	zp1 = z;
	zp2 = min(z0)*ones(1,length(x));
	zpatch = [min(z0) zp1 zp2(imax:-1:1)];

	subplot('position',[.15 .10 .75 .50])

%subplot(2,1,2) is Horizontal Distance versus Landmass Elevation

plot(xpatch,zpatch,'k');
patch(xpatch,zpatch,'k');

axis([0 xmaxplot min(z0) max(z0+(uprate*tmax))]);
xlabel(['Horizontal Distance (km)'],'fontsize',18,'fontname','times')
ylabel(['Elevation (m)'],'fontsize',18,'fontname','times')
set(gca,'fontsize',14,'fontname','times')

hold on
seanow = sea(j)*ones(1,ksea);
xnowplot = xplot(1:ksea);
plot(xnowplot,seanow,'b','linewidth',3.0)
hold off

zsave=[zsave;z];

edot
abs(time(j)-tmax)

M(:,numprint) = getframe(gcf); 

% if you dont want movies comment this out

end


end

zsend=zsave';
save zcliff.file zsend -ascii -tabs

time=time/1000;
erosave=[time;ero]';
save erocliff.file erosave -ascii -tabs

xsea=xsea/1000;
xsave=[xsea;time]';

%save xsea.file xsave -ascii -tabs
% if you want files to manipulate, un-comment this last line

% You can run  a movie animation by typing in the command window:
% clf
% movie(M,3,2) and it will run the movie 3 times at 2 frames/second