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