mlosch/interp_llc/addshelfice.m

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);

1	mlosch	1.1	clear all
2
3			new = 'input.180x120x23_shelfice';
4			input = 'input.180x120x23';
5			eostype = 'mdjwf';
6
7			nx = 180;
8			ny = 120;
9			nz = 23;
10			nt = 12;
11
12			load FMT
13			hn = mit_readfield(fullfile(input,'bathymetry.bin'),[nx ny],fmt);
14			hnz = mit_readfield(fullfile(input,'shelfice_bath.bin'),[nx ny],fmt);
15			zn = mit_readfield(fullfile(input,'shelfice_topo.bin'),[nx ny],fmt);
16
17			h = hn+hnz;
18			mit_writefield(fullfile(new,'bathymetry.bin'),hn,fmt);
19			mit_writefield(fullfile(new,'bathymetry.shice'),h,fmt);
20			mit_writefield(fullfile(new,'shelfice_topo.bin'),zn,fmt);
21
22			% create hydrographic fields
23			levt = mit_readfield(fullfile(input,'lev_t.bin'),[nx ny nz nt],fmt);
24			levs = mit_readfield(fullfile(input,'lev_s.bin'),[nx ny nz nt],fmt);
25			is = find(zn~=0);
26			[ix,iy] = find(zn~=0);
27			[t,s] = shelfice_hydrography(ix,iy,is,levt,levs);
28			mit_writefield(fullfile(new,'lev_t.shice'),t,fmt);
29			mit_writefield(fullfile(new,'lev_s.shice'),s,fmt);
30
31			% create geopotential anomaly
32			gravity = 9.81;
33			rho0 = 1035;
34			tol = 0;
35			si2dbar = 1e-4;
36			phiHydC = zeros(nz,length(ix));
37			phiHydF = zeros(nz+1,length(ix));
38			disp('compute geopotential anomaly')
39			load VGRID
40			for ks=1:length(ix)
41			t0 = squeeze(mean(t(ix(ks),iy(ks),:,:),4));
42			s0 = squeeze(mean(s(ix(ks),iy(ks),:,:),4));
43			% compute potential anomaly exactly as in code
44			% for that we need the correct density
45			rho = [];
46			p = -zc(:)gravityrho0*si2dbar;
47			dp = p;
48			tol1 = 1;
49			tol2 = 2;
50			kw = 0;
51			while tol1 > tol
52			kw = kw+1;
53			if strcmp(eostype,'mdjwf')
54			rho = [rho densmdjwf(s0,t0,p(:,end))];
55			else
56			error(['unknown eostype: ' eostype]);
57			end
58			p = [p -zc(:).rho(:,end)gravity*si2dbar];
59			dp = p(:,end)-p(:,end-1);
60			tol2 = tol1;
61			tol1 = sqrt(sum(dp.^2));
62			if tol1==tol2; break; end;
63			end
64			% now intergrate
65			drho = rho(:,end)-rho0;
66			for k=1:nz
67			drm = .5*dz(k);
68			if k==1; drm = zf(k)-zc(k); end
69			if k==nz;
70			drp = zc(k)-zf(k+1);
71			else
72			drp = .5*dz(k+1);
73			end
74			phiHydC(k,ks)=phiHydF(k,ks) + drmgravitydrho(k)/rho0;
75			phiHydF(k+1,ks)=phiHydC(k,ks) +drpgravitydrho(k)/rho0;
76			end
77			% find the appropriate level
78			zloc = zn(is(ks));
79			kl = max(find(zloc < zf));
80			ph(ks) = phiHydF(kl,ks);
81			end
82
83			pload = zeros(nx,ny);
84			for ks=1:length(ix)
85			pload(ix(ks),iy(ks)) = -ph(ks)*rho0;
86			end
87
88			mit_writefield(fullfile(new,['pload.' eostype]),pload,fmt);
89