%clear all

new = 'input.shelfice';
input = 'input';
eostype = 'mdjwf';

nx = 45;
ny = nx*18;
nz = 23;
nt = 12;

load MASKS
hf = msk;
load FMT
load HN
load ZN
% $$$ hn = mit_readfield(fullfile(input,'bathy_llc_p90.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,'bathy_llc_p90.bin'),mdsiocompact(hn),fmt);
mit_writefield(fullfile(new,'bathy_llc_p90.shice'),mdsiocompact(h),fmt);
mit_writefield(fullfile(new,'shelfice_topo.bin'),mdsiocompact(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'),mdsiocompact(t),fmt);
mit_writefield(fullfile(new,'lev_s.shice'),mdsiocompact(s),fmt);

% create hydrographic fields
levt = mdsiocompact(mit_readfield(fullfile(input,'lev_t.init'),[nx ny nz],fmt),0);
levs = mdsiocompact(mit_readfield(fullfile(input,'lev_s.init'),[nx ny nz],fmt),0);
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.init'),mdsiocompact(t),fmt);
mit_writefield(fullfile(new,'lev_s.shice.init'),mdsiocompact(s),fmt);

% create geopotential anomaly
gravity = 9.81;
rho0 = 1035;
tol0 = 0;
si2dbar = 1e-4;
disp('compute geopotential anomaly')
load VGRID
zg = zf;
dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
hFacMin = 0.1;
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   = abs(zc(:))*gravity*rho0*si2dbar;
  dp  = p;
  tol1 = 1;
  tol2 = 2;
  kp = 0;
  while  tol1 > tol0 
    kp = kp+1;
    p0 = p;
    if strcmp(eostype,'mdjwf')
      drho =  densmdjwf(s0,t0,p(:,end))-rho0;
    else
      error(['unknown eostype: ' eostype]);
    end
    phiHydF(1) = 0;
    for k=1:length(zc(:));
      phiHydC(k)   = phiHydF(k) + dzm(k)*gravity*drho(k)/rho0;
      phiHydF(k+1) = phiHydC(k) + dzp(k)*gravity*drho(k)/rho0;
    end
    p   = [p (gravity*rho0*abs(zc(:)) + phiHydC(:)*rho0)/gravity/rho0];
    dp  = p(:,end)-p(:,end-1);
    tol2 = tol1;
    tol1 = sqrt(sum(dp.^2));
    if tol1==tol2; break; end;
  end
  % find the appropriate level
  zloc = zn(is(ks));
  kl0 = max(find(abs(zg-hFacMin*zg)<=abs(zloc)));
  hfloc= squeeze(hf(ix(ks),iy(ks),:));
  kl = min(find(hfloc>0));
  if isempty(kl);
    kl = 0;
    ph(ks) = 0;
  else
    ph(ks) = phiHydF(kl);
  end
  disp(sprintf('kl0 = %u, kl = %u',kl0,kl));
end

pload = 0*hn;
for ks=1:length(ix)
  pload(ix(ks),iy(ks)) = -ph(ks)*rho0;
end

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