isomip/input_ad/gendata.m

9a037d41df Patr*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;
                0009 nz=30;
                0010 deltaZ = 30;
                0011
                0012 dlat = 0.1; dy=dlat;
                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 rms(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 fid=fopen('icetopo.exp1','w','b'); fwrite(fid,icetopo,acc);fclose(fid);
                0089 fid=fopen('pload.exp1','w','b'); fwrite(fid,-icetopo,acc);fclose(fid);
6a65fe38f9 Mart*0090
                0091 % After modifying the code in calc_phi_hyd.F on Apr26,2012 this is the
                0092 % consistent way of computing phi0surf. For this, we need the grid
                0093 % information (hFacC's). For convenience, it's taken from a previous model
                0094 % run.
                0095 %
                0096 % The way of computing phi0surf consistent with code prior to Apr26,2012
                0097 % is recovered by setting drloc*dphi=0
                0098 g=rdmnc('grid.*','HFacC');
                0099 msk=sum(g.HFacC,3); msk(msk>0)=1;
9a037d41df Patr*0100 phi0surf = zeros(nx,ny);
                0101 for ix=1:nx
                0102   for iy=1:ny
                0103     k=max(find(abs(zg)<abs(icetopo(ix,iy))));
                0104     if isempty(k)
                0105       k=0;
                0106     end
                0107     if k>0
6a65fe38f9 Mart*0108       kp1=min(k+1,nz);
                0109       drloc=1-g.HFacC(ix,iy,k);
                0110       %drloc=(abs(icetopo(ix,iy))-abs(zg(k)))/dz(k);
                0111       dphi = phiHydF(kp1)-phiHydF(k);
                0112       phi0surf(ix,iy) = (phiHydF(k)+drloc*dphi)*rhoConst*msk(ix,iy);
9a037d41df Patr*0113     end
                0114   end
                0115 end
                0116 fid=fopen(['phi0surf.exp1.' eos],'w','b'); fwrite(fid,phi0surf,acc);fclose(fid);
                0117