clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read in and reshape files
fid = fopen('pickup.OLD.data','r','b');
pickupold = fread(fid,inf,'float64');
fclose(fid);

fid = fopen('pload.OLD.jmd95z','r','b');
ploadold = fread(fid,inf,'real*8');
fclose(fid);

fid = fopen('shelftopo.OLDPIG.bin','r','b');
stpigold=fread(fid,inf,'real*8');
fclose(fid);

fid = fopen('shelftopo.OLDROUND.bin','r','b');
stroundold=fread(fid,inf,'real*8');
fclose(fid);

nx=1;
ny=200;
nz=100;

delz=10;

dz = delz*ones(1,nz);
zgp1 = [0,cumsum(dz)];
zc = .5*(zgp1(1:end-1)+zgp1(2:end));
zg = zgp1(1:end-1);
dz = diff(zgp1);

s3=nx*ny*nz;
s2=nx*ny;

olduvel=reshape(pickupold(1:s3),[nx ny nz]);
oldvvel=reshape(pickupold(1+(1*s3):2*s3),[nx ny nz]);
oldtheta=reshape(pickupold(1+2*s3:3*s3),[nx ny nz]);
oldsalt=reshape(pickupold(1+3*s3:4*s3),[nx ny nz]);
oldgunm1=reshape(pickupold(1+4*s3:5*s3),[nx ny nz]);
oldgvnm1=reshape(pickupold(1+5*s3:6*s3),[nx ny nz]);
oldetan=reshape(pickupold(1+6*s3:(6*s3)+(1*s2)),[nx ny]);
olddetahdt=reshape(pickupold(1+6*s3+(1*s2):(6*s3)+(2*s2)),[nx ny]);
oldetah=reshape(pickupold(1+6*s3+(2*s2):(6*s3)+(3*s2)),[nx ny]);

