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