--- MITgcm_contrib/mlosch/interp_llc/addshelfice.m	2007/05/04 14:53:00	1.2
+++ MITgcm_contrib/mlosch/interp_llc/addshelfice.m	2007/05/16 16:09:06	1.3
@@ -1,4 +1,4 @@
-clear all
+%clear all
 
 new = 'input.shelfice';
 input = 'input';
@@ -9,6 +9,8 @@
 nz = 23;
 nt = 12;
 
+load MASKS
+hf = msk;
 load FMT
 load HN
 load ZN
@@ -42,53 +44,56 @@
 % 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 = 0*hn;