%----------------------------------------------------------------------
%Outputs
Var.qplot=1;                            %Display graphics? (1=yes, 0=no)
Var.save=1;                             %Save results?  (1=yes, 0=no)

%----------------------------------------------------------------------
%Adjustable parameters
Var.phi=[0:1:85];                             %Latitudes bounding boxes
Var.phi2=(Var.phi(1:end-1)+Var.phi(2:end))/2; %mid-point latitude of each box
Var.num=length(Var.phi2);                     %number of EBM boxes
Var.dphi=diff(Var.phi2(1:2));                 %length of box in degrees
Var.dt=24*3600;                                %Time step for EBM, seconds 
Var.ice_step=500;                             %Integation time for icesheet, years
Var.t_begin=-2.1e6;                            %Start time
Var.t_end=-0e6;                                %End time

%----------------------------------------------------------------------
%EBM constants
Var.rho_i=910;                          %Density of ice
Var.rho_w=1000;                         %Density of water
Var.rho_mantle=3300;                    %Density of the mantle
Var.Cp= 2100;                           %Specific heat capacity, J kg^-1 K^-1
Var.Ki=4;                               %J/m K s; thermal conductivity of ice
Var.k= 0.03;                            %Mass absorption coefficient for water m^2 kg^-1
Var.g=9.8;                              %Acceleration due to gravity, m s^-2
Var.sigma_r=5.67e-8;                    %Stefan Boltzmann constant, W m^-2 K^-4
Var.Cd=0.0011;                          %Drag coefficient
Var.Lv=2.5e6;                           %Latent heat of vaporization, J kg^-1
Var.Lm=3.34e5;                          %Latent heat of melting, J kg^-1
Var.Ls=2.84e6;                          %Latent heat of sublimation, J kg^-1
Var.a=6.37e6;                           %Radius of the earth, m
Var.K=273.15;                           %Melting point

%EBM Parameters
Var.min_interior_thickness=10;          %Minimum thickness of interior layer
Var.Ci=10*Var.rho_i*Var.Cp;   
Var.Cs=5*Var.rho_i*Var.Cp;              %Heat capacity of the surface, 1m*900kg/m^3*2100J/(kg C)=4e6 J/(m^2 C)
Var.Ca=5000*1.5*1000;                   %Heat capacity of atmosphere, 5000m*1.5kg/m^3*1000J/(kg C)=7.5e6 J/(m^2 C)
Var.Ksensible=5;                        %Heat flux coefficient between ground and atmosphere
Var.Hm=5;                               %Height of middle of atmosphere layer
Var.Htm=2;                              %Top to middle atmospheric layer distance, km
Var.alpha_l=0.3;                        %land albedo
Var.alpha_i=0.8;                        %ice albedo
Var.A=0.2;                              %Absorption of atmosphere 
Var.R=0.3;                              %Reflection of atmosphere
Var.T=0.5;                              %Transmission of atmosphere 
Var.eps_a=0.85;                          %Longwave atmospheric emissivity
Var.P=1/(365*24*3600);                %m/s
Var.lr = 6.5;                             %Moist adiabatic lapse rate

%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%Ice-sheet parameters

Var.dt_ice=5;                           %time (years)
Var.n=3;                                %exponent for Glen's law
Var.Nx=length(Var.phi2)-1; 
Var.dphi=diff(Var.phi(1:2))*pi/180;    %length of box in degrees
Var.dx=Var.a*Var.dphi/pi;               %meters 

Var.a1=3.615e-13; % Pa^-3 s^-1; T<263
Var.a2=1.733e3;   % Pa^-3 s^-1; T>=263
Var.Q1=6.0e4;     % J/mol; T<263
Var.Q2=13.9e4;    % J/mol; T>=263
Var.RgasPB=8.314; % J/(mol K);
Var.T_ice=200;    %Kelvin
Var.Aice = Var.a1*exp(-Var.Q1./(Var.RgasPB*Var.T_ice)).*(Var.T_ice<263) + Var.a2*exp(-Var.Q2./(Var.RgasPB*Var.T_ice)).*(Var.T_ice>=263);
Var.Gam=2*(Var.rho_i*Var.g)^Var.n*Var.Aice;   % overall constant in continuity eqn
Var.Rice=0.5*Var.dt_ice*365*24*3600/(Var.dx*Var.dx);     %related to stability for continuity eqn
Var.Tb=5000;            %Bed depression timescale, years
Var.Heq=0;              %Equilibrium height of the ice-free surface, meters

%Sedimentalogical constants
Var.rho_s=2390;                  %saturated bulk sediment density, kg/m^3
Var.rho_b=3370;                  %bedrock density, kg/m^3
Var.phi_sed=22*pi/180;           %angle of internal friction 
Var.Do=(7.9e-7)/(365*24*3600);   %reference deformation rate in per second
Var.ns=1.25;                     %exponent
Var.bs=(Var.rho_s-Var.rho_w)*Var.g*tan(Var.phi_sed); %change in sediment shear strength with z.
Var.uo=3e9;                       %sediment reference viscostiy
Var.hs = 10;               



%-------------------------------------------------------------------
%Initialize state variables

%ebm
State.year=Var.t_begin;                       %Initial year
State.Ti=linspace(Var.K+25,Var.K-10,Var.num); %Temperature of the interior, K
State.Ta=linspace(Var.K,Var.K-50,Var.num);    %Temperature of the atmopshere, K
State.Ts=linspace(Var.K+35,Var.K-20,Var.num); %Temperature of the surface, K
State.Da=linspace(-50,50,Var.num);            %Atmospheric heat fluxS
State.M=zeros(1,Var.num); 

%ice model
State.hh=zeros(1,Var.num);                        %height of ice-sheet above bed
State.H=ones(1,Var.num).*Var.Heq;                 %Height of bed
State.Flow=zeros(1,Var.num);                      %Flow of ice, calculated from ice model