%% Input for ISSM runs
% The MATLAB code in this file provides the input required for running the ISSM 'synthetic glacier' experiment described in
% the paper "Atypical landslide induces speedup, advance, and long-term slowdown of a tidewater glacier".
% Some code based on the ISSM ISMIP exaple file. All credit to the original authors.
% Please contact me at vanwy048@umn.edu with any queries.

%Create synthetic topography

%Create a parabolic bed topography, with 3.5% slope, 100m resolution, 3km
%iwde, 10 km long
bed = ones(30,200).*(3.5:3.5:700);

for loop = 1:30 %Parabolic valley
    bed(loop,:) =  bed(loop,:)+abs(loop-15.5).^2*400./210.25;
end


%Calculate an ice thickness
thick = zeros(30,200);

for loop = 1:30 %Parabolic valley
    thick(loop,:) =  10+400-abs(loop-15.5).^2*400./210.25;
end

%
%% Create mesh

%Create a square mesh
md=model();
md=squaremesh(md,20000,3000,200,30);

%
%% Set all ice to grounded

% set the mask 
md=setmask(md,'','');

% plot the given mask #md.mask to locate the field (if you like)
plotmodel(md,'data',md.mask.groundedice_levelset);

%
%% Parametrize model and build geometry

%Geometry
disp('   Constructing Geometry');


%Define the geometry of the simulation #md.geometry
surface_g = bed + thick; %surface is bed + thickness
md.geometry.surface = surface_g(:);
md.geometry.base=bed(:);
%thickness is the difference between surface and base #md.geometry
%->
md.geometry.thickness=thick(:);
%plot the geometry to check it outbed + thick

plotmodel(md,'data',md.geometry.thickness);

disp('   Defining friction parameters');
%These parameters will not be used but need to be fixed #md.friction
%one friciton coefficient per node (md.mesh.numberofvertices,1)
md.friction.coefficient=25.0*ones(md.mesh.numberofvertices,1);

%one friciton exponent (p,q) per element
md.friction.p=ones(md.mesh.numberofelements,1);
md.friction.q=ones(md.mesh.numberofelements,1);

disp('   Construct ice rheological properties');
%The rheology parameters sit in the material section #md.materials
%B has one value per vertex
md.materials.rheology_B=cuffey(273)*ones(md.mesh.numberofvertices,1);
%n has one value per element
md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);

disp('   Set boundary conditions');
%Set the default boundary conditions for an ice-sheet
md=SetIceSheetBC(md);

%Name the run
md.miscellaneous.name = 'Reclus_theoretical';

%Set flow equation to SSA, to calculate initial velocity field
md=setflowequation(md,'SSA','all');

%
%% Set stress balance parameters


%Needs one value per vertice
md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);

%Apply boundary conditions
maxY=find(md.mesh.y==max(md.mesh.y));
minY=find(md.mesh.y==min(md.mesh.y));

md.stressbalance.spcvx(maxY)=0.0;

md.stressbalance.spcvy(minY)=0.0;

md.stressbalance.spcvy(maxY)=0.0;

md.stressbalance.spcvx(minY)=0.0;

%
%% Launch run stressbalance

%Choose model parameters
md.cluster=generic('name',oshostname(),'np',1);

md.verbose=verbose('convergence',true);
% Solve #help solve -> this launches the ISSM run
md=solve(md,'Stressbalance');

%
%% Extrude and set full stokes flow law

extrudecoeff = 2;
md.initialization.vx = repmat(md.results.StressbalanceSolution.Vx,1,1);
md.initialization.vy = repmat(md.results.StressbalanceSolution.Vy,1,1);
md.initialization.vz= md.initialization.vy*0;
md=extrude(md,extrudecoeff,1);

% Now, set to run the 3D HO model 
md=setflowequation(md,'HO','all');

%
%% Set stress balance parameters

md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);


maxY=find(md.mesh.y==max(md.mesh.y));
% create minX
%->
minY=find(md.mesh.y==min(md.mesh.y));

md.stressbalance.spcvx(maxY)=0.0;

md.stressbalance.spcvy(minY)=0.0;

md.stressbalance.spcvy(maxY)=0.0;

md.stressbalance.spcvx(minY)=0.0;

%
%% Launch run stressbalance

md.cluster=generic('name',oshostname(),'np',1);
% Set which control message you want to see #help verbose
%->
md.verbose=verbose('convergence',true);


%Solve a 'diagnostic' HO stress balance model first
md=solve(md,'Stressbalance');

%
%% Launch run transient

%Now we start doing a forward run with timesteps.

md.initialization.vx = md.results.StressbalanceSolution.Vx;
md.initialization.vy = md.results.StressbalanceSolution.Vy;
md.initialization.vz= md.results.StressbalanceSolution.Vz;



md.cluster=generic('name',oshostname(),'np',1);
% Set which control message you want to see #help verbose
%->
md.verbose=verbose('convergence',true);
% set the transient model to ignore the thermal model
% #md.transient
%->
md.transient.isthermal=0;
% define the timestepping scheme
% everything here should be provided in years #md.timestepping
% give the length of the time_step (4 years)
%->
md.timestepping.time_step=0.1;
% give final_time (20*4 years time_steps)
%->
md.timestepping.final_time = 100*md.timestepping.time_step;
% Solve #help solve
% we are solving a TransientSolution
%->
md=solve(md,'Transient');

