source: issm/trunk/test/Validation/PattynStokes/test1_iceshelf/Square.par@ 16

Last change on this file since 16 was 16, checked in by Mathieu Morlighem, 16 years ago

Added new tests form ice1

  • Property svn:executable set to *
File size: 3.6 KB
Line 
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.
120if strcmp(md.cluster,'yes')
121 ServerDisconnect;
122end
Note: See TracBrowser for help on using the repository browser.