1 |
dgoldberg |
1.1 |
%Verion of gendata.m modified by Vero |
2 |
|
|
%This is a matlab script that generates the input data |
3 |
|
|
|
4 |
|
|
|
5 |
|
|
% the configuation approximately the ISOMIP experiment no. 1 |
6 |
|
|
% require matlab functions for equation of state |
7 |
|
|
|
8 |
|
|
|
9 |
|
|
% Dimensions of grid |
10 |
|
|
nx=30; |
11 |
|
|
ny=50; |
12 |
|
|
nz=100; |
13 |
|
|
delZ = 10; |
14 |
|
|
delz=10; |
15 |
|
|
dz=10; |
16 |
|
|
|
17 |
|
|
hfacMin = 0.2; |
18 |
|
|
|
19 |
|
|
dlat = 0.125/4; dy=dlat; |
20 |
|
|
dlon = 0.125; dx=dlon; |
21 |
|
|
|
22 |
|
|
%eos = 'linear'; |
23 |
|
|
eos = 'jmd95z'; |
24 |
|
|
% eos = 'mdjwf'; |
25 |
|
|
|
26 |
|
|
acc = 'real*8'; |
27 |
|
|
|
28 |
|
|
long = [-105.5:dlon:-99.375]; |
29 |
|
|
lonc = long+dlon/2; |
30 |
|
|
latg = [-75.4457:dlat:-74.5068-dlat]; |
31 |
|
|
latc = latg+dlat/2; |
32 |
|
|
size(latc) |
33 |
|
|
|
34 |
|
|
% Nominal depth of model (meters) |
35 |
|
|
H = -1000; %water depth in the ice shelf cavity |
36 |
|
|
Hmin = -900; % deepest point of cavern |
37 |
|
|
Hmax = -100; % shallowest point of cavern |
38 |
|
|
dHdy = (Hmax-Hmin)/(max(latg)-min(latc)); %Slope of ice shelf: if denominator = nx, shelf will cover the whole domain |
39 |
|
|
|
40 |
|
|
bathy = ones(nx,ny)*H; %For flat bathymetry: bathy = ones(nx,ny)*H; |
41 |
|
|
bathy(1,:) = 0; |
42 |
|
|
bathy(nx,:) = 0; |
43 |
|
|
bathy(:,1) = 0; |
44 |
|
|
%bathy(:,ny) = 0; |
45 |
|
|
fid=fopen('bathymetry.pig.bin','w','b'); fwrite(fid,bathy,acc);fclose(fid); |
46 |
|
|
|
47 |
|
|
|
48 |
|
|
|
49 |
|
|
dz = delZ*ones(1,nz); |
50 |
|
|
zgp1 = [0,cumsum(dz)]; |
51 |
|
|
zc = .5*(zgp1(1:end-1)+zgp1(2:end)); |
52 |
|
|
zg = zgp1(1:end-1); |
53 |
|
|
dz = diff(zgp1); |
54 |
|
|
sprintf('delZ = %d * %7.6g,',nz,dz) |
55 |
|
|
|
56 |
|
|
|
57 |
|
|
T_sfc = -1.9; |
58 |
|
|
T_bottom = 1.1967; |
59 |
|
|
del_T = (T_bottom - T_sfc)/((nz-1)*delZ); |
60 |
|
|
|
61 |
|
|
for iz = 1:nz; |
62 |
|
|
tref(iz) = T_sfc + del_T*((iz-1)*delZ); |
63 |
|
|
end |
64 |
|
|
|
65 |
|
|
S_sfc = 34.2050; |
66 |
|
|
S_bottom = 34.6967; |
67 |
|
|
del_S = (S_bottom - S_sfc)/((nz-1)*delZ); |
68 |
|
|
|
69 |
|
|
for iz = 1:nz; |
70 |
|
|
sref(iz) = S_sfc + del_S*((iz-1)*delZ); |
71 |
|
|
end |
72 |
|
|
|
73 |
|
|
% Gravity |
74 |
|
|
gravity=9.81; |
75 |
|
|
rhoConst = 1030; |
76 |
|
|
% compute potential field underneath ice shelf |
77 |
|
|
talpha = 2e-4; |
78 |
|
|
sbeta = 7.4e-4; |
79 |
|
|
% tref = -1.9*ones(nz,1); |
80 |
|
|
t = tref; |
81 |
|
|
% sref = 34.4*ones(nz,1); |
82 |
|
|
s = sref; |
83 |
|
|
gravity = 9.81; |
84 |
|
|
k=1; |
85 |
|
|
dzm = abs([zg(1)-zc(1) .5*diff(zc)]); |
86 |
|
|
dzp = abs([.5*diff(zc) zc(end)-zg(end)]); |
87 |
|
|
p = abs(zc)*gravity*rhoConst*1e-4; |
88 |
|
|
dp = p; |
89 |
|
|
kp = 0; |
90 |
|
|
|
91 |
|
|
|
92 |
|
|
|
93 |
|
|
while rms(dp) > 1e-13 |
94 |
|
|
phiHydF(k) = 0; |
95 |
|
|
p0 = p; |
96 |
|
|
kp = kp+1 |
97 |
|
|
for k = 1:nz |
98 |
|
|
switch eos |
99 |
|
|
case 'linear' |
100 |
|
|
drho = rhoConst*(1-talpha*(t(k)-tref(k))+sbeta*(s(k)-sref(k)))-rhoConst; |
101 |
|
|
case 'jmd95z' |
102 |
|
|
drho = densjmd95(s(k),t(k),p(k))-rhoConst; |
103 |
|
|
case 'mdjwf' |
104 |
|
|
drho = densmdjwf(s(k),t(k),p(k))-rhoConst; |
105 |
|
|
otherwise |
106 |
|
|
error(sprintf('unknown EOS: %s',eos)) |
107 |
|
|
end |
108 |
|
|
phiHydC(k) = phiHydF(k) + dzm(k)*gravity*drho/rhoConst; |
109 |
|
|
phiHydF(k+1) = phiHydC(k) + dzp(k)*gravity*drho/rhoConst; |
110 |
|
|
end |
111 |
|
|
switch eos |
112 |
|
|
case 'mdjwf' |
113 |
|
|
p = (gravity*rhoConst*abs(zc) + phiHydC*rhoConst)/gravity/rhoConst; |
114 |
|
|
end |
115 |
|
|
dp = p-p0; |
116 |
|
|
end |
117 |
|
|
|
118 |
|
|
|
119 |
|
|
%Modify icetopo (shape of ice shelf cavity) |
120 |
|
|
|
121 |
|
|
% %icetopo = ones(ny,1)*min(Hmin + 2*dHdy*(lonc(nx)-long),Hmax); |
122 |
|
|
% B=min(Hmin + 2*dHdy*(lonc(ny)-long),Hmax); |
123 |
|
|
% B=fliplr(B); |
124 |
|
|
% |
125 |
|
|
% icetopo = ones(nx,1)*B; |
126 |
|
|
% %icetopo = icetopo'; |
127 |
|
|
% icetopo(:,101:end)=0; |
128 |
|
|
|
129 |
|
|
% use streamice generated thickness |
130 |
|
|
fid = fopen('h_init.box','r','b'); hinit=fread(fid,[dy dx],'real*8'); fclose(fid); |
131 |
|
|
%hinit(:,1:5)=550; |
132 |
|
|
icetopo = -917/1028*hinit; |
133 |
|
|
|
134 |
|
|
|
135 |
|
|
%Modify icetopo (shape of ice shelf cavity) |
136 |
|
|
|
137 |
|
|
icetopo = ones(ny,1)*min(Hmax,Hmin + 2*dHdy*(lonc(nx)-long)); icetopo = icetopo'; |
138 |
|
|
icetopo(1:25,:)=0; |
139 |
|
|
|
140 |
|
|
|
141 |
|
|
% use streamice generated thickness |
142 |
|
|
fid = fopen('h_init.box','r','b'); hinit=fread(fid,[30 50],'real*8'); fclose(fid); |
143 |
|
|
%hinit(:,1:5)=550; |
144 |
|
|
icetopo = -917/1028*hinit; |
145 |
|
|
|
146 |
|
|
|
147 |
|
|
% adjust topo so that no hfac is smaller than hfacMin |
148 |
|
|
|
149 |
|
|
for ix=1:nx |
150 |
|
|
for iy=1:ny |
151 |
|
|
k=max(find(abs(zg)<abs(icetopo(ix,iy)))); |
152 |
|
|
|
153 |
|
|
if(~isempty(k)) |
154 |
|
|
|
155 |
|
|
|
156 |
|
|
hfacTemp = (icetopo(ix,iy) - (-zg(k+1)))/delZ; |
157 |
|
|
|
158 |
|
|
if (hfacTemp < hfacMin) |
159 |
|
|
if (hfacTemp < hfacMin/2) |
160 |
|
|
hfacTemp = 0; |
161 |
|
|
else |
162 |
|
|
hfacTemp = hfacMin; |
163 |
|
|
end |
164 |
|
|
end |
165 |
|
|
|
166 |
|
|
else |
167 |
|
|
|
168 |
|
|
hfacTemp = 0; |
169 |
|
|
|
170 |
|
|
end |
171 |
|
|
|
172 |
|
|
icetopo(ix,iy) = icetopo(ix,iy) + hfacTemp; |
173 |
|
|
|
174 |
|
|
end |
175 |
|
|
end |
176 |
|
|
|
177 |
|
|
% phi anomaly (relative to hydrostatic with rho_const) at icetopo |
178 |
|
|
|
179 |
|
|
% phi anomaly (relative to hydrostatic with rho_const) at icetopo |
180 |
|
|
|
181 |
|
|
phi0surf = zeros(nx,ny); |
182 |
|
|
|
183 |
|
|
for ix=1:nx |
184 |
|
|
for iy=1:ny |
185 |
|
|
k=max(find(abs(zg)<abs(icetopo(ix,iy)))); |
186 |
|
|
if isempty(k) |
187 |
|
|
k=0; |
188 |
|
|
end |
189 |
|
|
if k>0 |
190 |
|
|
|
191 |
|
|
dr = -zg(k) - icetopo(ix,iy); |
192 |
|
|
|
193 |
|
|
if (dr>=delZ/2) |
194 |
|
|
phi0surf(ix,iy) = phiHydF(k) + (delZ-dr) * (phiHydC(k)-phiHydF(k))/(delZ/2); |
195 |
|
|
else |
196 |
|
|
phi0surf(ix,iy) = phiHydC(k) + (delZ/2-dr) * (phiHydF(k+1)-phiHydC(k))/(delZ/2); |
197 |
|
|
end |
198 |
|
|
|
199 |
|
|
end |
200 |
|
|
end |
201 |
|
|
end |
202 |
|
|
|
203 |
|
|
mass = phi0surf / gravity - rhoConst * icetopo; |
204 |
|
|
|
205 |
|
|
%use streamicegenerated thickness |
206 |
|
|
|
207 |
|
|
%mass = hinit * 917; |
208 |
|
|
|
209 |
|
|
fid = fopen('shelftopo.pig.bin','w','b'); fwrite(fid,icetopo,'real*8'); fclose(fid); |
210 |
|
|
fid = fopen('shelficemassinit.bin','w','b'); fwrite(fid,mass,'real*8'); fclose(fid); |
211 |
|
|
fid = fopen('pload.pig.jmd95z','w','b'); fwrite(fid,phi0surf,'real*8'); fclose(fid); |
212 |
|
|
|
213 |
|
|
etainit = zeros(size(phi0surf)); |
214 |
|
|
|
215 |
|
|
% new topography: icetopo rounded to the nearest k * deltaZ |
216 |
|
|
% eta_init set to make difference |
217 |
|
|
|
218 |
|
|
icetopo2 = icetopo; |
219 |
|
|
|
220 |
|
|
for ix=1:nx |
221 |
|
|
for iy=1:ny |
222 |
|
|
k=max(find(abs(zg)<abs(icetopo2(ix,iy)))); |
223 |
|
|
if isempty(k) |
224 |
|
|
k=0; |
225 |
|
|
else |
226 |
|
|
|
227 |
|
|
dr = 1-(-zg(k) - icetopo2(ix,iy))/delz; |
228 |
|
|
if (dr > .5) |
229 |
|
|
% bring Ro_surf *up* to closest grid face & make etainit negative |
230 |
|
|
% to compensate |
231 |
|
|
icetopo2(ix,iy) = -zg(k); |
232 |
|
|
etainit(ix,iy) = (dr-1)*delz; |
233 |
|
|
else |
234 |
|
|
% bring Ro_surf *down* to closest grid face & make etainit pos |
235 |
|
|
% to compensate |
236 |
|
|
icetopo2(ix,iy) = -zg(k+1); |
237 |
|
|
etainit(ix,iy) = (dr)*delz; |
238 |
|
|
end |
239 |
|
|
|
240 |
|
|
end |
241 |
|
|
end |
242 |
|
|
end |
243 |
|
|
|
244 |
|
|
fid = fopen('shelftopo.round.bin','w','b'); fwrite(fid,icetopo2,'real*8'); fclose(fid); |
245 |
|
|
fid = fopen('etainit.round.bin','w','b'); fwrite(fid,etainit,'real*8'); fclose(fid); |