3D_example/input/gendata.m

% This is a matlab script that generates the input data

% the configuation approximately the ISOMIP experiment no. 1
% require matlab functions for equation of state

% Dimensions of grid
nx=200; 
ny=150;
nz=500;
deltaZ = 1;

eos = 'jmd95z';

acc = 'real*8';

% Nominal depth of model (meters)
H = -500;

  bathy = ones(nx,ny)*H;
  bathy(end,:) = -499;
  bathy(:,1)=0;
  bathy(:,end)=0;
% fid=fopen('bathy.box','w','b'); fwrite(fid,bathy,acc);fclose(fid);

dz = deltaZ*ones(1,nz);
zgp1 = [0,cumsum(dz)];
zc = .5*(zgp1(1:end-1)+zgp1(2:end));
zg = zgp1(1:end-1);
dz = diff(zgp1);
sprintf('delZ = %d * %7.6g,',nz,dz)

% Gravity
gravity=9.81;
rhoConst = 1030;
% compute potential field underneath ice shelf
talpha = 2e-4;
sbeta  = 7.4e-4;
tref = -1.9*ones(nz,1);
t    = tref;
sref = 34*ones(nz,1);
s    = sref;
gravity = 9.81;
k=1;
dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
p = abs(zc)*gravity*rhoConst*1e-4;
dp = p;
kp = 0;
while std(dp) > 1e-13
  phiHydF(k) = 0;
  p0 = p;
  kp = kp+1
  for k = 1:nz
    switch eos
     case 'linear'
      drho = rhoConst*(1-talpha*(t(k)-tref(k))+sbeta*(s(k)-sref(k)))-rhoConst;
     case 'jmd95z'
      drho = densjmd95(s(k),t(k),p(k))-rhoConst;
     case 'mdjwf'
      drho = densmdjwf(s(k),t(k),p(k))-rhoConst;
     otherwise
      error(sprintf('unknown EOS: %s',eos))
    end
    phiHydC(k)   = phiHydF(k) + dzm(k)*gravity*drho/rhoConst;
    phiHydF(k+1) = phiHydC(k) + dzp(k)*gravity*drho/rhoConst;
  end
  switch eos
   case 'mdjwf'
    p = (gravity*rhoConst*abs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
  end
  dp = p-p0;
end

ceil(ny/2)
icetopo = zeros(nx,ny);
icetopo(end,:)= -499;
hole1=[75:76];
icetopo(end,hole1) = -498;
fid=fopen('icetopo.exp1','w','b'); fwrite(fid,icetopo,acc);fclose(fid);
fid=fopen('pload.exp1','w','b'); fwrite(fid,-icetopo,acc);fclose(fid);
phi0surf = zeros(nx,ny);
for ix=1:nx
  for iy=1:ny
    k=max(find(abs(zg)<abs(icetopo(ix,iy))));
    if isempty(k)
      k=0;
    end
    if k>0
      phi0surf(ix,iy) = phiHydF(k)*rhoConst;
    end
  end
end
fid=fopen(['phi0surf.exp1.' eos],'w','b'); fwrite(fid,phi0surf,acc);fclose(fid);

1	yunx	1.1	% This is a matlab script that generates the input data
2
3			% the configuation approximately the ISOMIP experiment no. 1
4			% require matlab functions for equation of state
5
6			% Dimensions of grid
7			nx=200;
8			ny=150;
9			nz=500;
10			deltaZ = 1;
11
12			eos = 'jmd95z';
13
14			acc = 'real*8';
15
16			% Nominal depth of model (meters)
17			H = -500;
18
19			bathy = ones(nx,ny)*H;
20			bathy(end,:) = -499;
21			bathy(:,1)=0;
22			bathy(:,end)=0;
23			% fid=fopen('bathy.box','w','b'); fwrite(fid,bathy,acc);fclose(fid);
24
25			dz = deltaZ*ones(1,nz);
26			zgp1 = [0,cumsum(dz)];
27			zc = .5*(zgp1(1:end-1)+zgp1(2:end));
28			zg = zgp1(1:end-1);
29			dz = diff(zgp1);
30			sprintf('delZ = %d * %7.6g,',nz,dz)
31
32			% Gravity
33			gravity=9.81;
34			rhoConst = 1030;
35			% compute potential field underneath ice shelf
36			talpha = 2e-4;
37			sbeta = 7.4e-4;
38			tref = -1.9*ones(nz,1);
39			t = tref;
40			sref = 34*ones(nz,1);
41			s = sref;
42			gravity = 9.81;
43			k=1;
44			dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
45			dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
46			p = abs(zc)gravityrhoConst*1e-4;
47			dp = p;
48			kp = 0;
49			while std(dp) > 1e-13
50			phiHydF(k) = 0;
51			p0 = p;
52			kp = kp+1
53			for k = 1:nz
54			switch eos
55			case 'linear'
56			drho = rhoConst(1-talpha(t(k)-tref(k))+sbeta*(s(k)-sref(k)))-rhoConst;
57			case 'jmd95z'
58			drho = densjmd95(s(k),t(k),p(k))-rhoConst;
59			case 'mdjwf'
60			drho = densmdjwf(s(k),t(k),p(k))-rhoConst;
61			otherwise
62			error(sprintf('unknown EOS: %s',eos))
63			end
64			phiHydC(k) = phiHydF(k) + dzm(k)gravitydrho/rhoConst;
65			phiHydF(k+1) = phiHydC(k) + dzp(k)gravitydrho/rhoConst;
66			end
67			switch eos
68			case 'mdjwf'
69			p = (gravityrhoConstabs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
70			end
71			dp = p-p0;
72			end
73
74			ceil(ny/2)
75			icetopo = zeros(nx,ny);
76			icetopo(end,:)= -499;
77			hole1=[75:76];
78			icetopo(end,hole1) = -498;
79			fid=fopen('icetopo.exp1','w','b'); fwrite(fid,icetopo,acc);fclose(fid);
80			fid=fopen('pload.exp1','w','b'); fwrite(fid,-icetopo,acc);fclose(fid);
81			phi0surf = zeros(nx,ny);
82			for ix=1:nx
83			for iy=1:ny
84			k=max(find(abs(zg)<abs(icetopo(ix,iy))));
85			if isempty(k)
86			k=0;
87			end
88			if k>0
89			phi0surf(ix,iy) = phiHydF(k)*rhoConst;
90			end
91			end
92			end
93			fid=fopen(['phi0surf.exp1.' eos],'w','b'); fwrite(fid,phi0surf,acc);fclose(fid);
94