--- MITgcm_contrib/mlosch/interp_llc/addshelfice.m	2007/05/03 21:07:20	1.1
+++ MITgcm_contrib/mlosch/interp_llc/addshelfice.m	2007/05/04 14:53:00	1.2
@@ -1,23 +1,25 @@
 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 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,8 +27,17 @@
 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;
@@ -80,10 +91,10 @@
   ph(ks) = phiHydF(kl,ks);
 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);