--- MITgcm_contrib/mlosch/interp_llc/addshelfice.m	2007/05/03 21:07:20	1.1
+++ MITgcm_contrib/mlosch/interp_llc/addshelfice.m	2007/05/16 16:09:06	1.3
@@ -1,23 +1,27 @@
-clear all
+%clear all
 
-new = 'input.180x120x23_shelfice';
-input = 'input.180x120x23';
+new = 'input.shelfice';
+input = 'input';
 eostype = 'mdjwf';
 
-nx = 180;
-ny = 120;
+nx = 45;
+ny = nx*18;
 nz = 23;
 nt = 12;
 
+load MASKS
+hf = msk;
 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);
+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,'bathymetry.bin'),hn,fmt);
-mit_writefield(fullfile(new,'bathymetry.shice'),h,fmt);
-mit_writefield(fullfile(new,'shelfice_topo.bin'),zn,fmt);
+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);
@@ -25,65 +29,77 @@
 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);
+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;
-tol = 0;
+tol0 = 0;
 si2dbar = 1e-4;
-phiHydC = zeros(nz,length(ix));
-phiHydF = zeros(nz+1,length(ix));
 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   = -zc(:)*gravity*rho0*si2dbar;
+  p   = abs(zc(:))*gravity*rho0*si2dbar;
   dp  = p;
   tol1 = 1;
   tol2 = 2;
-  kw = 0;
-  while  tol1 > tol 
-    kw = kw+1;
+  kp = 0;
+  while  tol1 > tol0 
+    kp = kp+1;
+    p0 = p;
     if strcmp(eostype,'mdjwf')
-      rho = [rho densmdjwf(s0,t0,p(:,end))];
+      drho =  densmdjwf(s0,t0,p(:,end))-rho0;
     else
       error(['unknown eostype: ' eostype]);
     end
-    p   = [p -zc(:).*rho(:,end)*gravity*si2dbar];
+    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
-  % 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);
+  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 = zeros(nx,ny);
+pload = 0*hn;
 for ks=1:length(ix)
   pload(ix(ks),iy(ks)) = -ph(ks)*rho0;
 end
 
-mit_writefield(fullfile(new,['pload.' eostype]),pload,fmt);
+mit_writefield(fullfile(new,['pload.' eostype]),mdsiocompact(pload),fmt);