[5176] | 1 | %Ok, start defining model parameters here
|
---|
| 2 | disp(' creating thickness');
|
---|
[8305] | 3 | md.thickness=10*ones(md.numberofnodes,1);
|
---|
| 4 | md.bed=zeros(md.numberofnodes,1);
|
---|
[5176] | 5 | md.surface=md.bed+md.thickness;
|
---|
| 6 |
|
---|
| 7 | disp(' creating drag');
|
---|
[9610] | 8 | md.friction.coefficient=20*ones(md.numberofnodes,1); %q=1. %no drag is specified in the analytical solution
|
---|
| 9 | md.friction.p=ones(md.numberofelements,1);
|
---|
| 10 | md.friction.q=ones(md.numberofelements,1);
|
---|
[5176] | 11 |
|
---|
| 12 | disp(' creating temperatures');
|
---|
| 13 | tmin=238.15; %K
|
---|
| 14 | st=1.67*10^-2/1000; %k/m;
|
---|
[5181] | 15 | radius=sqrt((md.x).^2+(md.y).^2);
|
---|
[9425] | 16 | md.temperature=(tmin+st*radius);
|
---|
[9611] | 17 | md.basalforcings.geothermalflux=4.2*10^-2*ones(md.numberofnodes,1);
|
---|
[5176] | 18 |
|
---|
| 19 | disp(' creating flow law paramter');
|
---|
[8305] | 20 | md.rheology_B=6.81*10^(7)*ones(md.numberofnodes,1); %to have the same B as the analytical solution
|
---|
[5176] | 21 | md.rheology_n=3*ones(md.numberofelements,1);
|
---|
| 22 |
|
---|
[8399] | 23 | disp(' creating surface mass balance');
|
---|
| 24 | smb_max=0.5; %m/yr
|
---|
[5176] | 25 | sb=10^-2/1000; %m/yr/m
|
---|
| 26 | rel=450*1000; %m
|
---|
[9607] | 27 | md.surfaceforcings.mass_balance=min(smb_max,sb*(rel-radius));
|
---|
[5176] | 28 |
|
---|
| 29 | disp(' creating velocities');
|
---|
| 30 | constant=0.3;
|
---|
| 31 | md.vx_obs=constant/2*md.x.*(md.thickness).^-1;
|
---|
| 32 | md.vy_obs=constant/2*md.y.*(md.thickness).^-1;
|
---|
| 33 | md.vel_obs=(sqrt((md.vx_obs).^2+(md.vy_obs).^2));
|
---|
[8305] | 34 | md.vx=zeros(md.numberofnodes,1);
|
---|
| 35 | md.vy=zeros(md.numberofnodes,1);
|
---|
| 36 | md.vz=zeros(md.numberofnodes,1);
|
---|
| 37 | md.pressure=zeros(md.numberofnodes,1);
|
---|
[5176] | 38 |
|
---|
| 39 | %Deal with boundary conditions:
|
---|
| 40 | disp(' boundary conditions for diagnostic model: ');
|
---|
| 41 | md=SetMarineIceSheetBC(md,'../Exp/RoundFrontEISMINT.exp');
|
---|
| 42 |
|
---|
| 43 | radius=sqrt((md.x).*md.x+(md.y).*md.y);
|
---|
| 44 | pos=find(radius==min(radius));
|
---|
| 45 | md.x(pos)=0; md.y(pos)=0; %the closest node to the center is changed to be exactly at the center
|
---|
| 46 |
|
---|
[8824] | 47 | md.spcvx(pos)=0;
|
---|
| 48 | md.spcvy(pos)=0;
|
---|
| 49 | md.spcvz(pos)=0;
|
---|
[5176] | 50 |
|
---|
| 51 | %parallel options
|
---|
| 52 | md.ndt=50000;
|
---|
| 53 |
|
---|
| 54 | %Constants
|
---|
| 55 | md.rho_ice=910;
|
---|
| 56 | md.thermalconductivity=2.1;
|
---|
| 57 | md.latentheat=3.35*10^5;
|
---|
[9597] | 58 | md.beta=8.66*10^-4/(md.rho_ice*md.constants.g); %conversion from K/m to K/Pa
|
---|
| 59 | md.constants.yts=31556926;
|
---|