shelfice_remeshing/input/gendata.m

%Verion of gendata.m modified by DNG
%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=90;
delz = 10;

hfacMin = 0.2;
%mwct = 3;

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


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)

%%%%%%%%% stratification %%%%%%%%%%%%%%%%

T_sfc = -1.9;
T_bottom = 1.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

%%%%%%%%%%% density %%%%%%%%%%%%%%%%%

% Gravity
gravity=9.81;
rhoConst = 1000;
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;

Rho = zeros(nz,1);

while rms(dp) > 1e-13
  phiHydF(k) = 0;
  p0 = p;
  kp = kp+1;
  for k = 1:nz
    switch eos
     case 'linear'

     case 'jmd95z'
      drho = densjmd95(sref(k),tref(k),p(k))-rhoConst;
     case 'mdjwf'
      drho = densmdjwf(sref(k),tref(k),p(k))-rhoConst;
     otherwise
      error(sprintf('unknown EOS: %s',eos))
    end
    Rho(k) = drho+rhoConst;
    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

shelficemass = binread('h0.bin',3,200) * 917;

%phi0surf = zeros(nx,ny);
topo = zeros(nx,ny);

for ix=1:nx
  for iy=1:ny
%    k=max(find(abs(zg)<abs(topo(ix,iy))));
%    if isempty(k)
%      k=0;
%    end
%    if k>0
      
%      dr = -zg(k) - topo(ix,iy);
      
%      if (dr<delz/2)
%          phi0surf(ix,iy) = phiHydF(k) + (delz/2-dr) * (phiHydC(k)-phiHydF(k))/(delz/2);
%      else
%          phi0surf(ix,iy) = phiHydC(k) + (dr-delz/2) * (phiHydF(k+1)-phiHydC(k))/(delz/2);
%      end
      
%    end

     mass = shelficemass (ix,iy);
     massFuncC = rhoConst * (phiHydC/gravity + zc);
     massFuncF = rhoConst * (phiHydF/gravity + zgp1);

     k = max (find ( massFuncF < mass ));
     if (isempty(k))
         k=0;
     end
     if (k>0)
     if (mass < massFuncC(k))
      topo(ix,iy) = -zg(k) - (mass-massFuncF(k)) * delz/2 / (massFuncC(k)-massFuncF(k));
     else
      topo(ix,iy) = -zc(k) - (mass-massFuncC(k)) * delz/2 / (massFuncF(k+1)-massFuncC(k));
     end
     end

  end
end

%mass = rhoConst * (phi0surf / gravity - topo);
%mass(:,1:100) = mass(:,1:100) + 917 * 10;
%topo(:,1:100) = bathy(:,1:100)+mwct;

etainit = zeros(size(topo));

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

icetopo2 = topo;

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)=0;
bathy(:,1)=0;

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);
fid = fopen('shelficemassinit.bin','w','b'); fwrite(fid,shelficemass,'real*8'); fclose(fid);
%fid = fopen('bathy_step.bin','w','b'); fwrite(fid,bathy,'real*8'); fclose(fid);