stroundold=reshape(stroundold,[nx ny]);
ploadold=reshape(ploadold,[nx ny]);
stpigold=reshape(stpigold,[nx ny]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Find where eta indicates a change in ice shelf topo above threshold.
%Change eta, iceshelf topo and r0surfto reflect change in mesh

splitthreshold=5;
mergethreshold=-5.6;


[splitx,splity] =find(oldetan>splitthreshold & stroundold<0);
[mergex,mergey] =find(oldetan<mergethreshold & stroundold<0);

%[splitx,splity] =find(oldetan>splitthreshold);
%[mergex,mergey] =find(oldetan>mergethreshold);

%Change single, not all cells >threshold
%etax=etax(3);
%etay=etay(3);


    
for i = 1:length(splitx)
    splitz(i)=((stroundold(splitx(i),splity(i)))/10*-1)+1;
end

for i = 1:length(mergex)
    mergez(i)=((stroundold(mergex(i),mergey(i)))/10*-1)+1;
end


OSPLITX=splitx;
OSPLITY=splity;
OSPLITZ=splitz;
NSPLITZ=splitz-1;

OMERGEX=mergex;
OMERGEY=mergey;
OMERGEZ=mergez;
NMERGEZ=mergez+1;

stroundnew=stroundold;
stpignew=stpigold;
stroundnew(splitx,splity)=stroundold(splitx,splity)+delz;
stpignew(splitx,splity)=stpigold(splitx,splity)+delz;
stroundnew(mergex,mergey)=stroundold(mergex,mergey)-delz;
stpignew(mergex,mergey)=stpigold(mergex,mergey)-delz;

newetan=oldetan;
newetan(splitx,splity)=oldetan(splitx,splity)-delz;
newetan(mergex,mergey)=oldetan(mergex,mergey)+delz;

newdetahdt=olddetahdt;
%%%%%%%%%%%%%%%%%%%%%%
%Etah=etan for now
newetah=newetan;

ploadnew=ploadold;
rhoConst = 1030;
talpha = 2e-4;
sbeta  = 7.4e-4;
gravity = 9.81;
k=1;
dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
p = abs(zc)*gravity*rhoConst*1e-4;
dp = p;
kp = 0;
eos = 'jmd95z';
T_sfc = -1.9;
T_bottom = 2;
del_T = (T_bottom - T_sfc)/(59*delz);

for iz = 1:nz;
    tref(iz) = T_sfc + del_T*((iz-30)*delz);
    if iz<=30;
        tref(iz)=-1.9;
    end
    if iz>=90
        tref(iz) =2;
    end
end

S_sfc = 34.2;
S_bottom = 34.7;
del_S = (S_bottom - S_sfc)/(59*delz);

for iz = 1:nz;
    
    
    sref(iz) = S_sfc + del_S*((iz-30)*delz);
    if iz<=30;
        sref(iz)=34.2;
    end
    if iz>=90
        sref(iz) =34.7;
    end
end

t    = tref;

s    = sref;


while rms(dp) > 1e-13
    phiHydF(k) = 0;
    p0 = p;
    kp = kp+1;
    for k = 1:nz
        switch eos
            case 'linear'
                drho = rhoConst*(1-talpha*(t(k)-tref(k))+sbeta*(s(k)-sref(k)))-rhoConst;
            case 'jmd95z'
                drho = densjmd95(s(k),t(k),p(k))-rhoConst;
            case 'mdjwf'
                drho = densmdjwf(s(k),t(k),p(k))-rhoConst;
            otherwise
                error(sprintf('unknown EOS: %s',eos))
        end
        phiHydC(k)   = phiHydF(k) + dzm(k)*gravity*drho/rhoConst;
        phiHydF(k+1) = phiHydC(k) + dzp(k)*gravity*drho/rhoConst;
    end
    switch eos
        case 'mdjwf'
            p = (gravity*rhoConst*abs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
    end
    dp = p-p0;
end


for ix=1:nx
    for iy=1:ny
        k=max(find(abs(zg)<abs(stpignew(ix,iy))));
        if isempty(k)
            k=0;
        end
        if k>0
            
            dr = -zg(k) - stpignew(ix,iy);
            
            if (dr>=delz/2)
                ploadnew(ix,iy) = phiHydF(k) + (delz-dr) * (phiHydC(k)-phiHydF(k))/(delz/2);
            else
                ploadnew(ix,iy) = phiHydC(k) + (delz/2-dr) * (phiHydF(k+1)-phiHydC(k))/(delz/2);
            end
            
        end
    end
end




newsalt=oldsalt;
newtheta=oldtheta;
newuvel=olduvel;
newvvel=oldvvel;
newgunm1=oldgunm1;
newgvnm1=oldgvnm1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Split cells, new cells take values of the cell they are split from

for i=1:length(OSPLITX)
    newsalt(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=oldsalt(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
    newtheta(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=oldtheta(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
    
    newuvel(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=olduvel(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
    %newuvel(OLDX(i)+1,OLDY(i),NEWZ(i))=olduvel(OLDX(i)+1,OLDY(i),OLDZ(i));
    newvvel(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=oldvvel(OSPLITX(i),OSPLITY(i),OSPLITZ(i));
    newvvel(OSPLITX(i),OSPLITY(i)+1,NSPLITZ(i))=oldvvel(OSPLITX(i),OSPLITY(i)+1,OSPLITZ(i));
    
    newgunm1(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=0;
    %newgunm1(OLDX(i)+1,OLDY(i),NEWZ(i))=0;
    newgvnm1(OSPLITX(i),OSPLITY(i),NSPLITZ(i))=0;
    newgvnm1(OSPLITX(i),OSPLITY(i)+1,NSPLITZ(i))=0;
    
end

for i=1:length(OMERGEX)
    newsalt(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((oldsalt(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*oldsalt(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    newtheta(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((oldtheta(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*oldtheta(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    
     newsalt(OMERGEX(i),OMERGEY(i),OMERGEZ(i))=0;
    newtheta(OMERGEX(i),OMERGEY(i),OMERGEZ(i))=0;
    
    newuvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((olduvel(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*olduvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    %newuvel(OMERGEX(i)+1,OMERGEY(i),NMERGEZ(i))=((olduvel(OMERGEX(i)+1,OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i)+1,OMERGEY(i))))...
    %    +(dz*olduvel(OMERGEX(i)+1,OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i)+1,OMERGEY(i)));
    newvvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=((oldvvel(OMERGEX(i),OMERGEY(i),OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*oldvvel(OMERGEX(i),OMERGEY(i),NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    newvvel(OMERGEX(i),OMERGEY(i)+1,NMERGEZ(i))=((oldvvel(OMERGEX(i),OMERGEY(i)+1,OMERGEZ(i))*(dz+oldetan(OMERGEX(i),OMERGEY(i))))...
        +(dz*oldvvel(OMERGEX(i),OMERGEY(i)+1,NMERGEZ(i))))/(2*dz+oldetan(OMERGEX(i),OMERGEY(i)));
    
    newgunm1(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=0;
    %newgunm1(OMERGEX(i)+1,OMERGEY(i),NMERGEZ(i))=0;
    newgvnm1(OMERGEX(i),OMERGEY(i),NMERGEZ(i))=0;
    newgvnm1(OMERGEX(i),OMERGEY(i)+1,NMERGEZ(i))=0;
    
end



% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %Write out pickup file
uvel=reshape(newuvel,[1 nx*ny*nz]);
vvel=reshape(newvvel,[1 nx*ny*nz]);
theta=reshape(newtheta,[1 nx*ny*nz]);
salt=reshape(newsalt,[1 nx*ny*nz]);
gunm1=reshape(newgunm1,[1 nx*ny*nz]);
gvnm1=reshape(newgvnm1,[1 nx*ny*nz]);
etan=reshape(newetan,[1 nx*ny]);
detahdt=reshape(newdetahdt,[1 nx*ny]);
etah=reshape(newetah,[1 nx*ny]);

pickupnew=[uvel,vvel,theta,salt,gunm1,gvnm1,etan,detahdt,etah];
pickupnew=pickupnew';

stroundnew=reshape(stroundnew,[nx*ny 1]);
ploadnew=reshape(ploadnew,[nx*ny 1]);


fid = fopen('pickup.0000004320.data','w','b'); fwrite(fid,pickupnew,'float64'); fclose(fid);
fid = fopen('pload.pig.jmd95z','w','b'); fwrite(fid,ploadnew,'real*8'); fclose(fid)
fid = fopen('shelftopo.round.bin','w','b'); fwrite(fid,stroundnew,'real*8'); fclose(fid)