clear all

new = 'input.180x120x23_shelfice';
input = 'input.180x120x23';
eostype = 'mdjwf';

nx = 180;
ny = 120;
nz = 23;
nt = 12;

load FMT
hn = mit_readfield(fullfile(input,'bathymetry.bin'),[nx ny],fmt);
hnz = mit_readfield(fullfile(input,'shelfice_bath.bin'),[nx ny],fmt);
zn = mit_readfield(fullfile(input,'shelfice_topo.bin'),[nx ny],fmt);

h = hn+hnz;
mit_writefield(fullfile(new,'bathymetry.bin'),hn,fmt);
mit_writefield(fullfile(new,'bathymetry.shice'),h,fmt);
mit_writefield(fullfile(new,'shelfice_topo.bin'),zn,fmt);

% create hydrographic fields
levt = mit_readfield(fullfile(input,'lev_t.bin'),[nx ny nz nt],fmt);
levs = mit_readfield(fullfile(input,'lev_s.bin'),[nx ny nz nt],fmt);
is = find(zn~=0);
[ix,iy] = find(zn~=0);
[t,s] = shelfice_hydrography(ix,iy,is,levt,levs);
mit_writefield(fullfile(new,'lev_t.shice'),t,fmt);
mit_writefield(fullfile(new,'lev_s.shice'),s,fmt);

% create geopotential anomaly
gravity = 9.81;
rho0 = 1035;
tol = 0;
si2dbar = 1e-4;
phiHydC = zeros(nz,length(ix));
phiHydF = zeros(nz+1,length(ix));
disp('compute geopotential anomaly')
load VGRID
for ks=1:length(ix)
  t0 = squeeze(mean(t(ix(ks),iy(ks),:,:),4)); 
  s0 = squeeze(mean(s(ix(ks),iy(ks),:,:),4));
  % compute potential anomaly exactly as in code
  % for that we need the correct density
  rho = [];
  p   = -zc(:)*gravity*rho0*si2dbar;
  dp  = p;
  tol1 = 1;
  tol2 = 2;
  kw = 0;
  while  tol1 > tol 
    kw = kw+1;
    if strcmp(eostype,'mdjwf')
      rho = [rho densmdjwf(s0,t0,p(:,end))];
    else
      error(['unknown eostype: ' eostype]);
    end
    p   = [p -zc(:).*rho(:,end)*gravity*si2dbar];
    dp  = p(:,end)-p(:,end-1);
    tol2 = tol1;
    tol1 = sqrt(sum(dp.^2));
    if tol1==tol2; break; end;
  end
  % now intergrate
  drho = rho(:,end)-rho0;
  for k=1:nz
    drm = .5*dz(k);
    if k==1; drm = zf(k)-zc(k); end
    if k==nz; 
      drp = zc(k)-zf(k+1); 
    else
      drp = .5*dz(k+1);
    end
    phiHydC(k,ks)=phiHydF(k,ks) + drm*gravity*drho(k)/rho0;
    phiHydF(k+1,ks)=phiHydC(k,ks) +drp*gravity*drho(k)/rho0;
  end
  % find the appropriate level
  zloc = zn(is(ks));
  kl = max(find(zloc < zf));
  ph(ks) = phiHydF(kl,ks);
end

pload = zeros(nx,ny);
for ks=1:length(ix)
  pload(ix(ks),iy(ks)) = -ph(ks)*rho0;
end

mit_writefield(fullfile(new,['pload.' eostype]),pload,fmt);