1	dgoldberg	1.1	%Verion of gendata.m modified by DNG
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=90;
13			delz = 10;
14
15			hfacMin = 0.2;
16			%mwct = 3;
17
18			dlat = 0.125/32; dy=dlat;
19			dlon = 0.125/4; dx=dlon;
20
21			%eos = 'linear';
22			eos = 'jmd95z';
23			% eos = 'mdjwf';
24
25			acc = 'real*8';
26
27			long = [-105.5:dlon:-105.5];
28			lonc = long+dlon/2;
29			latg = [-75.4457:dlat:-73.8809-dlat];
30			latc = latg+dlat/2;
31			size(latc);
32
33
34
35			dz = delz*ones(1,nz);
36			zgp1 = [0,cumsum(dz)];
37			zc = .5*(zgp1(1:end-1)+zgp1(2:end));
38			zg = zgp1(1:end-1);
39			dz = diff(zgp1);
40			% sprintf('delZ = %d * %7.6g,',nz,dz)
41
42			%%%%%%%%% stratification %%%%%%%%%%%%%%%%
43
44			T_sfc = -1.9;
45			T_bottom = 1.2;
46			del_T = (T_bottom - T_sfc)/(59*delz);
47
48			for iz = 1:nz;
49
50
51			tref(iz) = T_sfc + del_T((iz-30)delz);
52			if iz<=30;
53			tref(iz)=-1.9;
54			end
55			if iz>=90
56			tref(iz) =2;
57			end
58			end
59
60			S_sfc = 34.2;
61			S_bottom = 34.7;
62			del_S = (S_bottom - S_sfc)/(59*delz);
63
64			for iz = 1:nz;
65
66
67			sref(iz) = S_sfc + del_S((iz-30)delz);
68			if iz<=30;
69			sref(iz)=34.2;
70			end
71			if iz>=90
72			sref(iz) =34.7;
73			end
74			end
75
76			%%%%%%%%%%% density %%%%%%%%%%%%%%%%%
77
78			% Gravity
79			gravity=9.81;
80			rhoConst = 1000;
81			k=1;
82			dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
83			dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
84			p = abs(zc)gravityrhoConst*1e-4;
85			dp = p;
86			kp = 0;
87
88			Rho = zeros(nz,1);
89
90			while rms(dp) > 1e-13
91			phiHydF(k) = 0;
92			p0 = p;
93			kp = kp+1;
94			for k = 1:nz
95			switch eos
96			case 'linear'
97
98			case 'jmd95z'
99			drho = densjmd95(sref(k),tref(k),p(k))-rhoConst;
100			case 'mdjwf'
101			drho = densmdjwf(sref(k),tref(k),p(k))-rhoConst;
102			otherwise
103			error(sprintf('unknown EOS: %s',eos))
104			end
105			Rho(k) = drho+rhoConst;
106			phiHydC(k) = phiHydF(k) + dzm(k)gravitydrho/rhoConst;
107			phiHydF(k+1) = phiHydC(k) + dzp(k)gravitydrho/rhoConst;
108			end
109			switch eos
110			case 'mdjwf'
111			p = (gravityrhoConstabs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
112			end
113			dp = p-p0;
114			end
115
116			shelficemass = binread('h0.bin',3,200) * 917;
117
118			%phi0surf = zeros(nx,ny);
119			topo = zeros(nx,ny);
120
121			for ix=1:nx
122			for iy=1:ny
123			% k=max(find(abs(zg)<abs(topo(ix,iy))));
124			% if isempty(k)
125			% k=0;
126			% end
127			% if k>0
128
129			% dr = -zg(k) - topo(ix,iy);
130
131			% if (dr<delz/2)
132			% phi0surf(ix,iy) = phiHydF(k) + (delz/2-dr) * (phiHydC(k)-phiHydF(k))/(delz/2);
133			% else
134			% phi0surf(ix,iy) = phiHydC(k) + (dr-delz/2) * (phiHydF(k+1)-phiHydC(k))/(delz/2);
135			% end
136
137			% end
138
139			mass = shelficemass (ix,iy);
140			massFuncC = rhoConst * (phiHydC/gravity + zc);
141			massFuncF = rhoConst * (phiHydF/gravity + zgp1);
142
143			k = max (find ( massFuncF < mass ));
144			if (isempty(k))
145			k=0;
146			end
147			if (k>0)
148			if (mass < massFuncC(k))
149			topo(ix,iy) = -zg(k) - (mass-massFuncF(k)) * delz/2 / (massFuncC(k)-massFuncF(k));
150			else
151			topo(ix,iy) = -zc(k) - (mass-massFuncC(k)) * delz/2 / (massFuncF(k+1)-massFuncC(k));
152			end
153			end
154
155			end
156			end
157
158			%mass = rhoConst * (phi0surf / gravity - topo);
159			%mass(:,1:100) = mass(:,1:100) + 917 * 10;
160			%topo(:,1:100) = bathy(:,1:100)+mwct;
161
162			etainit = zeros(size(topo));
163
164			% new topography: icetopo rounded to the nearest k * deltaZ
165			% eta_init set to make difference
166
167			icetopo2 = topo;
168
169			for ix=1:nx
170			for iy=1:ny
171			k=max(find(abs(zg)<abs(icetopo2(ix,iy))));
172			if isempty(k)
173			k=0;
174			else
175
176			dr = 1-(-zg(k) - icetopo2(ix,iy))/delz;
177			if (dr > .25)
178			% bring Ro_surf up to closest grid face & make etainit negative
179			% to compensate
180			icetopo2(ix,iy) = -zg(k);
181			etainit(ix,iy) = (dr-1)*delz;
182			else
183			% bring Ro_surf down to closest grid face & make etainit pos
184			% to compensate
185			icetopo2(ix,iy) = -zg(k+1);
186			etainit(ix,iy) = (dr)*delz;
187			end
188
189			end
190			end
191			end
192
193			etainit(:,1)=0;
194			icetopo2(:,1)=0;
195			bathy(:,1)=0;
196
197			fid = fopen('shelftopo.round.bin','w','b'); fwrite(fid,icetopo2,'real*8'); fclose(fid);
198			fid = fopen('etainit.round.bin','w','b'); fwrite(fid,etainit,'real*8'); fclose(fid);
199			fid = fopen('shelficemassinit.bin','w','b'); fwrite(fid,shelficemass,'real*8'); fclose(fid);
200			%fid = fopen('bathy_step.bin','w','b'); fwrite(fid,bathy,'real*8'); fclose(fid);
201