isomip/input.icefront/gendata.m

200c9094e5 Yun *0001 % This is a matlab script that generates the input data
                0002
                0003 % the configuation approximately the ISOMIP experiment no. 1
                0004 % require matlab functions for equation of state
                0005
                0006 % Dimensions of grid
                0007 nx=50; nxi=20;
                0008 ny=100; nyi=[51:100];
                0009 nz=30;
                0010 deltaZ = 30;
                0011
                0012 dlat = 0.1; dy=6.4e6*dlat./pi;
                0013 dlon = 0.3; dx=dlon;
                0014
                0015 %eos = 'linear';
                0016 eos = 'jmd95z';
                0017 %eos = 'mdjwf';
                0018
                0019 acc = 'real*8';
                0020
                0021 long = [0:dlon:10-dlon];
                0022 lonc = long+dlon/2;
                0023 latg = [-80:dlat:-70-dlat];
                0024 latc = latg+dlat/2;
                0025
                0026 % Nominal depth of model (meters)
                0027 H = -900;
                0028 Hmin = -700; % deepest point of cavern
                0029 Hmax = -200; % shallowest point of cavern
                0030 dHdx = (Hmax-Hmin)/4;
                0031
                0032 bathy = ones(nx,ny)*H;
                0033 bathy(1,:) = 0;
                0034 bathy(:,1) = 0;
                0035 %fid=fopen('bathy.box','w','b'); fwrite(fid,bathy,acc);fclose(fid);
                0036
                0037
                0038 dz = deltaZ*ones(1,nz);
                0039 zgp1 = [0,cumsum(dz)];
                0040 zc = .5*(zgp1(1:end-1)+zgp1(2:end));
                0041 zg = zgp1(1:end-1);
                0042 dz = diff(zgp1);
                0043 sprintf('delZ = %d * %7.6g,',nz,dz)
                0044
                0045 % Gravity
                0046 gravity=9.81;
                0047 rhoConst = 1030;
                0048 % compute potential field underneath ice shelf
                0049 talpha = 2e-4;
                0050 sbeta  = 7.4e-4;
                0051 tref = -1.9*ones(nz,1);
                0052 t    = tref;
                0053 sref = 34.4*ones(nz,1);
                0054 s    = sref;
                0055 gravity = 9.81;
                0056 k=1;
                0057 dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
                0058 dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
                0059 p = abs(zc)*gravity*rhoConst*1e-4;
                0060 dp = p;
                0061 kp = 0;
                0062 while std(dp) > 1e-13
                0063   phiHydF(k) = 0;
                0064   p0 = p;
                0065   kp = kp+1
                0066   for k = 1:nz
                0067     switch eos
                0068      case 'linear'
                0069       drho = rhoConst*(1-talpha*(t(k)-tref(k))+sbeta*(s(k)-sref(k)))-rhoConst;
                0070      case 'jmd95z'
                0071       drho = densjmd95(s(k),t(k),p(k))-rhoConst;
                0072      case 'mdjwf'
                0073       drho = densmdjwf(s(k),t(k),p(k))-rhoConst;
                0074      otherwise
                0075       error(sprintf('unknown EOS: %s',eos))
                0076     end
                0077     phiHydC(k)   = phiHydF(k) + dzm(k)*gravity*drho/rhoConst;
                0078     phiHydF(k+1) = phiHydC(k) + dzp(k)*gravity*drho/rhoConst;
                0079   end
                0080   switch eos
                0081    case 'mdjwf'
                0082     p = (gravity*rhoConst*abs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
                0083   end
                0084   dp = p-p0;
                0085 end
                0086
                0087 icetopo = ones(nx,1)*min(Hmax,Hmin + dHdx*(latc-latg(1)));
                0088 icetopo(:,nyi)=0;
                0089 fid=fopen('icetopo.exp1','w','b'); fwrite(fid,icetopo,acc);fclose(fid);
                0090 frontdepth=zeros(nx,ny);frontdepth(:,nyi(1))=-icetopo(:,nyi(1)-1);
                0091 fid=fopen('frontdepth.xuyun','w','b'); fwrite(fid,frontdepth,acc);fclose(fid);
                0092 frontcircum=zeros(nx,ny);
                0093  frontcircum(:,nyi(1))=deltaZ./dy;
                0094 fid=fopen('frontcircum.xuyun','w','b'); fwrite(fid,frontcircum,acc);fclose(fid);
                0095
                0096 % After modifying the code in calc_phi_hyd.F on Apr26,2012 this is the
                0097 % consistent way of computing phi0surf. For this, we need the grid
                0098 % information (hFacC's). For convenience, it's taken from a previous model
                0099 % run.
                0100 %
                0101 % The way of computing phi0surf consistent with code prior to Apr26,2012
                0102 % is recovered by setting drloc*dphi=0
96266b3693 Mart*0103 hf=rdmds('../tr_run.icefront/hFacC');
                0104 msk=sum(hf,3); msk(msk>0)=1;
200c9094e5 Yun *0105 phi0surf = zeros(nx,ny);
                0106 for ix=1:nx
                0107   for iy=1:ny
                0108     k=max(find(abs(zg)<abs(icetopo(ix,iy))));
                0109     if isempty(k)
                0110       k=0;
                0111     end
                0112     if k>0
                0113       kp1=min(k+1,nz);
96266b3693 Mart*0114       drloc=1-hf(ix,iy,k);
200c9094e5 Yun *0115       %drloc=(abs(icetopo(ix,iy))-abs(zg(k)))/dz(k);
                0116       dphi = phiHydF(kp1)-phiHydF(k);
                0117       phi0surf(ix,iy) = (phiHydF(k)+drloc*dphi)*rhoConst*msk(ix,iy);
                0118     end
                0119   end
                0120 end
                0121 fid=fopen(['phi0surf.exp1.' eos],'w','b'); fwrite(fid,phi0surf,acc);fclose(fid);
                0122