1 |
|
---|
2 | %Ok, start defining model parameters here
|
---|
3 |
|
---|
4 | %material parameters
|
---|
5 | md.g=9.8;
|
---|
6 | md.rho_ice=917;
|
---|
7 | md.rho_water=1023;
|
---|
8 | di=md.rho_ice/md.rho_water;
|
---|
9 | md.yts=365*24*3600;
|
---|
10 | md.heatcapacity=2009;
|
---|
11 | md.thermalconductivity=2.2; %W/mK
|
---|
12 | md.beta=9.8*10^-8;
|
---|
13 |
|
---|
14 | %Solution parameters
|
---|
15 | %parallelization
|
---|
16 | md.cluster='none';
|
---|
17 | md.np=2;
|
---|
18 | md.time=1;
|
---|
19 | md.exclusive=0;
|
---|
20 |
|
---|
21 | %statics
|
---|
22 | md.eps_rel=0.01;
|
---|
23 | md.eps_abs=10;
|
---|
24 | md.lowmem=1;
|
---|
25 | if md.numberofgrids<1000000,
|
---|
26 | md.sparsity=.001;
|
---|
27 | else
|
---|
28 | md.sparsity=.0001;
|
---|
29 | end
|
---|
30 |
|
---|
31 | %dynamics
|
---|
32 | md.dt=1*md.yts; %1 year
|
---|
33 | md.ndt=md.dt*10;
|
---|
34 | md.artificial_diffusivity=1;
|
---|
35 |
|
---|
36 | %control
|
---|
37 | md.control_type={'drag'}; %'drag', 'B'
|
---|
38 | md.nsteps=5;
|
---|
39 | md.tolx=10^-4;
|
---|
40 | md.maxiter=20;
|
---|
41 | md.optscal=10;
|
---|
42 | md.fit='logarithmic'; %'absolute','relative','logarithmic'
|
---|
43 | md.meanvel=1000/md.yts; %1000 meters/year
|
---|
44 | md.epsvel=eps;
|
---|
45 |
|
---|
46 |
|
---|
47 | disp(' creating thickness');
|
---|
48 | hmin=300;
|
---|
49 | hmax=1000;
|
---|
50 | ymin=min(md.y);
|
---|
51 | ymax=max(md.y);
|
---|
52 | md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
|
---|
53 | md.firn_layer=10*ones(md.numberofgrids,1);
|
---|
54 | md.bed=-di*md.thickness;
|
---|
55 | md.surface=md.bed+md.thickness;
|
---|
56 |
|
---|
57 | disp(' creating velocities');
|
---|
58 | md.vx_obs=zeros(md.numberofgrids,1);
|
---|
59 | md.vy_obs=zeros(md.numberofgrids,1);
|
---|
60 | md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
|
---|
61 |
|
---|
62 | disp(' creating drag');
|
---|
63 | md.drag_type=2; %0 none 1 plastic 2 viscous
|
---|
64 | md.drag=200*ones(md.numberofgrids,1); %q=1.
|
---|
65 | %Take care of iceshelves: no basal drag
|
---|
66 | pos=find(md.elementoniceshelf);
|
---|
67 | md.drag(md.elements(pos,:))=0;
|
---|
68 | md.p=ones(md.numberofelements,1);
|
---|
69 | md.q=ones(md.numberofelements,1);
|
---|
70 |
|
---|
71 | disp(' creating temperature');
|
---|
72 | md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
|
---|
73 |
|
---|
74 | disp(' creating flow law paramter');
|
---|
75 | md.B=paterson(md.observed_temperature);
|
---|
76 | md.n=3*ones(md.numberofelements,1);
|
---|
77 |
|
---|
78 | disp(' creating accumulation rates');
|
---|
79 | md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
|
---|
80 | md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
|
---|
81 |
|
---|
82 | %Deal with boundary conditions:
|
---|
83 |
|
---|
84 | disp(' boundary conditions for diagnostic model: ');
|
---|
85 | %Build gridonicefront, array of boundary grids belonging to the icefront:
|
---|
86 | gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
|
---|
87 | gridonicefront=double(md.gridonboundary & gridinsideicefront);
|
---|
88 |
|
---|
89 | md.gridondirichlet_diag=zeros(md.numberofgrids,1);
|
---|
90 | pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
|
---|
91 | md.dirichletvalues_diag=zeros(md.numberofgrids,2);
|
---|
92 |
|
---|
93 | pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
|
---|
94 | md.segmentonneumann_diag=md.segments(pos,:);
|
---|
95 | md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
|
---|
96 |
|
---|
97 | disp(' boundary conditions for prognostic model: ');
|
---|
98 | md.gridondirichlet_prog=zeros(md.numberofgrids,1);
|
---|
99 | md.dirichletvalues_prog=zeros(md.numberofgrids,1);
|
---|
100 | pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
|
---|
101 | md.segmentonneumann_prog=md.segments(pos,:);
|
---|
102 | md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
|
---|
103 | md.neumannvalues_prog(:)=NaN; %free radiation
|
---|
104 |
|
---|
105 | pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
|
---|
106 | md.segmentonneumann_prog2=md.segments(pos,:);
|
---|
107 | md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
|
---|
108 | md.neumannvalues_prog2(:)=NaN; %free radiation
|
---|
109 |
|
---|
110 | disp(' boundary conditions for thermal model: ');
|
---|
111 | md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
|
---|
112 | md.dirichletvalues_thermal=md.observed_temperature;
|
---|
113 | md.geothermalflux=zeros(md.numberofgrids,1);
|
---|
114 | pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
|
---|
115 |
|
---|
116 |
|
---|
117 |
|
---|
118 |
|
---|
119 | % Some Cielo code, ignore.
|
---|
120 | if strcmp(md.cluster,'yes')
|
---|
121 | ServerDisconnect;
|
---|
122 | end
|
---|