CLEAN/input/gendataSHICE.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=30; 
ny=50;
nz=100;
delZ = 10;
delz=10;
dz=10;

hfacMin = 0.2;

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

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

acc = 'real*8';

long = [-105.5:dlon:-99.375];
lonc = long+dlon/2;
latg = [-75.4457:dlat:-74.5068-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 = -100; % shallowest point of cavern
dHdy = (Hmax-Hmin)/(max(latg)-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(nx,:) = 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 = 1.1967;
del_T = (T_bottom - T_sfc)/((nz-1)*delZ);

for iz = 1:nz;
    tref(iz) = T_sfc + del_T*((iz-1)*delZ);    
end

S_sfc = 34.2050;
S_bottom = 34.6967;
del_S = (S_bottom - S_sfc)/((nz-1)*delZ);

for iz = 1:nz;
    sref(iz) = S_sfc + del_S*((iz-1)*delZ);    
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*(lonc(ny)-long),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,[dy dx],'real*8'); fclose(fid);
%hinit(:,1:5)=550;
icetopo = -917/1028*hinit;


%Modify icetopo (shape of ice shelf cavity)

icetopo = ones(ny,1)*min(Hmax,Hmin + 2*dHdy*(lonc(nx)-long)); icetopo = icetopo';
icetopo(1:25,:)=0;


% use streamice generated thickness
fid = fopen('h_init.box','r','b'); hinit=fread(fid,[30 50],'real*8'); fclose(fid);
%hinit(:,1:5)=550;
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

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

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 > .5)
          % 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

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=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)gravityrhoConst*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)gravitydrho/rhoConst;
109			phiHydF(k+1) = phiHydC(k) + dzp(k)gravitydrho/rhoConst;
110			end
111			switch eos
112			case 'mdjwf'
113			p = (gravityrhoConstabs(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 + 2dHdy*(lonc(nx)-long),Hmax);
122			% B=min(Hmin + 2dHdy(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 + 2dHdy*(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);