Warning, /verification/shelfice_2d_remesh/input/gendata.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 9952f046 on 2019-12-17 17:59:00 UTC
9952f046d7 dngo*0001 % This is a matlab script that generates the input data
0002 % require matlab functions for equation of state
0003
0004 kwr=0;
0005 % Dimensions of grid
0006 nx=1;
0007 ny=200;
0008 nr=90;
0009 delz=10;
0010
0011 v0 = 2e3;
0012 h0 = 800;
0013
0014 hfacMin = 0.2;
0015
0016 dlat = 0.125/16; dy=dlat;
0017 dlon = 0.125/1; dx=dlon;
0018
0019 eos = 'jmd95z';
0020 prec = 'real*8';
0021
0022 long = 0;
0023 lonc = long+dlon/2;
0024 latg = -75.5+[0:ny-1]*dlat;
0025 latc = latg+dlat/2;
0026 zC=-delz*([1:nr]-0.5);
0027 zF=-delz*[0:nr];
0028 %size(latc)
0029
0030 % Gravity
0031 gravity= 9.81;
0032 rhoConst= 1030;
0033 rhoIce = 917;
0034
0035 % Nominal depth of model (meters)
0036 H = -900; %water depth in the ice shelf cavity
0037 Hmin = -600; % deepest point of cavern
0038 Hmax = -300; % shallowest point of cavern
0039 jEnd = ny*3/4; % where ice-shelf ends
0040 dHdy = (Hmax-Hmin)/dlat/(jEnd-2); %Slope of ice shelf
0041
0042 bathy = ones(nx,ny)*H; %For flat bathymetry: bathy = ones(nx,ny)*H;
0043 bathy(:,1) = 0;
0044
0045 j2=jEnd+1;
0046 hIce=Hmin+dHdy*[-1:ny-2]*dlat;
0047 hIce(1)=0; hIce(j2:ny)=0;
0048
0049 var=([1:ny]-2)/(jEnd-2);
0050 %var(2), var(jEnd)
0051 dMdt_fy=-cos(pi*var);
0052 dMdt_fy(1)=0; dMdt_fy(j2:ny)=0;
0053
0054 figure(1);clf;
0055 subplot(211)
0056 yax=[1:ny];
0057 plot(yax,hIce,'r-');
0058 hold on;
0059 plot(yax,bathy,'b-');
0060 hold off;
0061 grid
0062 title('ice-shelf & model depth [m]')
0063
0064 subplot(212)
0065 plot(yax,dMdt_fy,'b-');
0066 grid
0067 title('shape of ice-mass input')
0068
0069 namF='bathy_flat.bin';
0070 if kwr > 2,
0071 fprintf(' writing file: %s ...',namF);
0072 fid=fopen(namF,'w','b'); fwrite(fid,bathy,prec);fclose(fid);
0073 fprintf(' done\n');
0074 end
0075
0076 namF='shelficeTopo.Lin.bin';
0077 if kwr > 2,
0078 fprintf(' writing file: %s ...',namF);
0079 fid=fopen(namF,'w','b'); fwrite(fid,hIce,prec);fclose(fid);
0080 fprintf(' done\n');
0081 end
0082
0083 regMsk=ones(ny,1);
0084 regMsk(1)=0; regMsk(j2:ny)=2;
0085 namF='under_Ice_mask.bin';
0086 if kwr > 2,
0087 fprintf(' writing file: %s ...',namF);
0088 fid=fopen(namF,'w','b'); fwrite(fid,regMsk,prec);fclose(fid);
0089 fprintf(' done\n');
0090 end
0091
0092 %- rate of change due to ice-stream dynamics
0093 %rateDyn=rhoConst*0.1/86400; sfx='r01';
0094 rateDyn=rhoConst*0.1/3600; sfx='r02';
0095
0096 dMdt=rateDyn*dMdt_fy;
0097 namF=sprintf('%s.%s.%s','shelfice_dMdt',sfx,'bin');
0098 if kwr > 0,
0099 fprintf(' writing file: %s ...',namF);
0100 fid=fopen(namF,'w','b'); fwrite(fid,dMdt,prec);fclose(fid);
0101 fprintf(' done\n');
0102 end
0103
0104 dz = delz*ones(1,nr);
0105 zgp1 = [0,cumsum(dz)];
0106 zc = .5*(zgp1(1:end-1)+zgp1(2:end));
0107 zg = zgp1(1:end-1);
0108 dz = diff(zgp1);
0109 fprintf(' delZ = %d * %7.6g\n',nr,dz(1))
0110
0111 T_sfc = -1.9;
0112 T_bot = 2;
0113 del_T = (T_bot - T_sfc)/(59*delz);
0114 for k = 1:nr;
0115 tref(k) = T_sfc + del_T*((k-20)*delz);
0116 tref(k)= max(T_sfc,min(tref(k),T_bot));
0117 end
0118
0119 S_sfc = 34.2;
0120 S_bot = 34.7;
0121 del_S = (S_bot - S_sfc)/(59*delz);
0122 for k = 1:nr;
0123 sref(k) = S_sfc + del_S*((k-20)*delz);
0124 sref(k)= max(S_sfc,min(sref(k),S_bot));
0125 end
0126 pEOS=-rhoConst*gravity*zC; % in Pa
0127 pEOS=pEOS*1.e-4; % in dBar
0128 rhoAn=densjmd95(sref,tref,pEOS);
0129 rhoAn=rhoAn-rhoConst;
0130
0131 pF=-rhoConst*gravity*zF*1.e-4; % in dBar
0132 rhoUp=densjmd95(sref,tref,pF(2:end));
0133 rhoDw=densjmd95(sref,tref,pF(1:nr));
0134 dRho=rhoUp(1:nr-1)-rhoDw(2:nr);
0135 NSq=-gravity*dRho/delz/rhoConst;
0136
0137 mnV=min(NSq); MxV=max(NSq); Avr=mean(NSq);
0138 fprintf(' Stratif (N^2): min= %e , max= %e , Avr= %e\n',mnV,MxV,Avr);
0139
0140 zax=[1:nr];
0141
0142 v1=2.5e-2;
0143 var=1+nr-2*zax; var=var/(nr-1);
0144 vobc=v1*var;
0145
0146 figure(2);clf;
0147 subplot(131);
0148 plot(tref,-zax,'r-');
0149 grid
0150 title('tRef')
0151 subplot(132);
0152 plot(sref,-zax,'b-');
0153 grid
0154 title('sRef')
0155 subplot(133);
0156 %plot(vobc,-zax,'b-');
0157 plot(rhoAn,-zax,'b-');
0158 %plot(NSq*1.e+6,0.5-[2:nr],'r-');
0159 grid
0160 title('rhoAn')
0161
0162 if kwr > 2,
0163 namF='temp_obc.bin';
0164 fprintf(' writing file: %s ...',namF);
0165 fid=fopen(namF,'w','b'); fwrite(fid,tref,prec);fclose(fid);
0166 fprintf(' done\n');
0167 namF='salt_obc.bin';
0168 fprintf(' writing file: %s ...',namF);
0169 fid=fopen(namF,'w','b'); fwrite(fid,sref,prec);fclose(fid);
0170 fprintf(' done\n');
0171 namF='vVel_obc.bin';
0172 fprintf(' writing file: %s ...',namF);
0173 fid=fopen(namF,'w','b'); fwrite(fid,vobc,prec);fclose(fid);
0174 fprintf(' done\n');
0175 end
0176
0177 if kwr > 2,
0178 var=ones(ny,1)*tref; %size(var)
0179 namF='temp_ini.bin';
0180 fprintf(' writing file: %s ...',namF);
0181 fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
0182 fprintf(' done\n');
0183 %-
0184 var=ones(ny,1)*sref; %size(var)
0185 namF='salt_ini.bin';
0186 fprintf(' writing file: %s ...',namF);
0187 fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
0188 fprintf(' done\n');
0189 end
0190
0191 rhoAvr=rhoConst-1.345;
0192 fprintf(' convert Ice topo using rhoAve = %10.6f\n',rhoAvr);
0193 mIce0=-rhoAvr*hIce;
0194
0195 namF='shelficeMass.Lin.bin';
0196 if kwr > 1,
0197 fprintf(' writing file: %s ...',namF);
0198 fid=fopen(namF,'w','b'); fwrite(fid,mIce,prec);fclose(fid);
0199 fprintf(' done\n');
0200 end
0201
0202 return
