
%%%
% Reading and saving GlaDS results
%%%

%%% Choose glacier %%%
glacier    = 'Hmb';
jump12step = false; %to save one timestep every two timesteps

curdir    = pwd();
DirRes    = ''; %this should be the path to the directory with the GlaDS output mat files
DirSave   = [curdir,'/OutputGlaDS/'];
nMouls    = 30;
FileRes   = ['sal',glacier,'pc.mat'];
LocRes    = [DirRes,FileRes];
Saving    = true;
limitonN  = 0.06; %E2022

%%% Load %%%
md        = loadmodel(LocRes);
if(strcmp(glacier,'Ptm'))
    glacfullname = 'Petermann';
elseif(strcmp(glacier,'Hlm'))
    glacfullname = 'Helheim';
elseif(strcmp(glacier,'Jkb'))
    glacfullname = 'Jakobshavn';
elseif(strcmp(glacier,'Kgl'))
    glacfullname = 'Kangerlussuaq';
elseif(strcmp(glacier,'Upv'))
    glacfullname = 'Upernavik';
elseif(strcmp(glacier,'Hmb'))
    glacfullname = 'Humboldt';
elseif(strcmp(glacier,'Kgb'))
    glacfullname = 'Kogebugt';
elseif(strcmp(glacier,'Rsl'))
    glacfullname = 'Russell';
end

%%% Span of domain %%%
xx = md.mesh.x;
yy = md.mesh.y;

%%% Constants %%%
rhoi  = md.materials.rho_ice;
rhofw = md.materials.rho_freshwater;
grv   = md.constants.g;

%%% Ice thickness and velocity %%%
icethick = md.geometry.thickness; %stand-alone GlaDS
icevel   = md.initialization.vel; %stand-alone GlaDS
pice     = rhoi*grv*icethick;

%%% Extract hydrology results %%%
peff0      = [md.results.TransientSolution(:).EffectivePressure];
peff       = peff0;
% inds       = peff0<limitonN*pice;
% peff(inds) = limitonN*pice(inds);
pwater     = rhoi*grv*icethick-peff;
sdevpwater = std(pwater,0,2);
mnpwater   = mean(pwater,2);
phi        = [md.results.TransientSolution(:).HydraulicPotential];
sdevphi    = std(phi,0,2);
mnphi      = mean(phi,2);
sheetthick = [md.results.TransientSolution(:).HydrologySheetThickness];
timing = [md.results.TransientSolution(:).time];


ntot   = numel(timing);
[dummy,ntinput]  = size(md.hydrology.moulin_input);
tminput          = md.hydrology.moulin_input(end,:);
totmoulininput   = [];
totmoulinfluxrun = [];
totchannelarea   = [];
totchanneldisch  = [];
for(tt=[1:ntinput])
    totmoulininput(end+1) = sum(md.hydrology.moulin_input(1:end-1,tt));
end
for(tt=[1:ntot])
    ttrun = timing(tt);
    totmoulinfluxrun(end+1) = interp1(tminput,totmoulininput,ttrun,'linear',min(totmoulininput));
    totchannelarea(end+1)   = sum(md.results.TransientSolution(tt).ChannelArea); 
    totchanneldisch(end+1)  = sum(abs(md.results.TransientSolution(tt).ChannelDischarge)); 
end
   
if(Saving)
    if(jump12step)
       timing           = timing(1:2:end);
       totmoulinfluxrun = totmoulinfluxrun(1:2:end);
       phi              = phi(:,1:2:end);
    end
    outx = md.mesh.x;
    outy = md.mesh.y;
    % Get moulin locations %
    [dummy,ind] = max(totmoulininput); 
    vals        = md.hydrology.moulin_input(1:end-1,ind); 
    vrtmoul     = find(vals>0);
    locmouls    = zeros(md.mesh.numberofvertices,1);
    locmouls(vrtmoul) = 1;
    % Get spc vertices %
    isspc       = zeros(md.mesh.numberofvertices,1);
    isspc(find(isnan(md.hydrology.spcphi)==false)) = 1;
    % Get ocean vertices %
    isocn       = zeros(md.mesh.numberofvertices,1);
    isocn(find(md.mask.ocean_levelset>=0)) = 1;    
    % Netcdf %
    outpathnc   = [DirSave,FileRes(1:end-4),'.nc'];
    ncout       = netcdf.create(outpathnc,'NC_WRITE');
    % Define dimensions %
    space_dim    = netcdf.defDim(ncout,'dim_space',numel(outx));
    tm_dim       = netcdf.defDim(ncout,'dim_time',numel(timing));
    % Define variables %
    time_out     = netcdf.defVar(ncout,'time','NC_FLOAT',tm_dim);
    x_out      = netcdf.defVar(ncout,'x','NC_FLOAT',space_dim);
    y_out      = netcdf.defVar(ncout,'y','NC_FLOAT',space_dim);
    isspc_out  = netcdf.defVar(ncout,'isSPC','NC_FLOAT',space_dim);
    iceH_out   = netcdf.defVar(ncout,'IceThickness','NC_FLOAT',space_dim);
    icevel_out = netcdf.defVar(ncout,'IceVel','NC_FLOAT',space_dim);
    phi_out    = netcdf.defVar(ncout,'HydraulicPotential','NC_FLOAT',[tm_dim,space_dim]);
    mflux_out  = netcdf.defVar(ncout,'InfluxMoulins','NC_FLOAT',tm_dim);
    mlocs_out  = netcdf.defVar(ncout,'PresenceMoulin','NC_INT',space_dim);
    bed_out    = netcdf.defVar(ncout,'Bed','NC_FLOAT',space_dim);
    % Variable units %
    netcdf.putAtt(ncout,x_out,'units','m');
    netcdf.putAtt(ncout,y_out,'units','m');
    netcdf.putAtt(ncout,bed_out,'units','m');
    netcdf.putAtt(ncout,iceH_out,'units','m');
    netcdf.putAtt(ncout,mflux_out,'units','m3 s-1');
    netcdf.putAtt(ncout,isspc_out,'units','binary');
    netcdf.putAtt(ncout,mlocs_out,'units','binary');
    netcdf.putAtt(ncout,icevel_out,'units','m s-1');
    netcdf.putAtt(ncout,time_out,'units','year');
    netcdf.putAtt(ncout,phi_out,'units','Pa');
    netcdf.endDef(ncout);
    % Save varaibles %
    netcdf.putVar(ncout,x_out,md.mesh.x);
    netcdf.putVar(ncout,y_out,md.mesh.y);
    netcdf.putVar(ncout,isspc_out,isspc);
    netcdf.putVar(ncout,bed_out,md.geometry.bed);
    netcdf.putVar(ncout,time_out,timing);
    netcdf.putVar(ncout,iceH_out,icethick);
    netcdf.putVar(ncout,icevel_out,icevel);
    netcdf.putVar(ncout,phi_out,phi');
    netcdf.putVar(ncout,mflux_out,totmoulinfluxrun);
    netcdf.putVar(ncout,mlocs_out,locmouls);
    % Close %
    netcdf.close(ncout)
end

