DIG/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=400;
nz=210;
delz = 10;

hfacMin = 0.2;

dlat = 0.125/32; 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('bathysteep.box','r','b');
bathy = fread(fid,inf,'real*8');
fclose(fid);
bathy=reshape(bathy, [nx ny]);
fid=fopen('bathy.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('hinitsteep3.box','r','b'); hinit=fread(fid,[3 400],'real*8'); fclose(fid);
%hinit(:,1:5)=550;
% load hinit
%hinit(:,1)=1740;
icetopo = -917/1028*hinit;

      for i =1:3
       for j=1:400
     if icetopo(i,j)<bathy(i,j)
        icetopo(i,j)=bathy(i,j)+3;

    end 
end
      end

            for i =1:3
       for j=1:400
     if icetopo(i,j)-bathy(i,j)<3;
        icetopo(i,j)=bathy(i,j)+3;

    end 
end
end
% 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
      for i =1:3
       for j=1:400
     if icetopo(i,j)<bathy(i,j)
%        icetopo(i,j)=bathy(i,j)+6;

    end 
end
      end

                  for i =1:3
       for j=1:400
     if icetopo(i,j)-bathy(i,j)<10
%        icetopo(i,j)=bathy(i,j)+delz;

    end 
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)=-1000;
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 > .25)
          % 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

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