save HO10yr_steady md;

%
%% The run should now have reached steady state

%Run for 5 years with:
%1) Regular, no changes
%2) Regular, increased slipperiness.
%3 Regular, decreased slipperiness
%4) Edited topography
%5) Edited topography, increased slipperiness
%6 Edited topography, decreased sliperiness

%
%% 1 Regular

md = loadmodel('./HO10yr_steady');

md.initialization.vx = md.results.TransientSolution(end).Vx;
md.initialization.vy = md.results.TransientSolution(end).Vy;
md.initialization.vz= md.results.TransientSolution(end).Vz;

md.cluster=generic('name',oshostname(),'np',1);
% Set which control message you want to see #help verbose
%->
md.verbose=verbose('convergence',true);
% set the transient model to ignore the thermal model
% #md.transient
%->
md.transient.isthermal=0;
% define the timestepping scheme
% everything here should be provided in years #md.timestepping
% give the length of the time_step (4 years)
%->
md.timestepping.time_step=0.01;
% give final_time (20*4 years time_steps)
%->
md.timestepping.final_time = 500*md.timestepping.time_step;
% Solve #help solve
% we are solving a TransientSolution
%->
md=solve(md,'Transient');

save HO5yr_case1 md;

%
%% 2 Regular, increased slippyness

md = loadmodel('./HO10yr_steady');

md.initialization.vx = md.results.TransientSolution(end).Vx;
md.initialization.vy = md.results.TransientSolution(end).Vy;
md.initialization.vz= md.results.TransientSolution(end).Vz;

md.friction.coefficient=20.0*ones(md.mesh.numberofvertices,1);

md.cluster=generic('name',oshostname(),'np',1);
% Set which control message you want to see #help verbose
%->
md.verbose=verbose('convergence',true);
% set the transient model to ignore the thermal model
% #md.transient
%->
md.transient.isthermal=0;
% define the timestepping scheme
% everything here should be provided in years #md.timestepping
% give the length of the time_step (4 years)
%->
md.timestepping.time_step=0.01;
% give final_time (20*4 years time_steps)
%->
md.timestepping.final_time = 500*md.timestepping.time_step;
% Solve #help solve
% we are solving a TransientSolution
%->
md=solve(md,'Transient');

save HO5yr_case2 md;

%
%% 3 Regular, decreased slippyness

md = loadmodel('./HO10yr_steady');

md.initialization.vx = md.results.TransientSolution(end).Vx;
md.initialization.vy = md.results.TransientSolution(end).Vy;
md.initialization.vz= md.results.TransientSolution(end).Vz;

md.friction.coefficient=30.0*ones(md.mesh.numberofvertices,1);

md.cluster=generic('name',oshostname(),'np',1);
% Set which control message you want to see #help verbose
%->
md.verbose=verbose('convergence',true);
% set the transient model to ignore the thermal model
% #md.transient
%->
md.transient.isthermal=0;
% define the timestepping scheme
% everything here should be provided in years #md.timestepping
% give the length of the time_step (4 years)
%->
md.timestepping.time_step=0.01;
% give final_time (20*4 years time_steps)
%->
md.timestepping.final_time = 500*md.timestepping.time_step;
% Solve #help solve
% we are solving a TransientSolution
%->
md=solve(md,'Transient');

save HO5yr_case3 md;

%
%% 4 Regular

md = loadmodel('./HO10yr_steady');



md.initialization.vx = md.results.TransientSolution(end).Vx;
md.initialization.vy = md.results.TransientSolution(end).Vy;
md.initialization.vz= md.results.TransientSolution(end).Vz;

%Calculate an ice thickness
thick = zeros(30,200);

for loop = 1:30 %Parabolic valley
    thick(loop,:) =  10+400-abs(loop-15.5).^2*400./210.25;
end
bed = ones(30,200).*(3.5:3.5:700);
bed2 = bed;
%add landslide to the topography
for xloop = 1:size(bed,2)
    for yloop = 1:size(bed,1)
        if (abs(yloop-size(bed,1))+abs(xloop-110)) < 14
            t=2;
        end
        
        bed(yloop,xloop) = bed(yloop,xloop) + max(0,250-17.857*(abs(yloop-0)+abs(xloop-110)));
        
    end
end
%Define the geometry of the simulation #md.geometry
surface_g = bed + thick;
md.geometry.surface = repmat(surface_g(:),2,1);
%base is [surface-1000+500*sin(x*2*pi/L).*sin(y*2*pi/L)]
%L is the size of the side of the square #max(md.mesh.x)-min(md.mesh.x)
%->

md.geometry.base=repmat(bed(:),2,1);
%thickness is the difference between surface and base #md.geometry
%->
md.geometry.thickness=repmat(thick(:),2,1);

