%function IFace=getIface(Ny,Nx,NSURF,NLAY)
% let's obtain the interface/s elevation/s from the density values
% first we will obtain them for the area where the data is provided (Zuid
% Holland) and afterwards we will put them into our model. Our area of
% interest is not fully contained in the given area but that will be solved
% later. In the meanwhile the elevation at the points which are not
% contained will be given a value NaN

pth='Z:\tolsthoorn On My Mac\GRWMODELS\NHI_and_AGV\mocdens3d';

load([pth filesep 'mocdense.mat']);

moc=rmfield(moc,'DZ');
moc=rmfield(moc,'SFLAG');
clear bas bcf

%% Remove all -1 values which are north sea water

moc.CHLOR(moc.CHLOR==-1)=max(moc.CHLOR(:));

%% Remove inversions by copying higher chlorides of overlying layer + DCHLOR
% this is primitive but soit

DCHLOR=10;


L1=zeros(size(moc.CHLOR));
L2=zeros(size(moc.CHLOR));
k=0;
while 1
    k=k+1;
    L1(:)=0; L2(:)=0;
    I=(moc.CHLOR(:,:,2:end)<=moc.CHLOR(:,:,1:end-1));
    fprintf('iteration %d, number of inversion points: %d\n',k,sum(I(:)));
    if any(I(:)),
        L1(:,:,1:end-1)=I;
        L2(:,:,2:end  )=I;
        moc.CHLOR(logical(L2))=moc.CHLOR(logical(L1))+DCHLOR;
    else
        break;
    end
end
clear L1 L2 I

%%
IFACEVALS=[300 1000 10000];

IFACE=zeros(moc.NROW,moc.NCOL,length(IFACEVALS));
CHLOR=NaN(moc.NROW,moc.NCOL,moc.NLAY+2);
CHLOR(:,:,2:end-1)=moc.CHLOR;
CHLOR(:,:,1  )=min(moc.CHLOR(:))   -1;  % actual minimum=5
CHLOR(:,:,end)=max(moc.CHLOR(:))+1000;  % actual maximum = 18700

Z=NaN(size(CHLOR)); Z(:,:,2:end-1)=moc.ZM; Z(:,:,1)=1000; Z(:,:,end)=-1000;
for ix=1:moc.NCOL
    fprintf('.')
    if rem(ix,50)==0, fprintf('%d/%d\n',ix,moc.NCOL); end
    for iy=1:moc.NROW
        IFACE(iy,ix,:)=interp1(squeeze(CHLOR(iy,ix,:)),squeeze(Z(iy,ix,:)),IFACEVALS);
    end
end

figure
for i=1:size(IFACE,3)
    subplot(3,1,i)
    contour(xm,ym,IFACE(:,:,i));
    title(sprintf('Interface %.0f mg/L',IFACEVALS(i)));
end

%IZ=interp3(XM,YM,moc.CHLOR,ZM,xm,ym,500); % runs out of memory




% switch NSURF
%     
%     
%     
%     
%     case 1
%         zeta=NaN(mocm.NROW,mocm.NCOL,NSURF);
%         for i=1:mocm.NROW
%             for j=1:mocm.NCOL
%                 zeta(i,j,NSURF)=interp1(DENS(:,i,j),zm(:,i,j),10000); %The interface is interpolated with interp1
%             end
%         end
%     case 2
%         zeta=NaN(mocm.NROW,mocm.NCOL,NSURF);
%         for i=1:mocm.NROW
%             for j=1:mocm.NCOL
%                 zeta(i,j,1)=interp1(DENS(:,i,j),zm(:,i,j),10000);
%                 zeta(i,j,2)=interp1(DENS(:,i,j),zm(:,i,j),15000);
%             end
%         end
%     case 3
%         zeta=NaN(mocm.NROW,mocm.NCOL,NSURF);
%         for i=1:mocm.NROW
%             for j=1:mocm.NCOL
%                 zeta(i,j,1)=interp1(DENS(:,i,j),zm(:,i,j),10000);
%                 zeta(i,j,2)=interp1(DENS(:,i,j),zm(:,i,j),12500);
%                 zeta(i,j,3)=interp1(DENS(:,i,j),zm(:,i,j),15000);
%             end
%         end
% end
% 
% interface=NaN(Ny,Nx,NSURF);
% 
% interface(25:200,1:180,:)=zeta(1:176,221:end,:);
% 
% IFace=IfaceIS(interface,Ny,Nx,NLAY,NSURF,TOP,BOTTOM);

fclose('all');
