function MISMIPplus(resolution,schoof)

% MISMIP+ experiment with Kori-ULB v0.9
% Details in Cornford et al (2020)

name=[int2str(resolution),'km_schoof',int2str(schoof)];

%% Initial ice sheet creation

ctr.runmode=1;

if resolution==2
    ctr.delta=2e3;
    ctr.imax=23; % need number of cells + 2
    ctr.jmax=352; % need number of cells + 1
elseif resolution==1
    ctr.delta=1e3;
    ctr.imax=43; % need number of cells + 2
    ctr.jmax=702; % need number of cells + 1
elseif resolution==4
    ctr.delta=4e3;
    ctr.imax=13; % need number of cells + 2
    ctr.jmax=177; % need number of cells + 1
elseif resolution==8
    ctr.delta=8e3;
    ctr.imax=8; % need number of cells + 2
    ctr.jmax=89; % need number of cells + 1
end
ctr.m=3;
ctr.mismip=1; % indicates that BCs need to be applied
ctr.SSA=1;
ctr.shelf=1;
ctr.Ao=2e-17; 
ctr.shelftune=1;
ctr.Asin=zeros(ctr.imax,ctr.jmax)+1e-6; 
ctr.schoof=schoof;

if ctr.schoof==0
    ctr.dt=0.5;
    ctr.nsteps=16001;
else
    ctr.dt=0.1; 
    ctr.nsteps=80001; 
end

Li=(ctr.imax-2)*ctr.delta;
Lj=(ctr.jmax-2)*ctr.delta;
[X,Y]=meshgrid(-ctr.delta:ctr.delta:Lj,-ctr.delta:ctr.delta:Li);
B=Bx(X)+By(Y,4.0e3,5.0e2,24.0e3);
B=max(B,-720.0);
% periodic BCs
B(1,:)=B(3,:);
B(ctr.imax,:)=B(ctr.imax-2,:);
H=zeros(ctr.imax,ctr.jmax)+10;
Mb=zeros(ctr.imax,ctr.jmax)+0.3;
Ts=zeros(ctr.imax,ctr.jmax)-5.0;
save('MISMIPplusin','B','H','Mb','Ts');

KoriModelAll('MISMIPplusin',['mismiptestin1_',name],ctr);

%% Refining grounding line position on upward bed slope: 

ctr.timeslice=1;
ctr.snapshot=10;
% value of the rate factor A is adapted to obtain similar steady state
% positions on the retrograde slope
if schoof==0
    if resolution==1
        ctr.Ao=4.5e-17; 
    elseif resolution==2
        ctr.Ao=5e-17;
    elseif resolution==4
        ctr.Ao=5e-17;
    elseif resolution==8
        ctr.Ao=11e-17;
    end
else
    ctr.Ao=6e-17; 
    if resolution==4
        ctr.Ao=7.5e-17;
    elseif resolution==8
        ctr.Ao=6e-17;
    end
end
if ctr.schoof==0
    ctr.dt=0.5;
    ctr.nsteps=16001;
else
    ctr.dt=0.1;
    ctr.nsteps=80001;
end
KoriModelAll(['mismiptestin1_',name],['mismiptest1_',name],ctr);

%% Ice0: steady-state check for 200 yr

ctr.timeslice=1;
ctr.snapshot=20;
if ctr.schoof==0
    ctr.dt=0.5;
    ctr.nsteps=401;
else
    ctr.dt=0.1;
    ctr.nsteps=2001;
end
KoriModelAll(['mismiptest1_',name],['mismiptest1_ice0_',name],ctr);

%% Ice1r experiment for 100 year with melting

ctr.snapshot=10;
if ctr.schoof==0
    ctr.dt=0.5;
    ctr.nsteps=201;
else
    ctr.dt=0.1;
    ctr.nsteps=1001;
end
ctr.meltfunc=10;
ctr.meltfac=1;
ctr.BetaIter=ctr.nsteps;
KoriModelAll(['mismiptest1_',name],['mismiptest1r_',name],ctr);

%% Ice1rr experiment for 200 year with melting

ctr.snapshot=20;
if ctr.schoof==0
    ctr.dt=0.5;
    ctr.nsteps=401;
else
    ctr.dt=0.1;
    ctr.nsteps=2001;
end
KoriModelAll(['mismiptest1_',name],['mismiptest1rr_',name],ctr);

%% Ice1ra experiment starting from Ice1r without melt

ctr.snapshot=10;
if ctr.schoof==0
    ctr.dt=0.5;
    ctr.nsteps=201;
else
    ctr.dt=0.1;
    ctr.nsteps=1001;
end
ctr.meltfunc=0;
KoriModelAll(['mismiptest1r_',name],['mismiptest1ra_',name],ctr);

%% Plot results

load(['mismiptest1_',name]);
Li=(ctr.imax-2)*ctr.delta/1e3;
Lj=(ctr.jmax-2)*ctr.delta/1e3;
[Xo,Yo] = meshgrid(-ctr.delta/1e3:ctr.delta/1e3:Lj,-ctr.delta/1e3:ctr.delta/1e3:Li);
figure;
contour(Xo,-Yo,MASK,[1],'k','linewidth',2);
hold on;
load(['mismiptest1r_',name]);
contour(Xo,-Yo,MASK,[1],'r','linewidth',2);
load(['mismiptest1rr_',name]);
contour(Xo,-Yo,MASK,[1],'g','linewidth',2);
load(['mismiptest1ra_',name]);
contour(Xo,-Yo,MASK,[1],'b','linewidth',2);
xlim([300 550]);
xlabel('x (km)');
ylabel('y (km)');
grid on;
legend('Ice0','Ice1r','Ice1rr','Ice1ra');

figure;
% Ice0
time=[0:10:200];
for i=1:21
    load(['mismiptest1_ice0_',name,'_',num2str(i-1,'%03i')])
    z=find(glMASK(1,:)==2);
    plot(time(i),X(1,z)/1000,'.','color','k'); hold on
end
% 1rr
time=[0:10:200];
for i=1:21
    load(['mismiptest1rr_',name,'_',num2str(i-1,'%03i')])
    z=find(glMASK(1,:)==2);
    plot(time(i),X(1,z)/1000,'.','color','g'); hold on
end
% 1r
time=0:10:100;
for i=1:11
    load(['mismiptest1r_',name,'_',num2str(i-1,'%03i')])
    z=find(glMASK(1,:)==2);
    plot(time(i),X(1,z)/1000,'.','color','b'); hold on
end
% 1ra
time=[100:10:200];
for i=1:11
    load(['mismiptest1ra_',name,'_',num2str(i-1,'%03i')])
    z=find(glMASK(1,:)==2);
    plot(time(i),X(1,z)/1000,'.','color','r'); hold on
end
ylim([320 460])


end


function [Bxx]=Bx(x)
    B0=-150;
    B2=-728.8;
    B4=343.91;
    B6=-50.57;
    xx=x./300.0e3;
    xx2=xx.^2;
    xx4=xx2.^2;
    xx6=xx4.*xx2;
    Bxx=B0+B2*xx2+B4*xx4+B6*xx6;
end

function [Byy]=By(y,fc,dc,wc)
    Byy=dc./(1.0+exp(-2.0*(y-wc)/fc)) + dc./(1.0+exp(2.0*(y+wc)./fc));
end