md.cluster=generic('name',oshostname(),'np',1);
% Set which control message you want to see #help verbose
%->
md.verbose=verbose('convergence',true);
% set the transient model to ignore the thermal model
% #md.transient
%->
md.transient.isthermal=0;
% define the timestepping scheme
% everything here should be provided in years #md.timestepping
% give the length of the time_step (4 years)
%->
md.timestepping.time_step=0.01;
% give final_time (20*4 years time_steps)
%->
md.timestepping.final_time = 500*md.timestepping.time_step;
% Solve #help solve
% we are solving a TransientSolution
%->
md=solve(md,'Transient');

save HO5yr_case4 md;

%
%% 5 Regular, increased slippyness

md = loadmodel('./HO10yr_steady');

md.initialization.vx = md.results.TransientSolution(end).Vx;
md.initialization.vy = md.results.TransientSolution(end).Vy;
md.initialization.vz= md.results.TransientSolution(end).Vz;

md.friction.coefficient=20.0*ones(md.mesh.numberofvertices,1);

%Calculate an ice thickness
thick = zeros(30,200);

for loop = 1:30 %Parabolic valley
    thick(loop,:) =  10+400-abs(loop-15.5).^2*400./210.25;
end
bed = ones(30,200).*(3.5:3.5:700);
%add landslide to the topography
for xloop = 1:size(bed,2)
    for yloop = 1:size(bed,1)
        if (abs(yloop-size(bed,1))+abs(xloop-110)) < 14
            t=2;
        end
        
        bed(yloop,xloop) = bed(yloop,xloop) + max(0,250-17.857*(abs(yloop-0)+abs(xloop-110)));
        
    end
end
%Define the geometry of the simulation #md.geometry
%surface is [-x*tan(0.5*pi/180)] #md.mesh
%->
surface_g = bed + thick;
md.geometry.surface = repmat(surface_g(:),2,1);
%base is [surface-1000+500*sin(x*2*pi/L).*sin(y*2*pi/L)]
%L is the size of the side of the square #max(md.mesh.x)-min(md.mesh.x)
%->

md.geometry.base=repmat(bed(:),2,1);
%thickness is the difference between surface and base #md.geometry
%->
md.geometry.thickness=repmat(thick(:),2,1);

md.cluster=generic('name',oshostname(),'np',1);
% Set which control message you want to see #help verbose
%->
md.verbose=verbose('convergence',true);
% set the transient model to ignore the thermal model
% #md.transient
%->
md.transient.isthermal=0;
% define the timestepping scheme
% everything here should be provided in years #md.timestepping
% give the length of the time_step (4 years)
%->
md.timestepping.time_step=0.01;
% give final_time (20*4 years time_steps)
%->
md.timestepping.final_time = 500*md.timestepping.time_step;
% Solve #help solve
% we are solving a TransientSolution
%->
md=solve(md,'Transient');

save HO5yr_case5 md;

%
%% 6 Regular, decreased slippyness

md = loadmodel('./HO10yr_steady');

md.initialization.vx = md.results.TransientSolution(end).Vx;
md.initialization.vy = md.results.TransientSolution(end).Vy;
md.initialization.vz= md.results.TransientSolution(end).Vz;

md.friction.coefficient=30.0*ones(md.mesh.numberofvertices,1);

%Calculate an ice thickness
thick = zeros(30,200);

for loop = 1:30 %Parabolic valley
    thick(loop,:) =  10+400-abs(loop-15.5).^2*400./210.25;
end
bed = ones(30,200).*(3.5:3.5:700);
%add landslide to the topography
for xloop = 1:size(bed,2)
    for yloop = 1:size(bed,1)
        if (abs(yloop-size(bed,1))+abs(xloop-110)) < 14
            t=2;
        end
        
        bed(yloop,xloop) = bed(yloop,xloop) + max(0,250-17.857*(abs(yloop-0)+abs(xloop-110)));
        
    end
end
%Define the geometry of the simulation #md.geometry
%surface is [-x*tan(0.5*pi/180)] #md.mesh
%->
surface_g = bed + thick;
md.geometry.surface = repmat(surface_g(:),2,1);
%base is [surface-1000+500*sin(x*2*pi/L).*sin(y*2*pi/L)]
%L is the size of the side of the square #max(md.mesh.x)-min(md.mesh.x)
%->

md.geometry.base=repmat(bed(:),2,1);
%thickness is the difference between surface and base #md.geometry
%->
md.geometry.thickness=repmat(thick(:),2,1);

md.cluster=generic('name',oshostname(),'np',1);
% Set which control message you want to see #help verbose
%->
md.verbose=verbose('convergence',true);
% set the transient model to ignore the thermal model
% #md.transient
%->
md.transient.isthermal=0;
% define the timestepping scheme
% everything here should be provided in years #md.timestepping
% give the length of the time_step (4 years)
%->
md.timestepping.time_step=0.01;
% give final_time (20*4 years time_steps)
%->
md.timestepping.final_time = 500*md.timestepping.time_step;
% Solve #help solve
% we are solving a TransientSolution
%->
md=solve(md,'Transient');

save HO5yr_case6 md;