function [ZETA,xGrMD,yGrMD,xm,ym]=SWI_getIface(Chlorides,mocdensematfile)
% [ZETA,xGrMd,yGrMd,xm,ym]=SWI_getIface(Chlorides,mocdensematfile)
% the xGrMd etc come from the stored and loaded mocdensematfile
% ZETA are surfaces interpolated from the MOCDENSE input given in
% mocedensematifle. This file was generated by readMOCDENSE
% Chlorides are the clorides for which the IFaces are desired, these
% interaces are interpolated here using teh specified chlorieds and the
% mocdensematfile
% xGr,yGr,xm,ym are grid and cell coordinates stored in the mocdensefile
% TO 090814


% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later


fprintf('SWI_getIFace now running ...\n')

%% Loading chloride data
fprintf('loading %s\n',mocdensematfile);
load(mocdensematfile);

%% save memory be removing unwanted items in the matfile
moc=rmfield(moc,'DZ');
moc=rmfield(moc,'SFLAG');
clear bas bcf

%% Remove all -1 values which are north sea water
maxchlor=max(moc.CHLOR(:));

xx =  1.0e+005 *[0.8776  1.0569 1.2959 1.2109 0.9787 0.8569];
yy =  1.0e+005 *[4.3319  4.3535 4.7015 4.8337 4.8256 4.4668];
[XM,YM]=meshgrid(xm,ym);
IN=zeros(size(moc.CHLOR));
in=inpolygon(XM,YM,xx,yy);
IN(:,:,1)=in;

fprintf('replacing -1 values by max %g mg/L\n',maxchlor);
moc.CHLOR(~IN & moc.CHLOR<0)=maxchlor;

%% Squeeze down inversions by copying higher chlorides of overlying layer + DCHLOR
% this is primitive but soit, but so be it. We can't have inversions in SWI
% and we cannot interpolate without a monotonic dpeth increasing chloried
% concentration
fprintf('sqeezing downward out inversions, this takes some time ...\n');

DCHLOR=10;

L1=zeros(size(moc.CHLOR));  % to store logicals
L2=zeros(size(moc.CHLOR));  % to shore logicals
k=0;
while 1
    k=k+1;
    L1(:)=0; L2(:)=0;
    I=(moc.CHLOR(:,:,2:end)<=moc.CHLOR(:,:,1:end-1));   % which cell had chloride lower than upper neighor?
    fprintf('iteration %d, number of inversion points: %d\n',k,sum(I(:)));
    if any(I(:)),   % as long as there are any cells with chloride inversions
        L1(:,:,1:end-1)=I;   % note which cells have a lower neighbor with lower concentration
        L2(:,:,2:end  )=I;   % note which cells have an upper neighbor with higher concentration
        moc.CHLOR(logical(L2))=moc.CHLOR(logical(L1))+DCHLOR; % replace lower by higher + DCHLOR
    else
        break;  % repeat loop until all inversions have been squeezed out
    end
end

%% clear unnesessary memory usage
clear L1 L2 I

%%  interpolate vertically to get the z-position of all wanted interfaces
%Chlorides=[300 1000 10000];
fprintf('Interpolating to get the z=positions of all desired interfaces ...\n');

%% Add an upper layer with minimum chlored
CHLOR=NaN(moc.NROW,moc.NCOL,moc.NLAY+2); % allocate to store all chlorides + an upper an a lower face outside the model
                                         % to make sure that we will always obtain valid interpolations and no NaNs
CHLOR(:,:,2:end-1)=moc.CHLOR;            % insert existing chlorides in cell centers
CHLOR(:,:,1  )=min(moc.CHLOR(:))   -1;  % add chlorides to artificial upper layer (actual minimum=5)
CHLOR(:,:,end)=max(moc.CHLOR(:))+1000;  % add chlorides to artificial lower layer (actual maximum = 18700)

%% Make z-coordinates with an extra above and one below the actual model
% (1000 m up and 1000 m below), to make sure we always get leag
% z-coordinaties from our interpolation

Z=NaN(size(CHLOR)); Z(:,:,2:end-1)=moc.ZM; Z(:,:,1)=1000; Z(:,:,end)=-1000;

%% Use interp1 for every row-col to interpolate the interface positions vertically
% The resuls is a set of 3D interafaces. SWI still needs intersection with
% the layers, which will be done inside write SWI and need not be stored
% intermediately

ZETA=zeros(moc.NROW,moc.NCOL,length(Chlorides)); % allocate
for ix=1:moc.NCOL
    fprintf('.')
    if rem(ix,50)==0, fprintf('%d/%d\n',ix,moc.NCOL); end
    for iy=1:moc.NROW
        ZETA(iy,ix,:)=interp1(squeeze(CHLOR(iy,ix,:)),squeeze(Z(iy,ix,:)),Chlorides);
    end
end

%% Show the results just to see if we did the job right

fprintf('Showing resulting interfaces as contours ...\n'); 
figure
NPLN=size(ZETA,3);  % SWI manual uses NPLN for the number of surfaces/interfaces
for ipln=1:NPLN
    subplot(NPLN,1,ipln)
    contour(xm,ym,ZETA(:,:,ipln));
    title(sprintf('Interface %.0f mg/L',Chlorides(ipln)));
end
%% Close file and done

fclose('all');

fprintf('... done (thanks for your patience\n');
