CLEAN/input/gendata_vero.m

%Verion of gendata.m modified by Vero
%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=3; 
ny=200;
nz=100;
delz = 10;

hfacMin = 0.2;

dlat = 0.125/16; dy=dlat;
dlon = 0.125/4; dx=dlon;

%eos = 'linear';
eos = 'jmd95z';
% eos = 'mdjwf';

acc = 'real*8';

long = [-105.5:dlon:-105.5];
lonc = long+dlon/2;
latg = [-75.4457:dlat:-73.8809-dlat];
latc = latg+dlat/2;
size(latc)

% Nominal depth of model (meters)
H = -1000;                                      %water depth in the ice shelf cavity
Hmin = -900; % deepest point of cavern          
Hmax = -300; % shallowest point of cavern
dHdy = (Hmax-Hmin)/(max(latc)-min(latc));                               %Slope of ice shelf: if denominator = nx, shelf will cover the whole domain

bathy = ones(nx,ny)*H;                          %For flat bathymetry: bathy = ones(nx,ny)*H;
%bathy(1,:) = 0;
%bathy(2,:) = 0;
%bathy(nx,:) = 0; 
%bathy(nx-1,:) = 0; 
bathy(:,1) = 0;
%bathy(:,ny) = 0;
fid=fopen('bathymetry.pig.bin','w','b'); fwrite(fid,bathy,acc);fclose(fid);


dz = delz*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)


T_sfc = -1.9;
T_bottom = 2;
del_T = (T_bottom - T_sfc)/(59*delz);

for iz = 1:nz;
    
    
    tref(iz) = T_sfc + del_T*((iz-30)*delz); 
    if iz<=30;
        tref(iz)=-1.9;
    end
    if iz>=90
        tref(iz) =2;
    end
end

S_sfc = 34.2;
S_bottom = 34.7;
del_S = (S_bottom - S_sfc)/(59*delz);

for iz = 1:nz;
    
    
    sref(iz) = S_sfc + del_S*((iz-30)*delz); 
    if iz<=30;
        sref(iz)=34.2;
    end
    if iz>=90
        sref(iz) =34.7;
    end
end

% 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.4*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 rms(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


%Modify icetopo (shape of ice shelf cavity)

%icetopo = ones(ny,1)*min(Hmin + 2*dHdy*(lonc(nx)-long),Hmax); 
B=min(Hmin + 2*dHdy*(latc(ny)-latg),Hmax);
B=fliplr(B);

icetopo = ones(nx,1)*B; 
%icetopo = icetopo';
icetopo(:,101:end)=0;

% use streamice generated thickness
 fid = fopen('h_init.box','r','b'); hinit=fread(fid,[3 200],'real*8'); fclose(fid);
%hinit(:,1:5)=550;
% load hinit
icetopo = -917/1028*hinit;


% adjust topo so that no hfac is smaller than hfacMin

for ix=1:nx
  for iy=1:ny
    k=max(find(abs(zg)<abs(icetopo(ix,iy))));
    
    if(~isempty(k))


    hfacTemp = (icetopo(ix,iy) - (-zg(k+1)))/delz;
    
    if (hfacTemp < hfacMin)
        if (hfacTemp < hfacMin/2)
            hfacTemp = 0;
        else
            hfacTemp = hfacMin;
        end
    end
    
    else
        
        hfacTemp = 0;
        
    end
    
    icetopo(ix,iy) = icetopo(ix,iy) + hfacTemp;
    
  end
end

% phi anomaly (relative to hydrostatic with rho_const) at icetopo

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
      
      dr = -zg(k) - icetopo(ix,iy);
      
      if (dr>=delz/2)
          phi0surf(ix,iy) = phiHydF(k) + (delz-dr) * (phiHydC(k)-phiHydF(k))/(delz/2);
      else
          phi0surf(ix,iy) = phiHydC(k) + (delz/2-dr) * (phiHydF(k+1)-phiHydC(k))/(delz/2);
      end
      
    end
  end
end

mass = phi0surf / gravity - rhoConst * icetopo;

%use streamicegenerated thickness

%mass = hinit * 917;
icetopo(:,1)=-750;
fid = fopen('shelftopo.pig.bin','w','b'); fwrite(fid,icetopo,'real*8'); fclose(fid);
fid = fopen('shelficemassinit.bin','w','b'); fwrite(fid,mass,'real*8'); fclose(fid);
fid = fopen('pload.pig.jmd95z','w','b'); fwrite(fid,phi0surf,'real*8'); fclose(fid);

etainit = zeros(size(phi0surf));

% new topography: icetopo rounded to the nearest k * deltaZ
%                 eta_init set to make difference

icetopo2 = icetopo;

for ix=1:nx
  for iy=1:ny
    k=max(find(abs(zg)<abs(icetopo2(ix,iy))));
    if isempty(k)
      k=0;
    else
      
      dr = 1-(-zg(k) - icetopo2(ix,iy))/delz;
      if (dr > .4)
          % bring Ro_surf *up* to closest grid face & make etainit negative
          % to compensate
          icetopo2(ix,iy) = -zg(k);
          etainit(ix,iy) = (dr-1)*delz;
      else
          % bring Ro_surf *down* to closest grid face & make etainit pos
          % to compensate
          icetopo2(ix,iy) = -zg(k+1);
          etainit(ix,iy) = (dr)*delz;
      end
       
    end
  end
end


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