

%Make sure to addpath the paths required for ISSM here

% Cluster options (example for PACE)
optionloadonly       = 0;
outname0             = 'salKgb';
codepath             = ''; %this should be the path to the ISSM bin/ directory
execpath             = ''; %this should be the path to the execution directory
loginname            = ''; %this should be your login name
coresrequested       = 10;
memrequested         = 2; %GB
timerequested        = 400*60; %minutes

%%% Choose glacier %%%
glacier = 'KogeBugt';

curdir    = pwd();
DirInit   = [curdir,'/InputsGlaDS/'];
DirRunoff = [curdir,'/InputsGlaDS/'];
DirMouls  = [curdir,'/InputsGlaDS/'];
nMouls    = 30;

%%% Run params %%%
deltat          = 2*1/(365.25*24);
outputfreq      = 36;
ntotyr          = 40;
endtime         = ntotyr;
reqoutputs      = {};

%%% Params for hydrology %%%
englVoidRatio   = 1e-3; %W2013,dF2018
bumpHeight      = 0.1; %W2013,dF2018 [m]
cavitySpace     = 2; %W2013,dF2018 [m]
sheetK          = 0.005; %dF2018 [m7/4 kg-1/2] (note unit errors in E2022supp)
channelK        = 0.1; %W2013,dF2018 [m3/2 kg-1/2]
widthShCh       = 2.0; %W2013,dF2018 [m] 
initWaterCol    = 0.03; %E2022 [m] (from runme file)
limitonN        = 0.06; %E2022

%%% Load %%%
fileInit  = [DirInit,'init',glacier,'04Mar2023.mat'];
fileRnf   = [DirRunoff,'runoffTot',glacier,'19702009monthlydEBM.csv'];
fileMouls = [DirMouls,'locMoulins',int2str(nMouls),glacier,'04Mar2023.csv'];
md        = loadmodel(fileInit);
% Runoff input time series %
datarnf   = csvread(fileRnf);
tmrnf     = datarnf(1,:)-datarnf(1,1); %starting timing of runoff at 0
totrnfinp = datarnf(2,:);
% Moulins locations %
datamouls = csvread(fileMouls);
moulsxx   = datamouls(:,1);
moulsyy   = datamouls(:,2);
iimoul    = zeros(md.mesh.numberofvertices,1);
flbnd     = md.mesh.vertexonboundary; %vertices on boundary of the domain
for(ii=[1:nMouls])
    dists              = sqrt((md.mesh.x-moulsxx(ii)).^2+(md.mesh.y-moulsyy(ii)).^2);
    dists(iimoul==1)   = 1e20; %in order to find closest gridpoint that does not have a moulin yet
    dists(flbnd==1)    = 1e20; %exclude vertices on boundary
    [dummy,vrt]        = min(dists);
    iimoul(vrt)        = 1;
end

%%% Set up GLADS %%%
md.hydrology                       = hydrologyglads();
% Geometry params %
md.hydrology.englacial_void_ratio  = englVoidRatio;
md.hydrology.bump_height           = bumpHeight*ones(md.mesh.numberofvertices,1);
md.hydrology.cavity_spacing        = cavitySpace;
% Sheet params %
md.hydrology.sheet_conductivity    = sheetK*ones(md.mesh.numberofvertices,1);
md.initialization.watercolumn      = initWaterCol*ones(md.mesh.numberofvertices,1);
% Channel params %
md.hydrology.ischannels            = 1;
md.hydrology.channel_sheet_width   = widthShCh;
md.hydrology.channel_conductivity  = channelK*ones(md.mesh.numberofvertices,1);
% Basal melt params %
md.hydrology.melt_flag                    = 1;
md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
% Moulins params %
moulininput    = zeros(md.mesh.numberofvertices+1,numel(tmrnf));
for(tt=[1:numel(tmrnf)])
    indivinput = totrnfinp(tt)/nMouls;
    moulininput(1:end-1,tt) = indivinput*iimoul;
    moulininput(end,tt)     = tmrnf(tt);
end
md.hydrology.moulin_input = moulininput;
% Initialize hydraulic potential %
% note: GlaDS runs start with Pw=(1/2)*Pice (Ehrenfeucht et al., 2022)
rhoi                                  = md.materials.rho_ice;
rhofw                                 = md.materials.rho_freshwater;
grv                                   = md.constants.g;
md.initialization.hydraulic_potential = (1/2)*(rhoi*grv*md.geometry.thickness);
% Boundary conditions %
md.hydrology.neumannflux          = zeros(md.mesh.numberofelements,1);
md.hydrology.spcphi               = NaN(md.mesh.numberofvertices,1);
indsCloseGL                       = find(md.mask.ocean_levelset<=50); 
md.hydrology.spcphi(indsCloseGL)  = rhoi*grv*md.geometry.thickness(indsCloseGL)+rhofw*grv*md.geometry.base(indsCloseGL);

%%% Additional params %%%
md.initialization.pressure = md.results.StressbalanceSolution.Pressure;

%%% Fix spc at some nodes (13Feb20223) %%%
if(strcmp(glacier,'Kogebugt'))
   inds = [131,132,176,177,947,1766,2033,2925,3123,3567,3594,3666,4473,4852,4892];   
   md.hydrology.spcphi(inds) = rhoi*grv*md.geometry.thickness(inds)+rhofw*grv*md.geometry.base(inds);
end

%%% Set up transient %%%
md.inversion.iscontrol         = 0;
md.transient                   = deactivateall(md.transient);
md.transient.ishydrology       = 1;
md.settings.output_frequency   = outputfreq;
md.timestepping.start_time     = 0;
md.timestepping.time_step      = deltat;
md.timestepping.final_time     = endtime;
md.transient.requested_outputs = reqoutputs;
md.results                     = {}; %remove old results

%%% Solving %%%
md.settings.waitonlock = 0;
md.cluster = pace('np',coresrequested,'login',loginname,'codepath',codepath,'executionpath',execpath,'mem',memrequested,'time',timerequested);
md.miscellaneous.name = outname0; %Run name
md.verbose.solution   = true;
md.toolkits.DefaultAnalysis=bcgslbjacobioptions(); %from runme_GlaDS of E2022 
md.toolkits.RecoveryAnalysis=mumpsoptions(); %from runme_GlaDS of E2022
md = solve(md,'Transient','runtimename',false,'loadonly',optionloadonly);
if(optionloadonly==1)
    md = loadresultsfromcluster(md);
    savename = [execpath,outname0,'.mat'];
    save(savename,'md','-v7.3') %save md
    %rmdir([execpath,outname0],'s'); %clean-up
end


if(isfile('s1init.mat'))
    delete 's1init.mat';
end

