clear all
close all
clc

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

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

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

fid = fopen('shelftopo.OLDROUND.bin','r','b');
shelftopoold=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]);

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

threshold=5;
[etax,etay] =find(oldetan>threshold);

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

for i = 1:size(etax,1)
    etaz(i)=((shelftopoold(etax(i),etay(i)))/10*-1)+1;
end
     
OLDX=etax;
OLDY=etay;
OLDZ=etaz;
NEWZ=etaz-1;

shelftoponew=shelftopoold;
pignew=pigold;
shelftoponew(etax,etay)=shelftopoold(etax,etay)+delz;
pignew(etax,etay)=pigold(etax,etay)+delz;

newetan=oldetan;
newetan(etax,etay)=oldetan(etax,etay)-delz;

newdetahdt=olddetahdt;
%newetah=oldetah;
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.2050;
S_bottom = 34.6967;
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(pignew(ix,iy))));
    if isempty(k)
    k=0;
    end
    if k>0
    
      dr = -zg(k) - pignew(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;



for i=1:size(OLDX)
    newsalt(OLDX(i),OLDY(i),NEWZ(i))=oldsalt(OLDX(i),OLDY(i),OLDZ(i));
    newtheta(OLDX(i),OLDY(i),NEWZ(i))=oldtheta(OLDX(i),OLDY(i),OLDZ(i));
  
    newuvel(OLDX(i),OLDY(i),NEWZ(i))=olduvel(OLDX(i),OLDY(i),OLDZ(i));
    %newuvel(OLDX(i)+1,OLDY(i),NEWZ(i))=olduvel(OLDX(i)+1,OLDY(i),OLDZ(i));
    newvvel(OLDX(i),OLDY(i),NEWZ(i))=oldvvel(OLDX(i),OLDY(i),OLDZ(i));    
    newvvel(OLDX(i),OLDY(i)+1,NEWZ(i))=oldvvel(OLDX(i),OLDY(i)+1,OLDZ(i)); 
    
    newgunm1(OLDX(i),OLDY(i),NEWZ(i))=0;
    %newgunm1(OLDX(i)+1,OLDY(i),NEWZ(i))=0;
    newgvnm1(OLDX(i),OLDY(i),NEWZ(i))=0;
    newgvnm1(OLDX(i),OLDY(i)+1,NEWZ(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';
  
   shelftoponew=reshape(shelftoponew,[nx*ny 1]);
   ploadnew=reshape(ploadnew,[nx*ny 1]);


 fid = fopen('pickup.0000025920.data','w','b'); fwrite(fid,pickupnew,'real*8'); 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,shelftoponew,'real*8'); fclose(fid)



%